cpow(3) implementations.

[Peter Jeremy]

On 2012-Aug-07 18:34:59 -0700, Steve Kargl <sgk at troutmask.apl.washington.edu> wrote:
>On Wed, Aug 08, 2012 at 10:03:57AM +1000, Peter Jeremy wrote:
>>   cpow(z, w) = exp(a*u - b*v) * cis(a*v + b*u)
...
>>   exp(a*u - b*v) = exp(a*u) / exp(b*v)
>>                  = pow(a, hypot(c, d)) / exp(b*v)

Oops, that last line is wrong and should be:
   exp(a*u - b*v) = pow(hypot(c, d), a) / exp(b*v)

>I cheated and peaked at NetBSD's implementation, which is
>based on Moshier's cephes implementation.

The pseudocode didn't make sense to me (I couldn't see the 2nd
argument anywhere) so I had a peek at both NetBSD and cephes-2.8
and the code is identical in both.  I believe the pseudo code
should be:

cpow(z,w):
x = creal(w)
y = cimag(w)
r = cabs(z)
t = carg(z)

if y == 0
   if x == 0
      return 0
   else
      return pow(r,x)*(cos(t*x)+I*sin(t*x))
else
   return pow(r,x)*exp(-t*y)*(cos(t*x+y*log(r))+I*sin(t*x+y*log(r)))

Note that the final line looks the same as my final line (modulo
algebraic errors).  The rest is just special casing.

On 2012-Aug-07 21:58:39 -0500, Stephen Montgomery-Smith <stephen at missouri.edu> wrote:
>I wouldn't regard errors in a*u-b*v as catastrophic cancellation.  This 
>is because exp(0) = 1.  So if the error in computing a*u-b*v is approx 
>DBL_EPSILON, and a*u-b*v approx zero, even though the relative error in 
>computing a*u-b*v is going to large (perhaps even infinite), 
>nevertheless the error in exp(a*u-b*v) may still be bounded by 1 or 2 ulp.

OK.

>I would be pleased if cpow(x,y) made special provisions for when y is a 
>small integer, and so, for example, cpow(z,2) was computed as z*z = 
>(x+y)*(x-y) + 2x*y*I.

Thu square-and-multiply approach makes sense.  Where do you draw the
line on "small"?  At some point the rounding errors will exceed the
trig approach.

>Other than that, if your cpow produced cpow(-1,0.5) = I + 1e-16, I 
>wouldn't be shocked at all, and I would find this kind of error totally 
>acceptable.

I like the idea of also special-casing "simple" rational powers.  This
offers things like cpow(-1.,0.5) = I and cpow(-8.,1./3.) = -2

Thank you both for your input.

-- 
Peter Jeremy
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 196 bytes
Desc: not available
Url : http://lists.freebsd.org/pipermail/freebsd-numerics/attachments/20120808/02b7571d/attachment.pgp