Skip site navigation (1)Skip section navigation (2)
Date:      Tue, 07 Aug 2012 21:58:39 -0500
From:      Stephen Montgomery-Smith <stephen@missouri.edu>
To:        freebsd-numerics@freebsd.org
Subject:   Re: cpow(3) implementations.
Message-ID:  <5021D5DF.7010709@missouri.edu>
In-Reply-To: <20120808000357.GA11128@server.rulingia.com>
References:  <20120808000357.GA11128@server.rulingia.com>

next in thread | previous in thread | raw e-mail | index | archive | help
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.



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5021D5DF.7010709>