cpow(3) implementations.

Bruce Evans brde at optusnet.com.au
Thu Aug 9 21:34:00 UTC 2012


On Wed, 8 Aug 2012, Stephen Montgomery-Smith wrote:

> On 08/08/2012 02:49 PM, Stephen Montgomery-Smith wrote:
>> On 08/08/2012 01:11 PM, Bruce Evans wrote:
>>> On Tue, 7 Aug 2012, Stephen Montgomery-Smith wrote:
> ...
>>>> 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.
>
> I tested it.  It performs very poorly.

x^2 - y^2 in z^2 may have large cancelation, and then squaring again
may have further large cancelation, so the problem doesn't really change

> I wrote a cpow function, just for integers.  cpow(x+I*y,n) does very poorly 
> when x/y approx tan(Pi/(2n)), that is, x+I*y is parallel to an nth root of 
> -1.  Only the cases n=0,1,2,-1,-2 do well.
>
> As a mathematician, I wouldn't expect any better.  In fact I would be 
> pleasantly surprised if it got the case n=2 correct.

I think that is "as a numerical analyst" :-).  The magic of analytic
functions involves delicate cancelations.

Bruce


More information about the freebsd-numerics mailing list