Date: Fri, 10 Aug 2012 07:33:55 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-numerics@freebsd.org Subject: Re: cpow(3) implementations. Message-ID: <20120810072650.E4653@besplex.bde.org> In-Reply-To: <50232C2A.8040704@missouri.edu> References: <20120808000357.GA11128@server.rulingia.com> <5021D5DF.7010709@missouri.edu> <20120809033123.V4514@besplex.bde.org> <5022C2B2.2030803@missouri.edu> <50232C2A.8040704@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120810072650.E4653>