double precision in LL tests

[drew]

Hi all,

In order to learn more of the math behind Lucas-Lehmer testing, I wrote my own C code implementation of the algorithm. It's working, but I think I may be losing some precision somewhere in the calculations, and I'm wondering what I should expect.

Through trial-and-error, I've discovered that the computation will succeed if the largest base of the fft array (prior to weighting) is 2^20 or less, and will fail if it is 2^22 or more. I haven't encountered a case, yet, that required 2^21, so the performace there is still unknown.

What I'm wondering is what amount of precision should I expect? Double-precision floating point format uses a 53-bit mantissa. I'd expect that to be cut in half because of the squaring, and then subtract perhaps a couple of bits to account for additions and rounding errors. That should still leave things at 24 bits or so, mayble 23, as far as I can tell.

Are my expectations unreasonable, or am I losing some precision in my code?

Thanks,
Drew

[jasonp]

Originally Posted by drew

What I'm wondering is what amount of precision should I expect? Double-precision floating point format uses a 53-bit mantissa. I'd expect that to be cut in half because of the squaring, and then subtract perhaps a couple of bits to account for additions and rounding errors. That should still leave things at 24 bits or so, mayble 23, as far as I can tell.

Are my expectations unreasonable, or am I losing some precision in my code?

Your experience mirrors that of other implementations; for latter-day size exponents I think 20-21 bits per array element is the best you can hope for, even with balanced digit representation.

jasonp

[drew]

Originally Posted by jasonp

Your experience mirrors that of other implementations; for latter-day size exponents I think 20-21 bits per array element is the best you can hope for, even with balanced digit representation.

jasonp

Thanks.

[ewmayer]

A more precise estimate, which also shows how the accuracy decreases (quite slowly though) with increasing FFT length can be found here - two out of three authors of that paper have posted to this thread, so hey, it must be true. ;)

[drew]

Originally Posted by ewmayer

A more precise estimate, which also shows how the accuracy decreases (quite slowly though) with increasing FFT length can be found here - two out of three authors of that paper have posted to this thread, so hey, it must be true. ;)

That's why I came here to ask.

FWIW, I switched to a balanced-digit representation, and now it works through 2^22, at least when testing the primes M11213, M21701, M44497 and M86243.

I now understand the cause of that limit...the depth of the FFT requires more additions of values, and I had been thinking in terms of decimal, not binary digits when I said it should only account for a couple of digits of precision. I expect that the loss of precision should grow *much* more slowly beyond the numbers I've already tested (up to the 34th Mersenne prime).

My next optimization will be to convert to a real-valued FFT, hopefully cutting the FFT size for most exponents in half, but since the data will be twice as dense, I expect to lose that bit I gained going to the balanced-digit representation.

Thanks,
Drew