Date: Wed, 8 Aug 2012 19:24:45 +1000 From: Peter Jeremy <peter@rulingia.com> To: freebsd-numerics@freebsd.org Subject: Re: cpow(3) implementations. Message-ID: <20120808092445.GB65056@server.rulingia.com> In-Reply-To: <5021D5DF.7010709@missouri.edu> References: <20120808000357.GA11128@server.rulingia.com> <5021D5DF.7010709@missouri.edu>
index | next in thread | previous in thread | raw e-mail
[-- Attachment #1 --]
On 2012-Aug-07 18:34:59 -0700, Steve Kargl <sgk@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@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
[-- Attachment #2 --]
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.19 (FreeBSD)
iEYEARECAAYFAlAiMF0ACgkQ/opHv/APuIebdACfeAT7FOE2T0hKyq4Pe9cAkvuc
xUsAn3rJsudDpsacUnPZbDIAN6wT6/IX
=J5iy
-----END PGP SIGNATURE-----
help
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120808092445.GB65056>