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