Fix LP64 libm for 128-bit long doubles
* reworked amd64/_fpmath.h and arm64/_fpmath.h to support 128-bit long
doubles.
* improved tests to cover long double cases
* made modfl an alias for LP32
Tests pass on x86, x86_64, arm, arm64 and mips.
Bug: 12921273
Change-Id: Ibe39acde57972447a8950fa45b1e702acc68ebeb
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acoshl.c b/libm/upstream-freebsd/lib/msun/src/e_acoshl.c
new file mode 100644
index 0000000..59faeb0
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/e_acoshl.c
@@ -0,0 +1,89 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_acosh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_acosh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */
+#if LDBL_MANT_DIG == 64
+#define EXP_LARGE 34
+#elif LDBL_MANT_DIG == 113
+#define EXP_LARGE 58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const double
+one = 1.0;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define ln2 u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+acoshl(long double x)
+{
+ long double t;
+ int16_t hx;
+
+ ENTERI();
+ GET_LDBL_EXPSIGN(hx, x);
+ if (hx < 0x3fff) { /* x < 1, or misnormal */
+ RETURNI((x-x)/(x-x));
+ } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */
+ if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */
+ RETURNI(x+x);
+ } else
+ RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */
+ } else if (hx == 0x3fff && x == 1) {
+ RETURNI(0.0); /* acosh(1) = 0 */
+ } else if (hx >= 0x4000) { /* LARGE > x >= 2, or misnormal */
+ t=x*x;
+ RETURNI(logl(2.0*x-one/(x+sqrtl(t-one))));
+ } else { /* 1<x<2 */
+ t = x-one;
+ RETURNI(log1pl(t+sqrtl(2.0*t+t*t)));
+ }
+}
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atanhl.c b/libm/upstream-freebsd/lib/msun/src/e_atanhl.c
new file mode 100644
index 0000000..a888426
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/e_atanhl.c
@@ -0,0 +1,74 @@
+/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_atanh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_atanh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define EXP_TINY -34
+#elif LDBL_MANT_DIG == 113
+#define EXP_TINY -58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const double one = 1.0, huge = 1e300;
+static const double zero = 0.0;
+
+long double
+atanhl(long double x)
+{
+ long double t;
+ uint16_t hx, ix;
+
+ ENTERI();
+ GET_LDBL_EXPSIGN(hx, x);
+ ix = hx & 0x7fff;
+ if (ix >= 0x3fff) /* |x| >= 1, or NaN or misnormal */
+ RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
+ if (ix < BIAS + EXP_TINY && (huge + x) > zero)
+ RETURNI(x); /* x is tiny */
+ SET_LDBL_EXPSIGN(x, ix);
+ if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */
+ t = x+x;
+ t = 0.5*log1pl(t+t*x/(one-x));
+ } else
+ t = 0.5*log1pl((x+x)/(one-x));
+ RETURNI((hx & 0x8000) == 0 ? t : -t);
+}
diff --git a/libm/upstream-freebsd/lib/msun/src/imprecise.c b/libm/upstream-freebsd/lib/msun/src/imprecise.c
new file mode 100644
index 0000000..a7503bf
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/imprecise.c
@@ -0,0 +1,69 @@
+/*-
+ * Copyright (c) 2013 David Chisnall
+ * 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.
+ *
+ * $FreeBSD$
+ */
+
+#include <float.h>
+#include <math.h>
+
+/*
+ * If long double is not the same size as double, then these will lose
+ * precision and we should emit a warning whenever something links against
+ * them.
+ */
+#if (LDBL_MANT_DIG > 53)
+#define WARN_IMPRECISE(x) \
+ __warn_references(x, # x " has lower than advertised precision");
+#else
+#define WARN_IMPRECISE(x)
+#endif
+/*
+ * Declare the functions as weak variants so that other libraries providing
+ * real versions can override them.
+ */
+#define DECLARE_WEAK(x)\
+ __weak_reference(imprecise_## x, x);\
+ WARN_IMPRECISE(x)
+
+long double
+imprecise_powl(long double x, long double y)
+{
+
+ return pow(x, y);
+}
+DECLARE_WEAK(powl);
+
+#define DECLARE_IMPRECISE(f) \
+ long double imprecise_ ## f ## l(long double v) { return f(v); }\
+ DECLARE_WEAK(f ## l)
+
+DECLARE_IMPRECISE(cosh);
+DECLARE_IMPRECISE(erfc);
+DECLARE_IMPRECISE(erf);
+DECLARE_IMPRECISE(lgamma);
+DECLARE_IMPRECISE(sinh);
+DECLARE_IMPRECISE(tanh);
+DECLARE_IMPRECISE(tgamma);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_asinhl.c b/libm/upstream-freebsd/lib/msun/src/s_asinhl.c
new file mode 100644
index 0000000..ba28f59
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/s_asinhl.c
@@ -0,0 +1,91 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)s_asinh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_asinh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
+/* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define EXP_LARGE 34
+#define EXP_TINY -34
+#elif LDBL_MANT_DIG == 113
+#define EXP_LARGE 58
+#define EXP_TINY -58
+#else
+#error "Unsupported long double format"
+#endif
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const double
+one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+huge= 1.00000000000000000000e+300;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define ln2 u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+asinhl(long double x)
+{
+ long double t, w;
+ uint16_t hx, ix;
+
+ ENTERI();
+ GET_LDBL_EXPSIGN(hx, x);
+ ix = hx & 0x7fff;
+ if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */
+ if (ix < BIAS + EXP_TINY) { /* |x| < TINY, or misnormal */
+ if (huge + x > one) RETURNI(x); /* return x inexact except 0 */
+ }
+ if (ix >= BIAS + EXP_LARGE) { /* |x| >= LARGE, or misnormal */
+ w = logl(fabsl(x))+ln2;
+ } else if (ix >= 0x4000) { /* LARGE > |x| >= 2.0, or misnormal */
+ t = fabsl(x);
+ w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
+ } else { /* 2.0 > |x| >= TINY, or misnormal */
+ t = x*x;
+ w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
+ }
+ RETURNI((hx & 0x8000) == 0 ? w : -w);
+}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_remquol.c b/libm/upstream-freebsd/lib/msun/src/s_remquol.c
index 712651c..8899293 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_remquol.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_remquol.c
@@ -5,7 +5,7 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
diff --git a/libm/upstream-freebsd/lib/msun/src/s_rintl.c b/libm/upstream-freebsd/lib/msun/src/s_rintl.c
index b43df89..e55b40e 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_rintl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_rintl.c
@@ -29,6 +29,7 @@
#include <float.h>
#include <math.h>
+#include <stdint.h>
#include "fpmath.h"