From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 8 02:58:41 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 7C675106566C for ; Wed, 8 Aug 2012 02:58:41 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id 2A6288FC0A for ; Wed, 8 Aug 2012 02:58:41 +0000 (UTC) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q782wdOY030787 for ; Tue, 7 Aug 2012 21:58:40 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <5021D5DF.7010709@missouri.edu> Date: Tue, 07 Aug 2012 21:58:39 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: freebsd-numerics@freebsd.org References: <20120808000357.GA11128@server.rulingia.com> In-Reply-To: <20120808000357.GA11128@server.rulingia.com> Content-Type: text/plain; charset=UTF-8; format=flowed Content-Transfer-Encoding: 7bit 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 02:58:41 -0000 On 08/07/2012 07:03 PM, 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) I wouldn't regard errors in a*u-b*v as catastrophic cancellation. This is because exp(0) = 1. So if the error in computing a*u-b*v is approx DBL_EPSILON, and a*u-b*v approx zero, even though the relative error in computing a*u-b*v is going to large (perhaps even infinite), nevertheless the error in exp(a*u-b*v) may still be bounded by 1 or 2 ulp. More generally, as a mathematician, I would be far more concerned that cpow(z,w) return accurate answers when w is real, and especially when w is a small integer. Real life applications of cpow(z,y) when w is not real are very few and far between. I would be pleased if cpow(x,y) made special provisions for when y is a small integer, and so, for example, cpow(z,2) was computed as z*z = (x+y)*(x-y) + 2x*y*I. For cpow(z,3), you are going to have a hard time avoiding large relative errors when x^2 = 3y^2 (i.e. z is parallel to a cube root of 1). Frankly I see that as somewhat unavoidable. Nevertheless, if you pumped up cpow(z,w) so that when w was a relatively small integer, that it broke w into its base 2 expansion, and then multiplied lots of terms of the form cpow(z,2^n), each of which was computed by n repetitions of cpow(z,2), and then for negative integers by taking the reciprocal of the whole thing, and then in all other cases simply use cexp(w*clog(z)), I would be very happy. Other than that, if your cpow produced cpow(-1,0.5) = I + 1e-16, I wouldn't be shocked at all, and I would find this kind of error totally acceptable.