Mixed floating/integer transform for faster LL testing

[ewmayer]

David Willmore (hi, David!) wrote (in the Software forum):

if I remember right, Ernst Mayer is still working on the dual floating/integer LL tester program. If I understand his approach correctly, he's doing the 'hyperthreading' manually--by interleaving the code for the two tests.

I think he expects some speed improvement. Is that just due to the processors he's targeting? AXP and such?

The short answer to the latter question is, I expect a speed improvement by making use of the integer execution units, which are largely idle when the machine is doing a floating-point FFT.

But it's far from being as simple as "Duh - just slap a modular transform in there next to the FFT" - long-winded explication follows. Math-o-phobes, beware! It's not as bad as a typical posting to the Number Theory mailing list, but it does require some patience and basic algebra.

* * *

Indeed, I am still planning on implementing a mixed floating-integer transform, but currently my top priority is finishing work on a releasable version of a C implementation of my Mlucas code. There's not terribly much work remaining to be done on the basic algorithmic source code, but my for-pay job doesn't leave me oodles of time to work on this.

The idea behind the mixed floating/integer (more precisely, floating/ modular) algorithm is this: do 2 transforms side-by-side. One is the usual approximate-arithmetic floating FFT, the other is an exact modular transform over some suitably chosen ring of integers. For those of you not familiar with commutative rings (and I am far from an expert in this area myself), the most important things to realize here are the following:

1) A ring defined by a real-integer prime modulus q has group order q-1; a ring defined by gaussian-integer (i.e. a complex pair a + b*I, with a and b separately modded) modulus q has group order q^2-1;

2) For length-N complex transform, we need N roots of unity. For complex FFT this is never a problem (since the complex numbers admit roots of unity of any desired order), but in the modular case the modulus must be chosen which admits a desired variety of runlengths (see (4) below). Typically we want a variety of runlengths N = k*2^n, where k is a small odd number and n is decently large, say at least 20 or so. In order to be able to split each interval between adjacent power-of-2 runlengths into 4 equal-length pieces (like Prime95 does), we need k = 1,3,5,7; in order to be able to split each interval into 8 equal-length pieces (like Mlucas does), we need k = 1,3,5,7,9,11,13,15.

3) For a Crandall/Fagin-style Mersenne-mod DWT to be built on top of the FFT, we also need N roots of 2. This latter requirement wrecks many schemes based on cleverly-found primes having the needed roots of unity.

4) For a modular transform of length N to be possible, N must divide the group order. For example, Nick Craig-Wood's modular code for the StrongArm CPU uses the prime q = 2^64 - 2^32 + 1, which defines a ring having multiplicative group order q-1 = 2^64 - 2^32 = 2^32*(3.5.17.257.65537). Thus, this prime admits power-of-two components of the runlength as large as 2^32, as well as the small odd primes k = 3 and 5, but not k = 7,9,11,13,15. This modulus also lends itself to fast hardware operations, another desirable property.

* * *

The implementation I am planning on will do modular arithmetic using the ring of Gaussian integers GF(q^2), which is a fancy term for gaussian (complex) integers modulo q. Here, q = 2^p-1 itself will be a Mersenne prime, so this will be an interesting recursive use of Mersenne primes - we use arithmetic modulo a small one to help us look for a big one. These integers have many nice properties - modding is very fast, and their lower-order power-of-2 roots of unity look just like those for the floating FFT, i.e. the modular transform will look very much like the floating FFT. What about the needed roots of 1 and 2?

Roots of 2 turn out to be easy for Mersenne-mod arithmetic, since (2^j)^N == 2^(j*N) == 2^(j*N mod p) modulo q. That is, for 2^j to be an Nth root of 2 modulo q, we need j*N == 1 modulo p. Since p is prime, arithmetic modulo p defines a field, which means we are guaranteed to be able to find such a multiplicative inverse of N modulo p. (Example: for p = 61 and runlength N = 3*2^18 = 768K, we have that j*N mod p == j*(N mod p) mod p = j*3 mod p, and it's easy to see that j = 41 is the multiplicative inverse of 3 mod p, so 2^41 is our primitive Nth root of 2.) In other words, all the roots of 2 needed for the DWT are themselves powers of 2, so multiplies by them can be effected simply via shift, mask and add operations.

How about roots of unity? For q = 2^p-1, the group order is q^2-1 = (q+1)*(q-1) = 2^(p+1)*[2^(p-1)-1]. That means all powers of two up to 2^(p+1) are allowed in the transform length, as well as all factors of the composite 2^(p-1)-1. For p > 2, his latter integer is guaranteed to always have 3 as a factor (i.e. we can always at least split each interval between adjacent power-of-2 runlengths into 2 equal-length pieces); the full factorization for some of the smaller Mersennes is:

p       p-1       factors of 2^(p-1)-1
------  ------    --------------------
3       2         3
5       22        3*5
7       2*3       32*7
13      22*3      32*5*7*13
17      24        3*5*17*257
19      2*32      33*7*19*73
31      2*3*5     32*7*11*31*151*331
61      22*3*5    32*52*7*11*13*31*41*61*151*331*1321
89      23*11     3*5*17*23*89*353*397*683*2113*P10
107     2*53      3*107*6361*69431*20394401*P14
127     2*32*7    33*72*19*43*73*127*337*5419*92737*649657*P11

Note that M61 is the smallest modulus that allows all of 3,5,7,9,11,13,15; the next-smaller modulus allowing all of these odd k's in the runlength is M2281, which isn't exactly a convenient modulus for hardware arithmetic. The two most-tempting moduli with respect to fast hardware arithmetic are M31 and M61. M31 needs double-width integer products 62 bits wide, which can be done on a large variety of hardware and fit nicely in a 64-bit integer. OTOH, M31 doesn't allow all the odd radices we'd like. M61 allows the best variety of transform lengths, but needs a double-wide integer product. This brings up a key issue: the speed of hardware integer multiply is the key determinant of whether this approach has a chance of competing with an all-floating code. That's because unlike floating multiply, where most popular RISC and Intel CPU s have similar capabilities (typically one pipelined FMUL per clock cycle, or for the pre-P4 Intel CPUs, every other cycle), their integer multiply capabilities differ vastly. For example, the Sparc (up to and including the latest Ultra-3 chips) is extremely good at floating-point and simple integer operations (add, sub, shift, and, or, xor, not), but is simply dismal at integer multiply - it uses a 2-bit-at-a-time Booth algorithm, so it needs over 30 cycles to generate the lower half of a full 64x64-bit integer product. At the other end of the scale, the Alpha has always been good at IMUL - it has hardware instructions for getting both the upper and lower halves of a 64x64-bit integer product, for instance - but these are fully pipelined only beginning with the third-generation 21264 CPU (a.k.a. ev6). Since we still need 2 separate instructions to get a full 128-bit product, this has pipelined cycle count of 2, half as fast as FMUL, but most optimized FFT code is dominated by adds, so this should be OK. Intel IA64 (Itanium) also can generate a 128-bit integer product every 2 cycles. My initial implementation will use a gaussian integer transform with q = 2^61, and will target these latter two processors, since they are best-of-breed. Basically, if one can't get the hybrid transform competitively fast on these 2 CPUs, it's probably not worth trying on lesser hardware.
* * *

At this point, the astute reader who knows a bit about all-integer transforms may be asking, "But won't the floating FFT and modular transform simply return the same set of outputs? What's the use of that?" It is true that if we do either by itself, we must limit our input digit sizes so that each floating or integer output digit contains all of the bits of the corresponding exact convolution output digit. For 64-bit hardware arithmetic, that typically limits us to something around 20 bits per input digit, i.e. to do an LL test of 2^p-1 we need a transform length N ~= p/20. The key thing that doing a floating and modular transform side-by-side allows is MUCH LARGER INPUTS, because even if the bottom log2(q) bits of each floating output have been contaminated by roundoff error, we can use the modular output (which is exact modulo q) to exactly reconstruct the messed-up lower bits. How much larger can the inputs be? Well, convolution outputs scale as the square of the inputs, so e.g. for q = 2^61-1, our inputs can be roughly sqrt(q), i.e. 30.5 bits larger. That means a transform length less than 40% that an all-floating code would need to test the same number for primality. In other words, even if the parallel-transform scheme winds up needing twice as many clock cycles for a given transform length (which could happen due to less-than-optimal instruction scheduling, especially with HLL code), it will still be a small gain. In an ideal scenario, we would get all of the added integer instructions "for free" and the resulting code would run more than 2X as fast as all-floating. I expect the actual performance to lie somewhere between these two values. Of course, it could also happen that the compiler does a horrible job and the parallel scheme winds up being so slow it's not worth doing (that's what happened to my first prototype of this scheme, written in Fortran-90 code about 3 years ago, and lacking lots of algorithmic tricks i've learned since.) In that case one would need to start writing assembly - I hope to be able to get away with no more than a handful of key assembly-code macros, since I don't have the time to write reams of ASM for various platforms.

[Dresdenboy]

Replying to a 13 months old thread

I'm currently working on mixed FGT/FFT - starting from 2 sources available by carey bloodworth. First tests (of reconstructing digits) are very promising.

On the GIMPS Opteron I also did first tests of fast multiply chains with a throughput of 1 every 2 cycles. It requires quick calculation of the modulo right after the product is available that the multiplier can store the next product into rax:rdx. It is necessary to load something into rax before a imul to break the dependency chain (and to multiply 2 new values ).

The LEA instruction would have been fine for shifting one value and adding another but it doesn't set any flags.

I think that will be a different hybrid transform (at least soucecode wise). I'm open for exchanging experiences.

[LoKI.GuZ]

Let me take this opportunity to ask: i have been wondering on a similar optimization: doing a multiprime fermat style fft, using both the ALU and the FPU at the same time... this way i would be able to do 2 transforms "simultaneouly"... is that possible?? or is it completely non sense??

[ewmayer]

The 2nd edition of Crandall & Pomerance adds a nice little exercise regarding primitive roots in F^*_{q[sup]2[/sup]} (the multiplicative group of complex (Gaussian) integers modulo q), where q=2^p-1 is a Mersenne prime - the first edition only mentioned the well-known Creutzburg-Tasche closed-form formula (Theorem 9.5.20 in the 1st edition; 9.5.21 in the 2nd) yielding roots of any desired power-of-2 order in F^*_{q[sup]2[/sup]}, but I whined to Richard Crandall to say something about roots of the full order q²-1, which may have prompted him to add the following exercise, relating primitive roots in F^*_{q[sup]2[/sup]} to those (of order q-1) in the real multiplicative field F^*_q - note that I took the liberty of swapping p and q in the problem statement, to make the notation agree with that in my original post in this thread:

Exercise 9.61: Show that if a+bi is a primitive root of maximum order q²-1 in F^*_{q[sup]2[/sup]} (with q == 3 (mod 4), so that "i" exists), then a²+b² must be a primitive root of maximum order in F^*_q. Is the converse true?

Give some Mersenne primes q=2^p-1 for which 6+i is a primitive root in F^*_{q[sup]2[/sup]}.

Hint for proving the first part of the exercise: First prove a useful property of (a+bi)^q. Multiplying the result by a+bi yields the (q+1)th power - if a+bi is a primitive root (order q²-1) in F^*_{q[sup]2[/sup]}, it follows that (a+bi)^(q+1) is a primitive root in F^*_q.

I would add the following, as well:

If a+bi is a primitive root in F^*_{q[sup]2[/sup]}, what can you say about b+ai?

What special property do roots of odd order in F^*_{q[sup]2[/sup]} have, and how might this be algorithmically useful?

[ewmayer]

Revival of this long-dormant thread by way of a quiz question for our math-o-phile readers, related to the search for primitive roots in F^*_{q[sup]2[/sup]}, that is identifying roots of the full order q²-1:

True or False: To test primitivity of a candidate root R = a + bi it suffices to test whether R^(q[sup]2-1)/2[/sup] == -1 (mod q).

[R.D. Silverman]

Originally Posted by ewmayer

Revival of this long-dormant thread by way of a quiz question for our math-o-phile readers, related to the search for primitive roots in F^*_{q[sup]2[/sup]}, that is identifying roots of the full order q²-1:

True or False: To test primitivity of a candidate root R = a + bi it suffices to test whether R^(q[sup]2-1)/2[/sup] == -1 (mod q).

Should I offer a hint?

[ewmayer]

Originally Posted by R.D. Silverman

Should I offer a hint?

Feel free to comment as desired, in the "help steer the interested reader toward the answer" sense, preferably.