The following reply was made to PR bin/170206; it has been noted by GNATS.
From: Stephen Montgomery-Smith <step...@missouri.edu>
To: Bruce Evans <b...@optusnet.com.au>
Cc: freebsd-bugs@freebsd.org, freebsd-gnats-sub...@freebsd.org,
Stephen Montgomery-Smith <step...@freebsd.org>
Subject: Re: bin/170206: complex arcsinh, log, etc.
Date: Sat, 28 Jul 2012 10:46:30 -0500
This is a multi-part message in MIME format.
--------------060304040300070501030603
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
OK. This clog really seems to work.
x*x + y*y - 1 is computed with a ULP less than 0.8. The rest of the
errors seem to be due to the implementation of log1p. The ULP of the
final answer seems to be never bigger than a little over 2.
--------------060304040300070501030603
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 10
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. We compute ax*ax + ay*ay - 1 using
* long multiplication in base 2^16.
*/
/*
* Split ax and ay into fixed point numbers with 160 bits after the
* "decimal" point.
*/
for (i=-1; i<NR_BLOCKS; i++) {
xx[i+1] = get_bits(ax,16*i);
yy[i+1] = get_bits(ay,16*i);
}
/*
* Long multiplication. We use 64 bit integers instead of 32 bit
* because we might get slightly bigger than 32 bit numbers due to
* the additions. (But probably 36 bit integers would be more than
* enough.)
*/
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];
}
/* Subtract 1. */
zz[0]--;
/* Handle the carries. */
carry = 0;
for (i=NR_BLOCKS-1; i>=-1; i--) {
zz[i+1] += carry;
carry = zz[i+1] >> 16;
zz[i+1] &= 0xffff;
}
/*
* If the number is negative, compute the 1's complement. (We
* should compute the 2's complement, but the the error will be
* negligable.)
*/
if ((zz[0] & 0x8000) != 0) {
sign = 1;
for (i=-1; i<NR_BLOCKS; i++)
zz[i+1] = 0xffff & (~zz[i+1]);
} else
sign = 0;
/* Convert fixed point into floating point. */
hm1 = 0;
for (i=NR_BLOCKS-1; i>=-1; 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);
}
--------------060304040300070501030603--
_______________________________________________
freebsd-bugs@freebsd.org mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-bugs
To unsubscribe, send any mail to "freebsd-bugs-unsubscr...@freebsd.org"
