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>
next in thread | previous in thread | raw e-mail | index | archive | help
--5vNYLRcllDrimb99
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
Content-Transfer-Encoding: quoted-printable
On 2012-Aug-07 18:34:59 -0700, Steve Kargl <sgk@troutmask.apl.washington.ed=
u> wrote:
>On Wed, Aug 08, 2012 at 10:03:57AM +1000, Peter Jeremy wrote:
>> cpow(z, w) =3D exp(a*u - b*v) * cis(a*v + b*u)
=2E..
>> exp(a*u - b*v) =3D exp(a*u) / exp(b*v)
>> =3D pow(a, hypot(c, d)) / exp(b*v)
Oops, that last line is wrong and should be:
exp(a*u - b*v) =3D 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 =3D creal(w)
y =3D cimag(w)
r =3D cabs(z)
t =3D carg(z)
if y =3D=3D 0
if x =3D=3D 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.e=
du> wrote:
>I wouldn't regard errors in a*u-b*v as catastrophic cancellation. This=20
>is because exp(0) =3D 1. So if the error in computing a*u-b*v is approx=
=20
>DBL_EPSILON, and a*u-b*v approx zero, even though the relative error in=20
>computing a*u-b*v is going to large (perhaps even infinite),=20
>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=20
>small integer, and so, for example, cpow(z,2) was computed as z*z =3D=20
>(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) =3D I + 1e-16, I=20
>wouldn't be shocked at all, and I would find this kind of error totally=20
>acceptable.
I like the idea of also special-casing "simple" rational powers. This
offers things like cpow(-1.,0.5) =3D I and cpow(-8.,1./3.) =3D -2
Thank you both for your input.
--=20
Peter Jeremy
--5vNYLRcllDrimb99
Content-Type: application/pgp-signature
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.19 (FreeBSD)
iEYEARECAAYFAlAiMF0ACgkQ/opHv/APuIebdACfeAT7FOE2T0hKyq4Pe9cAkvuc
xUsAn3rJsudDpsacUnPZbDIAN6wT6/IX
=J5iy
-----END PGP SIGNATURE-----
--5vNYLRcllDrimb99--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120808092445.GB65056>