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
