Date: Sat, 15 Sep 2012 13:59:45 -0500 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: freebsd-numerics@freebsd.org Subject: Re: cexp error Message-ID: <5054D021.1070508@missouri.edu> In-Reply-To: <5054C72C.1070903@missouri.edu> 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> <5054C72C.1070903@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On 09/15/2012 01:21 PM, Stephen Montgomery-Smith wrote: > 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. I definitely screwed up the details in my presentation. For one thing, I forgot it is |g pi - f/2^n| where f and g are integers, that we want to see how small it is. But I suspect I got the rough idea correct.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5054D021.1070508>