From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 8 09:24:58 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 14690106564A for ; Wed, 8 Aug 2012 09:24:58 +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 9C5878FC08 for ; Wed, 8 Aug 2012 09:24:57 +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 q789Ot27050752 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Wed, 8 Aug 2012 19:24:56 +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 q789OkvT066394 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Wed, 8 Aug 2012 19:24:46 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q789OjF4066393 for freebsd-numerics@freebsd.org; Wed, 8 Aug 2012 19:24:45 +1000 (EST) (envelope-from peter) Date: Wed, 8 Aug 2012 19:24:45 +1000 From: Peter Jeremy To: freebsd-numerics@freebsd.org Message-ID: <20120808092445.GB65056@server.rulingia.com> References: <20120808000357.GA11128@server.rulingia.com> <5021D5DF.7010709@missouri.edu> MIME-Version: 1.0 Content-Type: multipart/signed; micalg=pgp-sha1; protocol="application/pgp-signature"; boundary="5vNYLRcllDrimb99" Content-Disposition: inline In-Reply-To: <5021D5DF.7010709@missouri.edu> X-PGP-Key: http://www.rulingia.com/keys/peter.pgp User-Agent: Mutt/1.5.21 (2010-09-15) Subject: Re: 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 09:24:58 -0000 --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 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 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--