From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 8 01:35:54 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 D17EE106564A for ; Wed, 8 Aug 2012 01:35:54 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) by mx1.freebsd.org (Postfix) with ESMTP id A7B658FC0A for ; Wed, 8 Aug 2012 01:35:54 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.5/8.14.5) with ESMTP id q781Z0Pf019395; Tue, 7 Aug 2012 18:35:00 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.5/8.14.5/Submit) id q781Yx77019394; Tue, 7 Aug 2012 18:34:59 -0700 (PDT) (envelope-from sgk) Date: Tue, 7 Aug 2012 18:34:59 -0700 From: Steve Kargl To: Peter Jeremy Message-ID: <20120808013459.GA19280@troutmask.apl.washington.edu> References: <20120808000357.GA11128@server.rulingia.com> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20120808000357.GA11128@server.rulingia.com> User-Agent: Mutt/1.4.2.3i Cc: freebsd-numerics@freebsd.org 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 01:35:54 -0000 On Wed, Aug 08, 2012 at 10:03:57AM +1000, Peter Jeremy wrote: > 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) = 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 = a + I*b > z = c + I*d > cis(r) = cos(r) + I*sin(r) > t = u + I*v = clog(c + I*d) > = log(hypot(c, d)) + I*atan2(d, c) > > cpow(z, w) = cexp(w * clog(z)) > = cpow(c + I*d, a + I*b) > = cexp((a + I*b) * clog(c + I*d)) > = cexp((a + I*b) * (u + I*v)) > = cexp((a*u - b*v) + I*(a*v + b*u)) > = 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) = exp(a*u) / exp(b*v) > = 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 ??? 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? I cheated and peaked at NetBSD's implementation, which is based on Moshier's cephes implementation. Your approach not only has to worry about catastrophic cancellation, but may need to switch between exp() and expm1() and log() and log1p() to achieve accuracy near the problematic areas for these functions. Anyway, Moshier converts z in polar coordinates and use pow(), exp(), sin(), and cos(). psuedo-code: x = creal(z) y = cimag(z) r = hypot(x,y) 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)+I*sin(t*x+y)) I don't if this is any better. -- Steve