From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 15 15:20:00 2012 Return-Path: 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 1CD791065675 for ; Wed, 15 Aug 2012 15:19:59 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail04.syd.optusnet.com.au (mail04.syd.optusnet.com.au [211.29.132.185]) by mx1.freebsd.org (Postfix) with ESMTP id 3DEBD8FC18 for ; Wed, 15 Aug 2012 15:19:58 +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 mail04.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q7FFJoS5021129 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Thu, 16 Aug 2012 01:19:51 +1000 Date: Thu, 16 Aug 2012 01:19:50 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <502A820C.6060804@missouri.edu> Message-ID: <20120816010912.Q1751@besplex.bde.org> References: <5017111E.6060003@missouri.edu> <501C361D.4010807@missouri.edu> <20120804165555.X1231@besplex.bde.org> <501D51D7.1020101@missouri.edu> <20120805030609.R3101@besplex.bde.org> <501D9C36.2040207@missouri.edu> <20120805175106.X3574@besplex.bde.org> <501EC015.3000808@missouri.edu> <20120805191954.GA50379@troutmask.apl.washington.edu> <20120807205725.GA10572@server.rulingia.com> <20120809025220.N4114@besplex.bde.org> <5027F07E.9060409@missouri.edu> <20120814003614.H3692@besplex.bde.org> <50295887.2010608@missouri.edu> <20120814055931.Q4897@besplex.bde.org> <50297468.20902@missouri.edu> <20120814173931.V934@besplex.bde.org> <502A820C.6060804@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-numerics@freebsd.org Subject: Re: Complex arg-trig functions 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: Wed, 15 Aug 2012 15:20:00 -0000 On Tue, 14 Aug 2012, Stephen Montgomery-Smith wrote: > On 08/14/2012 05:09 AM, Bruce Evans wrote: > > When x is incredibly large (close to DBL_MAX), a mathematician would consider > x to represent all the numbers between x-x*DBL_EPSILON to x+x*DBL_EPSILON > (approximately), or more precisely, all the numbers that are within 0.5 ULP > of x. > > So as a Mathematician I would prefer to think of sin(close_to_DBL_MAX) as > undefined. (Although as a programmer, I would hate it if it spat out NaN - I > would prefer the meaningless answer.) Isn't it both as a programmer? Mathematicians don't stop at DBL_MAX, but go to real infinity and then large cardinals :-). Old libraries had a TLOSS error for this. Even fdlibm had this. We seem to have lost a bit by removing this together with historical cruft. Grep shows the following lines matching TLOSS in fdlibm-5.3: % fdlibm.h: * set X_TLOSS = pi*2**52, which is possibly defined in % fdlibm.h:#define X_TLOSS 1.41484755040568800000e+16 % fdlibm.h:#define TLOSS 5 % k_standard.c: * 34-- j0(|x|>X_TLOSS) % k_standard.c: * 35-- y0(x>X_TLOSS) % k_standard.c: * 36-- j1(|x|>X_TLOSS) % k_standard.c: * 37-- y1(x>X_TLOSS) % k_standard.c: * 38-- jn(|x|>X_TLOSS, n) % k_standard.c: * 39-- yn(x>X_TLOSS, n) % k_standard.c: /* j0(|x|>X_TLOSS) */ % k_standard.c: exc.type = TLOSS; % k_standard.c: (void) WRITE2(": TLOSS error\n", 14); % ... Though sin() doesn't lose precision at large or not so large zeros, Bessel functions probably do. I think no one knows where many of their zeros near DBL_MAX are, since they are not evenly spaced like multiples of pi are, so every one needs an individual calculation. pari become incredibly slow at just finding them after the first few hundred, at least using the following perhaps too simple script: for (i = 0, 999, print(solve(x=i*Pi, (i+1)*Pi, besselj(0,x)))) This gets slower and slower as i increases and takes about 1 second per zero starting at i = 900. But after pari takes a few minutes to generate 1000 zeros, FreeBSD libm j0 and j0f pass checks of them all in 1.5 msec. The check is that there is a sign change on each side of the supposed zero. So FreeBSD j0 and j0f get at least 1 bit right (the sign bit) for the first 1000 zeros. > I see code like "if (y!=y) return (y+y)". Does "y+y" quieten the NaNs as > well as "y+0.0"? Is my code compliant in this regard? Yes, y+y works fine and is probably most efficient if there is only 1 NaN to return, since it doesn't require an extra instruction to load 0. >> I only try re-using >> variables near the start of a function (maybe x = creal(z); use(x); >> x = fabs(x)), to try to stop the compiler making so many copies of x. >> This somtimes works. (The typical pessimization avoided by this is >> when x passed on the stack. gcc likes to copy it to another place >> on the stack, using pessimal methods, and then never use the original >> copy. This is good for debugging, but otherwise not very good.). > > What do you mean by "pessimization"? Too much copying of data, without understanding of the memory hierarchy. On some arches, there is a large penalty for loads that don't match stores exactly. gcc-4.2 doesn't understand this, and sometimes enlarges the problem by unnecessary copying that creates mismatches where there were none in the original copies. Bruce