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/105991c8d5eb
branches: trunk
changeset: 347862:105991c8d5eb
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 422c63be8d56 -r 105991c8d5eb 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 422c63be8d56 -r 105991c8d5eb 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