From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 8 00:04:13 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 54164106566C for ; Wed, 8 Aug 2012 00:04:13 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id DA0088FC0A for ; Wed, 8 Aug 2012 00:04:12 +0000 (UTC) Received: from server.rulingia.com (c220-239-247-45.belrs5.nsw.optusnet.com.au [220.239.247.45]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q78045DW049406 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Wed, 8 Aug 2012 10:04:06 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q7803v7n034688 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Wed, 8 Aug 2012 10:03:57 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7803vGk034686 for freebsd-numerics@freebsd.org; Wed, 8 Aug 2012 10:03:57 +1000 (EST) (envelope-from peter) Date: Wed, 8 Aug 2012 10:03:57 +1000 From: Peter Jeremy To: freebsd-numerics@freebsd.org Message-ID: <20120808000357.GA11128@server.rulingia.com> MIME-Version: 1.0 Content-Type: multipart/signed; micalg=pgp-sha1; protocol="application/pgp-signature"; boundary="0OAP2g/MAC+5xKAE" Content-Disposition: inline X-PGP-Key: http://www.rulingia.com/keys/peter.pgp User-Agent: Mutt/1.5.21 (2010-09-15) Subject: cpow(3) implementations. X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 08 Aug 2012 00:04:13 -0000 --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--