What's a nice method to factor gaussian integers?

I already have prime factorization (for integers), but now I want开发者_StackOverflow中文版 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 i² = -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) = a² + b² .

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 a²+b² congruent to 1 mod 4 ,

a, where a is a prime congruent to 3 mod 4 , with norm a²

---

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, let q = (1+i). (Note that q = (1-i) would work equally well, since they only differ by a unit factor.)

If p = 3 mod 4, q = p. But the norm of q is p², 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 = a² + b², 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 k² = -1 mod p. Then k²+1 = 0 mod p, which is equivalent to saying that p divides k²+1 in the integers (and therefore also the Gaussian integers). In the Gaussian integers, k²+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 p₁, p₂ ... p_n.

For each remaining factor p:
   if p=2, u = (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 = 2² mod p = 4. 4 is congruent to -1 mod 5. Success! k = 2¹ = 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, 2⁸ = 1 mod 17, no good. Try n=3 instead. 3⁸ = 16 mod 17 = -1 mod 17. Success! So k = 3⁴ = 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, 2²⁶ = 52 mod 53 = -1 mod 53. Then k = 2¹³ 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?

---

Use floating point for the real and imaginary components if you want full single cell integer accuracy, and define gsub, gmul and a special division gdivr with rounded coefficients, not floored. That's because the Pollard rho factorization method needs gcd via Euclid's algorithm, whith a slightly modified gmodulo:

gmodulo((x,y),(x',y'))=gsub((x,y),gmul((x',y'),gdivr((x,y),(x',y'))))

Pollard rho

def poly((a,b),(x,y))=gmodulo(gsub(gmul((a,b),(a,b)),(1,0)),(x,y))

input (x,y),(a,b) % (x,y) is the Gaussian number to be factorized
(c,d)<-(a,b)
do 
   (a,b)=poly((a,b),(x,y))
   (c,d)=poly(poly((c,d),(x,y)),(x,y))
   (e,f)=ggcd((x,y),gsub((a,b),(c,d)))
   if (e,f)=(x,y) then return (x,y) % failure, try other (a,b)
until e^2+f^2>1
return (e,f)

A normal start value is a=1, b=0.

I have used this method programmed in Forth on my blog http://forthmath.blogspot.se

For safety, use rounded values in all calculations while using floating points for integers.