From owner-freebsd-numerics@FreeBSD.ORG Thu Aug 9 21:34:00 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 6C274106564A for ; Thu, 9 Aug 2012 21:34:00 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail03.syd.optusnet.com.au (mail03.syd.optusnet.com.au [211.29.132.184]) by mx1.freebsd.org (Postfix) with ESMTP id DE2C08FC08 for ; Thu, 9 Aug 2012 21:33:59 +0000 (UTC) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail03.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q79LXuvx028045 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Fri, 10 Aug 2012 07:33:57 +1000 Date: Fri, 10 Aug 2012 07:33:55 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <50232C2A.8040704@missouri.edu> Message-ID: <20120810072650.E4653@besplex.bde.org> References: <20120808000357.GA11128@server.rulingia.com> <5021D5DF.7010709@missouri.edu> <20120809033123.V4514@besplex.bde.org> <5022C2B2.2030803@missouri.edu> <50232C2A.8040704@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-numerics@freebsd.org Subject: Re: cpow(3) implementations. X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 09 Aug 2012 21:34:00 -0000 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