[PATCH] Evaluate pow(x,n) at compile-time

[Geoff Keating]

Roger Sayle <roger@www.eyesopen.com> writes:

> Hi Geoff,
> > The glibc pow() function is specified as "Given two IEEE double
> > machine numbers y,x it computes the correctly rounded (to nearest)
> > value of x^y." by which it means that it always gives the correct
> > result.  By comparison, your routine may under some circumstances give
> > results that differ from the correct result by 1ulp; it's easy to generate
> > an example for 128-bit floats, harder for 64-bit ones.
> >
> > I'm now worried about the sqrt() routine.  It's an IEEE requirement
> > that sqrt() always give the correct value.  Can you be sure that
> > there's no case in which your sqrt will produce an incorrect value?
> 
> 
> You should probably take a look at Gaston Gonnet's WWW site on "Error
> Analysis of System Mathematical Functions" at ETH, Zurich:
> http://www.inf.ethz.ch/personal/gonnet/FPAccuracy/Analysis.html
> 
> Not only does he show 0.56164 ulps error in glibc's pow implementation
> on i686-pc-linux-gnu, but also a 0.50024 ulp error in sqrt.  On exactly
> the same benchark cases, GCC's real_sqrt gets a perfect 0.5 ulp score!
> 
> http://gcc.gnu.org/ml/gcc-patches/2002-12/msg00471.html
> 
> 
> [Still don't trust me? :>]

That's on x86, where glibc uses the hardware, and suffers due to
various design/use problems with it.  The functions I referred to are
used on sane machines like sparc and powerpc.

While the paper is interesting, you can't prove correctness by running
a few thousand testcases.  If you wanted to prove correctness by
testing, you have to test *every* input value, all 2^64 of them, or
offer some argument that you only need to check a subset.

-- 
- Geoffrey Keating <geoffk@geoffk.org>