Re: New insns for the s390 backend (2)

[Bradley Lucier <lucier at math dot purdue dot edu>]

On Wednesday, August 27, 2003, at 04:43  PM, Roger Sayle wrote:

> I did however manage to find the following URL of interest:
> http://www.active-web.cc/html/research/sine/sin-cos.txt
> All of the algorithms for precomputing FFT sin/cos tables given
> on that page appear to be unaffected by -ffast-math.

Roger, look at the accuracy requirements for sin/cos in Percival's 
paper as I applied it to my code---each entry can be off by four 
rounding errors, max, in tables with tens of millions of values.  Your 
reference is not as complete as

http://www.dattalo.com/technical/theory/sinewave.html

The error for method 4 here grows linearly; the one for method 5 grows 
quadratically.  If you want to use the best method I know of

sin(a + da) = cos(da) * sin(a) + sin(da) * cos(a)
            = sin(a) + ((cos(da)-1) * sin(a) + sin(da) * cos(a))
            = sin(a) + ((-2*sin^2(da/2)) * sin(a) + sin(da) * cos(a))

and similarly for cos, you can make a probabilistic argument that the 
error grows like the square root of the number of iterations. So, *if* 
you have 80 bit intermediate values you *probably* have < 4 rounding 
errors at double precision with tens of millions of table entries (and 
on PPC, SPARC, ..., you don't have that extra precision).  But that 
doesn't cut it for the number theory folks---they (I) need guaranteed < 
4 rounding error accuracy, so I just do library calls.  On x86 I could 
use the above iteration in blocks of 2048 table entries *if* I could be 
assured that gcc wouldn't spill any temporaries to doubles, but then 
I'd still have jumps in values at the boundaries of the blocks, and I 
just gave up with that.

> Then just to double check, I've just built and run "make check"
> on GNU gmp-4.1.2 as configured by default and then again with
> CC="gcc -ffast-math".  No failures/differences.  This is especially
> relevant to your concerns, as mpn/generic/mul_fft.c contains code to
> perform bignum multiplications using Karatsuba, 3-way Toom-Cook,
> and Fermat FFT.

I don't believe that gmp uses floating-point FFTs at all, see the code 
in mpn/generic/mul_fft.c, they use FFTs in other fields.  My impression 
of the gmp folks is that they don't trust floating-point arithmetic, 
numerical analysis (or analysts ;-), or mathematical libraries, which 
is fine, but I'm not surprised to find out that -ffast-math doesn't 
affect the test results.

I think this is getting way afield of the topic of gcc-patches; perhaps 
we should take any more discussion off list.

Brad