cexp error

[Stephen Montgomery-Smith]

On 09/15/2012 07:00 AM, Bruce Evans wrote:
> On Fri, 14 Sep 2012, Steve Kargl wrote:
>
>> On Fri, Sep 14, 2012 at 07:26:55PM -0500, Stephen Montgomery-Smith wrote:
>>> On 09/14/2012 01:05 AM, Bruce Evans wrote:
>>>
>>>> My tests also determined the exact minimum for all multiples up to the
>>>> thresholds for "large" multiples in e_rem_pio2.c, except for ld128,
>>>> to verify that there are enough bits in the special approximations for
>>>> Pi/2 there, except for ld128.  The maximum multiples handled there are
>>>> 2**28*Pi/2, except for ld128 they are 2**45*Pi/2.  2**45 is too many
>>>> to check exhaustively.
>>>
>>> But presumably one could prove a result to this effect using continued
>>> fractions, right?
>
> Yes.  I'm not sure if the theory can produce a good bound for a relatively
> small range of x though.
>
>> See attached.  It shows a method for determining the number
>> of needed bits.
>
> It doesn't give the details for continued fractions, but supplies the bound
> for 53 mantissa bits (2**-62 -> 61 extra bits), and shows that this was
> known in 1992, and describes the 1992 fdlibm better/differently than I
> described the current implementation.  The 73 (75?) extra that I remembered
> must be for 64 mantissa bits.
>
> I have Kahan's program which does searches using the continued fraction
> method (das sent it to me in 2008), but don't really understand it, and
> in particular don't know how to hack it to give the bound for the 2**45
> range.  Saved results show that I didn't finish checking with it in 2008
> either.

I am just making a guess here.  But I am thinking that Kahan's program 
might be based around the following type of argument.  First, find 
integers p, q and r so that
|pi - p/q| < 1/r
The idea is that r should be much bigger than q.  There is a theorem 
that says you can find a p, q and r so that r=q^2.  In any case, I 
believe continued fractions is the way to find these kinds of numbers.

Now suppose that n = the number of bits in the mantissa.  Since non-zero 
multiples of pi/2 are bigger than 1, we are asking the question: how 
small can
|pi - f/2^n|
be for any integer f.  How use the inequalities:

|pi - f/2^n|
 >= |p/q - f/2^n| - 1/r
= |p*2^n - f*q|/(q*2^n) - 1/r.

If you have picked r bigger than q*2^n, which we know is possible, then 
it is just a case of seeing if |p*2^n - f*q| can equal zero.  Since we 
may assume that p and q have no common factors, this can only happen if 
f=p*2^m and q=2^(n-m).

Probably the argument used by Kahan or das@ is sharper than this.