Getting the last bit right is a challenge, see the end of http://www.research.ibm.com/journal/rd/341/ibmrd3401N.pdf
Henry Rich > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of Roger Hui > Sent: Wednesday, December 06, 2006 4:36 PM > To: General forum > Subject: [Jgeneral] exp(y) > > I don't understand something here: > > y=: 100 > n=: <.0.5*y% ^.2 > (2^n)*^y-n*^.2 > 2.68812e43 > ^ y > 2.68812e43 > y-n*^.2 > 50.0934 > > But 50.09 requires more than "relatively few terms" > of the Taylor series, and (at least in this case) > y-n*^.2 does not involve subtracting nearly equal > quantities. > > If I were implementing from scratch ^y where y is a > 64-bit IEEE floating point number, I would do it as > follows: > > a. Pre-compute the table T=: ^i.710 . > > b. If 0>y, ^y is %^-y . > > c. If y>:710, ^y is _ > > d. For y in the range from 0 to 710, ^y is > ((<.y){T) * PS y-<.y > where PS is the power series for ^ . Since y-<.y > is less than 1 and %!18 is 1.56192e_16, at most 19 > terms are required for 16 digits of precision. > > > > ----- Original Message ----- > From: John Randall <[EMAIL PROTECTED]> > Date: Tuesday, December 5, 2006 4:20 pm > Subject: Re: [Jgeneral] Bug? Sin of pi. > > > The exponentional function can be calculated accurately for large > > arguments using relatively few terms of the Taylor series > by using the > > identity > > > > ^y-:(2^n)*^(y-n*^.2) > > > > where n=<.0.5 y% ^.2 . > > > > This is the method used by Microsoft (and other) libm libraries. It > > has the advantage that multiplication by 2^n can be done by adding n > > to the exponent. > > > > Unfortunately, (y-n*^.2) involves subtracting nearly equal > quantities, > > which leads to loss of precision. In particular, ^.2 has to be > > calculated to greater precision than the result. This is > commonly not > > done carefully, and results are not reliable. I leave it to the > > reader to show that ^100 in J (and most languages) does not > have full > > precision. > > > ---------------------------------------------------------------------- > For information about J forums see > http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
