Date: Thu, 6 Sep 2012 23:50:58 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-numerics@freebsd.org Subject: Re: Complex arg-trig functions Message-ID: <20120906221028.O1542@besplex.bde.org> In-Reply-To: <502C0CF8.8040003@missouri.edu> 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> <50295F5C.6010800@missouri.edu> <20120814072946.S5260@besplex.bde.org> <50297CA5.5010900@missouri.edu> <50297E43.7090309@missouri.edu> <20120814201105.T934@besplex.bde.org> <502A780B.2010106@missouri.edu> <20120815223631.N1751@besplex.bde.org> <502C0CF8.8040003@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, 15 Aug 2012, Stephen Montgomery-Smith wrote:
> On 08/15/2012 08:35 AM, Bruce Evans wrote:
>> An up-front check may still be simpler, and gives more control. In
>> csqrt*(), I needed an explicit check and special expressions to get
>> uniform behaviour.
>
> I still like it the way I have it. There is a definite logic in the way infs
> and nans come out of casinh, etc.
Changed now :-).
> There is only one place I disagree with C99:
> catanh(1) = Inf + 0*I
>
> I think mpc gets it correct:
> atanh(1) = Inf + nan*I
C99 is trying be compatible with the real function to a fault. Real
atanh(1) must be Inf, and any imaginary part would make catanh(1)
different. C99 requires similar beahaviour for many other functions.
>> I added this to the NaN mixing in catan[h]*(),
>> and now all my tests pass:
>>
>> % diff -c2 catrig.c~ catrig.c
>> % *** catrig.c~ Sun Aug 12 17:29:18 2012
>> % --- catrig.c Wed Aug 15 11:57:02 2012
>> % ***************
>> % *** 605,609 ****
>> % */
>> % if (ISNAN(x) || ISNAN(y))
>> % ! return (cpack(x+y, x+y));
>> % % /* (x,inf) and (inf,y) and (inf,inf) -> (0,PI/2) */
>> % --- 609,613 ----
>> % */
>> % if (ISNAN(x) || ISNAN(y))
>> % ! return (cpack((x+0.0L)+(y+0), (x+0.0L)+(y+0)));
>> % % /* (x,inf) and (inf,y) and (inf,inf) -> (0,PI/2) */
>>
>> Use this expression in all precisions.
>
> Would this work?
>
> if (ISNAN(x) || ISNAN(y))
> return (cpack((x+x)+(y+y), (x+x)+(y+y)));
No, since this would be evaluated by a different ALU than my version,
if the ALU for long doubles is different from the ALU for doubles, as
happens on amd64. However, I could use y+y instead of y+0 in some
cases. The point of using y+0 is to use the same code and not change
the result much if y happens to be finite (only -0 gets changed, to 0).
ALUs:
(x+0.0L)+(y+0):
x+0.0L:
long double + long double, after promotion of x. Always in the i387.
y+0 (or y+y):
typeof(y) + typeof(y), after promotion of x. In the i387 on i386,
but in SSE on amd64. y+0 ends up quasi-promoted on i386.
final addition:
long double + long double, after promotion of y+0. Always in the i387.
Always long double result. Returning this then converts to double.
(x+x)+(y+y):
x+x:
as above for y+y. Only in the i387 on i386.
y+y:
as above
final addition:
still only in the i387 on i386. So the result of casin() on 2 NaNs is
quite different from that of casinl() on the same 2 NaNs (initially),
since the i387 ALU has different semantics for NaNs than the SSE ALU
despite the ALUs sharing the FPU hardware. Only quieting of NaNs works
perfectly in this version.
>> I forgot to comment it. Adding 0 quietens signaling NaNs before mixing
>> NaNs. I should have tried y+y. Adding 0.0L promotes part of the
>> expression to long double together with quietening signaling NaNs.
>> The rest of the expression is promoted to match. I should try the
>> old way again: of (long double)x+x.
Another explanation. Without testing it, I think it has the same speed
with the x+0[.0L] versions, but is slightly better for the x+x version,
since x+x can be done without another load for 0 in SSE only. But the
extra load is likely to free on x86 since it can be hidden in scheduling
latencies.
>> ...
>> % Index: ../s_csqrt.c
>> % ===================================================================
>> % RCS file: /home/ncvs/src/lib/msun/src/s_csqrt.c,v
>> % retrieving revision 1.4
>> % diff -u -2 -r1.4 s_csqrt.c
>> % --- ../s_csqrt.c 8 Aug 2008 00:15:16 -0000 1.4
>> % +++ ../s_csqrt.c 14 Aug 2012 20:34:07 -0000
>> % @@ -34,14 +34,5 @@
>> % #include "math_private.h"
>> % % -/*
>> % - * gcc doesn't implement complex multiplication or division correctly,
>> % - * so we need to handle infinities specially. We turn on this pragma to
>> % - * notify conforming c99 compilers that the fast-but-incorrect code that
>> % - * gcc generates is acceptable, since the special cases have already
>> been
>> % - * handled.
>> % - */
>> % -#pragma STDC CX_LIMITED_RANGE ON
>>
>> Remove this. There was only 1 complex expression, and it depended on the
>> negation of this pragma to work. Since gcc doesn't support this pragma,
>> the expression only worked accidentally when it was optimized.
>
> I removed it. (I copied it verbatim from csqrt without really understanding
> it.)
Good.
> The part that follows - is this all referencing csqrt?
Depends on whether you clipped an "Index" header :-).
>>
>> % -
>> % -/* We risk spurious overflow for components >= DBL_MAX / (1 +
>> sqrt(2)). */
>> % +/* For avoiding overflow for components >= DBL_MAX / (1 + sqrt(2)). */
>> % #define THRESH 0x1.a827999fcef32p+1022
>>
>> .............. snip
Looks like csqrt.
I didn't check that THRESH is actually correct. It seems to have more
bits than are strictly necessary.
>> This is like a fix in clog(). hypot() handles denormals OK, but
>> necessarily loses accuracy when it returns a denormal result, so
>> the expression (a + hypot(a, b)) is more inaccurate than necessary.
>
> Which code is being referenced here? I use expressions like this catrig.
> Although I think when I use it, I am somewhat certain that neither a nor b
> are denormal.
Now "This" is the change in csqrt() in the snipped quote and is being
compared to related code in clog() that wasn't quoted recently in this
thread.
Snipped quote:
% + scale = 2;
% } else {
% - scale = 0;
% + scale = 1;
% + }
% +
% + /* Scale to reduce inaccuracies when both components are denormal. */
% + if (fabs(a) <= 0x1p-1023 && fabs(b) <= 0x1p-1023) {
% + a *= 0x1p54;
% + b *= 0x1p54;
% + scale = 0x1p-27;
Related code in clog():
% /* Reduce inaccuracies and avoid underflow when ax is denormal. */
% if (kx <= MIN_EXP - 2)
% return (cpack(log(hypot(x * 0x1p1023, y * 0x1p1023)) +
% (MIN_EXP - 2) * ln2_lo + (MIN_EXP - 2) * ln2_hi, v));
%
% /* Avoid remaining underflows (when ax is small but not denormal). */
% if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
% return (cpack(log(hypot(x, y)), v));
This code in catrig.c seems to be related:
@ /*
@ * sum_squares(x,y) = x*x + y*y.
@ * Assumes x*x and y*y will not overflow.
@ * Assumes x and y are finite.
@ * Assumes fabs(x) >= DBL_EPSILON.
@ */
@ inline static double
@ sum_squares(double x, double y)
@ {
@ /*
@ * The following code seems to do better than
@ * return (hypot(x, y) * hypot(x, y));
@ */
@
@ if (fabs(y) < MIN_4TH_ROOT)
@ if ((int)y==0) /* raise inexact */
@ return (x*x);
@ return (x*x + y*y);
@ }
but it isn't really -- it is more like the case of a large difference
of exponents, so x*x+y*y reduces to x*x (at a smaller exponent difference
than for hypot(x, y) reducing to |x|). hypot() sets inexact by returning
|x|+|y| in that case, but the above avoids using y*y since it would give
spurious underflow.
Hmm, MIN_4TH_ROOT seems to be logically wrong here, and physically
wrong for float precision. In double precision, MIN_4TH_ROOT/DBL_EPSILON
is ~0.5e-59, but in float precision, MIN_4TH_ROOT/FLT_EPSILON is only
~1e-2, so when x = FLT_EPSILON, y = MIN_4TH_ROOT is very far from giving
a discardable y*y.
The logically correct threshold is something like y/x < 2**-MANT_DIG/2
or y < 2**-MANT_DIG/2 * EPSILON to avoid the division. clog*() uses
essentially y/x < 2**-MANT_DIG for some reason (it actually uses a fuzzy
check based on exponent differences: kx - ky > MANT_DIG; MANT_DIG/2
would be a little too small here, but doubling it seems wasteful).
hypot() uses essentially y/x < 2**-60, where the magic 60 is DBL_MANT_DIG
plys a safety margin. hypotf() uses a magic 30. The non-magic integer
exponent thresholds work very well in clog*(), but depend on having the
exponents handy for efficiency, which may be possible via optimization
for an inline function like the above but takes more source code.
I don't see how the above avoids cancelation problems like the ones in
clog(). You never subtract 1 from sum_squares(), but there are open-coded
expressions like '(1-ax)*(1+ax) - ay*ay' to reduce the cancelation error.
I found that for clog() it was worth using extra precision out to about
|z| = 4 (or |z| < 1/4) to reduce the cancelation error. But somehow my
tests didn't show large errors. Perhaps because I'm happy with < 4 ulps
instead of < 1.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120906221028.O1542>