| Sure. For each trig function with 1 argument, there is a corresponding
2-argument function; for, e.g., sin, this is defined by sin(x,cycle) =
sin(x*2pi/cycle). That is, the 2-arg version has an explicit cycle (always
greater than 0) which needn't be 2pi.
Now, given two *binary* floating point numbers x and cycle, it is fairly simple
to show that if x > 0, x mod cycle can be computed exactly (no roundoff error)
as a floating point number. So, if we calculate (x mod cycle)/cycle, we have a
very accurate result which can be fed to a routine which calculates sin(x,1) by
power series or some other method. And this routine will be accurate over the
entire range of floating point numbers, which is somewhat surprising: for high
exponent values, consecutive machine numbers are millions of cycles apart. (By
"accurate" I mean that the error is small relative to the result, i.e., the sin,
not merely small relative to the argument.)
The exercise I wanted to solve was to find a method which was similarly accurate
for the 1-argument case, where the cycle is 2pi exactly. This involves
calculating (x/2pi) mod 1 if I am reusing the routine that finds sin(x,1). It
should now be at least partly evident why the long bit string is useful: if I
want to find x*(1/2pi) mod 1, and x = 2^e*f where f is the fraction, then I can
find this by "shifting left" 1/2pi by about e-m bits, where m is the mantissa
length, and ignoring bits that "shift out" on the left (because they can only
contribute integral multiples of the cycle to the computation).
|
| It can be shown that x mod cycle is a machine number also in the nonbinary case,
but I haven't found an algorithm for computing the exact result that is
obviously right. Exercise: find a reasonable algorithm, obvious or not, that
calculates it.
|
| Here's one, just for the record. Stolen, I mean, cut-and-pasted from:
http://www.ast.univie.ac.at/~wasi/PI/stuart_lyster/bb.htm
"
More recently, Borwein & Borwein found another quartically converging
(and better - 4x better) algorithm.
This finds 1 / PI:
a(0) = 6 - (4 * square root (2))
y(0) = (square root (2)) - 1
y(n+1) = [1 - (fourth root( 1 - y^4))] / [1 + (fourth root( 1 - y^4))]
a(n+1) = a(n)*( 1 + y(n+1))^4 - 2^(2n+3) * y(n+1) * (1 + y(n+1) +
y(n+1)^2)
A(n) converges to 1/PI.
"
Have fun, and don't forget to post the code.
John
|
| Re .2:
> ... The exercise I wanted to solve was to find a method which was
> similarly accurate for the 1-argument case, where the cycle is 2pi
> exactly. This involves calculating (x/2pi) mod 1 if I am reusing
> the routine that finds sin(x,1).
Sometime around 1982 Mary Payne and I published an article in SIGNUM on how
to do exactly what you are attempting. I believe you can find copies of the
article by searching the Corporate Research groups archieves of published
reports for
"Argument Reduction for Trigonometric Functions"
The technique you are trying to develop is used in all of the the VAX and Alpha
trig functions. The actual number of bits of 1/pi required is roughly
equal to the size of the largest argument you intend to reduce. For VAX and
Alpha libraries, the largest value is on the order of 2^16385, so we keep 1/pi
store to roughly 17000 bits.
In order to generate the 1/pi bits, we (the math library group) use a home
grown multi-precision package. I don't remember exactly how long it takes to
generate the points anymore, but I'm pretty sure it's under an hour.
Several other vendors, including SUN and IBM use the technique described in the
SIGNUM paper for their trig functions.
|