Questions on coding my FFT&LLR implementation

[fivemack]

You want Colin Percival's paper "Rapid multiplication modulo the sum and difference of highly composite numbers" (http://www.daemonology.net/papers/fft.pdf)

[Prime95]

Originally Posted by fivemack

You want Colin Percival's paper "Rapid multiplication modulo the sum and difference of highly composite numbers" (http://www.daemonology.net/papers/fft.pdf)

gwnum uses a different method where k does not have to be highly composite. I think the SoB people wrote about the method recently.

[nuggetprime]

is gwnum's method faster? where is it explained(link)?

[philmoore]

Yes, George made a clever (in my opinion) innovation that improved on Colin Percival's algorithm. The paper describing it has been held up for a few months in the referee stage, but I will look to see if I can put together a short description that might interest readers of this thread. Give me a couple days, as I am preparing for the start of fall term next Monday.

[nuggetprime]

Thanks for your work philmoore.
I might code,with this, a very fast assembler free proggie for Riesel,Proth,N+1,N-1,Prp and other tests.

Regards,
nuggetprime

[philmoore]

I take no credit, but I did summarize work that George described in an email. Here is a summary:

George expanded on the work of Percival to adapt his routines to do computations modulo numbers of the form k*2^n +/- c, thus making his routines of particular value to Seventeen or Bust and Rieselsieve. Treatment of the +/-c term essentially follows the scheme of Percival, but the k*2^n term is treated in a way that allows use of smaller FFT sizes for some exponents than those required by Percival's method. This new scheme involves performing a Mersenne-like DWT on 2^(n+log k) +/- c (All logs denote logarithms in base 2) with weights ranging from 1 to 2 rather than using Percival's weights for k in the FFT word which range from 1 to k, saving approximately log k bits of weighting data, which become 2*log k bits in the inverse FFT word after point-wise squaring. The complication of this scheme is that the ''wrap around'' data is now divided by k, which is dealt with by multiplying each result word by k before rounding to an integer. Carry out of the top word is done very carefully, but is not a major problem. The procedure creates a result that has been multiplied by k, which is circumvented by dealing with numbers in (x / k) mod k*2^n +/- c format. Therefore, to square x, the transform of (x / k) is computed and squared component-wise, after which the inverse transform is taken, giving a result x^2 / k still in the same format, easily adapted to a long series of squarings such as is performed in a probable prime or Proth test. At conclusion of the computation, a final multiplication by k converts the result back into the desired form. The net savings in each inverse FFT word amounts to the equivalent of (log k) / 2 bits in the original input FFT word, allowing in many cases for the use of smaller FFT sizes for the k values under investigation.

[Xyzzy]

Where George works:

[jrk]

At least he has a decent taste in video games. :)

[nuggetprime]

Is the code I attatched acceptably fast for a tabled dft of 8000 unsigned short elements?
Remember, it's doing DFT,not FFT. And it's not doing inverse DFT.
I'm pretty sure it's not fast enough. If you also think so, can you tell me what can be improved? That would be a great help

--nuggetprime

[retina]

Originally Posted by Xyzzy

Where George works:

I think that image is a fake. Do you think we are all stupid? Nobody, and I mean nobody, has a workbench as clean as that in the real word. Sheesh! What's next? You'll be telling us that the food is less than 3 months old!

[nuggetprime]

The prev.code is totally broken. Sorry. Will attatch a new one probably tomorrow.

--Nuggetprime