Source-Changes-HG archive

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index][Old Index]

[src/trunk]: src/lib/libm Add the complex trig functions from FreeBSD



details:   https://anonhg.NetBSD.org/src/rev/1afcd5cb1a73
branches:  trunk
changeset: 818004:1afcd5cb1a73
user:      christos <christos%NetBSD.org@localhost>
date:      Mon Sep 19 22:05:05 2016 +0000

description:
Add the complex trig functions from FreeBSD

diffstat:

 lib/libm/complex/Makefile.inc |   23 +-
 lib/libm/complex/catrig.c     |  653 ++++++++++++++++++++++++++++++++++++++++++
 lib/libm/complex/catrigf.c    |  406 ++++++++++++++++++++++++++
 lib/libm/complex/catrigl.c    |  443 ++++++++++++++++++++++++++++
 lib/libm/src/math_private.h   |   53 +++-
 5 files changed, 1571 insertions(+), 7 deletions(-)

diffs (truncated from 1619 to 300 lines):

diff -r 17c8e3af909f -r 1afcd5cb1a73 lib/libm/complex/Makefile.inc
--- a/lib/libm/complex/Makefile.inc     Mon Sep 19 20:46:55 2016 +0000
+++ b/lib/libm/complex/Makefile.inc     Mon Sep 19 22:05:05 2016 +0000
@@ -1,14 +1,25 @@
-# $NetBSD: Makefile.inc,v 1.8 2014/10/10 12:43:07 christos Exp $
+# $NetBSD: Makefile.inc,v 1.9 2016/09/19 22:05:05 christos Exp $
 
 .PATH: ${.CURDIR}/complex
 
-COMPLEX_SRCS = cabs.c cacos.c cacosh.c carg.c casin.c casinh.c catan.c \
+COMPLEX_SRCS = cabs.c carg.c \
        ccos.c ccosh.c cephes_subr.c cexp.c clog.c conj.c cpow.c cproj.c \
-       cimag.c creal.c csin.c csinh.c csqrt.c ctan.c ctanh.c catanh.c
+       cimag.c creal.c csin.c csinh.c csqrt.c ctan.c ctanh.c \
+       catrig.c
+CATRIG_SRCS = cacos.c cacosh.c casin.c casinh.c catan.c catanh.c
 
+CPPFLAGS+=-I${.CURDIR}/src
 .for i in ${COMPLEX_SRCS}
 SRCS+= $i ${i:S/.c/f.c/} ${i:S/.c/l.c/}
-MAN+= ${i:Ncephes_*:S/.c/.3/}
-MLINKS+= ${i:Ncephes_*:S/.c/.3/} ${i:Ncephes_*:S/.c/f.3/}
-MLINKS+= ${i:Ncephes_*:S/.c/.3/} ${i:Ncephes_*:S/.c/l.3/}
+MAN+= ${i:Ncatrig*:Ncephes_*:S/.c/.3/}
+MLINKS+= ${i:Ncatrig*:Ncephes_*:S/.c/.3/} ${i:Ncatrig*:Ncephes_*:S/.c/f.3/}
+MLINKS+= ${i:Ncatrig*:Ncephes_*:S/.c/.3/} ${i:Ncatrig*:Ncephes_*:S/.c/l.3/}
 .endfor
+
+.for i in ${CATRIG_SRCS}
+MAN+= ${i:S/.c/.3/}
+MLINKS+= ${i:S/.c/.3/} ${i:S/.c/f.3/}
+MLINKS+= ${i:S/.c/.3/} ${i:S/.c/l.3/}
+.endfor
+
+
diff -r 17c8e3af909f -r 1afcd5cb1a73 lib/libm/complex/catrig.c
--- /dev/null   Thu Jan 01 00:00:00 1970 +0000
+++ b/lib/libm/complex/catrig.c Mon Sep 19 22:05:05 2016 +0000
@@ -0,0 +1,653 @@
+/*     $NetBSD: catrig.c,v 1.1 2016/09/19 22:05:05 christos Exp $      */
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <stephen%FreeBSD.ORG@localhost>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 275819 2014-12-16 09:21:56Z ed $");
+#endif
+__RCSID("$NetBSD: catrig.c,v 1.1 2016/09/19 22:05:05 christos Exp $");
+
+#include "namespace.h"
+#ifdef __weak_alias
+__weak_alias(casin, _casin)
+#endif
+#ifdef __weak_alias
+__weak_alias(catan, _catan)
+#endif
+
+#include <complex.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+
+
+#undef isinf
+#define isinf(x)       (fabs(x) == INFINITY)
+#undef isnan
+#define isnan(x)       ((x) != (x))
+#define        raise_inexact() do { volatile float junk __unused = /*LINTED*/1 + tiny; } while(/*CONSTCOND*/0)
+#undef signbit
+#define signbit(x)     (__builtin_signbit(x))
+
+/* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */
+static const double
+A_crossover =          10, /* Hull et al suggest 1.5, but 10 works better */
+B_crossover =          0.6417,                 /* suggested by Hull et al */
+FOUR_SQRT_MIN =                0x1p-509,               /* >= 4 * sqrt(DBL_MIN) */
+QUARTER_SQRT_MAX =     0x1p509,                /* <= sqrt(DBL_MAX) / 4 */
+m_e =                  2.7182818284590452e0,   /*  0x15bf0a8b145769.0p-51 */
+m_ln2 =                        6.9314718055994531e-1,  /*  0x162e42fefa39ef.0p-53 */
+pio2_hi =              1.5707963267948966e0,   /*  0x1921fb54442d18.0p-52 */
+RECIP_EPSILON =                1 / DBL_EPSILON,
+SQRT_3_EPSILON =       2.5809568279517849e-8,  /*  0x1bb67ae8584caa.0p-78 */
+SQRT_6_EPSILON =       3.6500241499888571e-8,  /*  0x13988e1409212e.0p-77 */
+SQRT_MIN =             0x1p-511;               /* >= sqrt(DBL_MIN) */
+
+static const volatile double
+pio2_lo =              6.1232339957367659e-17; /*  0x11a62633145c07.0p-106 */
+static const volatile float
+tiny =                 0x1p-100; 
+
+static double complex clog_for_large_values(double complex z);
+
+/*
+ * Testing indicates that all these functions are accurate up to 4 ULP.
+ * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh.
+ * The functions catan(h) are a little under 2 times slower than atanh.
+ *
+ * The code for casinh, casin, cacos, and cacosh comes first.  The code is
+ * rather complicated, and the four functions are highly interdependent.
+ *
+ * The code for catanh and catan comes at the end.  It is much simpler than
+ * the other functions, and the code for these can be disconnected from the
+ * rest of the code.
+ */
+
+/*
+ *                     ================================
+ *                     | casinh, casin, cacos, cacosh |
+ *                     ================================
+ */
+
+/*
+ * The algorithm is very close to that in "Implementing the complex arcsine
+ * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
+ * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
+ * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
+ * http://dl.acm.org/citation.cfm?id=275324.
+ *
+ * Throughout we use the convention z = x + I*y.
+ *
+ * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
+ * where
+ * A = (|z+I| + |z-I|) / 2
+ * B = (|z+I| - |z-I|) / 2 = y/A
+ *
+ * These formulas become numerically unstable:
+ *   (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
+ *       is, Re(casinh(z)) is close to 0);
+ *   (b) for Im(casinh(z)) when z is close to either of the intervals
+ *       [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
+ *       close to PI/2).
+ *
+ * These numerical problems are overcome by defining
+ * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
+ * Then if A < A_crossover, we use
+ *   log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
+ *   A-1 = f(x, 1+y) + f(x, 1-y)
+ * and if B > B_crossover, we use
+ *   asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
+ *   A-y = f(x, y+1) + f(x, y-1)
+ * where without loss of generality we have assumed that x and y are
+ * non-negative.
+ *
+ * Much of the difficulty comes because the intermediate computations may
+ * produce overflows or underflows.  This is dealt with in the paper by Hull
+ * et al by using exception handling.  We do this by detecting when
+ * computations risk underflow or overflow.  The hardest part is handling the
+ * underflows when computing f(a, b).
+ *
+ * Note that the function f(a, b) does not appear explicitly in the paper by
+ * Hull et al, but the idea may be found on pages 308 and 309.  Introducing the
+ * function f(a, b) allows us to concentrate many of the clever tricks in this
+ * paper into one function.
+ */
+
+/*
+ * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
+ * Pass hypot(a, b) as the third argument.
+ */
+static inline double
+f(double a, double b, double hypot_a_b)
+{
+       if (b < 0)
+               return ((hypot_a_b - b) / 2);
+       if (b == 0)
+               return (a / 2);
+       return (a * a / (hypot_a_b + b) / 2);
+}
+
+/*
+ * All the hard work is contained in this function.
+ * x and y are assumed positive or zero, and less than RECIP_EPSILON.
+ * Upon return:
+ * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
+ * B_is_usable is set to 1 if the value of B is usable.
+ * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
+ * If returning sqrt_A2my2 has potential to result in an underflow, it is
+ * rescaled, and new_y is similarly rescaled.
+ */
+static inline void
+do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
+    double *sqrt_A2my2, double *new_y)
+{
+       double R, S, A; /* A, B, R, and S are as in Hull et al. */
+       double Am1, Amy; /* A-1, A-y. */
+
+       R = hypot(x, y + 1);            /* |z+I| */
+       S = hypot(x, y - 1);            /* |z-I| */
+
+       /* A = (|z+I| + |z-I|) / 2 */
+       A = (R + S) / 2;
+       /*
+        * Mathematically A >= 1.  There is a small chance that this will not
+        * be so because of rounding errors.  So we will make certain it is
+        * so.
+        */
+       if (A < 1)
+               A = 1;
+
+       if (A < A_crossover) {
+               /*
+                * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
+                * rx = log1p(Am1 + sqrt(Am1*(A+1)))
+                */
+               if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
+                       /*
+                        * fp is of order x^2, and fm = x/2.
+                        * A = 1 (inexactly).
+                        */
+                       *rx = sqrt(x);
+               } else if (x >= DBL_EPSILON * fabs(y - 1)) {
+                       /*
+                        * Underflow will not occur because
+                        * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
+                        */
+                       Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
+                       *rx = log1p(Am1 + sqrt(Am1 * (A + 1)));
+               } else if (y < 1) {
+                       /*
+                        * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
+                        * A = 1 (inexactly).
+                        */
+                       *rx = x / sqrt((1 - y) * (1 + y));
+               } else {                /* if (y > 1) */
+                       /*
+                        * A-1 = y-1 (inexactly).
+                        */
+                       *rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));
+               }
+       } else {
+               *rx = log(A + sqrt(A * A - 1));
+       }
+
+       *new_y = y;
+
+       if (y < FOUR_SQRT_MIN) {
+               /*
+                * Avoid a possible underflow caused by y/A.  For casinh this
+                * would be legitimate, but will be picked up by invoking atan2
+                * later on.  For cacos this would not be legitimate.
+                */
+               *B_is_usable = 0;
+               *sqrt_A2my2 = A * (2 / DBL_EPSILON);
+               *new_y = y * (2 / DBL_EPSILON);
+               return;
+       }
+
+       /* B = (|z+I| - |z-I|) / 2 = y/A */
+       *B = y / A;
+       *B_is_usable = 1;
+
+       if (*B > B_crossover) {
+               *B_is_usable = 0;
+               /*
+                * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
+                * sqrt_A2my2 = sqrt(Amy*(A+y))
+                */
+               if (y == 1 && x < DBL_EPSILON / 128) {
+                       /*
+                        * fp is of order x^2, and fm = x/2.
+                        * A = 1 (inexactly).
+                        */
+                       *sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);
+               } else if (x >= DBL_EPSILON * fabs(y - 1)) {
+                       /*
+                        * Underflow will not occur because
+                        * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
+                        * and
+                        * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
+                        */
+                       Amy = f(x, y + 1, R) + f(x, y - 1, S);
+                       *sqrt_A2my2 = sqrt(Amy * (A + y));
+               } else if (y > 1) {
+                       /*
+                        * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and



Home | Main Index | Thread Index | Old Index