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>