Date: Sat, 28 Jul 2012 05:40:10 GMT From: Stephen Montgomery-Smith <stephen@missouri.edu> To: freebsd-bugs@FreeBSD.org Subject: Re: bin/170206: complex arcsinh, log, etc. Message-ID: <201207280540.q6S5eAh2074533@freefall.freebsd.org>
next in thread | raw e-mail | index | archive | help
The following reply was made to PR bin/170206; it has been noted by GNATS.
From: Stephen Montgomery-Smith <stephen@missouri.edu>
To: Bruce Evans <brde@optusnet.com.au>
Cc: Stephen Montgomery-Smith <stephen@FreeBSD.org>,
FreeBSD-gnats-submit@FreeBSD.org, freebsd-bugs@FreeBSD.org
Subject: Re: bin/170206: complex arcsinh, log, etc.
Date: Sat, 28 Jul 2012 00:36:36 -0500
This is a multi-part message in MIME format.
--------------090903080003030709020200
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
Yes, everywhere I said "double precision" I meant "doubled precision."
I think the papers by Hull et al were perfectly happy with a ULP of
around 4.
I have been trying to do a little better, but like you I am noticing
that log1p isn't that good either.
I have tried some other things. I am attaching this example which gets
a ULP a little over 2. I simulate high precision arithmetic by
expanding everything out into integers. I certainly didn't aim for a
speedy program.
--------------090903080003030709020200
Content-Type: text/x-csrc;
name="cplex.c"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
filename="cplex.c"
#include <stdio.h>
#include <string.h>
#include <complex.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
/* Get binary digits -d through -d-16. Assume x > 0 */
uint32_t get_bits(double x, int d) {
uint32_t hi, lo;
int e;
if (x == 0) return 0;
e = d+ilogb(x)-4;
EXTRACT_WORDS(hi, lo, x);
hi &= 0x000fffff;
hi |= 0x00100000;
if (e <= -32) return 0;
if (e <= 0) {
hi >>= -e;
return hi & 0xffff;
}
if (e < 32) {
hi <<= e;
lo >>= (32-e);
return (hi | lo) & 0xffff;
}
if (e == 32)
return lo & 0xffff;
if (e <= 63) {
lo <<= (e-32);
return lo & 0xffff;
}
return 0;
}
#define NR_BLOCKS 8
double complex
clog(double complex z)
{
double x, y;
double ax, ay, t, hm1;
uint64_t xx[NR_BLOCKS+1], yy[NR_BLOCKS+1];
uint64_t zz[NR_BLOCKS+1];
uint64_t carry;
int sign;
int i, j;
x = creal(z);
y = cimag(z);
/* Handle NaNs using the general formula to mix them right. */
if (x != x || y != y)
return (cpack(log(hypot(x, y)), atan2(y, x)));
ax = fabs(x);
ay = fabs(y);
if (ax < ay) {
t = ax;
ax = ay;
ay = t;
}
/*
* To avoid unnecessary overflow, if x or y are very large, divide x
* and y by M_E, and then add 1 to the logarithm. This depends on
* M_E being larger than sqrt(2).
* There is a potential loss of accuracy caused by dividing by M_E,
* but this case should happen extremely rarely.
*/
if (ay > 5e307)
return (cpack(log(hypot(x / M_E, y / M_E)) + 1, atan2(y, x)));
if (ax == 1) {
if (ay < 1e-150)
return (cpack((ay * 0.5) * ay, atan2(y, x)));
return (cpack(log1p(ay * ay) * 0.5, atan2(y, x)));
}
/*
* Because atan2 and hypot conform to C99, this also covers all the
* edge cases when x or y are 0 or infinite.
*/
if (ax < 1e-50 || ay < 1e-50 || ax > 1e50 || ay > 1e50)
return (cpack(log(hypot(x, y)), atan2(y, x)));
/*
* From this point on, we don't need to worry about underflow or
* overflow in calculating ax*ax or ay*ay.
*/
/* Some easy cases. */
if (ax*ax + ay*ay <= 0.1 || ax*ax + ay*ay >= 10)
return (cpack(log(ax*ax + ay*ay) * 0.5, atan2(y, x)));
/*
* Take extra care so that ULP of real part is small if hypot(x,y) is
* moderately close to 1.
*/
for (i=-1; i<NR_BLOCKS; i++) {
xx[i+1] = get_bits(ax,16*i);
yy[i+1] = get_bits(ay,16*i);
}
memset(zz,0,sizeof(zz));
for (i=-1; i<NR_BLOCKS; i++)
for (j=-1; j<NR_BLOCKS && i+j+1 < NR_BLOCKS; j++) {
zz[i+j+2] += xx[i+1]*xx[j+1];
zz[i+j+2] += yy[i+1]*yy[j+1];
}
zz[0]--;
carry = 0;
for (i=NR_BLOCKS-1; i>=-1; i--) {
zz[i+1] += carry;
carry = zz[i+1] >> 16;
zz[i+1] &= 0xffff;
}
if ((zz[0] & 0x8000) != 0) {
sign = 1;
for (i=-1; i<NR_BLOCKS; i++)
zz[i+1] = 0xffff & (~zz[i+1]);
} else
sign = 0;
hm1 = 0;
for (i=-1; i<NR_BLOCKS; i++)
hm1 += zz[i+1] * exp2(16*(-1-i));
if (sign == 1) hm1 = -hm1;
return (cpack(0.5 * log1p(hm1), atan2(y, x)));
}
float complex
clogf(float complex z)
{
return clog(z);
}
--------------090903080003030709020200--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201207280540.q6S5eAh2074533>