I already have prime factorization (for integers), but now I want to implement it for gaussian integers but how should I do it? thanks!
This turned out to be a bit verbose, but I hope it fully answers your question...
A Gaussian integer is a complex number of the form
G = a+bi
where i2 = -1, and a and b are integers.
The Gaussian integers form a unique factorization domain. Some of them act as units (e.g. 1, -1, i, and -i), some as primes (e.g. 1 + i), and the rest composite, that can be decomposed as a product of primes and units that is unique, aside from the order of factors and the presence of a set of units whose product is 1.
The norm of such a number G is defined as an integer: norm(G) = a2 + b2 .
It can be shown that the norm is a multiplicative property, that is:
norm(I*J) = norm(I)*norm(J)
So if you want to factor a Gaussian integer G, you could take advantage of the fact that any Gaussian integer I that divides G must satisfy the property that norm(I) divides norm(G), and you know how to find the factors of norm(G).
The primes of the Gaussian integers fall into three categories:
1 +/- i , with norm 2,
a +/- bi, with prime norm a2+b2 congruent to 1 mod 4 ,
a, where a is a prime congruent to 3 mod 4 , with norm a2
Now to turn this into an algorithm...if you want to factor a Gaussian integer G, you can find its norm N, and then factor that into prime integers. Then we work our way down this list, peeling off prime factors p of N that correspond to prime Gaussian factors q of our original number G.
There are only three cases to consider, and two of them are trivial.
If p = 2, try u = (1+i) . If G/u has remainder 0, q=u, else q=(1-i).
If p = 3 mod 4, q = p. But the norm of q is p2, so we can strike another factor of p from the list of remaining factors of norm(G).
The p = 1 mod 4 case is the only one that's a little tricky to deal with. It's equivalent to the problem of expressing p as the sum of two squares: if p = a2 + b2, then a+bi and a-bi form a conjugate pair of Gaussian primes with norm p, and one of them will be the factor we're looking for.
But with a little number theory, it turns out not to be too difficult. Consider the integers mod p. Suppose we can find an integer k such that k2 = -1 mod p. Then k2+1 = 0 mod p, which is equivalent to saying that p divides k2+1 in the integers (and therefore also the Gaussian integers). In the Gaussian integers, k2+1 factors into (k+i)(k-i). p divides the product, but does not divide either factor. Therefore, p has a nontrivial GCD with each of the factors (k+i) and (k-i), and that GCD or its conjugate is the factor we're looking for!
But how do we find such an integer k? Let n be some integer between 2 and p-1 inclusive. Calculate n(p-1)/2 mod p -- this value will be either 1 or -1. If -1, then k = n(p-1)/4, otherwise try a different n. Nearly half the possible values of n will give us a square root of -1 mod p, so it won't take many guesses to find a value of k that works.
To find the Gaussian primes with norm p, just use Euclid's algorithm (slightly modified to work with Gaussian integers) to compute the GCD of (p, k+i). That gives one trial divisor. If it evenly divides the Gaussian integer we're trying to factor (remainder = 0), we're done. Otherwise, its conjugate is the desired factor.
Euclid's GCD algorithm for Gaussian integers is almost identical to that for normal integers. Each iteration consists of a trial division with remainder. If we're looking for gcd(a,b),
q = floor(a/b), remainder = a - q*b, and if the remainder is nonzero we return gcd(b,remainder).
In the integers, if we get a fraction as the quotient, we round it toward zero.
In the Gaussian integers, if the real or imaginary parts of the quotient are fractions, they get rounded to the nearest integer. Besides that, the
algorithms are identical.
So the algorithm for factoring a Gaussian integer G looks something like this:
Calculate norm(G), then factor norm(G) into primes p1, p2 ... pn.
For each remaining factor p:
if p=2, u = (1 + i).
If G/u has remainder 0, q=u, else q=(1-i)
strike p from the list of remaining primes.
else if p mod 4 = 3, q = p, and strike 2 copies of p from the list of primes.
else find k such that k^2 = -1 mod p, then u = gcd(p, k+i)
if G/u has remainder 0, q = u
else q = conjugate(u)
strike p from the list of remaining primes.
Add q to the list of Gaussian factors.
Replace G with G/q.
endfor
At the end of this procedure, G is a unit with norm 1. But it's not necessarily 1 -- it could be -1, i, or -i, in which case add G to the list of factors, to make the signs come out right when you multiply all the factors to see if the product matches the original value of G.
Here's a worked example: factor G = 361 - 1767i over the Gaussian integers. norm(G) = 3252610 = 2 * 5 * 17 * 19 * 19 * 53
Considering 2, we try q = (1+i), and find G/q = -703 - 1064i with remainder 0.
G <= G/q = -703 - 1064i
Considering 5, we see it is congruent to 1 mod 4. We need to find a good k. Trying n=2, n(p-1)/2 mod p = 22 mod p = 4. 4 is congruent to -1 mod 5. Success! k = 21 = 2. u = gcd(5, 2+i) which works out to be 2+i. G/u = -494 - 285i, with remainder 0, so we find q = 2+i.
G <= G/q = -494 - 285i
Considering 17, it is also congruent to 1 mod 4, so we need to find another k mod 17. Trying n=2, 28 = 1 mod 17, no good. Try n=3 instead. 38 = 16 mod 17 = -1 mod 17. Success! So k = 34 = 13 mod 17. gcd(17, 13+i) = u = 4-i, G/u = -99 -96i with remainder -2. No good, so try conjugate(u) = 4+i. G/u = -133 - 38i with remainder 0. Another factor!
G <= G/(4+i) = -133 - 38i
Considering 19, it is congruent to 3 mod 4. So our next factor is 19, and we strike the second copy of 19 from the list.
G <= G/19 = -7 - 2i
Considering 53, it is congruent to 1 mod 4. Again with the k process... Trying n=2, 226 = 52 mod 53 = -1 mod 53. Then k = 213 mod 53 = 30. gcd(53,30+i) = u = -7 - 2i. That's identical to G, so the final quotient G/(-7-2i) = 1, and there are no factors of -1, i, or -i to worry about.
We have found factors (1+i)(2+i)(4+i)(19+0i)(-7-2i). And if you multiply that out (left as an exercise for the reader...), lo and behold, the product is 361-1767i, which is the number we started with.
Ain't number theory sweet?