cpow(3) implementations.

[Peter Jeremy]

The C99 standard (at least WG14/N1256) is quite liberal as regards
cpow() and specifically allows a conforming implementation to just do:
  cpow(z, w) = cexp(w * clog(z))

The downside of this approach is that log() inherently loses precision
by pushing (most of) the exponent bits into the fraction, displacing
original fraction bits.  I've therefore been looking at how to
implement cpow(3) with a precision similar to pow(3).  The following
are some thoughts, together with questions.

In the following:
  w = a + I*b
  z = c + I*d
  cis(r) = cos(r) + I*sin(r)
  t = u + I*v = clog(c + I*d)
              = log(hypot(c, d)) + I*atan2(d, c)

  cpow(z, w) = cexp(w * clog(z))
             = cpow(c + I*d, a + I*b)
             = cexp((a + I*b) * clog(c + I*d))
             = cexp((a + I*b) * (u + I*v))
             = cexp((a*u - b*v) + I*(a*v + b*u))
             = exp(a*u - b*v) * cis(a*v + b*u)

Unfortunately, either or both of (a*u - b*v) and (a*v + b*u) are
potentially subject to catastrophic cancellation.  However:
  exp(a*u - b*v) = exp(a*u) / exp(b*v)
                 = pow(a, hypot(c, d)) / exp(b*v)
Ignoring overflow/underflow issues, does this avoid the cancellation
issues associated with the subtraction?  (I realise there are still
issues when a ≤ 0 due to the pow(3) definition).

Since cis() is periodic, it's also possible to range-reduce each term
of (a*v + b*u) independently - but that just fiddles with the issue
since it just adds another layer of precision loss.  Can anyone
suggest an alternative approach?

-- 
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/5accea0a/attachment.pgp