Implementing cpow(3)
Stephen Montgomery-Smith
stephen at missouri.edu
Fri Oct 5 21:48:54 UTC 2012
On 10/05/2012 04:32 PM, Peter Jeremy wrote:
> Looking at "rth = rw * th + iw * lrad;": This is used solely as an
> argument for sin() or cos() - both of which have a period of 2π,
> therefore rth can be evaluated modulo 2π, as can each term in the
> expression. Assuming lrad is split as (lrad0 + lrad1) and mod_2pi()
> performs a high-precision modulo 2π:
> rth = mod_2pi(rw) * th + mod_2pi(iw) * (mod_2pi(lrad0) + mod_2pi(lrad1));
> This should minimise the catastrophic cancellation.
I just wanted to comment that you seem to be using that
mod_2pi(a*b) = mod_2pi(a) * mod_2pi(b) + 2*PI*n,
for some integer n. But this is generally only true if a and b are
integers.
If I were writing cpow(complex z, complex w), my code would look at the
special case when w is a real integer, and handle it separately from the
other cases.
When z and w are general complex numbers (or even when z is complex and
w is merely restricted to be real), I don't see how to avoid
catastrophic cancellation.
More information about the freebsd-numerics
mailing list