From owner-freebsd-numerics@FreeBSD.ORG  Fri Oct  5 21:48:54 2012
Return-Path: <owner-freebsd-numerics@FreeBSD.ORG>
Delivered-To: freebsd-numerics@freebsd.org
Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52])
	by hub.freebsd.org (Postfix) with ESMTP id BF89A106566B
	for <freebsd-numerics@freebsd.org>;
	Fri,  5 Oct 2012 21:48:54 +0000 (UTC)
	(envelope-from stephen@missouri.edu)
Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu
	[128.206.184.213])
	by mx1.freebsd.org (Postfix) with ESMTP id 894998FC0C
	for <freebsd-numerics@freebsd.org>;
	Fri,  5 Oct 2012 21:48:54 +0000 (UTC)
Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213])
	by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id
	q95LmkMj027265
	for <freebsd-numerics@freebsd.org>; Fri, 5 Oct 2012 16:48:47 -0500 (CDT)
	(envelope-from stephen@missouri.edu)
Message-ID: <506F55BE.7000304@missouri.edu>
Date: Fri, 05 Oct 2012 16:48:46 -0500
From: Stephen Montgomery-Smith <stephen@missouri.edu>
User-Agent: Mozilla/5.0 (X11; Linux i686;
	rv:15.0) Gecko/20120912 Thunderbird/15.0.1
MIME-Version: 1.0
To: freebsd-numerics@freebsd.org
References: <20121005213217.GA90440@vps.rulingia.com>
In-Reply-To: <20121005213217.GA90440@vps.rulingia.com>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
Subject: Re: Implementing cpow(3)
X-BeenThere: freebsd-numerics@freebsd.org
X-Mailman-Version: 2.1.5
Precedence: list
List-Id: "Discussions of high quality implementation of libm functions."
	<freebsd-numerics.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-numerics>, 
	<mailto:freebsd-numerics-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-numerics>
List-Post: <mailto:freebsd-numerics@freebsd.org>
List-Help: <mailto:freebsd-numerics-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-numerics>, 
	<mailto:freebsd-numerics-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Fri, 05 Oct 2012 21:48:54 -0000

On 10/05/2012 04:32 PM, Peter Jeremy wrote:

> Looking at "rth = rw * th + iw * lrad;": This is used solely as an
> argument for sin() or cos() - both of which have a period of 2π,
> therefore rth can be evaluated modulo 2π, as can each term in the
> expression.  Assuming lrad is split as (lrad0 + lrad1) and mod_2pi()
> performs a high-precision modulo 2π:
>    rth = mod_2pi(rw) * th + mod_2pi(iw) * (mod_2pi(lrad0) + mod_2pi(lrad1));
> This should minimise the catastrophic cancellation.

I just wanted to comment that you seem to be using that

mod_2pi(a*b) = mod_2pi(a) * mod_2pi(b) + 2*PI*n,

for some integer n.  But this is generally only true if a and b are 
integers.

If I were writing cpow(complex z, complex w), my code would look at the 
special case when w is a real integer, and handle it separately from the 
other cases.

When z and w are general complex numbers (or even when z is complex and 
w is merely restricted to be real), I don't see how to avoid 
catastrophic cancellation.