Skip site navigation (1)Skip section navigation (2)
Date:      Sat, 15 Sep 2012 13:21:32 -0500
From:      Stephen Montgomery-Smith <stephen@missouri.edu>
To:        Bruce Evans <brde@optusnet.com.au>
Cc:        freebsd-numerics@freebsd.org, Steve Kargl <sgk@troutmask.apl.washington.edu>
Subject:   Re: cexp error
Message-ID:  <5054C72C.1070903@missouri.edu>
In-Reply-To: <20120915210927.Q2023@besplex.bde.org>
References:  <5048D00B.8010401@missouri.edu> <504D3CCD.2050006@missouri.edu> <504FF726.9060001@missouri.edu> <20120912191556.F1078@besplex.bde.org> <20120912225847.J1771@besplex.bde.org> <50511B40.3070009@missouri.edu> <20120913204808.T1964@besplex.bde.org> <5052A923.9030806@missouri.edu> <20120914143553.X870@besplex.bde.org> <5053CB4F.3030709@missouri.edu> <20120915003408.GA70269@troutmask.apl.washington.edu> <20120915210927.Q2023@besplex.bde.org>

next in thread | previous in thread | raw e-mail | index | archive | help
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.



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5054C72C.1070903>