From owner-freebsd-bugs@FreeBSD.ORG Sat Jul 28 15:50:12 2012 Return-Path: Delivered-To: freebsd-bugs@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 30127106566B for ; Sat, 28 Jul 2012 15:50:12 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 1A10D8FC12 for ; Sat, 28 Jul 2012 15:50:12 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.5/8.14.5) with ESMTP id q6SFoBwg073015 for ; Sat, 28 Jul 2012 15:50:11 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.5/8.14.5/Submit) id q6SFoBdC073014; Sat, 28 Jul 2012 15:50:11 GMT (envelope-from gnats) Date: Sat, 28 Jul 2012 15:50:11 GMT Message-Id: <201207281550.q6SFoBdC073014@freefall.freebsd.org> To: freebsd-bugs@FreeBSD.org From: Stephen Montgomery-Smith Cc: Subject: Re: bin/170206: complex arcsinh, log, etc. X-BeenThere: freebsd-bugs@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: Stephen Montgomery-Smith List-Id: Bug reports List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 28 Jul 2012 15:50:12 -0000 The following reply was made to PR bin/170206; it has been noted by GNATS. From: Stephen Montgomery-Smith To: Bruce Evans Cc: freebsd-bugs@freebsd.org, FreeBSD-gnats-submit@freebsd.org, Stephen Montgomery-Smith 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 #include #include #include #include #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=-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=-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--