bin/170206: complex arcsinh, log, etc.

[Stephen Montgomery-Smith]

The following reply was made to PR bin/170206; it has been noted by GNATS.

From: Stephen Montgomery-Smith <stephen at missouri.edu>
To: Bruce Evans <brde at optusnet.com.au>
Cc: freebsd-bugs at freebsd.org, FreeBSD-gnats-submit at freebsd.org,
        Stephen Montgomery-Smith <stephen at 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--