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>