Complex arg-trig functions
Bruce Evans
brde at optusnet.com.au
Mon Aug 13 16:57:50 UTC 2012
On Sun, 12 Aug 2012, Stephen Montgomery-Smith wrote:
> Having brooded over the code for too many weeks, I now think I have finished
> my complex arg-trig functions. I have also written versions for float and
> long. So I am ready to have the code reviewed.
>
> http://people.freebsd.org/~stephen/
I finally tested a version of this. I only did simple comparisons (float
vs double and double vs long double). The results look promising after
fixing a few bugs:
% amd64 float prec, on 2**12 * 2**12 args:
% rcacos:max_er = 0x58460841 2.7585, avg_er = 0.317, #>=1:0.5 = 29084:255712
% rcacosh:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
% rcasin:max_er = 0x631b8183 3.0971, avg_er = 0.209, #>=1:0.5 = 38388:382508
% rcasinh:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
% rcatan:max_er = 0x51d7c47a 2.5576, avg_er = 0.290, #>=1:0.5 = 52984:318084
% rcatanh:max_er = 0x3693fccd4f 436.6246, avg_er = 0.223, #>=1:0.5 = 213124:1438280
% rclog: max_er = 0x26dfae4d 1.2148, avg_er = 0.247, #>=1:0.5 = 184:92244
% icacos:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
% icacosh:max_er = 0x3690000000 436.5000, avg_er = 0.317, #>=1:0.5 = 29104:255732
% icasin:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
% icasinh:max_er = 0x631b8183 3.0971, avg_er = 0.209, #>=1:0.5 = 38388:382508
% icatan:max_er = 0x3693fccd4f 436.6246, avg_er = 0.223, #>=1:0.5 = 213124:1438280
% icatanh:max_er = 0x51d7c47a 2.5576, avg_er = 0.290, #>=1:0.5 = 52984:318084
% iclog: max_er = 0x1fc2b4f5 0.9925, avg_er = 0.302, #>=1:0.5 = 0:349830
%
% amd64 double prec, on 2**12 x 2**12 args:
% rcacos:max_er = 0x1b5a 3.4189, avg_er = 0.228, #>=1:0.5 = 2394:125988
% rcacosh:max_er = 0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741860
% rcasin:max_er = 0x15c5 2.7212, avg_er = 0.113, #>=1:0.5 = 33296:99152
% rcasinh:max_er = 0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741796
% rcatan:max_er = 0x8000000000000000 4503599627370496.0000, avg_er = 268435456.212, #>=1:0.5 = 4681:81365
% rcatanh:max_er = 0x8000000000000000 4503599627370496.0000, avg_er = 268435456.047, #>=1:0.5 = 428997:691341
% rclog: max_er = 0x704 0.8770, avg_er = 0.250, #>=1:0.5 = 0:20152
% icacos:max_er = 0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741860
% icacosh:max_er = 0x1b5a 3.4189, avg_er = 0.228, #>=1:0.5 = 2394:125988
% icasin:max_er = 0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741796
% icasinh:max_er = 0x15c5 2.7212, avg_er = 0.113, #>=1:0.5 = 33296:99152
% icatan:max_er = 0x8000000000000000 4503599627370496.0000, avg_er = 268435456.047, #>=1:0.5 = 428997:691341
% icatanh:max_er = 0x8000000000000000 4503599627370496.0000, avg_er = 268435456.212, #>=1:0.5 = 4681:81365
% iclog: max_er = 0x6f4 0.8691, avg_er = 0.213, #>=1:0.5 = 0:181032
rfoo is the real part of foo, etc.
2**12 x 2**12 args is not enough. Some bugs showed up only with slightly
more args, but after fixing some bugs the behaviour doesn't seem to
depend so much on the number.
The error of 436.6246 upls has been turning up a lot. It was for a bug
in my test program for a result of 0 when the correctly rounded result is
the smallest denormal, but I thought I fixed it.
The error of 0x8000000000000000 for double precision *catan* is a sign
mismatch.
% i386 float prec, on 2**12 x 2**12 args:
% rcacos:max_er = 0x42cd19c6 2.0875, avg_er = 0.314, #>=1:0.5 = 3854:215116
% rcacosh:max_er = 0x3170e232 1.5450, avg_er = 0.254, #>=1:0.5 = 23008:3245028
% rcasin:max_er = 0x55adc0df 2.6775, avg_er = 0.208, #>=1:0.5 = 34304:353980
% rcasinh:max_er = 0x3170e232 1.5450, avg_er = 0.254, #>=1:0.5 = 23008:3245028
% rcatan:max_er = 0x3c4078ec 1.8829, avg_er = 0.284, #>=1:0.5 = 13260:190836
% rcatanh:max_er = 0x3693fccd4f 436.6246, avg_er = 0.186, #>=1:0.5 = 4796:421616
% rclog: max_er = 0x25830853 1.1722, avg_er = 0.246, #>=1:0.5 = 120:24892
% icacos:max_er = 0x3170e232 1.5450, avg_er = 0.254, #>=1:0.5 = 23008:3245028
% icacosh:max_er = 0x3690000000 436.5000, avg_er = 0.315, #>=1:0.5 = 3874:215136
% icasin:max_er = 0x3170e232 1.5450, avg_er = 0.254, #>=1:0.5 = 23008:3245028
% icasinh:max_er = 0x55adc0df 2.6775, avg_er = 0.208, #>=1:0.5 = 34304:353980
% icatan:max_er = 0x3693fccd4f 436.6246, avg_er = 0.186, #>=1:0.5 = 4796:421616
% icatanh:max_er = 0x3c4078ec 1.8829, avg_er = 0.284, #>=1:0.5 = 13260:190836
% iclog: max_er = 0x1fc2b4f5 0.9925, avg_er = 0.302, #>=1:0.5 = 0:338712
%
% i386 double prec, on 2**12 x 2**12 args:
% rcacos:max_er = 0x11e8 2.2383, avg_er = 0.165, #>=1:0.5 = 248:111850
% rcacosh:max_er = 0xb02 1.3760, avg_er = 0.256, #>=1:0.5 = 104:2715312
% rcasin:max_er = 0x13ce 2.4756, avg_er = 0.112, #>=1:0.5 = 5616:95060
% rcasinh:max_er = 0xb02 1.3760, avg_er = 0.256, #>=1:0.5 = 104:2715312
% rcatan:max_er = 0x9ed 1.2407, avg_er = 0.015, #>=1:0.5 = 4084:48920
% rcatanh:max_er = 0xb17 1.3862, avg_er = 0.014, #>=1:0.5 = 56:77456
% rclog: max_er = 0x704 0.8770, avg_er = 0.250, #>=1:0.5 = 0:20112
% icacos:max_er = 0xb02 1.3760, avg_er = 0.256, #>=1:0.5 = 104:2715312
% icacosh:max_er = 0x11e8 2.2383, avg_er = 0.165, #>=1:0.5 = 248:111850
% icasin:max_er = 0xb02 1.3760, avg_er = 0.256, #>=1:0.5 = 104:2715312
% icasinh:max_er = 0x13ce 2.4756, avg_er = 0.112, #>=1:0.5 = 5616:95060
% icatan:max_er = 0xb17 1.3862, avg_er = 0.014, #>=1:0.5 = 56:77456
% icatanh:max_er = 0x9ed 1.2407, avg_er = 0.015, #>=1:0.5 = 4084:48920
% iclog: max_er = 0x6f4 0.8691, avg_er = 0.213, #>=1:0.5 = 0:181032
Note that i386 doesn't have the sign errors. However, i386 with -O0 has
the sign errors for at least float precision in the same places that
amd64 with -O has them for double precision.
> The long versions require a logl and a log1pl, which I faked using mpfr.
>
> The float versions are more complicated, because FLT_EPSILON is too close to
> the 4th root of FLT_MIN. It is simpler to make the float versions wrappers
> for the double versions. But I wrote the float versions anyway, just in case
> some purist insists that the wrapper approach is morally wrong.
There are negative reasons to have the float versions unless they are not
wrappers. The reasons to have non-wrappers are to test the algorithm and
run faster.
Fixes needed for the above test results:
@ diff -c2 catrig.c~ catrig.c
@ *** catrig.c~ Sun Aug 12 17:29:18 2012
@ --- catrig.c Mon Aug 13 12:07:09 2012
@ ***************
@ *** 265,269 ****
@ return;
@ }
@ ! if (y < MIN_4TH_ROOT) {
@ /*
@ * Avoid a possible underflow caused by y/A. For casinh this
@ --- 265,269 ----
@ return;
@ }
@ ! if (!ISNAN(x) && y < MIN_4TH_ROOT) {
@ /*
@ * Avoid a possible underflow caused by y/A. For casinh this
This stops the result for NaN x depending on the size of y relative
to MIN_4TH_ROOT. MIN_4TH_ROOT varies with the precision, but the result
shouldn't vary with the precision. The result should probably be the
original NaN x quieted (unless y is also NaN). This happens naturally
if you do a computation with x (except -x). However, symmetry often
requires forcing or flipping signs, and we force or flip the sign even
for NaNs. This change gives consistent signs for NaN results, by taking
the same path in all precisions.
@ ***************
@ *** 408,416 ****
@
@ if (ISFINITE(bx) && ISFINITE(by) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ ! if (huge+x+y>one) { /* raise inexact flag */
@ ! w = clog_for_large_values(z) + M_LN2;
Addition clobbers the sign of -0 in the imaginary part in some cases.
I think it should clobber it in all cases, since it does for non-complex
doubles (-0.0 + +0.0 is +0.0). Clobbering was observed in the following
cases:
- on i386 (-march=athlon64) with -O0, presumably because the addition to
the imaginary part was not optimized away then
- on amd64 with -O0 and -O, presumably because:
- same as i386 for -O0
- with -O (no -march), float complex and double complex are represented
as a vector in an SSE2 register, and it is natural to add both the
components of the vector, and this may be optimal for at least float
complex with the default -march. I didn't check what happens for
double complex.
@ if (sy == 0)
@ ! return (cpack(cimag(w), -creal(w)));
@ ! return (cpack(-cimag(w), creal(w)));
The sign of creal(cacos()) is always 1, but this makes it +- the sign
of atan2(x, y).
@ }
@ }
@ --- 408,420 ----
@
@ if (ISFINITE(bx) && ISFINITE(by) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ ! /* XXX following can also raise overflow */
Just note this error.
@ ! if (huge+x+y>one) { /* raise inexact */
Start removing ' flag'.
@ ! w = clog_for_large_values(z);
@ ! /* Can't add M_LN2 to w since it should clobber -0*I. */
@ ! rx = fabs(cimag(w));
@ ! ry = creal(w) + M_LN2;
@ if (sy == 0)
@ ! ry = -ry;
@ ! return (cpack(rx, ry));
@ }
@ }
Fix the above bugs.
@ ***************
@ *** 482,486 ****
@ * but this case should happen extremely rarely.
@ */
@ ! if (ay > 0.5*DBL_MAX)
@ return (cpack(log(hypot(x / M_E, y / M_E)) + 1, atan2(y, x)));
@
@ --- 486,490 ----
@ * but this case should happen extremely rarely.
@ */
@ ! if (ax > 0.5*DBL_MAX)
@ return (cpack(log(hypot(x / M_E, y / M_E)) + 1, atan2(y, x)));
@
The smallest of ax and ay must be compared. I think I broke this in an
early version of clog(). Replacing the whole function by clog() made
little difference once this was fixed.
@ diff -c2 catrigf.c~ catrigf.c
@ *** catrigf.c~ Sun Aug 12 17:00:52 2012
@ --- catrigf.c Mon Aug 13 14:14:42 2012
@ ***************
@ *** 138,142 ****
@ return;
@ }
@ ! if (y < MIN_4TH_ROOT) {
@ *B_is_usable = 0;
@ if ((int)y==0) /* raise inexact flag */
@ --- 138,142 ----
@ return;
@ }
@ ! if (!isnan(x) && y < MIN_4TH_ROOT) {
@ *B_is_usable = 0;
@ if ((int)y==0) /* raise inexact flag */
@ ***************
@ *** 233,240 ****
@ if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ if (huge+x+y>one) { /* raise inexact flag */
@ ! w = clog_for_large_values(z) + M_LN2;
@ if (sy == 0)
@ ! return (cpackf(cimagf(w), -crealf(w)));
@ ! return (cpackf(-cimagf(w), crealf(w)));
@ }
@ }
@ --- 233,242 ----
@ if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ if (huge+x+y>one) { /* raise inexact flag */
@ ! w = clog_for_large_values(z);
@ ! rx = fabsf(cimagf(w));
@ ! ry = crealf(w) + M_LN2;
@ if (sy == 0)
@ ! ry = -ry;
@ ! return (cpackf(rx, ry));
@ }
@ }
@ ***************
@ *** 290,294 ****
@ }
@
@ ! if (ay > 0.5*FLT_MAX)
@ return (cpackf(logf(hypotf(x / M_E, y / M_E)) + 1, atan2f(y, x)));
@
@ --- 292,296 ----
@ }
@
@ ! if (ax > 0.5*FLT_MAX)
@ return (cpackf(logf(hypotf(x / M_E, y / M_E)) + 1, atan2f(y, x)));
@
@ diff -c2 catrigl.c~ catrigl.c
@ *** catrigl.c~ Sun Aug 12 06:54:46 2012
@ --- catrigl.c Mon Aug 13 12:08:21 2012
@ ***************
@ *** 119,123 ****
@ return;
@ }
@ ! if (y < MIN_4TH_ROOT) {
@ *B_is_usable = 0;
@ if ((int)y==0) /* raise inexact flag */
@ --- 119,123 ----
@ return;
@ }
@ ! if (!isnan(x) && y < MIN_4TH_ROOT) {
@ *B_is_usable = 0;
@ if ((int)y==0) /* raise inexact flag */
@ ***************
@ *** 207,214 ****
@ if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ if (huge+x+y>one) { /* raise inexact flag */
@ ! w = clog_for_large_values(z) + L_LN2;
@ if (sy == 0)
@ ! return (cpackl(cimagl(w), -creall(w)));
@ ! return (cpackl(-cimagl(w), creall(w)));
@ }
@ }
@ --- 207,216 ----
@ if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ if (huge+x+y>one) { /* raise inexact flag */
@ ! w = clog_for_large_values(z);
@ ! rx = fabsl(cimagl(w));
@ ! ry = creall(w) + M_LN2;
@ if (sy == 0)
@ ! ry = -ry;
@ ! return (cpackl(rx, ry));
@ }
@ }
@ ***************
@ *** 264,268 ****
@ }
@
@ ! if (ay > 0.5*LDBL_MAX)
@ return (cpackl(logl(hypotl(x / L_E, y / L_E)) + 1, atan2l(y, x)));
@
@ --- 266,270 ----
@ }
@
@ ! if (ax > 0.5*LDBL_MAX)
@ return (cpackl(logl(hypotl(x / L_E, y / L_E)) + 1, atan2l(y, x)));
@
Some other bugs that I noticed:
catrigf.c:
% ...
% #define MAX_4TH_ROOT 1e9 /* approx pow(FLT_MAX,0.25) */
% #define MIN_4TH_ROOT 1e-9 /* approx pow(FLT_MIN,0.25) */
% #define SQRT_EPSILON_100 1e-6 /* approx (sqrt(FLT_EPSILON)/100) */
% /* 100 is to "play it safe." */
% #define RECIP_SQRT_EPSILON_100 1e6 /* 1/SQRT_EPSILON_100 */
% #define EPSILON_100 1e-9 /* approx FLT_EPSILON/100 */
% #define RECIP_EPSILON_100 1e9 /* 1/EPSILON_100 */
These aren't float constants. Mixing them with floats may give unwanted
extra precision and slowness.
Most of these should be hex constants, so that they are easier to
understand. In clog() I originaly used magic hex constants, but but
almost everything was a power of 2 related to FOO_MAX or FOO_MANT_DIG,
as above, but it is difficult to translate from those to hex constants
(FLT_EPSILON / 128 could be done using token pasting as
((0x1p # FLT_MIN_EXP) * 2 / 128), but square and fourth roots are harder.
After changing the classification using bits, all the comparisions with
thresholds became comparisons of exponents, and it is now trivial to
take roots by dividing the expondent. ld128 clogl() now works with
identical code to ld80 clogl(), based on macros in <float.h>, and the
differences for other precisions are minor (I could hide them all using
more macros).
% ...
% static const float
% one = 1.00000000000000000000e+00,
% huge= 1.00000000000000000000e+30;
Using variables instead of macros avoids having to put F suffixes on
all float constants.
% ...
% inline static float
% f(float a, float b, float hypot_a_b)
% {
% if (b < 0) return (0.5 * (hypot_a_b - b));
% if (b == 0) return (0.5*a);
% return (0.5 * a*a / (hypot_a_b + b));
% }
There a squillions of 0.5's that should be 0.5F's and deserve being named
`half' more than 1 deserves being named 'one'.
catrigl.c:
% #define MAX_4TH_ROOT 1e1200L /* approx pow(LDBL_MAX,0.25) */
% #define MIN_4TH_ROOT 1e-1200L /* approx pow(LDBL_MIN,0.25) */
% #define SQRT_EPSILON_100 1e-13L /* approx (sqrtl(LDBL_EPSILON)/100) */
% #define RECIP_SQRT_EPSILON_100 1e13L /* 1/SQRT_EPSILON_100 */
% #define EPSILON_100 1e-24L /* approx LDBL_EPSILON/100 */
% #define RECIP_EPSILON_100 1e24L /* 1/EPSILON_100 */
Now there is no need for a suffix on the smaller constants. You don't
use one for the 0.5's later.
The ones related to epsilon look wrong for ld128. LDBL_EPSILON =
2**-(LDBL_MANT_DIG - 1) is way smaller when LDBL_MANT_DIG is 113 than
when it is 64.
I'm fairly happy with the attached clog() now. All refinements that
I tried lately give either less accuracy, lower speed or larger code.
Ones that don't quite work are:
- write log(ax*ax + ay*ay) * 0.5 = log(ax) + log1p(1 + (ay/ax)**2) * 0.5.
This avoids complications when ax is larger or ay is small (after
ensuring than ay/ax is not small), but tends to be slower and a
bit less accurate.
- write ax = 2**k * bx; ay = 2**k * by; log(ax*ax + ay*ay) * 0.5 =
k * log(2) + log(bx*bx + by*by) * 0.5. This gives marginally more
accuracy, but tends to be slower and a bit more complicated. One
complication is that bx*bx + by*by may need another scaling by
a factor of 2 or 1/2 to make it between sqrt(2)/2 and sqrt(2).
Many changes in this version:
% /* Avoid overflow. */
% if (kx >= MAX_EXP - 1)
% return (cpack(log(hypot(x * 0x1p-1022, y * 0x1p-1022)) +
% (MAX_EXP - 2) * ln2_lo + (MAX_EXP - 2) * ln2_hi, v));
% if (kx >= (MAX_EXP - 1) / 2)
% return (cpack(log(hypot(x, y)), v));
1. Use integer exponenents kx and ky for all classifications.
2. Get all the thresholds right.
3. Scale down to near 1 here, instead of just by a factor of 2. This
gives some extra accuracy for no runtime cost.
%
% /* 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));
4. Scale up to near 1 here, instead of just by a factor of 2. This
gives considerable extra accuracy for no runtime cost.
%
% /* 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));
5. Don't scale here. The correct scale factor is ~ 2**-kx for both cases,
but that has more overheads. So we only scale up when necessary to
avoid losing considerable accuracy (for denormal ax), and let log()
do the scaling for other cases.
% /* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
% t = (double)(ax * (0x1p27 + 1));
% axh = (double)(ax - t) + t;
% axl = ax - axh;
% ax2h = ax * ax;
% ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
% t = (double)(ay * (0x1p27 + 1));
% ayh = (double)(ay - t) + t;
% ayl = ay - ayh;
% ay2h = ay * ay;
% ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
6. Additive splitting like I had for clogl() before doesn't work unless
ax is near 1. Use the correct Veldkamp algorithm after googling
Dekker's algorithms. ax has type double_t, so casting to double
avoids compiler bugfeatures at some cost.
% sh -= 1;
% norm(sh, sl);
% norm(ax2l, ay2l);
% /* Briggs-Kahan algorithm (except we discard the final low term): */
% norm(sh, ax2l);
% norm(sl, ay2l);
% t = ax2l + sl;
% normF(sh, t);
% return (cpack(log1p(ay2l + t + sh) * 0.5, v));
7. Use a full modernised Dekker's algorithm for the critical case of
|z| very near 1 (and for many cases not so near 1). There are 2
Veldkamp splittings, 4 norm() (= 2sum) calls and 2 normF() (= 2sumF)
calls on the way here. norm() takes 6 additions and normF() takes 3,
so there are 30 additions just for the norm*()'s. But this is only
moderately slow (about the same speed as extra branches to avoid
doing it all so much).
Bruce
-------------- next part --------------
#include <complex.h>
#include <float.h>
#include "fpmath.h"
#include "local.h"
#include "math.h"
#include "math_private.h"
/* 2sum algorithm. */
#undef norm
#define norm(a, b) do { \
__typeof(a) __s, __w; \
\
__w = (a) + (b); \
__s = __w - (a); \
(b) = ((a) - (__w - __s)) + ((b) - __s); \
(a) = __w; \
} while (0)
/* 2sumF algorithm. */
#define normF(a, b) do { \
__typeof(a) __w; \
\
__w = (a) + (b); \
(b) = ((a) - __w) + (b); \
(a) = __w; \
} while (0)
#define MANT_DIG DBL_MANT_DIG
#define MAX_EXP DBL_MAX_EXP
#define MIN_EXP DBL_MIN_EXP
static const double
ln2_hi = 6.9314718055829871e-1, /* 0x162e42fefa0000.0p-53 */
ln2_lo = 1.6465949582897082e-12; /* 0x1cf79abc9e3b3a.0p-92 */
double complex
clog(double complex z)
{
double_t ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl, sh, sl, t;
double x, y, v;
uint32_t hax, hay;
int kx, ky;
x = creal(z);
y = cimag(z);
v = atan2(y, x);
ax = fabs(x);
ay = fabs(y);
if (ax < ay) {
t = ax;
ax = ay;
ay = t;
}
GET_HIGH_WORD(hax, ax);
kx = (hax >> 20) - 1023;
GET_HIGH_WORD(hay, ay);
ky = (hay >> 20) - 1023;
/* Handle NaNs and Infs using the general formula. */
if (kx == MAX_EXP || ky == MAX_EXP)
return (cpack(log(hypot(x, y)), v));
/* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
if (ax == 1) {
if (ky < (MIN_EXP - 1) / 2)
return (cpack((ay * 0.5) * ay, v));
return (cpack(log1p(ay * ay) * 0.5, v));
}
/* Avoid underflow when ax is not small. Also handle zero args. */
if (kx - ky > MANT_DIG || ay == 0)
return (cpack(log(ax), v));
/* Avoid overflow. */
if (kx >= MAX_EXP - 1)
return (cpack(log(hypot(x * 0x1p-1022, y * 0x1p-1022)) +
(MAX_EXP - 2) * ln2_lo + (MAX_EXP - 2) * ln2_hi, v));
if (kx >= (MAX_EXP - 1) / 2)
return (cpack(log(hypot(x, y)), v));
/* 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));
/* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
t = (double)(ax * (0x1p27 + 1));
axh = (double)(ax - t) + t;
axl = ax - axh;
ax2h = ax * ax;
ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
t = (double)(ay * (0x1p27 + 1));
ayh = (double)(ay - t) + t;
ayl = ay - ayh;
ay2h = ay * ay;
ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
/*
* When log(|z|) is far from 1, accuracy in calculating the sum
* of the squares is not very important since log() reduces
* inaccuracies. We depended on this to use the general
* formula when log(|z|) is very far from 1. When log(|z|) is
* moderately far from 1, we go through the extra-precision
* calculations to reduce branches and gain a little accuracy.
*
* When |z| is near 1, we subtract 1 and use log1p() and don't
* leave it to log() to subtract 1, since we gain at least 1 bit
* of accuracy in this way.
*
* When |z| is very near 1, subtracting 1 can cancel almost
* 3*MANT_DIG bits. We arrange that subtracting 1 is exact in
* doubled precision, and then do the rest of the calculation
* in sloppy doubled precision. Although large cancelations
* often lose lots of accuracy, here the final result is exact
* in doubled precision if the large calculation occurs (because
* then it is exact in tripled precision and the cancelation
* removes enough bits to fit in doubled precision). Thus the
* result is accurate in sloppy doubled precision, and the only
* significant loss of accuracy is when it is summed and passed
* to log1p().
*/
sh = ax2h;
sl = ay2h;
normF(sh, sl);
if (sh < 0.5 || sh >= 3)
return (cpack(log(ay2l + ax2l + sl + sh) * 0.5, v));
sh -= 1;
norm(sh, sl);
norm(ax2l, ay2l);
/* Briggs-Kahan algorithm (except we discard the final low term): */
norm(sh, ax2l);
norm(sl, ay2l);
t = ax2l + sl;
normF(sh, t);
return (cpack(log1p(ay2l + t + sh) * 0.5, v));
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(clog, clogl);
#endif
#undef MANT_DIG
#define MANT_DIG FLT_MANT_DIG
#undef MAX_EXP
#define MAX_EXP FLT_MAX_EXP
#undef MIN_EXP
#define MIN_EXP FLT_MIN_EXP
static const float
ln2f_hi = 6.9314575195e-1, /* 0xb17200.0p-24 */
ln2f_lo = 1.4286067653e-6; /* 0xbfbe8e.0p-43 */
float complex
clogf(float complex z)
{
float_t ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl, sh, sl, t;
float x, y, v;
uint32_t hax, hay;
int kx, ky;
x = crealf(z);
y = cimagf(z);
v = atan2f(y, x);
ax = fabsf(x);
ay = fabsf(y);
if (ax < ay) {
t = ax;
ax = ay;
ay = t;
}
GET_FLOAT_WORD(hax, ax);
kx = (hax >> 23) - 127;
GET_FLOAT_WORD(hay, ay);
ky = (hay >> 23) - 127;
/* Handle NaNs and Infs using the general formula. */
if (kx == MAX_EXP || ky == MAX_EXP)
return (cpackf(logf(hypotf(x, y)), v));
/* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
if (hax == 0x3f800000) {
if (ky < (MIN_EXP - 1) / 2)
return (cpackf((ay * 0.5F) * ay, v));
return (cpackf(log1pf(ay * ay) * 0.5F, v));
}
/* Avoid underflow when ax is not small. Also handle zero args. */
if (kx - ky > MANT_DIG || hay == 0)
return (cpackf(logf(ax), v));
/* Avoid overflow. */
if (kx >= MAX_EXP - 1)
return (cpackf(logf(hypotf(x * 0x1p-126F, y * 0x1p-126F)) +
(MAX_EXP - 2) * ln2_lo + (MAX_EXP - 2) * ln2_hi, v));
if (kx >= (MAX_EXP - 1) / 2)
return (cpackf(logf(hypotf(x, y)), v));
/* Reduce inaccuracies and avoid underflow when ax is denormal. */
if (kx <= MIN_EXP - 2)
return (cpackf(logf(hypotf(x * 0x1p127F, y * 0x1p127F)) +
(MIN_EXP - 2) * ln2f_lo + (MIN_EXP - 2) * ln2f_hi, v));
/* Avoid remaining underflows (when ax is small but not denormal). */
if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
return (cpackf(logf(hypotf(x, y)), v));
/* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
t = (float)(ax * (0x1p12F + 1));
axh = (float)(ax - t) + t;
axl = ax - axh;
ax2h = ax * ax;
ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
t = (float)(ay * (0x1p12F + 1));
ayh = (float)(ay - t) + t;
ayl = ay - ayh;
ay2h = ay * ay;
ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
/*
* When log(|z|) is far from 1, accuracy in calculating the sum
* of the squares is not very important since log() reduces
* inaccuracies. We depended on this to use the general
* formula when log(|z|) is very far from 1. When log(|z|) is
* moderately far from 1, we go through the extra-precision
* calculations to reduce branches and gain a little accuracy.
*
* When |z| is near 1, we subtract 1 and use log1p() and don't
* leave it to log() to subtract 1, since we gain at least 1 bit
* of accuracy in this way.
*
* When |z| is very near 1, subtracting 1 can cancel almost
* 3*MANT_DIG bits. We arrange that subtracting 1 is exact in
* doubled precision, and then do the rest of the calculation
* in sloppy doubled precision. Although large cancelations
* often lose lots of accuracy, here the final result is exact
* in doubled precision if the large calculation occurs (because
* then it is exact in tripled precision and the cancelation
* removes enough bits to fit in doubled precision). Thus the
* result is accurate in sloppy doubled precision, and the only
* significant loss of accuracy is when it is summed and passed
* to log1p().
*/
sh = ax2h;
sl = ay2h;
normF(sh, sl);
if (sh < 0.5F || sh >= 3)
return (cpackf(logf(ay2l + ax2l + sl + sh) * 0.5F, v));
sh -= 1;
norm(sh, sl);
norm(ax2l, ay2l);
/* Briggs-Kahan algorithm (except we discard the final low term): */
norm(sh, ax2l);
norm(sl, ay2l);
t = ax2l + sl;
normF(sh, t);
return (cpackf(log1pf(ay2l + t + sh) * 0.5F, v));
}
#undef MANT_DIG
#define MANT_DIG LDBL_MANT_DIG
#undef MAX_EXP
#define MAX_EXP LDBL_MAX_EXP
#undef MIN_EXP
#define MIN_EXP LDBL_MIN_EXP
long double complex
clogl(long double complex z)
{
long double ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl;
long double sh, sl, t;
long double x, y, v;
uint16_t hax, hay;
int kx, ky;
x = creall(z);
y = cimagl(z);
v = atan2l(y, x);
ax = fabsl(x);
ay = fabsl(y);
if (ax < ay) {
t = ax;
ax = ay;
ay = t;
}
GET_LDBL_EXPSIGN(hax, ax);
kx = hax - 16383;
GET_LDBL_EXPSIGN(hay, ay);
ky = hay - 16383;
/* Handle NaNs and Infs using the general formula. */
if (kx == MAX_EXP || ky == MAX_EXP)
return (cpackl(logl(hypotl(x, y)), v));
/* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
if (ax == 1) {
if (ky < (MIN_EXP - 1) / 2)
return (cpackl((ay * 0.5) * ay, v));
return (cpackl(log1pl(ay * ay) * 0.5, v));
}
/* Avoid underflow when ax is not small. Also handle zero args. */
if (kx - ky > MANT_DIG || ay == 0)
return (cpackl(logl(ax), v));
/* Avoid overflow. */
if (kx >= MAX_EXP - 1)
return (cpackl(logl(hypotl(x * 0x1p-16382L, y * 0x1p-16382L)) +
(MAX_EXP - 2) * ln2_lo + (MAX_EXP - 2) * ln2_hi, v));
if (kx >= (MAX_EXP - 1) / 2)
return (cpackl(logl(hypotl(x, y)), v));
/* Reduce inaccuracies and avoid underflow when ax is denormal. */
if (kx <= MIN_EXP - 2)
return (cpackl(logl(hypotl(x * 0x1p16383L, y * 0x1p16383L)) +
(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 (cpackl(logl(hypotl(x, y)), v));
/* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
t = (long double)(ax * (0x1p32 + 1));
axh = (long double)(ax - t) + t;
axl = ax - axh;
ax2h = ax * ax;
ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
t = (long double)(ay * (0x1p32 + 1));
ayh = (long double)(ay - t) + t;
ayl = ay - ayh;
ay2h = ay * ay;
ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
/*
* When log(|z|) is far from 1, accuracy in calculating the sum
* of the squares is not very important since log() reduces
* inaccuracies. We depended on this to use the general
* formula when log(|z|) is very far from 1. When log(|z|) is
* moderately far from 1, we go through the extra-precision
* calculations to reduce branches and gain a little accuracy.
*
* When |z| is near 1, we subtract 1 and use log1p() and don't
* leave it to log() to subtract 1, since we gain at least 1 bit
* of accuracy in this way.
*
* When |z| is very near 1, subtracting 1 can cancel almost
* 3*MANT_DIG bits. We arrange that subtracting 1 is exact in
* doubled precision, and then do the rest of the calculation
* in sloppy doubled precision. Although large cancelations
* often lose lots of accuracy, here the final result is exact
* in doubled precision if the large calculation occurs (because
* then it is exact in tripled precision and the cancelation
* removes enough bits to fit in doubled precision). Thus the
* result is accurate in sloppy doubled precision, and the only
* significant loss of accuracy is when it is summed and passed
* to log1p().
*/
sh = ax2h;
sl = ay2h;
normF(sh, sl);
if (sh < 0.5 || sh >= 3)
return (cpackl(logl(ay2l + ax2l + sl + sh) * 0.5, v));
sh -= 1;
norm(sh, sl);
norm(ax2l, ay2l);
/* Briggs-Kahan algorithm (except we discard the final low term): */
norm(sh, ax2l);
norm(sl, ay2l);
t = ax2l + sl;
normF(sh, t);
return (cpackl(log1pl(ay2l + t + sh) * 0.5, v));
}
#define GEN(func, type, ltype, part, lpart) \
type \
lpart ## c ## func ## ltype(type x, type y) \
{ \
type complex z; \
type complex w; \
\
z = cpack ## ltype(x, y); \
w = c ## func ## ltype(z); \
return (c ## part ## ltype(w)); \
}
#define GEN3(func) \
GEN(func, double, , real, r) \
GEN(func, double, , imag, i) \
GEN(func, float, f, real, r) \
GEN(func, float, f, imag, i) \
GEN(func, long double, l, real, r) \
GEN(func, long double, l, imag, i)
GEN3(acos)
GEN3(acosh)
GEN3(asin)
GEN3(asinh)
GEN3(atan)
GEN3(atanh)
GEN3(log)
-------------- next part --------------
diff -c2 catrig.c~ catrig.c
*** catrig.c~ Sun Aug 12 17:29:18 2012
--- catrig.c Mon Aug 13 12:07:09 2012
***************
*** 265,269 ****
return;
}
! if (y < MIN_4TH_ROOT) {
/*
* Avoid a possible underflow caused by y/A. For casinh this
--- 265,269 ----
return;
}
! if (!ISNAN(x) && y < MIN_4TH_ROOT) {
/*
* Avoid a possible underflow caused by y/A. For casinh this
***************
*** 408,416 ****
if (ISFINITE(bx) && ISFINITE(by) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
! if (huge+x+y>one) { /* raise inexact flag */
! w = clog_for_large_values(z) + M_LN2;
if (sy == 0)
! return (cpack(cimag(w), -creal(w)));
! return (cpack(-cimag(w), creal(w)));
}
}
--- 408,420 ----
if (ISFINITE(bx) && ISFINITE(by) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
! /* XXX following can also raise overflow */
! if (huge+x+y>one) { /* raise inexact */
! w = clog_for_large_values(z);
! /* Can't add M_LN2 to w since it should clobber -0*I. */
! rx = fabs(cimag(w));
! ry = creal(w) + M_LN2;
if (sy == 0)
! ry = -ry;
! return (cpack(rx, ry));
}
}
***************
*** 482,486 ****
* but this case should happen extremely rarely.
*/
! if (ay > 0.5*DBL_MAX)
return (cpack(log(hypot(x / M_E, y / M_E)) + 1, atan2(y, x)));
--- 486,490 ----
* but this case should happen extremely rarely.
*/
! if (ax > 0.5*DBL_MAX)
return (cpack(log(hypot(x / M_E, y / M_E)) + 1, atan2(y, x)));
diff -c2 catrigf.c~ catrigf.c
*** catrigf.c~ Sun Aug 12 17:00:52 2012
--- catrigf.c Mon Aug 13 14:14:42 2012
***************
*** 138,142 ****
return;
}
! if (y < MIN_4TH_ROOT) {
*B_is_usable = 0;
if ((int)y==0) /* raise inexact flag */
--- 138,142 ----
return;
}
! if (!isnan(x) && y < MIN_4TH_ROOT) {
*B_is_usable = 0;
if ((int)y==0) /* raise inexact flag */
***************
*** 233,240 ****
if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
if (huge+x+y>one) { /* raise inexact flag */
! w = clog_for_large_values(z) + M_LN2;
if (sy == 0)
! return (cpackf(cimagf(w), -crealf(w)));
! return (cpackf(-cimagf(w), crealf(w)));
}
}
--- 233,242 ----
if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
if (huge+x+y>one) { /* raise inexact flag */
! w = clog_for_large_values(z);
! rx = fabsf(cimagf(w));
! ry = crealf(w) + M_LN2;
if (sy == 0)
! ry = -ry;
! return (cpackf(rx, ry));
}
}
***************
*** 290,294 ****
}
! if (ay > 0.5*FLT_MAX)
return (cpackf(logf(hypotf(x / M_E, y / M_E)) + 1, atan2f(y, x)));
--- 292,296 ----
}
! if (ax > 0.5*FLT_MAX)
return (cpackf(logf(hypotf(x / M_E, y / M_E)) + 1, atan2f(y, x)));
diff -c2 catrigl.c~ catrigl.c
*** catrigl.c~ Sun Aug 12 06:54:46 2012
--- catrigl.c Mon Aug 13 12:08:21 2012
***************
*** 119,123 ****
return;
}
! if (y < MIN_4TH_ROOT) {
*B_is_usable = 0;
if ((int)y==0) /* raise inexact flag */
--- 119,123 ----
return;
}
! if (!isnan(x) && y < MIN_4TH_ROOT) {
*B_is_usable = 0;
if ((int)y==0) /* raise inexact flag */
***************
*** 207,214 ****
if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
if (huge+x+y>one) { /* raise inexact flag */
! w = clog_for_large_values(z) + L_LN2;
if (sy == 0)
! return (cpackl(cimagl(w), -creall(w)));
! return (cpackl(-cimagl(w), creall(w)));
}
}
--- 207,216 ----
if (isfinite(x) && isfinite(y) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
if (huge+x+y>one) { /* raise inexact flag */
! w = clog_for_large_values(z);
! rx = fabsl(cimagl(w));
! ry = creall(w) + M_LN2;
if (sy == 0)
! ry = -ry;
! return (cpackl(rx, ry));
}
}
***************
*** 264,268 ****
}
! if (ay > 0.5*LDBL_MAX)
return (cpackl(logl(hypotl(x / L_E, y / L_E)) + 1, atan2l(y, x)));
--- 266,270 ----
}
! if (ax > 0.5*LDBL_MAX)
return (cpackl(logl(hypotl(x / L_E, y / L_E)) + 1, atan2l(y, x)));
More information about the freebsd-numerics
mailing list