What?

If you're completely in the dark - browse around Yves Gallot's GFN search pages for more info.

How?

Firstly I must thank Jim Fougeron for so quickly implementing my new algorithm in his GFNSieve, which is ultra-fast x86 assembly. I spent many a month running that software. However, since the final 2x optimisation (see below), David Underbakke took the torch from Jim, and implemented the most up-to-date algorithm in his Athlon- and P4-optimised AthGFNSieve. This is the only software I use currently. (See the news section for a fuller story.)
Oh, I'm no longer tempted to race it on a big cluster of Alphas :-(!

The algorithm

No matter what people say, with all the trolls and cranks, Usenet is still an amazing place with some amazingly helpful people. I'm thinking of sci.math, which is basically the place that my algorithm was forged. By name, I have to thank Paul Pollack for a wonderfully lucid mail in response to my original post which answered every part of my question in infinite detail, and answered questions I'd not even asked yet, and to Robin Chapman too.

Here're the guts:

Factors of b^(2^n)+1 are of the form k.2^(n+1)+1. So for each prime of the form p=k.2^(n+1)+1 (from a sieve) we want to satisfy the equation

  b^(2^n)+1 == 0  (mod p)
  b^(2^n)   == -1 (mod p)

i.e. we are looking for the (2^n)-th roots of -1 (mod p). i.e. we are looking for the primitive (2^(n+1))-th roots of 1 (mod p), explicitly those that are (2^(n+1))-th roots but not (2^n)-th roots.

Now (2^i)-th roots of 1 are easy to find modulo numbers of the form k.2^(n+1)+1, because from Fermat's Little Theorem we have

  a^(p-1)       == 1 (mod p) 
  a^(k.2^(n+1)) == 1 (mod p)

  (a^k)^(2^(n+1)) == 1 (mod p)

i.e. /every/ a (mod p) can be used to generate a (2^(n+1))-th root, simply take a^k, and you have one.

However, not all such a's generate _primitive_ roots. Fortunately, 50% of them do, and it's trivial to find out which! We are of course still looking for the (2^n)-th roots of -1, i.e.

  (a^k)^(2^n) == -1 (mod p)

So we obviously want to check by calculating (a^k)^(2^n) (mod p), which can yield one of several values:

  (a^k)^(2^n) (mod p) =
    -1    : We have a primitive (2^(n+1))-th root of unity :-)
    1     : We have a non-primitive root of unity, try another a
    other : p wasn't prime!

In fact, (a^k)^(2^n) is a _strong pseudoprimality test_ for p. Every non-primitive root provides us with a Euler (or SS) PRP test, but the primitive root when found has actually caused a Strong (or MR) PRP test to be performed. This means that our prime generator is allowed to send us the occasional composite, as we will be able to detect that easily and move on to a different 'prime'!

Anyway, say we have c=(a^k) such that c^(2^n)==-1 (mod p) Then

  (c^3)^(2^n) == (c^(2^n))^3 == (-1)^3 == -1 (mod p)

In fact for all odd exponents this is true, i.e.

  (c^(2i+1))^(2^n) == -1 (mod p), for integer i

i.e. once we have generated _one_ primitive root, c1, using on average two trials, we can generate all (2^n) of them simply by multiplying by (c1^2), to generate (c1^3), (c1^5), ...

And those are _all_ of the primitive roots, and are the values, mod p, which we want to strike out from our sieve. (For small p, if c is a solution, then so is c+p, like a traditional modular sieve).

Now we halve the work-factor, courtesy of Bernhard Frey!
If we look at what is being calculated exactly half way through the loop above that calculates the 2^n primitive roots, we see that:

  c1^(2*(2^(n-1))+1) == c1^(2^n+1) == c1^(2^n)*c1 == -1*c1 == -c1

That means that we only need to caclulate _half_ of the primitive roots, and simply look at the additive inverse of the first 2^(n-1) of them in order to generate the other half of them.

The cost of the sieve technique is as follows: Assuming the small prime sieve returns a prime 1 time in s (s will be <2 hopefully, I'll assume s=2)

  Generation of primitive root = s * ~2 * (1.5 lg(k)) mulmods
  Generation of all roots      =          2^(n-1)     mulmods
  Total = 6.lg(k) + 2^(n-1)

(for 50 bit primes and large n (>12) 6.lg(k) is dwarfed by 2^(n-1).)

The cost of a traditional trial factoring technique is

  Trial division               = #candidates * (1.5 lg(k) + n) mulmods

(for 50 bit primes (1.5 lg(k) + n) > 60)

So it's clear that with a large number of candidates, the sieving technique is a win. Very roughly:

  Sieve time  < Trial division time
when 
  2^(n-1)     > #candidates * 60

Implementation issues of course can change this ratio hugely, of course, and the above should be taken with a pinch of salt. I'll put up some real data when I have it.

I appear to have a rough analysis for 131072 now. I will try to generalise for the other ranges with what I've learnt from doing that analysis.

Another hastily constructed page by Phil Carmody
Home / Maths / GFN Sieve / maths.html