Date: Thu, 26 Jul 2012 21:47:46 -0500 (CDT) From: Stephen Montgomery-Smith <stephen@FreeBSD.org> To: FreeBSD-gnats-submit@FreeBSD.org Subject: bin/170206: complex arcsinh, log, etc. Message-ID: <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> Resent-Message-ID: <201207270250.q6R2o9Rn056084@freefall.freebsd.org>
next in thread | raw e-mail | index | archive | help
>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 <stdio.h>
X#include <complex.h>
X#include <float.h>
X#include <math.h>
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 <complex.h>
X#include <float.h>
X#include <math.h>
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:
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201207270247.q6R2lkeR021134>