Possible optimization for GIMPS

[alpertron]

Please correct me if I'm wrong, because I have not programmed FFT.

The main loop of LL proceeds by squaring the number and subtracting 2 mod 2^p-1. The square is done by peforming two transforms, a convolution and an inverse transform using IBDWT.

I think that by using Fermat tests (in base 3 for example), where only exponentiations are needed, we could eliminate the need for transforms and inverse transforms. So only convolutions are done in the main loop which are a lot faster.

Then we could perform LL only for the pseudoprimes (a tiny fraction of the tests).

Is this approach correct?

[Dubslow]

I'm assuming you mean this: http://en.wikipedia.org/wiki/Fermat_primality_test

Note two things: First, you must take a^(x-1), where x itself is 2^p-1 for large p. Not a smart idea.

Second: How do you do an exponentiation?

Originally Posted by Fermat test article

Using fast algorithms for modular exponentiation, the running time of this algorithm is O(k × log² n × log log n × log log log n), where k is the number of times we test a random a, and n is the value we want to test for primality.

Note that we are testing n=2^p-1, so log n~p.
The LL test with FFT muls has complexity O(p² × log p × log log p), which is about the same as the same as the Fermat test once you substitute log n~p (and the LLT is deterministic).

You do exponentiation with lots of squaring and multiplication. After all, exponentiation is repeated multiplications by definition, and the LL test is already down to the relatively cheap p-2 squarings per test.

Originally Posted by exp by square

A brief analysis shows that such an algorithm uses log n squarings and at most log n multiplications.

Substituting log n ~ p, then we have p squarings plus some other multiplications, up to p of them. So p-2 squarings is really as low as we can get.

(I'm not confident of this last part, but it's possible that the LL test can be derived from the Fermat test with base 2 or something, and some of the specific algorithms for exponentiation.)

PS "only convolutions" is meaningless with respect to multiplication without doing a transform. It's like saying 6235*4865 == 6*4 2*8 3*6 5*5 == 24 16 18 25 == (with carrying) 25805, which is patently not the actual product of 30333275.

[alpertron]

My idea is to perform one transform at the beginning, then perform the convolutions for each loop and then the inverse transform at the end. I do not want to reduce the number of multiplications but the timing for each multiplication.

Since multiple convolutions will generate very big numbers, we could separate the exponent from the mantissa.

Edit: I do not know if this approach would work because of floating point approximation errors. Maybe we could do inverse and direct transforms every few iterations.

[axn]

Originally Posted by alpertron

Is this approach correct?

Short answer: no.

Long answer: You should program DFT (which is mathematically same as FFT, but much less efficient) and observe how things work.