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