From owner-freebsd-bugs@FreeBSD.ORG Sat Jul 28 05:40:10 2012 Return-Path: Delivered-To: freebsd-bugs@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id A1433106564A for ; Sat, 28 Jul 2012 05:40:10 +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 8B6F78FC0A for ; Sat, 28 Jul 2012 05:40:10 +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 q6S5eAov074534 for ; Sat, 28 Jul 2012 05:40:10 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.5/8.14.5/Submit) id q6S5eAh2074533; Sat, 28 Jul 2012 05:40:10 GMT (envelope-from gnats) Date: Sat, 28 Jul 2012 05:40:10 GMT Message-Id: <201207280540.q6S5eAh2074533@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 05:40:10 -0000 The following reply was made to PR bin/170206; it has been noted by GNATS. From: Stephen Montgomery-Smith To: Bruce Evans Cc: Stephen Montgomery-Smith , 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 #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 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=-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