Complex arg-trig functions

[Bruce Evans]

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 <values.h>
% 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