Skip site navigation (1)Skip section navigation (2)
Date:      Wed, 08 Aug 2012 14:49:06 -0500
From:      Stephen Montgomery-Smith <stephen@missouri.edu>
To:        Bruce Evans <brde@optusnet.com.au>
Cc:        freebsd-numerics@freebsd.org
Subject:   Re: cpow(3) implementations.
Message-ID:  <5022C2B2.2030803@missouri.edu>
In-Reply-To: <20120809033123.V4514@besplex.bde.org>
References:  <20120808000357.GA11128@server.rulingia.com> <5021D5DF.7010709@missouri.edu> <20120809033123.V4514@besplex.bde.org>

next in thread | previous in thread | raw e-mail | index | archive | help
On 08/08/2012 01:11 PM, Bruce Evans wrote:
> On Tue, 7 Aug 2012, Stephen Montgomery-Smith wrote:
>
>> 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.
>
> Indeed.  It's sort of the opposite of catastrophic cancelation.  We
> might have a*u-b*v off by a factor of 2, with value DBL_EPSILON when
> it should be 2*DBL_EPSILON.  That is about 2**52 ulps in double
> precision.  But then adding 1 turns it into an error of just 1 ulp.
>
> Even more so for the cos() part of cis().  Now the small term with the
> large error gets squared before adding 1.
>
> But for the sin() part of cis(), the final error is the one that the
> cancelation gave.  It may be small absolute for cis(), but on
> multiplication by large exp() it will become large absolute (though
> no larger relative, and always small relative to the cos() part which
> is about 1).
>
>> 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.
>
> This is easy when z is also real, by using pow().  Then there is the
> issue of continuity -- how close is cpow(x + I*eps, y) to pow(x, y)?
> Similarly for cpow(x, y + I*eps).
>
>> 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.
>
> x^2 - y^2 is easier than x^2 + y^2 - 1 for clog().
>
>> 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.
>
> x^2 - 3y^2 is still easier than for clog().  But go to z^4, and the real
> part is x^4 - 6 x^2 y^2 + y^4, so it's much harder than clog().  Seems
> to need quadrupled precision for the sum of the 4th powers, then a trick
> to avoid maybe octupled precision :-).

How about doing the real part of z^4 as
((x-y)*(x+y)-2*x*y)*((x-y)*(x+y)+2*x*y) ?
That is what my suggested process would do.




Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5022C2B2.2030803>