Other mersenne progs, why so slow?

[Cyclamen Persicum]

I'm not a Gimps member, because the
probability to be a winner is zero, zero, module zero and
and zero periodically...
I have tried other progs available in Internet, but they as slow as continent drift.
When I wrote my own, but it still 10 or 100 times slower.
I don't know why and what is the secret of Prime95?
It seems the 39th number was discovered without any SSE
for only 1.5 month. But my prog requier 41 year to check it!!!
I use two FFT for squaring, and modular reduction
which worth nothing for 11111111...b.

[cheesehead]

The "secret" of Prime95 is twofold -- hardware and software.

Intel Pentium CPUs execute floating-point (FP) instructions more rapidly and with higher precision than most other manufacturers' CPUs. They can pipeline as fast as one FP instruction completion per clock cycle, whereas many others can execute FP no faster than two, three, or more cycles per instruction completed.

In addition, Intel FP units have 80-bit internal registers, longer than the 64-bit registers common on other brands. With proper programming to keep computations inside those internal registers as long as possible before the eventual conversion to a 64-bit external result, Intel CPUs can produce 16 more bits of numerical result per FP instruction than other brands.

George Woltman has intensively studied the details of Intel CPUs and worked for years to maximize the throughput of Prime95 software, programming the most speed-critical modules in assembler language in order to carefully interleave machine instructions so as to keep the various parts of the CPU as busy as possible churning out results with minimum wasted time.

(I'm quite envious of George's work because I've done similar software optimization for speed, in the past, on older computers, but by the time I discovered what George was up to with Prime95, he was waaay ahead of where I'd have started. :) )

[ewmayer]

Originally Posted by cheesehead

The "secret" of Prime95 is twofold -- hardware and software.

Intel Pentium CPUs execute floating-point (FP) instructions more rapidly and with higher precision than most other manufacturers' CPUs. They can pipeline as fast as one FP instruction completion per clock cycle, whereas many others can execute FP no faster than two, three, or more cycles per instruction completed.

Pretty much all RISC CPUs (Alpha, MIPS, Sparc, PowerPPC) can also do one pipelined add and one mul per cycle. No advantage for Intel there, except that it's by far the most common CPU out there, hence making it worthwhile for someone with the time and skill to write very highly optimized assembly code in large quantities for. Because of the x86's CISC nature and dumb stack-based instruction pipeline, it seems necessary to write ASM to get anywhere near the theoretical maximum throughput out of the chip (especially with FPU-intensive, data-nonlocal algorithms like the FFT), but once that's done, the throughput is as good as most any high-end RISC chip can do. The advantage of the RISCs is that they can get close to the maximal performance starting with high-level code, i.e. you don't need to spend a lot of time writing ASM to get good performance out of them.

In fact George's code is not better on a per-cycle basis than other top-quality LL test codes (Guillermo Valor's Glucas and my Mlucas) - on an Alpha ev6 I get even better per-cycle performance than George's code does on the P4 - but P4's are way cheaper and more numerous, and kudos to George for getting such good performance out of a chip with smaller caches and lower bandwidth to main memory. Also, the speed demons at Intel have been able to ratchet up chip speeds much faster than the RISC manufacturers.

In addition, Intel FP units have 80-bit internal registers, longer than the 64-bit registers common on other brands. With proper programming to keep computations inside those internal registers as long as possible before the eventual conversion to a 64-bit external result, Intel CPUs can produce 16 more bits of numerical result per FP instruction than other brands.

Yes, but because of the x86's small register set and the fact that the FFT's we do are so large, this gives only a small accuracy enhancement - perhaps 1-2% over an IEEE64-compliant CPU.

George Woltman has intensively studied the details of Intel CPUs and worked for years to maximize the throughput of Prime95 software, programming the most speed-critical modules in assembler language in order to carefully interleave machine instructions so as to keep the various parts of the CPU as busy as possible churning out results with minimum wasted time.

No argument there - George has squeezed remarkable performance out of these cheap commodity CPUs.

(I'm quite envious of George's work because I've done similar software optimization for speed, in the past, on older computers, but by the time I discovered what George was up to with Prime95, he was waaay ahead of where I'd have started. :) )

Yes, but you surely learned lots of new stuff along the way, so it has its rewards even if your code doesn't discover the next Mersenne prime. :)

[cperciva]

Originally Posted by ewmayer

Yes, but because of the x86's small register set and the fact that the FFT's we do are so large, this gives only a small accuracy enhancement - perhaps 1-2% over an IEEE64-compliant CPU.

Are you sure? It looks to me like it should give an advantage of at least 1.5 bits.

[Prime95]

Originally Posted by cperciva

Are you sure? It looks to me like it should give an advantage of at least 1.5 bits.

Ignoring SSE2 which is standard 64-bit FP, the x86 code does all loads and stores in the 64-bit format. The reason for this is that the 80-bit load and store operations are very slow. Thus, the FFT code does not benefit greatly from the 80-bit internal format.

[cperciva]

Originally Posted by Prime95

Ignoring SSE2 which is standard 64-bit FP, the x86 code does all loads and stores in the 64-bit format. The reason for this is that the 80-bit load and store operations are very slow. Thus, the FFT code does not benefit greatly from the 80-bit internal format.

I don't know exactly how your FFT is implemented, but assuming that you calculate (a,b) <- (a+bw,a-bw) entirely in registers, the maximum error drops from B^2*N*(n*3.943... + O(1))*eps to B^2*N*(n*1.707... + O(1))*eps. I can't be so precise about the average error, but it should be similar.

The critical point is that complex multiplication and addition introduce up to sqrt(5)+1 ulp of error, while loading and storing introduce only 1 ulp. (Unavoidable inaccuracies in the trig data contributes the remaining sqrt(1/2) ulp.)

[cperciva]

(If anyone hasn't already read it, my paper on FFT error bounds has *finally* been published in Math. Comp.)

[Prime95]

Originally Posted by cperciva

I don't know exactly how your FFT is implemented, but assuming that you calculate (a,b) <- (a+bw,a-bw) entirely in registers, the maximum error drops from B^2*N*(n*3.943... + O(1))*eps to B^2*N*(n*1.707... + O(1))*eps. I can't be so precise about the average error, but it should be similar.

Obviously, I've not studied the error accumulation as much as you have. Prime95 uses a radix 4 implementation. Correct me if I'm wrong: When using 80-bit FP, you can pretty much assume all the error is introduced by the store instruction that occurs every 2 "levels" of the FFT. Thus, during a 1M FFT squaring each FFT data value is stored 21 times (10 in the forward FFT, 10 in the inverse FFT, and 1 in the carry propagation).
When using 64-bit FP, each radix-4 "unit" does a complex multiply (4 muls, 2 adds) and two complex adds. That is 2 muls and 3 adds for each data value. In the 1M FFT squaring that means there are about 105 opportunities for error to accumulate.

Are you saying this 5 to 1 ratio should result in an extra 1.5 bits of precision? Does that match the 1% to 2% larger exponents the x86
FFTs support?

[cperciva]

It's somewhat more complicated than that; the (not perfectly accurate) trig data contributes some error, and the error introduced in complex multiplication can be bounded somewhat more closely than a naive operation count would indicate.

But to answer your second question, in the 20M range you're packing about 20 bits per double, so reducing the error by one bit would mean a change of 2.5% in the size of numbers for which you can use a given FFT length.

That said, all my experience is with worst-case bounds; in practice the errors are much smaller; indeed, I could construct values for which the prime95 multiplication would fail -- but it doesn't matter since double-checking would find those anyway.

[ewmayer]

Originally Posted by cperciva

(If anyone hasn't already read it, my paper on FFT error bounds has *finally* been published in Math. Comp.)

Colin, could you kindly provide a link to your paper?

Now that you mention it, I see that Crandall, Papadopoulos' and my paper on F24 has also appeared, at long last - just over 3 years since we actually submitted the initial draft:

http://www.ams.org/journal-getitem?pii=S0025-5718-02-01479-5

The LaTeX manuscript can be gotten for free from

http://hogranch.com/mayer/F24.rev3.tex

Math. Comp. is very slow, and (at least based on our experience, where one reviewer didn't do anything with the manuscript for over a year, without MC doing anything to remind them of their responsibility as a reviewer) has poor timeliness-of-review controls. But grousing about tardiness aside, that paper also has some words about errors in floating-point FFTs, with an approach that is more heuristic in nature and more geared to most-likely error behavior than Colin's. The error estimates thus derived seem to match real-world experience quite closely.

[cperciva]

http://www.ams.org/journal-getitem?pii=S0025-5718-02-01419-9

Or, for those of you without a subscription, a preprint is at http://www.sfu.ca/~cperciva/mpaper.ps

I refer to your F24 paper in a section discussing previous error bound work -- Crandall sent me a copy a couple years ago.