From owner-freebsd-bugs@FreeBSD.ORG Fri Jul 27 02:50:09 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 866FF1065672 for ; Fri, 27 Jul 2012 02:50:09 +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 5AF3D8FC15 for ; Fri, 27 Jul 2012 02:50:09 +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 q6R2o9FR056085 for ; Fri, 27 Jul 2012 02:50:09 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.5/8.14.5/Submit) id q6R2o9Rn056084; Fri, 27 Jul 2012 02:50:09 GMT (envelope-from gnats) Resent-Date: Fri, 27 Jul 2012 02:50:09 GMT Resent-Message-Id: <201207270250.q6R2o9Rn056084@freefall.freebsd.org> Resent-From: FreeBSD-gnats-submit@FreeBSD.org (GNATS Filer) Resent-To: freebsd-bugs@FreeBSD.org Resent-Reply-To: FreeBSD-gnats-submit@FreeBSD.org, Stephen Montgomery-Smith Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 94722106564A for ; Fri, 27 Jul 2012 02:47:52 +0000 (UTC) (envelope-from stephen@wilberforce.math.missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id 6C2B58FC0C for ; Fri, 27 Jul 2012 02:47:52 +0000 (UTC) Received: from wilberforce.math.missouri.edu (localhost [127.0.0.1]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q6R2lkBj021155 for ; Thu, 26 Jul 2012 21:47:46 -0500 (CDT) (envelope-from stephen@wilberforce.math.missouri.edu) Received: (from stephen@localhost) by wilberforce.math.missouri.edu (8.14.5/8.14.5/Submit) id q6R2lkeR021134; Thu, 26 Jul 2012 21:47:46 -0500 (CDT) (envelope-from stephen) Message-Id: <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> Date: Thu, 26 Jul 2012 21:47:46 -0500 (CDT) From: Stephen Montgomery-Smith To: FreeBSD-gnats-submit@FreeBSD.org X-Send-Pr-Version: 3.113 Cc: Subject: 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: Fri, 27 Jul 2012 02:50:09 -0000 >Number: 170206 >Category: bin >Synopsis: complex arcsinh, log, etc. >Confidential: no >Severity: non-critical >Priority: low >Responsible: freebsd-bugs >State: open >Quarter: >Keywords: >Date-Required: >Class: change-request >Submitter-Id: current-users >Arrival-Date: Fri Jul 27 02:50:08 UTC 2012 >Closed-Date: >Last-Modified: >Originator: Stephen Montgomery-Smith >Release: FreeBSD 8.3-STABLE amd64 >Organization: >Environment: System: FreeBSD wilberforce 8.3-STABLE FreeBSD 8.3-STABLE #0: Tue Jun 12 20:29:45 CDT 2012 stephen@wilberforce:/usr/obj/usr/src/sys/GENERIC amd64 >Description: Implement casin(h), cacos(h), catan(h), clog functions. >How-To-Repeat: >Fix: These algorithms seem to have relative errors no worse than 4ULP for the arc-trig functions, and no worse than 5ULP for clog. # This is a shell archive. Save it in a file, remove anything before # this line, and then unpack it by entering "sh file". Note, it may # create directories; files and directories will be owned by you and # have default permissions. # # This archive contains: # # cplex.c # catrig.c # echo x - cplex.c sed 's/^X//' >cplex.c << 'ab710d5798014bb1e9398bb65fa5894f' X#include X#include X#include X#include X X#include "math_private.h" X X/* round down to 18 = 54/3 bits */ Xstatic double trim(double x) { X uint32_t hi; X X GET_HIGH_WORD(hi, x); X INSERT_WORDS(x, hi &0xfffffff8, 0); X return x; X} X Xdouble complex Xclog(double complex z) X{ X double x, y; X double ax, ay, x0, y0, x1, y1, x2, y2, t, hm1; X double val[12]; X int i, sorted; X X x = creal(z); X y = cimag(z); X X /* Handle NaNs using the general formula to mix them right. */ X if (x != x || y != y) X return (cpack(log(hypot(x, y)), atan2(y, x))); X X ax = fabs(x); X ay = fabs(y); X if (ax < ay) { X t = ax; X ax = ay; X ay = t; X } X X /* X * To avoid unnecessary overflow, if x or y are very large, divide x X * and y by M_E, and then add 1 to the logarithm. This depends on X * M_E being larger than sqrt(2). X * There is a potential loss of accuracy caused by dividing by M_E, X * but this case should happen extremely rarely. X */ X if (ay > 5e307) X return (cpack(log(hypot(x / M_E, y / M_E)) + 1, atan2(y, x))); X X if (ax == 1) { X if (ay < 1e-150) X return (cpack((ay * 0.5) * ay, atan2(y, x))); X return (cpack(log1p(ay * ay) * 0.5, atan2(y, x))); X } X X /* X * Because atan2 and hypot conform to C99, this also covers all the X * edge cases when x or y are 0 or infinite. X */ X if (ax < 1e-50 || ay < 1e-50 || ax > 1e50 || ay > 1e50) X return (cpack(log(hypot(x, y)), atan2(y, x))); X X /* X * From this point on, we don't need to worry about underflow or X * overflow in calculating ax*ax or ay*ay. X */ X X /* Some easy cases. */ X X if (ax >= 1) X return (cpack(log1p((ax-1)*(ax+1) + ay*ay) * 0.5, atan2(y, x))); X X if (ax*ax + ay*ay <= 0.7) X return (cpack(log(ax*ax + ay*ay) * 0.5, atan2(y, x))); X X /* X * Take extra care so that ULP of real part is small if hypot(x,y) is X * moderately close to 1. X */ X X x0 = trim(ax); X ax = ax-x0; X x1 = trim(ax); X x2 = ax-x1; X y0 = trim(ay); X ay = ay-y0; X y1 = trim(ay); X y2 = ay-y1; X X val[0] = x0*x0; X val[1] = y0*y0; X val[2] = 2*x0*x1; X val[3] = 2*y0*y1; X val[4] = x1*x1; X val[5] = y1*y1; X val[6] = 2*x0*x2; X val[7] = 2*y0*y2; X val[8] = 2*x1*x2; X val[9] = 2*y1*y2; X val[10] = x2*x2; X val[11] = y2*y2; X X /* Bubble sort. */ X do { X sorted = 1; X for (i=0;i<11;i++) { X if (val[i] < val[i+1]) { X sorted = 0; X t = val[i]; X val[i] = val[i+1]; X val[i+1] = t; X } X } X } while (!sorted); X X hm1 = -1; X for (i=0;i<12;i++) hm1 += val[i]; X return (cpack(0.5 * log1p(hm1), atan2(y, x))); X} X Xfloat complex Xclogf(float complex z) X{ X return clog(z); X} X ab710d5798014bb1e9398bb65fa5894f echo x - catrig.c sed 's/^X//' >catrig.c << 'e37ddaa44b334e25d827d6a69ee351aa' X#include X#include X#include X X#include "math_private.h" X X/* X * gcc doesn't implement complex multiplication or division correctly, X * so we need to handle infinities specially. We turn on this pragma to X * notify conforming c99 compilers that the fast-but-incorrect code that X * gcc generates is acceptable, since the special cases have already been X * handled. X */ X#pragma STDC CX_LIMITED_RANGE ON X Xcomplex double clog(complex double z); X Xstatic const double Xone = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ Xhuge= 1.00000000000000000000e+300; X X/* X * Testing indicates that all these functions are accurate up to 4 ULP. X */ X X/* X * The algorithm is very close to that in "Implementing the complex arcsine X * and arccosine functions using exception handling" by T. E. Hull, X * Thomas F. Fairgrieve, and Ping Tak Peter Tang, published in ACM X * Transactions on Mathematical Software, Volume 23 Issue 3, 1997, Pages X * 299-335, http://dl.acm.org/citation.cfm?id=275324 X * X * casinh(x+iy) = sign(x)*log(A+sqrt(A*A-1)) + sign(y)*I*asin(B) X * where X * A = 0.5(|z+I| + |z-I|) = f(x,1+y) + f(x,1-y) + 1 X * B = 0.5(|z+I| - |z-I|) X * z = x+I*y X * f(x,y) = 0.5*(hypot(x,y)-y) X * We also use X * asin(B) = atan2(sqrt(A*A-y*y),y) X * A-y = f(x,y+1) + f(x,y-1). X * X * Much of the difficulty comes because computing f(x,y) may produce X * underflows. X */ X X/* X * Returns 0.5*(hypot(x,y)-y). It assumes x is positive, and that y does X * not satisfy the inequalities 0 < fabs(y) < 1e-20. X * If reporting the answer risks an underflow, the underflow flag is set, X * and it returns 0.5*(hypot(x,y)-y)/x/x. X */ Xstatic double f(double x, double y, int *underflow) { X if (x==0) { X *underflow = 0; X if (y > 0) X return 0; X return -y; X } X if (y==0) { X *underflow = 0; X return 0.5*x; X } X if (x < 1e-100 && x < y) { X *underflow = 1; X return 0.5/(hypot(x,y)+y); X } X if (x < y) { X *underflow = 0; X return 0.5*x*x/(hypot(x,y)+y); X } X *underflow = 0; X return 0.5*(hypot(x,y)-y); X} X X/* X * All the hard work is contained in this function. X * Upon return: X * rx = Re(casinh(x+I*y)) X * B_good is set to 1 if the value of B is usable. X * If B_good is set to 0, A2my2 = A*A-y*y. X */ Xstatic void do_hard_work(double x, double y, double *rx, int *B_good, double *B, double *A2my2) X{ X double R, S, A, fp, fm; X int fpuf, fmuf; X X R = hypot(x,y+1); X S = hypot(x,y-1); X A = 0.5*(R + S); X X if (A < 10) { X fp = f(x,1+y,&fpuf); X fm = f(x,1-y,&fmuf); X if (fpuf == 1 && fmuf == 1) { X if (huge+x>one) /* set inexact flag. */ X *rx = log1p(x*sqrt((fp+fm)*(A+1))); X } else if (fmuf == 1) { X /* Overflow not possible because fp < 1e50 and x > 1e-100. X Underflow not possible because either fm=0 or fm X approximately bigger than 1e-200. */ X if (huge+x>one) /* set inexact flag. */ X *rx = log1p(fp+sqrt(x)*sqrt((fp/x+fm*x)*(A+1))); X } else if (fpuf == 1) { X /* Similar arguments against over/underflow. */ X if (huge+x>one) /* set inexact flag. */ X *rx = log1p(fm+sqrt(x)*sqrt((fm/x+fp*x)*(A+1))); X } else { X *rx = log1p(fp + fm + sqrt((fp+fm)*(A+1))); X } X } else X *rx = log(A + sqrt(A*A-1)); X X *B = y/A; /* = 0.5*(R - S) */ X *B_good = 1; X X if (*B > 0.5) { X *B_good = 0; X fp = f(x,y+1,&fpuf); X fm = f(x,y-1,&fmuf); X if (fpuf == 1 && fmuf == 1) X *A2my2 =x*sqrt((A+y)*(fp+fm)); X else if (fmuf == 1) X /* Overflow not possible because fp < 1e50 and x > 1e-100. X Underflow not possible because either fm=0 or fm X approximately bigger than 1e-200. */ X *A2my2 = sqrt(x)*sqrt((A+y)*(fp/x+fm*x)); X else if (fpuf == 1) X /* Similar arguments against over/underflow. */ X *A2my2 = sqrt(x)*sqrt((A+y)*(fm/x+fp*x)); X else X *A2my2 = sqrt((A+y)*(fp+fm)); X } X} X Xdouble complex Xcasinh(double complex z) X{ X double x, y, rx, ry, B, A2my2; X int sx, sy; X int B_good; X X x = creal(z); X y = cimag(z); X sx = signbit(x); X sy = signbit(y); X x = fabs(x); X y = fabs(y); X X if (cabs(z) > 1e20) { X if (huge+x>one) { /* set inexact flag. */ X if (sx == 0) return clog(2*z); X if (sx == 1) return -clog(-2*z); X } X } X X if (cabs(z) < 1e-20) X if (huge+x>one) /* set inexact flag. */ X return z; X X do_hard_work(x, y, &rx, &B_good, &B, &A2my2); X if (B_good) X ry = asin(B); X else X ry = atan2(y,A2my2); X X if (sx == 1) rx = -rx; X if (sy == 1) ry = -ry; X X return cpack(rx,ry); X} X X/* X * casin(z) = reverse(casinh(reverse(z))) X * where reverse(x+I*y) = y+x*I = I*conj(x+I*y). X */ X Xdouble complex Xcasin(double complex z) X{ X complex result; X X result = casinh(cpack(cimag(z),creal(z))); X return cpack(cimag(result),creal(result)); X} X X/* X * cacos(z) = PI/2 - casin(z) X * but do the computation carefully so cacos(z) is accurate when z is X * close to 1. X */ X Xdouble complex Xcacos(double complex z) X{ X double x, y, rx, ry, B, A2my2; X int sx, sy; X int B_good; X complex w; X X x = creal(z); X y = cimag(z); X sx = signbit(x); X sy = signbit(y); X x = fabs(x); X y = fabs(y); X X if (cabs(z) > 1e20) { X if (huge+x>one) { /* set inexact flag. */ X w = clog(2*z); X if (signbit(cimag(w)) == 0) X return cpack(cimag(w),-creal(w)); X return cpack(-cimag(w),creal(w)); X } X } X X if (cabs(z) < 1e-10) X if (huge+x>one) /* set inexact flag. */ X return cpack(M_PI_2-creal(z),-cimag(z)); X X do_hard_work(y, x, &ry, &B_good, &B, &A2my2); X if (B_good) { X if (sx==0) X rx = acos(B); X else X rx = acos(-B); X } else { X if (sx==0) X rx = atan2(A2my2,x); X else X rx = atan2(A2my2,-x); X } X X if (sy==0) ry = -ry; X X return cpack(rx,ry); X} X X/* X * cacosh(z) = I*cacos(z) or -I*cacos(z) X * where the sign is chosen so Re(cacosh(z)) >= 0. X */ X Xdouble complex Xcacosh(double complex z) X{ X complex double w; X X w = cacos(z); X if (signbit(cimag(w)) == 0) X return cpack(cimag(w),-creal(w)); X else X return cpack(-cimag(w),creal(w)); X} X X/* X * catanh(z) = 0.5 * log((z+1)/(z-1)) X * = 0.25 * log(|z+1|^2/|z-1|^2) + 0.5 * I * atan2(2y, (1-x*x-y*y)) X */ X Xdouble complex Xcatanh(double complex z) X{ X double x, y, rx, ry, hp, hm; X X x = creal(z); X y = cimag(z); X X if (cabs(z) < 1e-20) X if (huge+x>one) /* set inexact flag. */ X return z; X X if (cabs(z) > 1e20) X if (huge+x>one) { /* set inexact flag. */ X if (signbit(x) == 0) X return cpack(0,M_PI_2); X return cpack(0,-M_PI_2); X } X X if (fabs(y) < 1e-100) { X if (huge+x>one) { /* set inexact flag. */ X hp = (x+1)*(x+1); X hm = (x-1)*(x-1); X } X } else { X hp = (x+1)*(x+1)+y*y; /* |z+1|^2 */ X hm = (x-1)*(x-1)+y*y; /* |z-1|^2 */ X } X X if (hp < 0.5 || hm < 0.5) X rx = 0.25*(log(hp/hm)); X else if (x > 0) X rx = 0.25*log1p(4*x/hm); X else X rx = -0.25*log1p(-4*x/hp); X X if (x==1 || x==-1) { X if (signbit(y) == 0) X ry = atan2(2, -y)/2; X else X ry = atan2(-2, y)/2; X } else if (fabs(y) < 1e-100) { X if (huge+x>one) /* set inexact flag. */ X ry = atan2(2*y, (1-x)*(1+x))/2; X } else X ry = atan2(2*y, (1-x)*(1+x)-y*y)/2; X X return cpack(rx,ry); X} X X/* X * catan(z) = reverse(catanh(reverse(z))) X * where reverse(x+I*y) = y+x*I = I*conj(x+I*y). X */ X Xdouble complex Xcatan(double complex z) X{ X complex result; X X result = catanh(cpack(cimag(z),creal(z))); X return cpack(cimag(result),creal(result)); X} e37ddaa44b334e25d827d6a69ee351aa exit >Release-Note: >Audit-Trail: >Unformatted: