Date: Wed, 8 Aug 2012 10:03:57 +1000 From: Peter Jeremy <peter@rulingia.com> To: freebsd-numerics@freebsd.org Subject: cpow(3) implementations. Message-ID: <20120808000357.GA11128@server.rulingia.com>
next in thread | raw e-mail | index | archive | help
--0OAP2g/MAC+5xKAE
Content-Type: text/plain; charset=utf-8
Content-Disposition: inline
Content-Transfer-Encoding: quoted-printable
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) =3D 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 =3D a + I*b
z =3D c + I*d
cis(r) =3D cos(r) + I*sin(r)
t =3D u + I*v =3D clog(c + I*d)
=3D log(hypot(c, d)) + I*atan2(d, c)
cpow(z, w) =3D cexp(w * clog(z))
=3D cpow(c + I*d, a + I*b)
=3D cexp((a + I*b) * clog(c + I*d))
=3D cexp((a + I*b) * (u + I*v))
=3D cexp((a*u - b*v) + I*(a*v + b*u))
=3D 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) =3D exp(a*u) / exp(b*v)
=3D 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 =E2=89=A4 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?
--=20
Peter Jeremy
--0OAP2g/MAC+5xKAE
Content-Type: application/pgp-signature
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.19 (FreeBSD)
iEYEARECAAYFAlAhrO0ACgkQ/opHv/APuIc18gCdEa+5IIDC+yucx8uyX/69ZjUk
CsgAn25wsZLxw32By+J2sRTpbOguZoEj
=+T1K
-----END PGP SIGNATURE-----
--0OAP2g/MAC+5xKAE--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120808000357.GA11128>