Skip site navigation (1)Skip section navigation (2)
Date:      Tue, 7 Aug 2012 18:34:59 -0700
From:      Steve Kargl <sgk@troutmask.apl.washington.edu>
To:        Peter Jeremy <peter@rulingia.com>
Cc:        freebsd-numerics@freebsd.org
Subject:   Re: cpow(3) implementations.
Message-ID:  <20120808013459.GA19280@troutmask.apl.washington.edu>
In-Reply-To: <20120808000357.GA11128@server.rulingia.com>
References:  <20120808000357.GA11128@server.rulingia.com>

next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, Aug 08, 2012 at 10:03:57AM +1000, Peter Jeremy wrote:
> 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?

I cheated and peaked at NetBSD's implementation, which is
based on Moshier's cephes implementation.  Your approach
not only has to worry about catastrophic cancellation, but
may need to switch between exp() and expm1() and log() and
log1p() to achieve accuracy near the problematic areas for
these functions.  Anyway, Moshier converts z in polar coordinates
and use pow(), exp(), sin(), and cos().

psuedo-code:

x = creal(z)
y = cimag(z)
r = hypot(x,y)
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)+I*sin(t*x+y))

I don't if this is any better.  

-- 
Steve



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120808013459.GA19280>