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/Android.mk b/libm/Android.mk
index d06088b..4a826eb 100644
--- a/libm/Android.mk
+++ b/libm/Android.mk
@@ -179,42 +179,55 @@
libm_common_src_files += \
fake_long_double.c \
- upstream-freebsd/lib/msun/src/s_modfl.c \
-# TODO: on Android, "long double" is "double".
-# upstream-freebsd/lib/msun/src/e_acosl.c \
-# upstream-freebsd/lib/msun/src/e_asinl.c \
-# upstream-freebsd/lib/msun/src/e_atan2l.c \
-# upstream-freebsd/lib/msun/src/e_fmodl.c \
-# upstream-freebsd/lib/msun/src/e_hypotl.c \
-# upstream-freebsd/lib/msun/src/e_remainderl.c \
-# upstream-freebsd/lib/msun/src/e_sqrtl.c \
-# upstream-freebsd/lib/msun/src/s_atanl.c \
-# upstream-freebsd/lib/msun/src/s_cbrtl.c \
-# upstream-freebsd/lib/msun/src/s_ceill.c \
-# upstream-freebsd/lib/msun/src/s_copysignl.c \
-# upstream-freebsd/lib/msun/src/s_cosl.c \
-# upstream-freebsd/lib/msun/src/s_fabsl.c \
-# upstream-freebsd/lib/msun/src/s_floorl.c \
-# upstream-freebsd/lib/msun/src/s_fmal.c \
-# upstream-freebsd/lib/msun/src/s_fmaxl.c \
-# upstream-freebsd/lib/msun/src/s_fminl.c \
-# upstream-freebsd/lib/msun/src/s_frexpl.c \
-# upstream-freebsd/lib/msun/src/s_ilogbl.c \
-# upstream-freebsd/lib/msun/src/s_llrintl.c \
-# upstream-freebsd/lib/msun/src/s_llroundl.c \
-# upstream-freebsd/lib/msun/src/s_logbl.c \
-# upstream-freebsd/lib/msun/src/s_lrintl.c \
-# upstream-freebsd/lib/msun/src/s_lroundl.c \
-# upstream-freebsd/lib/msun/src/s_nextafterl.c \
-# upstream-freebsd/lib/msun/src/s_nexttoward.c \
-# upstream-freebsd/lib/msun/src/s_remquol.c \
-# upstream-freebsd/lib/msun/src/s_rintl.c \
-# upstream-freebsd/lib/msun/src/s_roundl.c \
-# upstream-freebsd/lib/msun/src/s_scalbnl.c \
-# upstream-freebsd/lib/msun/src/s_sinl.c \
-# upstream-freebsd/lib/msun/src/s_tanl.c \
-# upstream-freebsd/lib/msun/src/s_truncl.c \
+libm_ld_src_files = \
+ upstream-freebsd/lib/msun/src/e_acosl.c \
+ upstream-freebsd/lib/msun/src/e_acoshl.c \
+ upstream-freebsd/lib/msun/src/e_asinl.c \
+ upstream-freebsd/lib/msun/src/e_atan2l.c \
+ upstream-freebsd/lib/msun/src/e_atanhl.c \
+ upstream-freebsd/lib/msun/src/e_fmodl.c \
+ upstream-freebsd/lib/msun/src/e_hypotl.c \
+ upstream-freebsd/lib/msun/src/e_remainderl.c \
+ upstream-freebsd/lib/msun/src/e_sqrtl.c \
+ upstream-freebsd/lib/msun/src/s_asinhl.c \
+ upstream-freebsd/lib/msun/src/s_atanl.c \
+ upstream-freebsd/lib/msun/src/s_cbrtl.c \
+ upstream-freebsd/lib/msun/src/s_ceill.c \
+ upstream-freebsd/lib/msun/src/s_copysignl.c \
+ upstream-freebsd/lib/msun/src/s_cosl.c \
+ upstream-freebsd/lib/msun/src/s_fabsl.c \
+ upstream-freebsd/lib/msun/src/s_floorl.c \
+ upstream-freebsd/lib/msun/src/s_fmal.c \
+ upstream-freebsd/lib/msun/src/s_fmaxl.c \
+ upstream-freebsd/lib/msun/src/s_fminl.c \
+ upstream-freebsd/lib/msun/src/s_modfl.c \
+ upstream-freebsd/lib/msun/src/s_frexpl.c \
+ upstream-freebsd/lib/msun/src/s_ilogbl.c \
+ upstream-freebsd/lib/msun/src/s_llrintl.c \
+ upstream-freebsd/lib/msun/src/s_llroundl.c \
+ upstream-freebsd/lib/msun/src/s_logbl.c \
+ upstream-freebsd/lib/msun/src/s_lrintl.c \
+ upstream-freebsd/lib/msun/src/s_lroundl.c \
+ upstream-freebsd/lib/msun/src/s_nextafterl.c \
+ upstream-freebsd/lib/msun/src/s_nexttoward.c \
+ upstream-freebsd/lib/msun/src/s_remquol.c \
+ upstream-freebsd/lib/msun/src/s_rintl.c \
+ upstream-freebsd/lib/msun/src/s_roundl.c \
+ upstream-freebsd/lib/msun/src/s_scalbnl.c \
+ upstream-freebsd/lib/msun/src/s_sinl.c \
+ upstream-freebsd/lib/msun/src/s_tanl.c \
+ upstream-freebsd/lib/msun/src/s_truncl.c \
+
+libm_ld_src_files += \
+ upstream-freebsd/lib/msun/ld128/invtrig.c \
+ upstream-freebsd/lib/msun/ld128/k_cosl.c \
+ upstream-freebsd/lib/msun/ld128/k_sinl.c \
+ upstream-freebsd/lib/msun/ld128/k_tanl.c \
+ upstream-freebsd/lib/msun/ld128/s_exp2l.c \
+ upstream-freebsd/lib/msun/ld128/s_expl.c \
+ upstream-freebsd/lib/msun/ld128/s_logl.c \
+ upstream-freebsd/lib/msun/ld128/s_nanl.c \
# TODO: re-enable i387/e_sqrtf.S for x86, and maybe others.
@@ -224,6 +237,7 @@
-include $(LOCAL_PATH)/freebsd-compat.h \
libm_common_includes := $(LOCAL_PATH)/upstream-freebsd/lib/msun/src/
+libm_ld_includes := $(LOCAL_PATH)/upstream-freebsd/lib/msun/ld128/
#
# libm.a for target.
@@ -241,22 +255,23 @@
LOCAL_C_INCLUDES_arm := $(LOCAL_PATH)/arm
LOCAL_SRC_FILES_arm := arm/fenv.c
-LOCAL_C_INCLUDES_arm64 := $(LOCAL_PATH)/arm64
-LOCAL_SRC_FILES_arm64 := arm64/fenv.c
+LOCAL_C_INCLUDES_arm64 := $(LOCAL_PATH)/arm64 $(libm_ld_includes)
+LOCAL_SRC_FILES_arm64 := arm64/fenv.c $(libm_ld_src_files)
LOCAL_C_INCLUDES_x86 := $(LOCAL_PATH)/i386 $(LOCAL_PATH)/i387
LOCAL_SRC_FILES_x86 := i387/fenv.c
-LOCAL_C_INCLUDES_x86_64 := $(LOCAL_PATH)/amd64
-LOCAL_SRC_FILES_x86_64 := amd64/fenv.c
+LOCAL_CFLAGS_x86_64 := -include $(LOCAL_PATH)/fpmath.h
+LOCAL_C_INCLUDES_x86_64 := $(LOCAL_PATH)/amd64 $(libm_ld_includes)
+LOCAL_SRC_FILES_x86_64 := amd64/fenv.c $(libm_ld_src_files)
LOCAL_CFLAGS_mips := -fno-builtin-rintf -fno-builtin-rint
LOCAL_C_INCLUDES_mips := $(LOCAL_PATH)/mips
LOCAL_SRC_FILES_mips := mips/fenv.c
LOCAL_CFLAGS_mips64 := -fno-builtin-rintf -fno-builtin-rint
-LOCAL_C_INCLUDES_mips64 := $(LOCAL_PATH)/mips
-LOCAL_SRC_FILES_mips64 := mips/fenv.c
+LOCAL_C_INCLUDES_mips64 := $(LOCAL_PATH)/mips $(libm_ld_includes)
+LOCAL_SRC_FILES_mips64 := mips/fenv.c $(libm_ld_src_files)
include $(BUILD_STATIC_LIBRARY)
diff --git a/libm/amd64/_fpmath.h b/libm/amd64/_fpmath.h
index c2a7384..f8d5c85 100755
--- a/libm/amd64/_fpmath.h
+++ b/libm/amd64/_fpmath.h
@@ -26,30 +26,35 @@
* $FreeBSD$
*/
+// ANDROID changed
+// Android uses 128 bits long doubles for LP64, so the structure and macros
+// were reworked for the quad precision ieee representation.
+
union IEEEl2bits {
long double e;
struct {
- unsigned int manl :32;
- unsigned int manh :32;
- unsigned int exp :15;
+ unsigned long manl :64;
+ unsigned long manh :48;
+ unsigned int exp :15;
unsigned int sign :1;
- unsigned int junkl :16;
- unsigned int junkh :32;
} bits;
struct {
- unsigned long man :64;
+ unsigned long manl :64;
+ unsigned long manh :48;
unsigned int expsign :16;
- unsigned long junk :48;
} xbits;
};
-#define LDBL_NBIT 0x80000000
-#define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
+#define LDBL_NBIT 0
+#define LDBL_IMPLICIT_NBIT
+#define mask_nbit_l(u) ((void)0)
-#define LDBL_MANH_SIZE 32
-#define LDBL_MANL_SIZE 32
+#define LDBL_MANH_SIZE 48
+#define LDBL_MANL_SIZE 64
#define LDBL_TO_ARRAY32(u, a) do { \
(a)[0] = (uint32_t)(u).bits.manl; \
- (a)[1] = (uint32_t)(u).bits.manh; \
-} while (0)
+ (a)[1] = (uint32_t)((u).bits.manl >> 32); \
+ (a)[2] = (uint32_t)(u).bits.manh; \
+ (a)[3] = (uint32_t)((u).bits.manh >> 32); \
+} while(0)
diff --git a/libm/arm64/_fpmath.h b/libm/arm64/_fpmath.h
index a24632a..9f46640 100644
--- a/libm/arm64/_fpmath.h
+++ b/libm/arm64/_fpmath.h
@@ -26,31 +26,48 @@
* $FreeBSD: src/lib/libc/aarch64/_fpmath.h $
*/
+// ANDROID changed
+// Android uses 128 bits long doubles for LP64, so the structure and the macros
+// were reworked for the quad precision ieee representation.
+
union IEEEl2bits {
long double e;
struct {
#ifndef __AARCH64EB__
- unsigned int manl :32;
- unsigned int manh :20;
- unsigned int exp :11;
+ unsigned long manl :64;
+ unsigned long manh :48;
+ unsigned int exp :15;
unsigned int sign :1;
#else
- unsigned int sign :1;
- unsigned int exp :11;
- unsigned int manh :20;
- unsigned int manl :32;
+ unsigned int sign :1;
+ unsigned int exp :15;
+ unsigned long manh :48;
+ unsigned long manl :64;
#endif
} bits;
+ struct {
+#ifndef __AARCH64EB__
+ unsigned long manl :64;
+ unsigned long manh :48;
+ unsigned int expsign :16;
+#else
+ unsigned int expsign :16;
+ unsigned long manh :48;
+ unsigned long manl :64;
+#endif
+ } xbits;
};
#define LDBL_NBIT 0
#define LDBL_IMPLICIT_NBIT
#define mask_nbit_l(u) ((void)0)
-#define LDBL_MANH_SIZE 32
-#define LDBL_MANL_SIZE 32
+#define LDBL_MANH_SIZE 48
+#define LDBL_MANL_SIZE 64
#define LDBL_TO_ARRAY32(u, a) do { \
(a)[0] = (uint32_t)(u).bits.manl; \
- (a)[1] = (uint32_t)(u).bits.manh; \
+ (a)[1] = (uint32_t)((u).bits.manl >> 32); \
+ (a)[2] = (uint32_t)(u).bits.manh; \
+ (a)[3] = (uint32_t)((u).bits.manh >> 32); \
} while(0)
diff --git a/libm/fake_long_double.c b/libm/fake_long_double.c
index 756e35b..13adf2f 100644
--- a/libm/fake_long_double.c
+++ b/libm/fake_long_double.c
@@ -22,6 +22,18 @@
int (isnanf)(float a1) { return __isnanf(a1); }
+// FreeBSD falls back to the double variants of these functions as well.
+long double coshl(long double a1) { return cosh(a1); }
+long double erfcl(long double a1) { return erfc(a1); }
+long double erfl(long double a1) { return erf(a1); }
+long double lgammal(long double a1) { return lgamma(a1); }
+long double powl(long double a1, long double a2) { return pow(a1, a2); }
+long double sinhl(long double a1) { return sinh(a1); }
+long double tanhl(long double a1) { return tanh(a1); }
+long double tgammal(long double a1) { return tgamma(a1); }
+long double significandl(long double a1) { return significand(a1); }
+
+#ifndef __LP64__
/*
* The BSD "long double" functions are broken when sizeof(long double) == sizeof(double).
* Android works around those cases by replacing the broken functions with our own trivial stubs
@@ -45,10 +57,9 @@
long double cbrtl(long double a1) { return cbrt(a1); }
long double ceill(long double a1) { return ceil(a1); }
long double copysignl(long double a1, long double a2) { return copysign(a1, a2); }
-long double coshl(long double a1) { return cosh(a1); }
+
long double cosl(long double a1) { return cos(a1); }
-long double erfcl(long double a1) { return erfc(a1); }
-long double erfl(long double a1) { return erf(a1); }
+
long double exp2l(long double a1) { return exp2(a1); }
long double expl(long double a1) { return exp(a1); }
long double expm1l(long double a1) { return expm1(a1); }
@@ -62,7 +73,6 @@
long double hypotl(long double a1, long double a2) { return hypot(a1, a2); }
int ilogbl(long double a1) { return ilogb(a1); }
long double ldexpl(long double a1, int a2) { return ldexp(a1, a2); }
-long double lgammal(long double a1) { return lgamma(a1); }
long long llrintl(long double a1) { return llrint(a1); }
long double log10l(long double a1) { return log10(a1); }
long double log1pl(long double a1) { return log1p(a1); }
@@ -72,19 +82,17 @@
long lrintl(long double a1) { return lrint(a1); }
long long llroundl(long double a1) { return llround(a1); }
long lroundl(long double a1) { return lround(a1); }
+long double modfl(long double a1, long double* a2) { double i; double f = modf(a1, &i); *a2 = i; return f; }
long double nanl(const char* a1) { return nan(a1); }
long double nextafterl(long double a1, long double a2) { return nextafter(a1, a2); }
-long double powl(long double a1, long double a2) { return pow(a1, a2); }
long double remainderl(long double a1, long double a2) { return remainder(a1, a2); }
long double remquol(long double a1, long double a2, int* a3) { return remquo(a1, a2, a3); }
long double rintl(long double a1) { return rint(a1); }
long double roundl(long double a1) { return round(a1); }
long double scalbnl(long double a1, int a2) { return scalbn(a1, a2); }
-long double significandl(long double a1) { return significand(a1); }
-long double sinhl(long double a1) { return sinh(a1); }
long double sinl(long double a1) { return sin(a1); }
long double sqrtl(long double a1) { return sqrt(a1); }
-long double tanhl(long double a1) { return tanh(a1); }
long double tanl(long double a1) { return tan(a1); }
-long double tgammal(long double a1) { return tgamma(a1); }
long double truncl(long double a1) { return trunc(a1); }
+
+#endif // __LP64__
diff --git a/libm/fpmath.h b/libm/fpmath.h
index f3753c1..a7c612c 100644
--- a/libm/fpmath.h
+++ b/libm/fpmath.h
@@ -27,6 +27,9 @@
* $FreeBSD: src/lib/libc/include/fpmath.h,v 1.3 2005/02/06 03:23:31 das Exp $
*/
+#ifndef _FPMATH_
+#define _FPMATH_
+
#include <endian.h>
#include "_fpmath.h"
@@ -76,9 +79,15 @@
* Android works around those cases by replacing the broken functions with our own trivial stubs
* that call the regular "double" function.
*/
+#ifndef __LP64__
+
#define __fpclassifyl __broken__fpclassify
#define __isfinitel __broken__isfinitel
#define __isinfl __broken__isinfl
#define __isnanl __broken__isnanl
#define __isnormall __broken__isnormall
#define __signbitl __broken_signbitl
+
+#endif // __LP64__
+
+#endif // _FPMATH_
diff --git a/libm/upstream-freebsd/lib/msun/ld128/e_rem_pio2l.h b/libm/upstream-freebsd/lib/msun/ld128/e_rem_pio2l.h
new file mode 100644
index 0000000..078d0c3
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/e_rem_pio2l.h
@@ -0,0 +1,140 @@
+/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* ld128 version of __ieee754_rem_pio2l(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+/*
+ * XXX need to verify that nonzero integer multiples of pi/2 within the
+ * range get no closer to a long double than 2**-140, or that
+ * ilogb(x) + ilogb(min_delta) < 45 - -140.
+ */
+/*
+ * invpio2: 113 bits of 2/pi
+ * pio2_1: first 68 bits of pi/2
+ * pio2_1t: pi/2 - pio2_1
+ * pio2_2: second 68 bits of pi/2
+ * pio2_2t: pi/2 - (pio2_1+pio2_2)
+ * pio2_3: third 68 bits of pi/2
+ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 = 1.67772160000000000000e+07; /* 0x41700000, 0x00000000 */
+
+static const long double
+invpio2 = 6.3661977236758134307553505349005747e-01L, /* 0x145f306dc9c882a53f84eafa3ea6a.0p-113 */
+pio2_1 = 1.5707963267948966192292994253909555e+00L, /* 0x1921fb54442d18469800000000000.0p-112 */
+pio2_1t = 2.0222662487959507323996846200947577e-21L, /* 0x13198a2e03707344a4093822299f3.0p-181 */
+pio2_2 = 2.0222662487959507323994779168837751e-21L, /* 0x13198a2e03707344a400000000000.0p-181 */
+pio2_2t = 2.0670321098263988236496903051604844e-43L, /* 0x127044533e63a0105df531d89cd91.0p-254 */
+pio2_3 = 2.0670321098263988236499468110329591e-43L, /* 0x127044533e63a0105e00000000000.0p-254 */
+pio2_3t = -2.5650587247459238361625433492959285e-65L; /* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */
+
+static inline __always_inline int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+ union IEEEl2bits u,u1;
+ long double z,w,t,r,fn;
+ double tx[5],ty[3];
+ int64_t n;
+ int e0,ex,i,j,nx;
+ int16_t expsign;
+
+ u.e = x;
+ expsign = u.xbits.expsign;
+ ex = expsign & 0x7fff;
+ if (ex < BIAS + 45 || ex == BIAS + 45 &&
+ u.bits.manh < 0x921fb54442d1LL) {
+ /* |x| ~< 2^45*(pi/2), medium size */
+ /* Use a specialized rint() to get fn. Assume round-to-nearest. */
+ fn = x*invpio2+0x1.8p112;
+ fn = fn-0x1.8p112;
+#ifdef HAVE_EFFICIENT_I64RINT
+ n = i64rint(fn);
+#else
+ n = fn;
+#endif
+ r = x-fn*pio2_1;
+ w = fn*pio2_1t; /* 1st round good to 180 bit */
+ {
+ union IEEEl2bits u2;
+ int ex1;
+ j = ex;
+ y[0] = r-w;
+ u2.e = y[0];
+ ex1 = u2.xbits.expsign & 0x7fff;
+ i = j-ex1;
+ if(i>51) { /* 2nd iteration needed, good to 248 */
+ t = r;
+ w = fn*pio2_2;
+ r = t-w;
+ w = fn*pio2_2t-((t-r)-w);
+ y[0] = r-w;
+ u2.e = y[0];
+ ex1 = u2.xbits.expsign & 0x7fff;
+ i = j-ex1;
+ if(i>119) { /* 3rd iteration need, 316 bits acc */
+ t = r; /* will cover all possible cases */
+ w = fn*pio2_3;
+ r = t-w;
+ w = fn*pio2_3t-((t-r)-w);
+ y[0] = r-w;
+ }
+ }
+ }
+ y[1] = (r-y[0])-w;
+ return n;
+ }
+ /*
+ * all other (large) arguments
+ */
+ if(ex==0x7fff) { /* x is inf or NaN */
+ y[0]=y[1]=x-x; return 0;
+ }
+ /* set z = scalbn(|x|,ilogb(x)-23) */
+ u1.e = x;
+ e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
+ u1.xbits.expsign = ex - e0;
+ z = u1.e;
+ for(i=0;i<4;i++) {
+ tx[i] = (double)((int32_t)(z));
+ z = (z-tx[i])*two24;
+ }
+ tx[4] = z;
+ nx = 5;
+ while(tx[nx-1]==zero) nx--; /* skip zero term */
+ n = __kernel_rem_pio2(tx,ty,e0,nx,3);
+ t = (long double)ty[2] + ty[1];
+ r = t + ty[0];
+ w = ty[0] - (r - t);
+ if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+ y[0] = r; y[1] = w; return n;
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/invtrig.c b/libm/upstream-freebsd/lib/msun/ld128/invtrig.c
new file mode 100644
index 0000000..df67d16
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/invtrig.c
@@ -0,0 +1,100 @@
+/*-
+ * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * 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>
+__FBSDID("$FreeBSD$");
+
+#include "invtrig.h"
+
+/*
+ * asinl() and acosl()
+ */
+const long double
+pS0 = 1.66666666666666666666666666666700314e-01L,
+pS1 = -7.32816946414566252574527475428622708e-01L,
+pS2 = 1.34215708714992334609030036562143589e+00L,
+pS3 = -1.32483151677116409805070261790752040e+00L,
+pS4 = 7.61206183613632558824485341162121989e-01L,
+pS5 = -2.56165783329023486777386833928147375e-01L,
+pS6 = 4.80718586374448793411019434585413855e-02L,
+pS7 = -4.42523267167024279410230886239774718e-03L,
+pS8 = 1.44551535183911458253205638280410064e-04L,
+pS9 = -2.10558957916600254061591040482706179e-07L,
+qS1 = -4.84690167848739751544716485245697428e+00L,
+qS2 = 9.96619113536172610135016921140206980e+00L,
+qS3 = -1.13177895428973036660836798461641458e+01L,
+qS4 = 7.74004374389488266169304117714658761e+00L,
+qS5 = -3.25871986053534084709023539900339905e+00L,
+qS6 = 8.27830318881232209752469022352928864e-01L,
+qS7 = -1.18768052702942805423330715206348004e-01L,
+qS8 = 8.32600764660522313269101537926539470e-03L,
+qS9 = -1.99407384882605586705979504567947007e-04L;
+
+/*
+ * atanl()
+ */
+const long double atanhi[] = {
+ 4.63647609000806116214256231461214397e-01L,
+ 7.85398163397448309615660845819875699e-01L,
+ 9.82793723247329067985710611014666038e-01L,
+ 1.57079632679489661923132169163975140e+00L,
+};
+
+const long double atanlo[] = {
+ 4.89509642257333492668618435220297706e-36L,
+ 2.16795253253094525619926100651083806e-35L,
+ -2.31288434538183565909319952098066272e-35L,
+ 4.33590506506189051239852201302167613e-35L,
+};
+
+const long double aT[] = {
+ 3.33333333333333333333333333333333125e-01L,
+ -1.99999999999999999999999999999180430e-01L,
+ 1.42857142857142857142857142125269827e-01L,
+ -1.11111111111111111111110834490810169e-01L,
+ 9.09090909090909090908522355708623681e-02L,
+ -7.69230769230769230696553844935357021e-02L,
+ 6.66666666666666660390096773046256096e-02L,
+ -5.88235294117646671706582985209643694e-02L,
+ 5.26315789473666478515847092020327506e-02L,
+ -4.76190476189855517021024424991436144e-02L,
+ 4.34782608678695085948531993458097026e-02L,
+ -3.99999999632663469330634215991142368e-02L,
+ 3.70370363987423702891250829918659723e-02L,
+ -3.44827496515048090726669907612335954e-02L,
+ 3.22579620681420149871973710852268528e-02L,
+ -3.03020767654269261041647570626778067e-02L,
+ 2.85641979882534783223403715930946138e-02L,
+ -2.69824879726738568189929461383741323e-02L,
+ 2.54194698498808542954187110873675769e-02L,
+ -2.35083879708189059926183138130183215e-02L,
+ 2.04832358998165364349957325067131428e-02L,
+ -1.54489555488544397858507248612362957e-02L,
+ 8.64492360989278761493037861575248038e-03L,
+ -2.58521121597609872727919154569765469e-03L,
+};
+
+const long double pi_lo = 8.67181013012378102479704402604335225e-35L;
diff --git a/libm/upstream-freebsd/lib/msun/ld128/invtrig.h b/libm/upstream-freebsd/lib/msun/ld128/invtrig.h
new file mode 100644
index 0000000..12f598b
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/invtrig.h
@@ -0,0 +1,113 @@
+/*-
+ * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * 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 "fpmath.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+#define MANH_SIZE (LDBL_MANH_SIZE + 1)
+
+/* Approximation thresholds. */
+#define ASIN_LINEAR (BIAS - 56) /* 2**-56 */
+#define ACOS_CONST (BIAS - 113) /* 2**-113 */
+#define ATAN_CONST (BIAS + 113) /* 2**113 */
+#define ATAN_LINEAR (BIAS - 56) /* 2**-56 */
+
+/* 0.95 */
+#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
+
+/* Constants shared by the long double inverse trig functions. */
+#define pS0 _ItL_pS0
+#define pS1 _ItL_pS1
+#define pS2 _ItL_pS2
+#define pS3 _ItL_pS3
+#define pS4 _ItL_pS4
+#define pS5 _ItL_pS5
+#define pS6 _ItL_pS6
+#define pS7 _ItL_pS7
+#define pS8 _ItL_pS8
+#define pS9 _ItL_pS9
+#define qS1 _ItL_qS1
+#define qS2 _ItL_qS2
+#define qS3 _ItL_qS3
+#define qS4 _ItL_qS4
+#define qS5 _ItL_qS5
+#define qS6 _ItL_qS6
+#define qS7 _ItL_qS7
+#define qS8 _ItL_qS8
+#define qS9 _ItL_qS9
+#define atanhi _ItL_atanhi
+#define atanlo _ItL_atanlo
+#define aT _ItL_aT
+#define pi_lo _ItL_pi_lo
+
+#define pio2_hi atanhi[3]
+#define pio2_lo atanlo[3]
+#define pio4_hi atanhi[1]
+
+/* Constants shared by the long double inverse trig functions. */
+extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6, pS7, pS8, pS9;
+extern const long double qS1, qS2, qS3, qS4, qS5, qS6, qS7, qS8, qS9;
+extern const long double atanhi[], atanlo[], aT[];
+extern const long double pi_lo;
+
+static inline long double
+P(long double x)
+{
+
+ return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
+ (pS4 + x * (pS5 + x * (pS6 + x * (pS7 + x * (pS8 + x * \
+ pS9))))))))));
+}
+
+static inline long double
+Q(long double x)
+{
+
+ return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * \
+ (qS5 + x * (qS6 + x * (qS7 + x * (qS8 + x * qS9)))))))));
+}
+
+static inline long double
+T_even(long double x)
+{
+
+ return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
+ (aT[8] + x * (aT[10] + x * (aT[12] + x * (aT[14] + x * \
+ (aT[16] + x * (aT[18] + x * (aT[20] + x * aT[22])))))))))));
+}
+
+static inline long double
+T_odd(long double x)
+{
+
+ return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
+ (aT[9] + x * (aT[11] + x * (aT[13] + x * (aT[15] + x * \
+ (aT[17] + x * (aT[19] + x * (aT[21] + x * aT[23])))))))))));
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/k_cosl.c b/libm/upstream-freebsd/lib/msun/ld128/k_cosl.c
new file mode 100644
index 0000000..5f4aa37
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/k_cosl.c
@@ -0,0 +1,61 @@
+/* From: @(#)k_cos.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * 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$");
+
+/*
+ * ld128 version of k_cos.c. See ../src/k_cos.c for most comments.
+ */
+
+#include "math_private.h"
+
+/*
+ * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]:
+ * |cos(x) - c(x))| < 2**-122.0
+ *
+ * 113-bit precision requires more care than 64-bit precision, since
+ * simple methods give a minimax polynomial with coefficient for x^2
+ * that is 1 ulp below 0.5, but we want it to be precisely 0.5. See
+ * ../ld80/k_cosl.c for more details.
+ */
+static const double
+one = 1.0;
+
+static const long double
+C1 = 0.04166666666666666666666666666666658424671L,
+C2 = -0.001388888888888888888888888888863490893732L,
+C3 = 0.00002480158730158730158730158600795304914210L,
+C4 = -0.2755731922398589065255474947078934284324e-6L,
+C5 = 0.2087675698786809897659225313136400793948e-8L,
+C6 = -0.1147074559772972315817149986812031204775e-10L,
+C7 = 0.4779477332386808976875457937252120293400e-13L;
+
+static const double
+C8 = -0.1561920696721507929516718307820958119868e-15,
+C9 = 0.4110317413744594971475941557607804508039e-18,
+C10 = -0.8896592467191938803288521958313920156409e-21,
+C11 = 0.1601061435794535138244346256065192782581e-23;
+
+long double
+__kernel_cosl(long double x, long double y)
+{
+ long double hz,z,r,w;
+
+ z = x*x;
+ r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*(C7+
+ z*(C8+z*(C9+z*(C10+z*C11))))))))));
+ hz = 0.5*z;
+ w = one-hz;
+ return w + (((one-w)-hz) + (z*r-x*y));
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/k_sinl.c b/libm/upstream-freebsd/lib/msun/ld128/k_sinl.c
new file mode 100644
index 0000000..bd415c0
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/k_sinl.c
@@ -0,0 +1,59 @@
+/* From: @(#)k_sin.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * 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$");
+
+/*
+ * ld128 version of k_sin.c. See ../src/k_sin.c for most comments.
+ */
+
+#include "math_private.h"
+
+static const double
+half = 0.5;
+
+/*
+ * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37]
+ * |sin(x)/x - s(x)| < 2**-122.1
+ *
+ * See ../ld80/k_cosl.c for more details about the polynomial.
+ */
+static const long double
+S1 = -0.16666666666666666666666666666666666606732416116558L,
+S2 = 0.0083333333333333333333333333333331135404851288270047L,
+S3 = -0.00019841269841269841269841269839935785325638310428717L,
+S4 = 0.27557319223985890652557316053039946268333231205686e-5L,
+S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
+S6 = 0.16059043836821614596571832194524392581082444805729e-9L,
+S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
+S8 = 0.28114572543451292625024967174638477283187397621303e-14L;
+
+static const double
+S9 = -0.82206352458348947812512122163446202498005154296863e-17,
+S10 = 0.19572940011906109418080609928334380560135358385256e-19,
+S11 = -0.38680813379701966970673724299207480965452616911420e-22,
+S12 = 0.64038150078671872796678569586315881020659912139412e-25;
+
+long double
+__kernel_sinl(long double x, long double y, int iy)
+{
+ long double z,r,v;
+
+ z = x*x;
+ v = z*x;
+ r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*(S8+
+ z*(S9+z*(S10+z*(S11+z*S12)))))))));
+ if(iy==0) return x+v*(S1+z*r);
+ else return x-((z*(half*y-v*r)-y)-v*S1);
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/k_tanl.c b/libm/upstream-freebsd/lib/msun/ld128/k_tanl.c
new file mode 100644
index 0000000..d7ec0b9
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/k_tanl.c
@@ -0,0 +1,119 @@
+/* From: @(#)k_tan.c 1.5 04/04/22 SMI */
+
+/*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * ld128 version of k_tan.c. See ../src/k_tan.c for most comments.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37]
+ * |tan(x)/x - t(x)| < 2**-117.8 (XXX should be ~1e-37)
+ *
+ * See ../ld80/k_cosl.c for more details about the polynomial.
+ */
+static const long double
+T3 = 0x1.5555555555555555555555555553p-2L,
+T5 = 0x1.1111111111111111111111111eb5p-3L,
+T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L,
+T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L,
+T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L,
+T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L,
+T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L,
+T17 = 0x1.355824803674477dfcf726649efep-11L,
+T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L,
+T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L,
+T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L,
+T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L,
+T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L,
+T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L,
+T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L,
+T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L,
+T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L,
+T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L,
+pio4 = 0x1.921fb54442d18469898cc51701b8p-1L,
+pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L;
+
+static const double
+T39 = 0.000000028443389121318352, /* 0x1e8a7592977938.0p-78 */
+T41 = 0.000000011981013102001973, /* 0x19baa1b1223219.0p-79 */
+T43 = 0.0000000038303578044958070, /* 0x107385dfb24529.0p-80 */
+T45 = 0.0000000034664378216909893, /* 0x1dc6c702a05262.0p-81 */
+T47 = -0.0000000015090641701997785, /* -0x19ecef3569ebb6.0p-82 */
+T49 = 0.0000000029449552300483952, /* 0x194c0668da786a.0p-81 */
+T51 = -0.0000000022006995706097711, /* -0x12e763b8845268.0p-81 */
+T53 = 0.0000000015468200913196612, /* 0x1a92fc98c29554.0p-82 */
+T55 = -0.00000000061311613386849674, /* -0x151106cbc779a9.0p-83 */
+T57 = 1.4912469681508012e-10; /* 0x147edbdba6f43a.0p-85 */
+
+long double
+__kernel_tanl(long double x, long double y, int iy) {
+ long double z, r, v, w, s;
+ long double osign;
+ int i;
+
+ iy = (iy == 1 ? -1 : 1); /* XXX recover original interface */
+ osign = (x >= 0 ? 1.0 : -1.0); /* XXX slow, probably wrong for -0 */
+ if (fabsl(x) >= 0.67434) {
+ if (x < 0) {
+ x = -x;
+ y = -y;
+ }
+ z = pio4 - x;
+ w = pio4lo - y;
+ x = z + w;
+ y = 0.0;
+ i = 1;
+ } else
+ i = 0;
+ z = x * x;
+ w = z * z;
+ r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
+ w * (T25 + w * (T29 + w * (T33 +
+ w * (T37 + w * (T41 + w * (T45 + w * (T49 + w * (T53 +
+ w * T57))))))))))));
+ v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
+ w * (T27 + w * (T31 + w * (T35 +
+ w * (T39 + w * (T43 + w * (T47 + w * (T51 + w * T55))))))))))));
+ s = z * x;
+ r = y + z * (s * (r + v) + y);
+ r += T3 * s;
+ w = x + r;
+ if (i == 1) {
+ v = (long double) iy;
+ return osign *
+ (v - 2.0 * (x - (w * w / (w + v) - r)));
+ }
+ if (iy == 1)
+ return w;
+ else {
+ /*
+ * if allow error up to 2 ulp, simply return
+ * -1.0 / (x+r) here
+ */
+ /* compute -1.0 / (x+r) accurately */
+ long double a, t;
+ z = w;
+ z = z + 0x1p32 - 0x1p32;
+ v = r - (z - x); /* z+v = r+x */
+ t = a = -1.0 / w; /* a = -1.0/w */
+ t = t + 0x1p32 - 0x1p32;
+ s = 1.0 + t * z;
+ return t + a * (s + t * v);
+ }
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/s_exp2l.c b/libm/upstream-freebsd/lib/msun/ld128/s_exp2l.c
new file mode 100644
index 0000000..5ed514c
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/s_exp2l.c
@@ -0,0 +1,427 @@
+/*-
+ * Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG>
+ * 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>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#include <stdint.h>
+
+#include "fpmath.h"
+#include "math.h"
+
+#define TBLBITS 7
+#define TBLSIZE (1 << TBLBITS)
+
+#define BIAS (LDBL_MAX_EXP - 1)
+#define EXPMASK (BIAS + LDBL_MAX_EXP)
+
+static volatile long double
+ huge = 0x1p10000L,
+ twom10000 = 0x1p-10000L;
+
+static const long double
+ P1 = 0x1.62e42fefa39ef35793c7673007e6p-1L,
+ P2 = 0x1.ebfbdff82c58ea86f16b06ec9736p-3L,
+ P3 = 0x1.c6b08d704a0bf8b33a762bad3459p-5L,
+ P4 = 0x1.3b2ab6fba4e7729ccbbe0b4f3fc2p-7L,
+ P5 = 0x1.5d87fe78a67311071dee13fd11d9p-10L,
+ P6 = 0x1.430912f86c7876f4b663b23c5fe5p-13L;
+
+static const double
+ P7 = 0x1.ffcbfc588b041p-17,
+ P8 = 0x1.62c0223a5c7c7p-20,
+ P9 = 0x1.b52541ff59713p-24,
+ P10 = 0x1.e4cf56a391e22p-28,
+ redux = 0x1.8p112 / TBLSIZE;
+
+static const long double tbl[TBLSIZE] = {
+ 0x1.6a09e667f3bcc908b2fb1366dfeap-1L,
+ 0x1.6c012750bdabeed76a99800f4edep-1L,
+ 0x1.6dfb23c651a2ef220e2cbe1bc0d4p-1L,
+ 0x1.6ff7df9519483cf87e1b4f3e1e98p-1L,
+ 0x1.71f75e8ec5f73dd2370f2ef0b148p-1L,
+ 0x1.73f9a48a58173bd5c9a4e68ab074p-1L,
+ 0x1.75feb564267c8bf6e9aa33a489a8p-1L,
+ 0x1.780694fde5d3f619ae02808592a4p-1L,
+ 0x1.7a11473eb0186d7d51023f6ccb1ap-1L,
+ 0x1.7c1ed0130c1327c49334459378dep-1L,
+ 0x1.7e2f336cf4e62105d02ba1579756p-1L,
+ 0x1.80427543e1a11b60de67649a3842p-1L,
+ 0x1.82589994cce128acf88afab34928p-1L,
+ 0x1.8471a4623c7acce52f6b97c6444cp-1L,
+ 0x1.868d99b4492ec80e41d90ac2556ap-1L,
+ 0x1.88ac7d98a669966530bcdf2d4cc0p-1L,
+ 0x1.8ace5422aa0db5ba7c55a192c648p-1L,
+ 0x1.8cf3216b5448bef2aa1cd161c57ap-1L,
+ 0x1.8f1ae991577362b982745c72eddap-1L,
+ 0x1.9145b0b91ffc588a61b469f6b6a0p-1L,
+ 0x1.93737b0cdc5e4f4501c3f2540ae8p-1L,
+ 0x1.95a44cbc8520ee9b483695a0e7fep-1L,
+ 0x1.97d829fde4e4f8b9e920f91e8eb6p-1L,
+ 0x1.9a0f170ca07b9ba3109b8c467844p-1L,
+ 0x1.9c49182a3f0901c7c46b071f28dep-1L,
+ 0x1.9e86319e323231824ca78e64c462p-1L,
+ 0x1.a0c667b5de564b29ada8b8cabbacp-1L,
+ 0x1.a309bec4a2d3358c171f770db1f4p-1L,
+ 0x1.a5503b23e255c8b424491caf88ccp-1L,
+ 0x1.a799e1330b3586f2dfb2b158f31ep-1L,
+ 0x1.a9e6b5579fdbf43eb243bdff53a2p-1L,
+ 0x1.ac36bbfd3f379c0db966a3126988p-1L,
+ 0x1.ae89f995ad3ad5e8734d17731c80p-1L,
+ 0x1.b0e07298db66590842acdfc6fb4ep-1L,
+ 0x1.b33a2b84f15faf6bfd0e7bd941b0p-1L,
+ 0x1.b59728de559398e3881111648738p-1L,
+ 0x1.b7f76f2fb5e46eaa7b081ab53ff6p-1L,
+ 0x1.ba5b030a10649840cb3c6af5b74cp-1L,
+ 0x1.bcc1e904bc1d2247ba0f45b3d06cp-1L,
+ 0x1.bf2c25bd71e088408d7025190cd0p-1L,
+ 0x1.c199bdd85529c2220cb12a0916bap-1L,
+ 0x1.c40ab5fffd07a6d14df820f17deap-1L,
+ 0x1.c67f12e57d14b4a2137fd20f2a26p-1L,
+ 0x1.c8f6d9406e7b511acbc48805c3f6p-1L,
+ 0x1.cb720dcef90691503cbd1e949d0ap-1L,
+ 0x1.cdf0b555dc3f9c44f8958fac4f12p-1L,
+ 0x1.d072d4a07897b8d0f22f21a13792p-1L,
+ 0x1.d2f87080d89f18ade123989ea50ep-1L,
+ 0x1.d5818dcfba48725da05aeb66dff8p-1L,
+ 0x1.d80e316c98397bb84f9d048807a0p-1L,
+ 0x1.da9e603db3285708c01a5b6d480cp-1L,
+ 0x1.dd321f301b4604b695de3c0630c0p-1L,
+ 0x1.dfc97337b9b5eb968cac39ed284cp-1L,
+ 0x1.e264614f5a128a12761fa17adc74p-1L,
+ 0x1.e502ee78b3ff6273d130153992d0p-1L,
+ 0x1.e7a51fbc74c834b548b2832378a4p-1L,
+ 0x1.ea4afa2a490d9858f73a18f5dab4p-1L,
+ 0x1.ecf482d8e67f08db0312fb949d50p-1L,
+ 0x1.efa1bee615a27771fd21a92dabb6p-1L,
+ 0x1.f252b376bba974e8696fc3638f24p-1L,
+ 0x1.f50765b6e4540674f84b762861a6p-1L,
+ 0x1.f7bfdad9cbe138913b4bfe72bd78p-1L,
+ 0x1.fa7c1819e90d82e90a7e74b26360p-1L,
+ 0x1.fd3c22b8f71f10975ba4b32bd006p-1L,
+ 0x1.0000000000000000000000000000p+0L,
+ 0x1.0163da9fb33356d84a66ae336e98p+0L,
+ 0x1.02c9a3e778060ee6f7caca4f7a18p+0L,
+ 0x1.04315e86e7f84bd738f9a20da442p+0L,
+ 0x1.059b0d31585743ae7c548eb68c6ap+0L,
+ 0x1.0706b29ddf6ddc6dc403a9d87b1ep+0L,
+ 0x1.0874518759bc808c35f25d942856p+0L,
+ 0x1.09e3ecac6f3834521e060c584d5cp+0L,
+ 0x1.0b5586cf9890f6298b92b7184200p+0L,
+ 0x1.0cc922b7247f7407b705b893dbdep+0L,
+ 0x1.0e3ec32d3d1a2020742e4f8af794p+0L,
+ 0x1.0fb66affed31af232091dd8a169ep+0L,
+ 0x1.11301d0125b50a4ebbf1aed9321cp+0L,
+ 0x1.12abdc06c31cbfb92bad324d6f84p+0L,
+ 0x1.1429aaea92ddfb34101943b2588ep+0L,
+ 0x1.15a98c8a58e512480d573dd562aep+0L,
+ 0x1.172b83c7d517adcdf7c8c50eb162p+0L,
+ 0x1.18af9388c8de9bbbf70b9a3c269cp+0L,
+ 0x1.1a35beb6fcb753cb698f692d2038p+0L,
+ 0x1.1bbe084045cd39ab1e72b442810ep+0L,
+ 0x1.1d4873168b9aa7805b8028990be8p+0L,
+ 0x1.1ed5022fcd91cb8819ff61121fbep+0L,
+ 0x1.2063b88628cd63b8eeb0295093f6p+0L,
+ 0x1.21f49917ddc962552fd29294bc20p+0L,
+ 0x1.2387a6e75623866c1fadb1c159c0p+0L,
+ 0x1.251ce4fb2a63f3582ab7de9e9562p+0L,
+ 0x1.26b4565e27cdd257a673281d3068p+0L,
+ 0x1.284dfe1f5638096cf15cf03c9fa0p+0L,
+ 0x1.29e9df51fdee12c25d15f5a25022p+0L,
+ 0x1.2b87fd0dad98ffddea46538fca24p+0L,
+ 0x1.2d285a6e4030b40091d536d0733ep+0L,
+ 0x1.2ecafa93e2f5611ca0f45d5239a4p+0L,
+ 0x1.306fe0a31b7152de8d5a463063bep+0L,
+ 0x1.32170fc4cd8313539cf1c3009330p+0L,
+ 0x1.33c08b26416ff4c9c8610d96680ep+0L,
+ 0x1.356c55f929ff0c94623476373be4p+0L,
+ 0x1.371a7373aa9caa7145502f45452ap+0L,
+ 0x1.38cae6d05d86585a9cb0d9bed530p+0L,
+ 0x1.3a7db34e59ff6ea1bc9299e0a1fep+0L,
+ 0x1.3c32dc313a8e484001f228b58cf0p+0L,
+ 0x1.3dea64c12342235b41223e13d7eep+0L,
+ 0x1.3fa4504ac801ba0bf701aa417b9cp+0L,
+ 0x1.4160a21f72e29f84325b8f3dbacap+0L,
+ 0x1.431f5d950a896dc704439410b628p+0L,
+ 0x1.44e086061892d03136f409df0724p+0L,
+ 0x1.46a41ed1d005772512f459229f0ap+0L,
+ 0x1.486a2b5c13cd013c1a3b69062f26p+0L,
+ 0x1.4a32af0d7d3de672d8bcf46f99b4p+0L,
+ 0x1.4bfdad5362a271d4397afec42e36p+0L,
+ 0x1.4dcb299fddd0d63b36ef1a9e19dep+0L,
+ 0x1.4f9b2769d2ca6ad33d8b69aa0b8cp+0L,
+ 0x1.516daa2cf6641c112f52c84d6066p+0L,
+ 0x1.5342b569d4f81df0a83c49d86bf4p+0L,
+ 0x1.551a4ca5d920ec52ec620243540cp+0L,
+ 0x1.56f4736b527da66ecb004764e61ep+0L,
+ 0x1.58d12d497c7fd252bc2b7343d554p+0L,
+ 0x1.5ab07dd48542958c93015191e9a8p+0L,
+ 0x1.5c9268a5946b701c4b1b81697ed4p+0L,
+ 0x1.5e76f15ad21486e9be4c20399d12p+0L,
+ 0x1.605e1b976dc08b076f592a487066p+0L,
+ 0x1.6247eb03a5584b1f0fa06fd2d9eap+0L,
+ 0x1.6434634ccc31fc76f8714c4ee122p+0L,
+ 0x1.66238825522249127d9e29b92ea2p+0L,
+ 0x1.68155d44ca973081c57227b9f69ep+0L,
+};
+
+static const float eps[TBLSIZE] = {
+ -0x1.5c50p-101,
+ -0x1.5d00p-106,
+ 0x1.8e90p-102,
+ -0x1.5340p-103,
+ 0x1.1bd0p-102,
+ -0x1.4600p-105,
+ -0x1.7a40p-104,
+ 0x1.d590p-102,
+ -0x1.d590p-101,
+ 0x1.b100p-103,
+ -0x1.0d80p-105,
+ 0x1.6b00p-103,
+ -0x1.9f00p-105,
+ 0x1.c400p-103,
+ 0x1.e120p-103,
+ -0x1.c100p-104,
+ -0x1.9d20p-103,
+ 0x1.a800p-108,
+ 0x1.4c00p-106,
+ -0x1.9500p-106,
+ 0x1.6900p-105,
+ -0x1.29d0p-100,
+ 0x1.4c60p-103,
+ 0x1.13a0p-102,
+ -0x1.5b60p-103,
+ -0x1.1c40p-103,
+ 0x1.db80p-102,
+ 0x1.91a0p-102,
+ 0x1.dc00p-105,
+ 0x1.44c0p-104,
+ 0x1.9710p-102,
+ 0x1.8760p-103,
+ -0x1.a720p-103,
+ 0x1.ed20p-103,
+ -0x1.49c0p-102,
+ -0x1.e000p-111,
+ 0x1.86a0p-103,
+ 0x1.2b40p-103,
+ -0x1.b400p-108,
+ 0x1.1280p-99,
+ -0x1.02d8p-102,
+ -0x1.e3d0p-103,
+ -0x1.b080p-105,
+ -0x1.f100p-107,
+ -0x1.16c0p-105,
+ -0x1.1190p-103,
+ -0x1.a7d2p-100,
+ 0x1.3450p-103,
+ -0x1.67c0p-105,
+ 0x1.4b80p-104,
+ -0x1.c4e0p-103,
+ 0x1.6000p-108,
+ -0x1.3f60p-105,
+ 0x1.93f0p-104,
+ 0x1.5fe0p-105,
+ 0x1.6f80p-107,
+ -0x1.7600p-106,
+ 0x1.21e0p-106,
+ -0x1.3a40p-106,
+ -0x1.40c0p-104,
+ -0x1.9860p-105,
+ -0x1.5d40p-108,
+ -0x1.1d70p-106,
+ 0x1.2760p-105,
+ 0x0.0000p+0,
+ 0x1.21e2p-104,
+ -0x1.9520p-108,
+ -0x1.5720p-106,
+ -0x1.4810p-106,
+ -0x1.be00p-109,
+ 0x1.0080p-105,
+ -0x1.5780p-108,
+ -0x1.d460p-105,
+ -0x1.6140p-105,
+ 0x1.4630p-104,
+ 0x1.ad50p-103,
+ 0x1.82e0p-105,
+ 0x1.1d3cp-101,
+ 0x1.6100p-107,
+ 0x1.ec30p-104,
+ 0x1.f200p-108,
+ 0x1.0b40p-103,
+ 0x1.3660p-102,
+ 0x1.d9d0p-103,
+ -0x1.02d0p-102,
+ 0x1.b070p-103,
+ 0x1.b9c0p-104,
+ -0x1.01c0p-103,
+ -0x1.dfe0p-103,
+ 0x1.1b60p-104,
+ -0x1.ae94p-101,
+ -0x1.3340p-104,
+ 0x1.b3d8p-102,
+ -0x1.6e40p-105,
+ -0x1.3670p-103,
+ 0x1.c140p-104,
+ 0x1.1840p-101,
+ 0x1.1ab0p-102,
+ -0x1.a400p-104,
+ 0x1.1f00p-104,
+ -0x1.7180p-103,
+ 0x1.4ce0p-102,
+ 0x1.9200p-107,
+ -0x1.54c0p-103,
+ 0x1.1b80p-105,
+ -0x1.1828p-101,
+ 0x1.5720p-102,
+ -0x1.a060p-100,
+ 0x1.9160p-102,
+ 0x1.a280p-104,
+ 0x1.3400p-107,
+ 0x1.2b20p-102,
+ 0x1.7800p-108,
+ 0x1.cfd0p-101,
+ 0x1.2ef0p-102,
+ -0x1.2760p-99,
+ 0x1.b380p-104,
+ 0x1.0048p-101,
+ -0x1.60b0p-102,
+ 0x1.a1ccp-100,
+ -0x1.a640p-104,
+ -0x1.08a0p-101,
+ 0x1.7e60p-102,
+ 0x1.22c0p-103,
+ -0x1.7200p-106,
+ 0x1.f0f0p-102,
+ 0x1.eb4ep-99,
+ 0x1.c6e0p-103,
+};
+
+/*
+ * exp2l(x): compute the base 2 exponential of x
+ *
+ * Accuracy: Peak error < 0.502 ulp.
+ *
+ * Method: (accurate tables)
+ *
+ * Reduce x:
+ * x = 2**k + y, for integer k and |y| <= 1/2.
+ * Thus we have exp2(x) = 2**k * exp2(y).
+ *
+ * Reduce y:
+ * y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
+ * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
+ * with |z - eps[i]| <= 2**-8 + 2**-98 for the table used.
+ *
+ * We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
+ * a degree-10 minimax polynomial with maximum error under 2**-120.
+ * The values in exp2t[] and eps[] are chosen such that
+ * exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
+ * that exp2t[i] is accurate to 2**-122.
+ *
+ * Note that the range of i is +-TBLSIZE/2, so we actually index the tables
+ * by i0 = i + TBLSIZE/2.
+ *
+ * This method is due to Gal, with many details due to Gal and Bachelis:
+ *
+ * Gal, S. and Bachelis, B. An Accurate Elementary Mathematical Library
+ * for the IEEE Floating Point Standard. TOMS 17(1), 26-46 (1991).
+ */
+long double
+exp2l(long double x)
+{
+ union IEEEl2bits u, v;
+ long double r, t, twopk, twopkp10000, z;
+ uint32_t hx, ix, i0;
+ int k;
+
+ u.e = x;
+
+ /* Filter out exceptional cases. */
+ hx = u.xbits.expsign;
+ ix = hx & EXPMASK;
+ if (ix >= BIAS + 14) { /* |x| >= 16384 */
+ if (ix == BIAS + LDBL_MAX_EXP) {
+ if (u.xbits.manh != 0
+ || u.xbits.manl != 0
+ || (hx & 0x8000) == 0)
+ return (x + x); /* x is NaN or +Inf */
+ else
+ return (0.0); /* x is -Inf */
+ }
+ if (x >= 16384)
+ return (huge * huge); /* overflow */
+ if (x <= -16495)
+ return (twom10000 * twom10000); /* underflow */
+ } else if (ix <= BIAS - 115) { /* |x| < 0x1p-115 */
+ return (1.0 + x);
+ }
+
+ /*
+ * Reduce x, computing z, i0, and k. The low bits of x + redux
+ * contain the 16-bit integer part of the exponent (k) followed by
+ * TBLBITS fractional bits (i0). We use bit tricks to extract these
+ * as integers, then set z to the remainder.
+ *
+ * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
+ * Then the low-order word of x + redux is 0x000abc12,
+ * We split this into k = 0xabc and i0 = 0x12 (adjusted to
+ * index into the table), then we compute z = 0x0.003456p0.
+ *
+ * XXX If the exponent is negative, the computation of k depends on
+ * '>>' doing sign extension.
+ */
+ u.e = x + redux;
+ i0 = (u.bits.manl & 0xffffffff) + TBLSIZE / 2;
+ k = (int)i0 >> TBLBITS;
+ i0 = i0 & (TBLSIZE - 1);
+ u.e -= redux;
+ z = x - u.e;
+ v.xbits.manh = 0;
+ v.xbits.manl = 0;
+ if (k >= LDBL_MIN_EXP) {
+ v.xbits.expsign = LDBL_MAX_EXP - 1 + k;
+ twopk = v.e;
+ } else {
+ v.xbits.expsign = LDBL_MAX_EXP - 1 + k + 10000;
+ twopkp10000 = v.e;
+ }
+
+ /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
+ t = tbl[i0]; /* exp2t[i0] */
+ z -= eps[i0]; /* eps[i0] */
+ r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * (P5 + z * (P6
+ + z * (P7 + z * (P8 + z * (P9 + z * P10)))))))));
+
+ /* Scale by 2**k. */
+ if(k >= LDBL_MIN_EXP) {
+ if (k == LDBL_MAX_EXP)
+ return (r * 2.0 * 0x1p16383L);
+ return (r * twopk);
+ } else {
+ return (r * twopkp10000 * twom10000);
+ }
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/s_expl.c b/libm/upstream-freebsd/lib/msun/ld128/s_expl.c
new file mode 100644
index 0000000..176c932
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/s_expl.c
@@ -0,0 +1,494 @@
+/*-
+ * Copyright (c) 2009-2013 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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.
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * ld128 version of s_expl.c. See ../ld80/s_expl.c for most comments.
+ */
+
+#include <float.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#define INTERVALS 128
+#define LOG2_INTERVALS 7
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const long double
+huge = 0x1p10000L,
+twom10000 = 0x1p-10000L;
+/* XXX Prevent gcc from erroneously constant folding this: */
+static volatile const long double tiny = 0x1p-10000L;
+
+static const long double
+/* log(2**16384 - 0.5) rounded towards zero: */
+/* log(2**16384 - 0.5 + 1) rounded towards zero for expm1l() is the same: */
+o_threshold = 11356.523406294143949491931077970763428L,
+/* log(2**(-16381-64-1)) rounded towards zero: */
+u_threshold = -11433.462743336297878837243843452621503L;
+
+static const double
+/*
+ * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication). L1 must
+ * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
+ * bits zero so that multiplication of it by n is exact.
+ */
+INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */
+L2 = -1.0253670638894731e-29; /* -0x1.9ff0342542fc3p-97 */
+static const long double
+/* 0x1.62e42fefa39ef35793c768000000p-8 */
+L1 = 5.41521234812457272982212595914567508e-3L;
+
+static const long double
+/*
+ * Domain [-0.002708, 0.002708], range ~[-2.4021e-38, 2.4234e-38]:
+ * |exp(x) - p(x)| < 2**-124.9
+ * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
+ */
+A2 = 0.5,
+A3 = 1.66666666666666666666666666651085500e-1L,
+A4 = 4.16666666666666666666666666425885320e-2L,
+A5 = 8.33333333333333333334522877160175842e-3L,
+A6 = 1.38888888888888888889971139751596836e-3L;
+
+static const double
+A7 = 1.9841269841269471e-4,
+A8 = 2.4801587301585284e-5,
+A9 = 2.7557324277411234e-6,
+A10 = 2.7557333722375072e-7;
+
+static const struct {
+ /*
+ * hi must be rounded to at most 106 bits so that multiplication
+ * by r1 in expm1l() is exact, but it is rounded to 88 bits due to
+ * historical accidents.
+ */
+ long double hi;
+ long double lo;
+} tbl[INTERVALS] = {
+ 0x1p0L, 0x0p0L,
+ 0x1.0163da9fb33356d84a66aep0L, 0x3.36dcdfa4003ec04c360be2404078p-92L,
+ 0x1.02c9a3e778060ee6f7cacap0L, 0x4.f7a29bde93d70a2cabc5cb89ba10p-92L,
+ 0x1.04315e86e7f84bd738f9a2p0L, 0xd.a47e6ed040bb4bfc05af6455e9b8p-96L,
+ 0x1.059b0d31585743ae7c548ep0L, 0xb.68ca417fe53e3495f7df4baf84a0p-92L,
+ 0x1.0706b29ddf6ddc6dc403a8p0L, 0x1.d87b27ed07cb8b092ac75e311753p-88L,
+ 0x1.0874518759bc808c35f25cp0L, 0x1.9427fa2b041b2d6829d8993a0d01p-88L,
+ 0x1.09e3ecac6f3834521e060cp0L, 0x5.84d6b74ba2e023da730e7fccb758p-92L,
+ 0x1.0b5586cf9890f6298b92b6p0L, 0x1.1842a98364291408b3ceb0a2a2bbp-88L,
+ 0x1.0cc922b7247f7407b705b8p0L, 0x9.3dc5e8aac564e6fe2ef1d431fd98p-92L,
+ 0x1.0e3ec32d3d1a2020742e4ep0L, 0x1.8af6a552ac4b358b1129e9f966a4p-88L,
+ 0x1.0fb66affed31af232091dcp0L, 0x1.8a1426514e0b627bda694a400a27p-88L,
+ 0x1.11301d0125b50a4ebbf1aep0L, 0xd.9318ceac5cc47ab166ee57427178p-92L,
+ 0x1.12abdc06c31cbfb92bad32p0L, 0x4.d68e2f7270bdf7cedf94eb1cb818p-92L,
+ 0x1.1429aaea92ddfb34101942p0L, 0x1.b2586d01844b389bea7aedd221d4p-88L,
+ 0x1.15a98c8a58e512480d573cp0L, 0x1.d5613bf92a2b618ee31b376c2689p-88L,
+ 0x1.172b83c7d517adcdf7c8c4p0L, 0x1.0eb14a792035509ff7d758693f24p-88L,
+ 0x1.18af9388c8de9bbbf70b9ap0L, 0x3.c2505c97c0102e5f1211941d2840p-92L,
+ 0x1.1a35beb6fcb753cb698f68p0L, 0x1.2d1c835a6c30724d5cfae31b84e5p-88L,
+ 0x1.1bbe084045cd39ab1e72b4p0L, 0x4.27e35f9acb57e473915519a1b448p-92L,
+ 0x1.1d4873168b9aa7805b8028p0L, 0x9.90f07a98b42206e46166cf051d70p-92L,
+ 0x1.1ed5022fcd91cb8819ff60p0L, 0x1.121d1e504d36c47474c9b7de6067p-88L,
+ 0x1.2063b88628cd63b8eeb028p0L, 0x1.50929d0fc487d21c2b84004264dep-88L,
+ 0x1.21f49917ddc962552fd292p0L, 0x9.4bdb4b61ea62477caa1dce823ba0p-92L,
+ 0x1.2387a6e75623866c1fadb0p0L, 0x1.c15cb593b0328566902df69e4de2p-88L,
+ 0x1.251ce4fb2a63f3582ab7dep0L, 0x9.e94811a9c8afdcf796934bc652d0p-92L,
+ 0x1.26b4565e27cdd257a67328p0L, 0x1.d3b249dce4e9186ddd5ff44e6b08p-92L,
+ 0x1.284dfe1f5638096cf15cf0p0L, 0x3.ca0967fdaa2e52d7c8106f2e262cp-92L,
+ 0x1.29e9df51fdee12c25d15f4p0L, 0x1.a24aa3bca890ac08d203fed80a07p-88L,
+ 0x1.2b87fd0dad98ffddea4652p0L, 0x1.8fcab88442fdc3cb6de4519165edp-88L,
+ 0x1.2d285a6e4030b40091d536p0L, 0xd.075384589c1cd1b3e4018a6b1348p-92L,
+ 0x1.2ecafa93e2f5611ca0f45cp0L, 0x1.523833af611bdcda253c554cf278p-88L,
+ 0x1.306fe0a31b7152de8d5a46p0L, 0x3.05c85edecbc27343629f502f1af2p-92L,
+ 0x1.32170fc4cd8313539cf1c2p0L, 0x1.008f86dde3220ae17a005b6412bep-88L,
+ 0x1.33c08b26416ff4c9c8610cp0L, 0x1.96696bf95d1593039539d94d662bp-88L,
+ 0x1.356c55f929ff0c94623476p0L, 0x3.73af38d6d8d6f9506c9bbc93cbc0p-92L,
+ 0x1.371a7373aa9caa7145502ep0L, 0x1.4547987e3e12516bf9c699be432fp-88L,
+ 0x1.38cae6d05d86585a9cb0d8p0L, 0x1.bed0c853bd30a02790931eb2e8f0p-88L,
+ 0x1.3a7db34e59ff6ea1bc9298p0L, 0x1.e0a1d336163fe2f852ceeb134067p-88L,
+ 0x1.3c32dc313a8e484001f228p0L, 0xb.58f3775e06ab66353001fae9fca0p-92L,
+ 0x1.3dea64c12342235b41223ep0L, 0x1.3d773fba2cb82b8244267c54443fp-92L,
+ 0x1.3fa4504ac801ba0bf701aap0L, 0x4.1832fb8c1c8dbdff2c49909e6c60p-92L,
+ 0x1.4160a21f72e29f84325b8ep0L, 0x1.3db61fb352f0540e6ba05634413ep-88L,
+ 0x1.431f5d950a896dc7044394p0L, 0x1.0ccec81e24b0caff7581ef4127f7p-92L,
+ 0x1.44e086061892d03136f408p0L, 0x1.df019fbd4f3b48709b78591d5cb5p-88L,
+ 0x1.46a41ed1d005772512f458p0L, 0x1.229d97df404ff21f39c1b594d3a8p-88L,
+ 0x1.486a2b5c13cd013c1a3b68p0L, 0x1.062f03c3dd75ce8757f780e6ec99p-88L,
+ 0x1.4a32af0d7d3de672d8bcf4p0L, 0x6.f9586461db1d878b1d148bd3ccb8p-92L,
+ 0x1.4bfdad5362a271d4397afep0L, 0xc.42e20e0363ba2e159c579f82e4b0p-92L,
+ 0x1.4dcb299fddd0d63b36ef1ap0L, 0x9.e0cc484b25a5566d0bd5f58ad238p-92L,
+ 0x1.4f9b2769d2ca6ad33d8b68p0L, 0x1.aa073ee55e028497a329a7333dbap-88L,
+ 0x1.516daa2cf6641c112f52c8p0L, 0x4.d822190e718226177d7608d20038p-92L,
+ 0x1.5342b569d4f81df0a83c48p0L, 0x1.d86a63f4e672a3e429805b049465p-88L,
+ 0x1.551a4ca5d920ec52ec6202p0L, 0x4.34ca672645dc6c124d6619a87574p-92L,
+ 0x1.56f4736b527da66ecb0046p0L, 0x1.64eb3c00f2f5ab3d801d7cc7272dp-88L,
+ 0x1.58d12d497c7fd252bc2b72p0L, 0x1.43bcf2ec936a970d9cc266f0072fp-88L,
+ 0x1.5ab07dd48542958c930150p0L, 0x1.91eb345d88d7c81280e069fbdb63p-88L,
+ 0x1.5c9268a5946b701c4b1b80p0L, 0x1.6986a203d84e6a4a92f179e71889p-88L,
+ 0x1.5e76f15ad21486e9be4c20p0L, 0x3.99766a06548a05829e853bdb2b52p-92L,
+ 0x1.605e1b976dc08b076f592ap0L, 0x4.86e3b34ead1b4769df867b9c89ccp-92L,
+ 0x1.6247eb03a5584b1f0fa06ep0L, 0x1.d2da42bb1ceaf9f732275b8aef30p-88L,
+ 0x1.6434634ccc31fc76f8714cp0L, 0x4.ed9a4e41000307103a18cf7a6e08p-92L,
+ 0x1.66238825522249127d9e28p0L, 0x1.b8f314a337f4dc0a3adf1787ff74p-88L,
+ 0x1.68155d44ca973081c57226p0L, 0x1.b9f32706bfe4e627d809a85dcc66p-88L,
+ 0x1.6a09e667f3bcc908b2fb12p0L, 0x1.66ea957d3e3adec17512775099dap-88L,
+ 0x1.6c012750bdabeed76a9980p0L, 0xf.4f33fdeb8b0ecd831106f57b3d00p-96L,
+ 0x1.6dfb23c651a2ef220e2cbep0L, 0x1.bbaa834b3f11577ceefbe6c1c411p-92L,
+ 0x1.6ff7df9519483cf87e1b4ep0L, 0x1.3e213bff9b702d5aa477c12523cep-88L,
+ 0x1.71f75e8ec5f73dd2370f2ep0L, 0xf.0acd6cb434b562d9e8a20adda648p-92L,
+ 0x1.73f9a48a58173bd5c9a4e6p0L, 0x8.ab1182ae217f3a7681759553e840p-92L,
+ 0x1.75feb564267c8bf6e9aa32p0L, 0x1.a48b27071805e61a17b954a2dad8p-88L,
+ 0x1.780694fde5d3f619ae0280p0L, 0x8.58b2bb2bdcf86cd08e35fb04c0f0p-92L,
+ 0x1.7a11473eb0186d7d51023ep0L, 0x1.6cda1f5ef42b66977960531e821bp-88L,
+ 0x1.7c1ed0130c1327c4933444p0L, 0x1.937562b2dc933d44fc828efd4c9cp-88L,
+ 0x1.7e2f336cf4e62105d02ba0p0L, 0x1.5797e170a1427f8fcdf5f3906108p-88L,
+ 0x1.80427543e1a11b60de6764p0L, 0x9.a354ea706b8e4d8b718a672bf7c8p-92L,
+ 0x1.82589994cce128acf88afap0L, 0xb.34a010f6ad65cbbac0f532d39be0p-92L,
+ 0x1.8471a4623c7acce52f6b96p0L, 0x1.c64095370f51f48817914dd78665p-88L,
+ 0x1.868d99b4492ec80e41d90ap0L, 0xc.251707484d73f136fb5779656b70p-92L,
+ 0x1.88ac7d98a669966530bcdep0L, 0x1.2d4e9d61283ef385de170ab20f96p-88L,
+ 0x1.8ace5422aa0db5ba7c55a0p0L, 0x1.92c9bb3e6ed61f2733304a346d8fp-88L,
+ 0x1.8cf3216b5448bef2aa1cd0p0L, 0x1.61c55d84a9848f8c453b3ca8c946p-88L,
+ 0x1.8f1ae991577362b982745cp0L, 0x7.2ed804efc9b4ae1458ae946099d4p-92L,
+ 0x1.9145b0b91ffc588a61b468p0L, 0x1.f6b70e01c2a90229a4c4309ea719p-88L,
+ 0x1.93737b0cdc5e4f4501c3f2p0L, 0x5.40a22d2fc4af581b63e8326efe9cp-92L,
+ 0x1.95a44cbc8520ee9b483694p0L, 0x1.a0fc6f7c7d61b2b3a22a0eab2cadp-88L,
+ 0x1.97d829fde4e4f8b9e920f8p0L, 0x1.1e8bd7edb9d7144b6f6818084cc7p-88L,
+ 0x1.9a0f170ca07b9ba3109b8cp0L, 0x4.6737beb19e1eada6825d3c557428p-92L,
+ 0x1.9c49182a3f0901c7c46b06p0L, 0x1.1f2be58ddade50c217186c90b457p-88L,
+ 0x1.9e86319e323231824ca78ep0L, 0x6.4c6e010f92c082bbadfaf605cfd4p-92L,
+ 0x1.a0c667b5de564b29ada8b8p0L, 0xc.ab349aa0422a8da7d4512edac548p-92L,
+ 0x1.a309bec4a2d3358c171f76p0L, 0x1.0daad547fa22c26d168ea762d854p-88L,
+ 0x1.a5503b23e255c8b424491cp0L, 0xa.f87bc8050a405381703ef7caff50p-92L,
+ 0x1.a799e1330b3586f2dfb2b0p0L, 0x1.58f1a98796ce8908ae852236ca94p-88L,
+ 0x1.a9e6b5579fdbf43eb243bcp0L, 0x1.ff4c4c58b571cf465caf07b4b9f5p-88L,
+ 0x1.ac36bbfd3f379c0db966a2p0L, 0x1.1265fc73e480712d20f8597a8e7bp-88L,
+ 0x1.ae89f995ad3ad5e8734d16p0L, 0x1.73205a7fbc3ae675ea440b162d6cp-88L,
+ 0x1.b0e07298db66590842acdep0L, 0x1.c6f6ca0e5dcae2aafffa7a0554cbp-88L,
+ 0x1.b33a2b84f15faf6bfd0e7ap0L, 0x1.d947c2575781dbb49b1237c87b6ep-88L,
+ 0x1.b59728de559398e3881110p0L, 0x1.64873c7171fefc410416be0a6525p-88L,
+ 0x1.b7f76f2fb5e46eaa7b081ap0L, 0xb.53c5354c8903c356e4b625aacc28p-92L,
+ 0x1.ba5b030a10649840cb3c6ap0L, 0xf.5b47f297203757e1cc6eadc8bad0p-92L,
+ 0x1.bcc1e904bc1d2247ba0f44p0L, 0x1.b3d08cd0b20287092bd59be4ad98p-88L,
+ 0x1.bf2c25bd71e088408d7024p0L, 0x1.18e3449fa073b356766dfb568ff4p-88L,
+ 0x1.c199bdd85529c2220cb12ap0L, 0x9.1ba6679444964a36661240043970p-96L,
+ 0x1.c40ab5fffd07a6d14df820p0L, 0xf.1828a5366fd387a7bdd54cdf7300p-92L,
+ 0x1.c67f12e57d14b4a2137fd2p0L, 0xf.2b301dd9e6b151a6d1f9d5d5f520p-96L,
+ 0x1.c8f6d9406e7b511acbc488p0L, 0x5.c442ddb55820171f319d9e5076a8p-96L,
+ 0x1.cb720dcef90691503cbd1ep0L, 0x9.49db761d9559ac0cb6dd3ed599e0p-92L,
+ 0x1.cdf0b555dc3f9c44f8958ep0L, 0x1.ac51be515f8c58bdfb6f5740a3a4p-88L,
+ 0x1.d072d4a07897b8d0f22f20p0L, 0x1.a158e18fbbfc625f09f4cca40874p-88L,
+ 0x1.d2f87080d89f18ade12398p0L, 0x9.ea2025b4c56553f5cdee4c924728p-92L,
+ 0x1.d5818dcfba48725da05aeap0L, 0x1.66e0dca9f589f559c0876ff23830p-88L,
+ 0x1.d80e316c98397bb84f9d04p0L, 0x8.805f84bec614de269900ddf98d28p-92L,
+ 0x1.da9e603db3285708c01a5ap0L, 0x1.6d4c97f6246f0ec614ec95c99392p-88L,
+ 0x1.dd321f301b4604b695de3cp0L, 0x6.30a393215299e30d4fb73503c348p-96L,
+ 0x1.dfc97337b9b5eb968cac38p0L, 0x1.ed291b7225a944efd5bb5524b927p-88L,
+ 0x1.e264614f5a128a12761fa0p0L, 0x1.7ada6467e77f73bf65e04c95e29dp-88L,
+ 0x1.e502ee78b3ff6273d13014p0L, 0x1.3991e8f49659e1693be17ae1d2f9p-88L,
+ 0x1.e7a51fbc74c834b548b282p0L, 0x1.23786758a84f4956354634a416cep-88L,
+ 0x1.ea4afa2a490d9858f73a18p0L, 0xf.5db301f86dea20610ceee13eb7b8p-92L,
+ 0x1.ecf482d8e67f08db0312fap0L, 0x1.949cef462010bb4bc4ce72a900dfp-88L,
+ 0x1.efa1bee615a27771fd21a8p0L, 0x1.2dac1f6dd5d229ff68e46f27e3dfp-88L,
+ 0x1.f252b376bba974e8696fc2p0L, 0x1.6390d4c6ad5476b5162f40e1d9a9p-88L,
+ 0x1.f50765b6e4540674f84b76p0L, 0x2.862baff99000dfc4352ba29b8908p-92L,
+ 0x1.f7bfdad9cbe138913b4bfep0L, 0x7.2bd95c5ce7280fa4d2344a3f5618p-92L,
+ 0x1.fa7c1819e90d82e90a7e74p0L, 0xb.263c1dc060c36f7650b4c0f233a8p-92L,
+ 0x1.fd3c22b8f71f10975ba4b2p0L, 0x1.2bcf3a5e12d269d8ad7c1a4a8875p-88L
+};
+
+long double
+expl(long double x)
+{
+ union IEEEl2bits u, v;
+ long double q, r, r1, t, twopk, twopkp10000;
+ double dr, fn, r2;
+ int k, n, n2;
+ uint16_t hx, ix;
+
+ /* Filter out exceptional cases. */
+ u.e = x;
+ hx = u.xbits.expsign;
+ ix = hx & 0x7fff;
+ if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */
+ if (ix == BIAS + LDBL_MAX_EXP) {
+ if (hx & 0x8000) /* x is -Inf or -NaN */
+ return (-1 / x);
+ return (x + x); /* x is +Inf or +NaN */
+ }
+ if (x > o_threshold)
+ return (huge * huge);
+ if (x < u_threshold)
+ return (tiny * tiny);
+ } else if (ix < BIAS - 114) { /* |x| < 0x1p-114 */
+ return (1 + x); /* 1 with inexact iff x != 0 */
+ }
+
+ ENTERI();
+
+ /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
+ /* Use a specialized rint() to get fn. Assume round-to-nearest. */
+ /* XXX assume no extra precision for the additions, as for trig fns. */
+ /* XXX this set of comments is now quadruplicated. */
+ fn = (double)x * INV_L + 0x1.8p52 - 0x1.8p52;
+#if defined(HAVE_EFFICIENT_IRINT)
+ n = irint(fn);
+#else
+ n = (int)fn;
+#endif
+ n2 = (unsigned)n % INTERVALS;
+ k = n >> LOG2_INTERVALS;
+ r1 = x - fn * L1;
+ r2 = fn * -L2;
+ r = r1 + r2;
+
+ /* Prepare scale factors. */
+ /* XXX sparc64 multiplication is so slow that scalbnl() is faster. */
+ v.e = 1;
+ if (k >= LDBL_MIN_EXP) {
+ v.xbits.expsign = BIAS + k;
+ twopk = v.e;
+ } else {
+ v.xbits.expsign = BIAS + k + 10000;
+ twopkp10000 = v.e;
+ }
+
+ /* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
+ dr = r;
+ q = r2 + r * r * (A2 + r * (A3 + r * (A4 + r * (A5 + r * (A6 +
+ dr * (A7 + dr * (A8 + dr * (A9 + dr * A10))))))));
+ t = tbl[n2].lo + tbl[n2].hi;
+ t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+
+ /* Scale by 2**k. */
+ if (k >= LDBL_MIN_EXP) {
+ if (k == LDBL_MAX_EXP)
+ RETURNI(t * 2 * 0x1p16383L);
+ RETURNI(t * twopk);
+ } else {
+ RETURNI(t * twopkp10000 * twom10000);
+ }
+}
+
+/*
+ * Our T1 and T2 are chosen to be approximately the points where method
+ * A and method B have the same accuracy. Tang's T1 and T2 are the
+ * points where method A's accuracy changes by a full bit. For Tang,
+ * this drop in accuracy makes method A immediately less accurate than
+ * method B, but our larger INTERVALS makes method A 2 bits more
+ * accurate so it remains the most accurate method significantly
+ * closer to the origin despite losing the full bit in our extended
+ * range for it.
+ *
+ * Split the interval [T1, T2] into two intervals [T1, T3] and [T3, T2].
+ * Setting T3 to 0 would require the |x| < 0x1p-113 condition to appear
+ * in both subintervals, so set T3 = 2**-5, which places the condition
+ * into the [T1, T3] interval.
+ */
+static const double
+T1 = -0.1659, /* ~-30.625/128 * log(2) */
+T2 = 0.1659, /* ~30.625/128 * log(2) */
+T3 = 0.03125;
+
+/*
+ * Domain [-0.1659, 0.03125], range ~[2.9134e-44, 1.8404e-37]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-122.03
+ */
+static const long double
+C3 = 1.66666666666666666666666666666666667e-1L,
+C4 = 4.16666666666666666666666666666666645e-2L,
+C5 = 8.33333333333333333333333333333371638e-3L,
+C6 = 1.38888888888888888888888888891188658e-3L,
+C7 = 1.98412698412698412698412697235950394e-4L,
+C8 = 2.48015873015873015873015112487849040e-5L,
+C9 = 2.75573192239858906525606685484412005e-6L,
+C10 = 2.75573192239858906612966093057020362e-7L,
+C11 = 2.50521083854417203619031960151253944e-8L,
+C12 = 2.08767569878679576457272282566520649e-9L,
+C13 = 1.60590438367252471783548748824255707e-10L;
+
+static const double
+C14 = 1.1470745580491932e-11, /* 0x1.93974a81dae30p-37 */
+C15 = 7.6471620181090468e-13, /* 0x1.ae7f3820adab1p-41 */
+C16 = 4.7793721460260450e-14, /* 0x1.ae7cd18a18eacp-45 */
+C17 = 2.8074757356658877e-15, /* 0x1.949992a1937d9p-49 */
+C18 = 1.4760610323699476e-16; /* 0x1.545b43aabfbcdp-53 */
+
+/*
+ * Domain [0.03125, 0.1659], range ~[-2.7676e-37, -1.0367e-38]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-121.44
+ */
+static const long double
+D3 = 1.66666666666666666666666666666682245e-1L,
+D4 = 4.16666666666666666666666666634228324e-2L,
+D5 = 8.33333333333333333333333364022244481e-3L,
+D6 = 1.38888888888888888888887138722762072e-3L,
+D7 = 1.98412698412698412699085805424661471e-4L,
+D8 = 2.48015873015873015687993712101479612e-5L,
+D9 = 2.75573192239858944101036288338208042e-6L,
+D10 = 2.75573192239853161148064676533754048e-7L,
+D11 = 2.50521083855084570046480450935267433e-8L,
+D12 = 2.08767569819738524488686318024854942e-9L,
+D13 = 1.60590442297008495301927448122499313e-10L;
+
+static const double
+D14 = 1.1470726176204336e-11, /* 0x1.93971dc395d9ep-37 */
+D15 = 7.6478532249581686e-13, /* 0x1.ae892e3D16fcep-41 */
+D16 = 4.7628892832607741e-14, /* 0x1.ad00Dfe41feccp-45 */
+D17 = 3.0524857220358650e-15; /* 0x1.D7e8d886Df921p-49 */
+
+long double
+expm1l(long double x)
+{
+ union IEEEl2bits u, v;
+ long double hx2_hi, hx2_lo, q, r, r1, t, twomk, twopk, x_hi;
+ long double x_lo, x2;
+ double dr, dx, fn, r2;
+ int k, n, n2;
+ uint16_t hx, ix;
+
+ /* Filter out exceptional cases. */
+ u.e = x;
+ hx = u.xbits.expsign;
+ ix = hx & 0x7fff;
+ if (ix >= BIAS + 7) { /* |x| >= 128 or x is NaN */
+ if (ix == BIAS + LDBL_MAX_EXP) {
+ if (hx & 0x8000) /* x is -Inf or -NaN */
+ return (-1 / x - 1);
+ return (x + x); /* x is +Inf or +NaN */
+ }
+ if (x > o_threshold)
+ return (huge * huge);
+ /*
+ * expm1l() never underflows, but it must avoid
+ * unrepresentable large negative exponents. We used a
+ * much smaller threshold for large |x| above than in
+ * expl() so as to handle not so large negative exponents
+ * in the same way as large ones here.
+ */
+ if (hx & 0x8000) /* x <= -128 */
+ return (tiny - 1); /* good for x < -114ln2 - eps */
+ }
+
+ ENTERI();
+
+ if (T1 < x && x < T2) {
+ x2 = x * x;
+ dx = x;
+
+ if (x < T3) {
+ if (ix < BIAS - 113) { /* |x| < 0x1p-113 */
+ /* x (rounded) with inexact if x != 0: */
+ RETURNI(x == 0 ? x :
+ (0x1p200 * x + fabsl(x)) * 0x1p-200);
+ }
+ q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +
+ x * (C7 + x * (C8 + x * (C9 + x * (C10 +
+ x * (C11 + x * (C12 + x * (C13 +
+ dx * (C14 + dx * (C15 + dx * (C16 +
+ dx * (C17 + dx * C18))))))))))))));
+ } else {
+ q = x * x2 * D3 + x2 * x2 * (D4 + x * (D5 + x * (D6 +
+ x * (D7 + x * (D8 + x * (D9 + x * (D10 +
+ x * (D11 + x * (D12 + x * (D13 +
+ dx * (D14 + dx * (D15 + dx * (D16 +
+ dx * D17)))))))))))));
+ }
+
+ x_hi = (float)x;
+ x_lo = x - x_hi;
+ hx2_hi = x_hi * x_hi / 2;
+ hx2_lo = x_lo * (x + x_hi) / 2;
+ if (ix >= BIAS - 7)
+ RETURNI(hx2_lo + x_lo + q + (hx2_hi + x_hi));
+ else
+ RETURNI(hx2_lo + q + hx2_hi + x);
+ }
+
+ /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
+ /* Use a specialized rint() to get fn. Assume round-to-nearest. */
+ fn = (double)x * INV_L + 0x1.8p52 - 0x1.8p52;
+#if defined(HAVE_EFFICIENT_IRINT)
+ n = irint(fn);
+#else
+ n = (int)fn;
+#endif
+ n2 = (unsigned)n % INTERVALS;
+ k = n >> LOG2_INTERVALS;
+ r1 = x - fn * L1;
+ r2 = fn * -L2;
+ r = r1 + r2;
+
+ /* Prepare scale factor. */
+ v.e = 1;
+ v.xbits.expsign = BIAS + k;
+ twopk = v.e;
+
+ /*
+ * Evaluate lower terms of
+ * expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2).
+ */
+ dr = r;
+ q = r2 + r * r * (A2 + r * (A3 + r * (A4 + r * (A5 + r * (A6 +
+ dr * (A7 + dr * (A8 + dr * (A9 + dr * A10))))))));
+
+ t = tbl[n2].lo + tbl[n2].hi;
+
+ if (k == 0) {
+ t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 +
+ (tbl[n2].hi - 1);
+ RETURNI(t);
+ }
+ if (k == -1) {
+ t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 +
+ (tbl[n2].hi - 2);
+ RETURNI(t / 2);
+ }
+ if (k < -7) {
+ t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+ RETURNI(t * twopk - 1);
+ }
+ if (k > 2 * LDBL_MANT_DIG - 1) {
+ t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+ if (k == LDBL_MAX_EXP)
+ RETURNI(t * 2 * 0x1p16383L - 1);
+ RETURNI(t * twopk - 1);
+ }
+
+ v.xbits.expsign = BIAS - k;
+ twomk = v.e;
+
+ if (k > LDBL_MANT_DIG - 1)
+ t = tbl[n2].lo - twomk + t * (q + r1) + tbl[n2].hi;
+ else
+ t = tbl[n2].lo + t * (q + r1) + (tbl[n2].hi - twomk);
+ RETURNI(t * twopk);
+}
diff --git a/libm/upstream-freebsd/lib/msun/ld128/s_logl.c b/libm/upstream-freebsd/lib/msun/ld128/s_logl.c
new file mode 100644
index 0000000..391d623
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/s_logl.c
@@ -0,0 +1,737 @@
+/*-
+ * Copyright (c) 2007-2013 Bruce D. Evans
+ * 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 unmodified, 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 ``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 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>
+__FBSDID("$FreeBSD$");
+
+/**
+ * Implementation of the natural logarithm of x for 128-bit format.
+ *
+ * First decompose x into its base 2 representation:
+ *
+ * log(x) = log(X * 2**k), where X is in [1, 2)
+ * = log(X) + k * log(2).
+ *
+ * Let X = X_i + e, where X_i is the center of one of the intervals
+ * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256)
+ * and X is in this interval. Then
+ *
+ * log(X) = log(X_i + e)
+ * = log(X_i * (1 + e / X_i))
+ * = log(X_i) + log(1 + e / X_i).
+ *
+ * The values log(X_i) are tabulated below. Let d = e / X_i and use
+ *
+ * log(1 + d) = p(d)
+ *
+ * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of
+ * suitably high degree.
+ *
+ * To get sufficiently small roundoff errors, k * log(2), log(X_i), and
+ * sometimes (if |k| is not large) the first term in p(d) must be evaluated
+ * and added up in extra precision. Extra precision is not needed for the
+ * rest of p(d). In the worst case when k = 0 and log(X_i) is 0, the final
+ * error is controlled mainly by the error in the second term in p(d). The
+ * error in this term itself is at most 0.5 ulps from the d*d operation in
+ * it. The error in this term relative to the first term is thus at most
+ * 0.5 * |-0.5| * |d| < 1.0/1024 ulps. We aim for an accumulated error of
+ * at most twice this at the point of the final rounding step. Thus the
+ * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps. Exhaustive
+ * testing of a float variant of this function showed a maximum final error
+ * of 0.5008 ulps. Non-exhaustive testing of a double variant of this
+ * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256).
+ *
+ * We made the maximum of |d| (and thus the total relative error and the
+ * degree of p(d)) small by using a large number of intervals. Using
+ * centers of intervals instead of endpoints reduces this maximum by a
+ * factor of 2 for a given number of intervals. p(d) is special only
+ * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen
+ * naturally. The most accurate minimax polynomial of a given degree might
+ * be different, but then we wouldn't want it since we would have to do
+ * extra work to avoid roundoff error (especially for P0*d instead of d).
+ */
+
+#ifdef DEBUG
+#include <assert.h>
+#include <fenv.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#ifndef NO_STRUCT_RETURN
+#define STRUCT_RETURN
+#endif
+#include "math_private.h"
+
+#if !defined(NO_UTAB) && !defined(NO_UTABL)
+#define USE_UTAB
+#endif
+
+/*
+ * Domain [-0.005280, 0.004838], range ~[-1.1577e-37, 1.1582e-37]:
+ * |log(1 + d)/d - p(d)| < 2**-122.7
+ */
+static const long double
+P2 = -0.5L,
+P3 = 3.33333333333333333333333333333233795e-1L, /* 0x15555555555555555555555554d42.0p-114L */
+P4 = -2.49999999999999999999999999941139296e-1L, /* -0x1ffffffffffffffffffffffdab14e.0p-115L */
+P5 = 2.00000000000000000000000085468039943e-1L, /* 0x19999999999999999999a6d3567f4.0p-115L */
+P6 = -1.66666666666666666666696142372698408e-1L, /* -0x15555555555555555567267a58e13.0p-115L */
+P7 = 1.42857142857142857119522943477166120e-1L, /* 0x1249249249249248ed79a0ae434de.0p-115L */
+P8 = -1.24999999999999994863289015033581301e-1L; /* -0x1fffffffffffffa13e91765e46140.0p-116L */
+/* Double precision gives ~ 53 + log2(P9 * max(|d|)**8) ~= 120 bits. */
+static const double
+P9 = 1.1111111111111401e-1, /* 0x1c71c71c71c7ed.0p-56 */
+P10 = -1.0000000000040135e-1, /* -0x199999999a0a92.0p-56 */
+P11 = 9.0909090728136258e-2, /* 0x1745d173962111.0p-56 */
+P12 = -8.3333318851855284e-2, /* -0x1555551722c7a3.0p-56 */
+P13 = 7.6928634666404178e-2, /* 0x13b1985204a4ae.0p-56 */
+P14 = -7.1626810078462499e-2; /* -0x12562276cdc5d0.0p-56 */
+
+static volatile const double zero = 0;
+
+#define INTERVALS 128
+#define LOG2_INTERVALS 7
+#define TSIZE (INTERVALS + 1)
+#define G(i) (T[(i)].G)
+#define F_hi(i) (T[(i)].F_hi)
+#define F_lo(i) (T[(i)].F_lo)
+#define ln2_hi F_hi(TSIZE - 1)
+#define ln2_lo F_lo(TSIZE - 1)
+#define E(i) (U[(i)].E)
+#define H(i) (U[(i)].H)
+
+static const struct {
+ float G; /* 1/(1 + i/128) rounded to 8/9 bits */
+ float F_hi; /* log(1 / G_i) rounded (see below) */
+ /* The compiler will insert 8 bytes of padding here. */
+ long double F_lo; /* next 113 bits for log(1 / G_i) */
+} T[TSIZE] = {
+ /*
+ * ln2_hi and each F_hi(i) are rounded to a number of bits that
+ * makes F_hi(i) + dk*ln2_hi exact for all i and all dk.
+ *
+ * The last entry (for X just below 2) is used to define ln2_hi
+ * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly
+ * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1.
+ * This is needed for accuracy when x is just below 1. (To avoid
+ * special cases, such x are "reduced" strangely to X just below
+ * 2 and dk = -1, and then the exact cancellation is needed
+ * because any the error from any non-exactness would be too
+ * large).
+ *
+ * The relevant range of dk is [-16445, 16383]. The maximum number
+ * of bits in F_hi(i) that works is very dependent on i but has
+ * a minimum of 93. We only need about 12 bits in F_hi(i) for
+ * it to provide enough extra precision.
+ *
+ * We round F_hi(i) to 24 bits so that it can have type float,
+ * mainly to minimize the size of the table. Using all 24 bits
+ * in a float for it automatically satisfies the above constraints.
+ */
+ 0x800000.0p-23, 0, 0,
+ 0xfe0000.0p-24, 0x8080ac.0p-30, -0x14ee431dae6674afa0c4bfe16e8fd.0p-144L,
+ 0xfc0000.0p-24, 0x8102b3.0p-29, -0x1db29ee2d83717be918e1119642ab.0p-144L,
+ 0xfa0000.0p-24, 0xc24929.0p-29, 0x1191957d173697cf302cc9476f561.0p-143L,
+ 0xf80000.0p-24, 0x820aec.0p-28, 0x13ce8888e02e78eba9b1113bc1c18.0p-142L,
+ 0xf60000.0p-24, 0xa33577.0p-28, -0x17a4382ce6eb7bfa509bec8da5f22.0p-142L,
+ 0xf48000.0p-24, 0xbc42cb.0p-28, -0x172a21161a107674986dcdca6709c.0p-143L,
+ 0xf30000.0p-24, 0xd57797.0p-28, -0x1e09de07cb958897a3ea46e84abb3.0p-142L,
+ 0xf10000.0p-24, 0xf7518e.0p-28, 0x1ae1eec1b036c484993c549c4bf40.0p-151L,
+ 0xef0000.0p-24, 0x8cb9df.0p-27, -0x1d7355325d560d9e9ab3d6ebab580.0p-141L,
+ 0xed8000.0p-24, 0x999ec0.0p-27, -0x1f9f02d256d5037108f4ec21e48cd.0p-142L,
+ 0xec0000.0p-24, 0xa6988b.0p-27, -0x16fc0a9d12c17a70f7a684c596b12.0p-143L,
+ 0xea0000.0p-24, 0xb80698.0p-27, 0x15d581c1e8da99ded322fb08b8462.0p-141L,
+ 0xe80000.0p-24, 0xc99af3.0p-27, -0x1535b3ba8f150ae09996d7bb4653e.0p-143L,
+ 0xe70000.0p-24, 0xd273b2.0p-27, 0x163786f5251aefe0ded34c8318f52.0p-145L,
+ 0xe50000.0p-24, 0xe442c0.0p-27, 0x1bc4b2368e32d56699c1799a244d4.0p-144L,
+ 0xe38000.0p-24, 0xf1b83f.0p-27, 0x1c6090f684e6766abceccab1d7174.0p-141L,
+ 0xe20000.0p-24, 0xff448a.0p-27, -0x1890aa69ac9f4215f93936b709efb.0p-142L,
+ 0xe08000.0p-24, 0x8673f6.0p-26, 0x1b9985194b6affd511b534b72a28e.0p-140L,
+ 0xdf0000.0p-24, 0x8d515c.0p-26, -0x1dc08d61c6ef1d9b2ef7e68680598.0p-143L,
+ 0xdd8000.0p-24, 0x943a9e.0p-26, -0x1f72a2dac729b3f46662238a9425a.0p-142L,
+ 0xdc0000.0p-24, 0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9691aed4d5e3df94.0p-140L,
+ 0xda8000.0p-24, 0xa2315d.0p-26, -0x11b26121629c46c186384993e1c93.0p-142L,
+ 0xd90000.0p-24, 0xa93f2f.0p-26, 0x1286d633e8e5697dc6a402a56fce1.0p-141L,
+ 0xd78000.0p-24, 0xb05988.0p-26, 0x16128eba9367707ebfa540e45350c.0p-144L,
+ 0xd60000.0p-24, 0xb78094.0p-26, 0x16ead577390d31ef0f4c9d43f79b2.0p-140L,
+ 0xd50000.0p-24, 0xbc4c6c.0p-26, 0x151131ccf7c7b75e7d900b521c48d.0p-141L,
+ 0xd38000.0p-24, 0xc3890a.0p-26, -0x115e2cd714bd06508aeb00d2ae3e9.0p-140L,
+ 0xd20000.0p-24, 0xcad2d7.0p-26, -0x1847f406ebd3af80485c2f409633c.0p-142L,
+ 0xd10000.0p-24, 0xcfb620.0p-26, 0x1c2259904d686581799fbce0b5f19.0p-141L,
+ 0xcf8000.0p-24, 0xd71653.0p-26, 0x1ece57a8d5ae54f550444ecf8b995.0p-140L,
+ 0xce0000.0p-24, 0xde843a.0p-26, -0x1f109d4bc4595412b5d2517aaac13.0p-141L,
+ 0xcd0000.0p-24, 0xe37fde.0p-26, 0x1bc03dc271a74d3a85b5b43c0e727.0p-141L,
+ 0xcb8000.0p-24, 0xeb050c.0p-26, -0x1bf2badc0df841a71b79dd5645b46.0p-145L,
+ 0xca0000.0p-24, 0xf29878.0p-26, -0x18efededd89fbe0bcfbe6d6db9f66.0p-147L,
+ 0xc90000.0p-24, 0xf7ad6f.0p-26, 0x1373ff977baa6911c7bafcb4d84fb.0p-141L,
+ 0xc80000.0p-24, 0xfcc8e3.0p-26, 0x196766f2fb328337cc050c6d83b22.0p-140L,
+ 0xc68000.0p-24, 0x823f30.0p-25, 0x19bd076f7c434e5fcf1a212e2a91e.0p-139L,
+ 0xc58000.0p-24, 0x84d52c.0p-25, -0x1a327257af0f465e5ecab5f2a6f81.0p-139L,
+ 0xc40000.0p-24, 0x88bc74.0p-25, 0x113f23def19c5a0fe396f40f1dda9.0p-141L,
+ 0xc30000.0p-24, 0x8b5ae6.0p-25, 0x1759f6e6b37de945a049a962e66c6.0p-139L,
+ 0xc20000.0p-24, 0x8dfccb.0p-25, 0x1ad35ca6ed5147bdb6ddcaf59c425.0p-141L,
+ 0xc10000.0p-24, 0x90a22b.0p-25, 0x1a1d71a87deba46bae9827221dc98.0p-139L,
+ 0xbf8000.0p-24, 0x94a0d8.0p-25, -0x139e5210c2b730e28aba001a9b5e0.0p-140L,
+ 0xbe8000.0p-24, 0x974f16.0p-25, -0x18f6ebcff3ed72e23e13431adc4a5.0p-141L,
+ 0xbd8000.0p-24, 0x9a00f1.0p-25, -0x1aa268be39aab7148e8d80caa10b7.0p-139L,
+ 0xbc8000.0p-24, 0x9cb672.0p-25, -0x14c8815839c5663663d15faed7771.0p-139L,
+ 0xbb0000.0p-24, 0xa0cda1.0p-25, 0x1eaf46390dbb2438273918db7df5c.0p-141L,
+ 0xba0000.0p-24, 0xa38c6e.0p-25, 0x138e20d831f698298adddd7f32686.0p-141L,
+ 0xb90000.0p-24, 0xa64f05.0p-25, -0x1e8d3c41123615b147a5d47bc208f.0p-142L,
+ 0xb80000.0p-24, 0xa91570.0p-25, 0x1ce28f5f3840b263acb4351104631.0p-140L,
+ 0xb70000.0p-24, 0xabdfbb.0p-25, -0x186e5c0a42423457e22d8c650b355.0p-139L,
+ 0xb60000.0p-24, 0xaeadef.0p-25, -0x14d41a0b2a08a465dc513b13f567d.0p-143L,
+ 0xb50000.0p-24, 0xb18018.0p-25, 0x16755892770633947ffe651e7352f.0p-139L,
+ 0xb40000.0p-24, 0xb45642.0p-25, -0x16395ebe59b15228bfe8798d10ff0.0p-142L,
+ 0xb30000.0p-24, 0xb73077.0p-25, 0x1abc65c8595f088b61a335f5b688c.0p-140L,
+ 0xb20000.0p-24, 0xba0ec4.0p-25, -0x1273089d3dad88e7d353e9967d548.0p-139L,
+ 0xb10000.0p-24, 0xbcf133.0p-25, 0x10f9f67b1f4bbf45de06ecebfaf6d.0p-139L,
+ 0xb00000.0p-24, 0xbfd7d2.0p-25, -0x109fab904864092b34edda19a831e.0p-140L,
+ 0xaf0000.0p-24, 0xc2c2ac.0p-25, -0x1124680aa43333221d8a9b475a6ba.0p-139L,
+ 0xae8000.0p-24, 0xc439b3.0p-25, -0x1f360cc4710fbfe24b633f4e8d84d.0p-140L,
+ 0xad8000.0p-24, 0xc72afd.0p-25, -0x132d91f21d89c89c45003fc5d7807.0p-140L,
+ 0xac8000.0p-24, 0xca20a2.0p-25, -0x16bf9b4d1f8da8002f2449e174504.0p-139L,
+ 0xab8000.0p-24, 0xcd1aae.0p-25, 0x19deb5ce6a6a8717d5626e16acc7d.0p-141L,
+ 0xaa8000.0p-24, 0xd0192f.0p-25, 0x1a29fb48f7d3ca87dabf351aa41f4.0p-139L,
+ 0xaa0000.0p-24, 0xd19a20.0p-25, 0x1127d3c6457f9d79f51dcc73014c9.0p-141L,
+ 0xa90000.0p-24, 0xd49f6a.0p-25, -0x1ba930e486a0ac42d1bf9199188e7.0p-141L,
+ 0xa80000.0p-24, 0xd7a94b.0p-25, -0x1b6e645f31549dd1160bcc45c7e2c.0p-139L,
+ 0xa70000.0p-24, 0xdab7d0.0p-25, 0x1118a425494b610665377f15625b6.0p-140L,
+ 0xa68000.0p-24, 0xdc40d5.0p-25, 0x1966f24d29d3a2d1b2176010478be.0p-140L,
+ 0xa58000.0p-24, 0xdf566d.0p-25, -0x1d8e52eb2248f0c95dd83626d7333.0p-142L,
+ 0xa48000.0p-24, 0xe270ce.0p-25, -0x1ee370f96e6b67ccb006a5b9890ea.0p-140L,
+ 0xa40000.0p-24, 0xe3ffce.0p-25, 0x1d155324911f56db28da4d629d00a.0p-140L,
+ 0xa30000.0p-24, 0xe72179.0p-25, -0x1fe6e2f2f867d8f4d60c713346641.0p-140L,
+ 0xa20000.0p-24, 0xea4812.0p-25, 0x1b7be9add7f4d3b3d406b6cbf3ce5.0p-140L,
+ 0xa18000.0p-24, 0xebdd3d.0p-25, 0x1b3cfb3f7511dd73692609040ccc2.0p-139L,
+ 0xa08000.0p-24, 0xef0b5b.0p-25, -0x1220de1f7301901b8ad85c25afd09.0p-139L,
+ 0xa00000.0p-24, 0xf0a451.0p-25, -0x176364c9ac81cc8a4dfb804de6867.0p-140L,
+ 0x9f0000.0p-24, 0xf3da16.0p-25, 0x1eed6b9aafac8d42f78d3e65d3727.0p-141L,
+ 0x9e8000.0p-24, 0xf576e9.0p-25, 0x1d593218675af269647b783d88999.0p-139L,
+ 0x9d8000.0p-24, 0xf8b47c.0p-25, -0x13e8eb7da053e063714615f7cc91d.0p-144L,
+ 0x9d0000.0p-24, 0xfa553f.0p-25, 0x1c063259bcade02951686d5373aec.0p-139L,
+ 0x9c0000.0p-24, 0xfd9ac5.0p-25, 0x1ef491085fa3c1649349630531502.0p-139L,
+ 0x9b8000.0p-24, 0xff3f8c.0p-25, 0x1d607a7c2b8c5320619fb9433d841.0p-139L,
+ 0x9a8000.0p-24, 0x814697.0p-24, -0x12ad3817004f3f0bdff99f932b273.0p-138L,
+ 0x9a0000.0p-24, 0x821b06.0p-24, -0x189fc53117f9e54e78103a2bc1767.0p-141L,
+ 0x990000.0p-24, 0x83c5f8.0p-24, 0x14cf15a048907b7d7f47ddb45c5a3.0p-139L,
+ 0x988000.0p-24, 0x849c7d.0p-24, 0x1cbb1d35fb82873b04a9af1dd692c.0p-138L,
+ 0x978000.0p-24, 0x864ba6.0p-24, 0x1128639b814f9b9770d8cb6573540.0p-138L,
+ 0x970000.0p-24, 0x87244c.0p-24, 0x184733853300f002e836dfd47bd41.0p-139L,
+ 0x968000.0p-24, 0x87fdaa.0p-24, 0x109d23aef77dd5cd7cc94306fb3ff.0p-140L,
+ 0x958000.0p-24, 0x89b293.0p-24, -0x1a81ef367a59de2b41eeebd550702.0p-138L,
+ 0x950000.0p-24, 0x8a8e20.0p-24, -0x121ad3dbb2f45275c917a30df4ac9.0p-138L,
+ 0x948000.0p-24, 0x8b6a6a.0p-24, -0x1cfb981628af71a89df4e6df2e93b.0p-139L,
+ 0x938000.0p-24, 0x8d253a.0p-24, -0x1d21730ea76cfdec367828734cae5.0p-139L,
+ 0x930000.0p-24, 0x8e03c2.0p-24, 0x135cc00e566f76b87333891e0dec4.0p-138L,
+ 0x928000.0p-24, 0x8ee30d.0p-24, -0x10fcb5df257a263e3bf446c6e3f69.0p-140L,
+ 0x918000.0p-24, 0x90a3ee.0p-24, -0x16e171b15433d723a4c7380a448d8.0p-139L,
+ 0x910000.0p-24, 0x918587.0p-24, -0x1d050da07f3236f330972da2a7a87.0p-139L,
+ 0x908000.0p-24, 0x9267e7.0p-24, 0x1be03669a5268d21148c6002becd3.0p-139L,
+ 0x8f8000.0p-24, 0x942f04.0p-24, 0x10b28e0e26c336af90e00533323ba.0p-139L,
+ 0x8f0000.0p-24, 0x9513c3.0p-24, 0x1a1d820da57cf2f105a89060046aa.0p-138L,
+ 0x8e8000.0p-24, 0x95f950.0p-24, -0x19ef8f13ae3cf162409d8ea99d4c0.0p-139L,
+ 0x8e0000.0p-24, 0x96dfab.0p-24, -0x109e417a6e507b9dc10dac743ad7a.0p-138L,
+ 0x8d0000.0p-24, 0x98aed2.0p-24, 0x10d01a2c5b0e97c4990b23d9ac1f5.0p-139L,
+ 0x8c8000.0p-24, 0x9997a2.0p-24, -0x1d6a50d4b61ea74540bdd2aa99a42.0p-138L,
+ 0x8c0000.0p-24, 0x9a8145.0p-24, 0x1b3b190b83f9527e6aba8f2d783c1.0p-138L,
+ 0x8b8000.0p-24, 0x9b6bbf.0p-24, 0x13a69fad7e7abe7ba81c664c107e0.0p-138L,
+ 0x8b0000.0p-24, 0x9c5711.0p-24, -0x11cd12316f576aad348ae79867223.0p-138L,
+ 0x8a8000.0p-24, 0x9d433b.0p-24, 0x1c95c444b807a246726b304ccae56.0p-139L,
+ 0x898000.0p-24, 0x9f1e22.0p-24, -0x1b9c224ea698c2f9b47466d6123fe.0p-139L,
+ 0x890000.0p-24, 0xa00ce1.0p-24, 0x125ca93186cf0f38b4619a2483399.0p-141L,
+ 0x888000.0p-24, 0xa0fc80.0p-24, -0x1ee38a7bc228b3597043be78eaf49.0p-139L,
+ 0x880000.0p-24, 0xa1ed00.0p-24, -0x1a0db876613d204147dc69a07a649.0p-138L,
+ 0x878000.0p-24, 0xa2de62.0p-24, 0x193224e8516c008d3602a7b41c6e8.0p-139L,
+ 0x870000.0p-24, 0xa3d0a9.0p-24, 0x1fa28b4d2541aca7d5844606b2421.0p-139L,
+ 0x868000.0p-24, 0xa4c3d6.0p-24, 0x1c1b5760fb4571acbcfb03f16daf4.0p-138L,
+ 0x858000.0p-24, 0xa6acea.0p-24, 0x1fed5d0f65949c0a345ad743ae1ae.0p-140L,
+ 0x850000.0p-24, 0xa7a2d4.0p-24, 0x1ad270c9d749362382a7688479e24.0p-140L,
+ 0x848000.0p-24, 0xa899ab.0p-24, 0x199ff15ce532661ea9643a3a2d378.0p-139L,
+ 0x840000.0p-24, 0xa99171.0p-24, 0x1a19e15ccc45d257530a682b80490.0p-139L,
+ 0x838000.0p-24, 0xaa8a28.0p-24, -0x121a14ec532b35ba3e1f868fd0b5e.0p-140L,
+ 0x830000.0p-24, 0xab83d1.0p-24, 0x1aee319980bff3303dd481779df69.0p-139L,
+ 0x828000.0p-24, 0xac7e6f.0p-24, -0x18ffd9e3900345a85d2d86161742e.0p-140L,
+ 0x820000.0p-24, 0xad7a03.0p-24, -0x1e4db102ce29f79b026b64b42caa1.0p-140L,
+ 0x818000.0p-24, 0xae768f.0p-24, 0x17c35c55a04a82ab19f77652d977a.0p-141L,
+ 0x810000.0p-24, 0xaf7415.0p-24, 0x1448324047019b48d7b98c1cf7234.0p-138L,
+ 0x808000.0p-24, 0xb07298.0p-24, -0x1750ee3915a197e9c7359dd94152f.0p-138L,
+ 0x800000.0p-24, 0xb17218.0p-24, -0x105c610ca86c3898cff81a12a17e2.0p-141L,
+};
+
+#ifdef USE_UTAB
+static const struct {
+ float H; /* 1 + i/INTERVALS (exact) */
+ float E; /* H(i) * G(i) - 1 (exact) */
+} U[TSIZE] = {
+ 0x800000.0p-23, 0,
+ 0x810000.0p-23, -0x800000.0p-37,
+ 0x820000.0p-23, -0x800000.0p-35,
+ 0x830000.0p-23, -0x900000.0p-34,
+ 0x840000.0p-23, -0x800000.0p-33,
+ 0x850000.0p-23, -0xc80000.0p-33,
+ 0x860000.0p-23, -0xa00000.0p-36,
+ 0x870000.0p-23, 0x940000.0p-33,
+ 0x880000.0p-23, 0x800000.0p-35,
+ 0x890000.0p-23, -0xc80000.0p-34,
+ 0x8a0000.0p-23, 0xe00000.0p-36,
+ 0x8b0000.0p-23, 0x900000.0p-33,
+ 0x8c0000.0p-23, -0x800000.0p-35,
+ 0x8d0000.0p-23, -0xe00000.0p-33,
+ 0x8e0000.0p-23, 0x880000.0p-33,
+ 0x8f0000.0p-23, -0xa80000.0p-34,
+ 0x900000.0p-23, -0x800000.0p-35,
+ 0x910000.0p-23, 0x800000.0p-37,
+ 0x920000.0p-23, 0x900000.0p-35,
+ 0x930000.0p-23, 0xd00000.0p-35,
+ 0x940000.0p-23, 0xe00000.0p-35,
+ 0x950000.0p-23, 0xc00000.0p-35,
+ 0x960000.0p-23, 0xe00000.0p-36,
+ 0x970000.0p-23, -0x800000.0p-38,
+ 0x980000.0p-23, -0xc00000.0p-35,
+ 0x990000.0p-23, -0xd00000.0p-34,
+ 0x9a0000.0p-23, 0x880000.0p-33,
+ 0x9b0000.0p-23, 0xe80000.0p-35,
+ 0x9c0000.0p-23, -0x800000.0p-35,
+ 0x9d0000.0p-23, 0xb40000.0p-33,
+ 0x9e0000.0p-23, 0x880000.0p-34,
+ 0x9f0000.0p-23, -0xe00000.0p-35,
+ 0xa00000.0p-23, 0x800000.0p-33,
+ 0xa10000.0p-23, -0x900000.0p-36,
+ 0xa20000.0p-23, -0xb00000.0p-33,
+ 0xa30000.0p-23, -0xa00000.0p-36,
+ 0xa40000.0p-23, 0x800000.0p-33,
+ 0xa50000.0p-23, -0xf80000.0p-35,
+ 0xa60000.0p-23, 0x880000.0p-34,
+ 0xa70000.0p-23, -0x900000.0p-33,
+ 0xa80000.0p-23, -0x800000.0p-35,
+ 0xa90000.0p-23, 0x900000.0p-34,
+ 0xaa0000.0p-23, 0xa80000.0p-33,
+ 0xab0000.0p-23, -0xac0000.0p-34,
+ 0xac0000.0p-23, -0x800000.0p-37,
+ 0xad0000.0p-23, 0xf80000.0p-35,
+ 0xae0000.0p-23, 0xf80000.0p-34,
+ 0xaf0000.0p-23, -0xac0000.0p-33,
+ 0xb00000.0p-23, -0x800000.0p-33,
+ 0xb10000.0p-23, -0xb80000.0p-34,
+ 0xb20000.0p-23, -0x800000.0p-34,
+ 0xb30000.0p-23, -0xb00000.0p-35,
+ 0xb40000.0p-23, -0x800000.0p-35,
+ 0xb50000.0p-23, -0xe00000.0p-36,
+ 0xb60000.0p-23, -0x800000.0p-35,
+ 0xb70000.0p-23, -0xb00000.0p-35,
+ 0xb80000.0p-23, -0x800000.0p-34,
+ 0xb90000.0p-23, -0xb80000.0p-34,
+ 0xba0000.0p-23, -0x800000.0p-33,
+ 0xbb0000.0p-23, -0xac0000.0p-33,
+ 0xbc0000.0p-23, 0x980000.0p-33,
+ 0xbd0000.0p-23, 0xbc0000.0p-34,
+ 0xbe0000.0p-23, 0xe00000.0p-36,
+ 0xbf0000.0p-23, -0xb80000.0p-35,
+ 0xc00000.0p-23, -0x800000.0p-33,
+ 0xc10000.0p-23, 0xa80000.0p-33,
+ 0xc20000.0p-23, 0x900000.0p-34,
+ 0xc30000.0p-23, -0x800000.0p-35,
+ 0xc40000.0p-23, -0x900000.0p-33,
+ 0xc50000.0p-23, 0x820000.0p-33,
+ 0xc60000.0p-23, 0x800000.0p-38,
+ 0xc70000.0p-23, -0x820000.0p-33,
+ 0xc80000.0p-23, 0x800000.0p-33,
+ 0xc90000.0p-23, -0xa00000.0p-36,
+ 0xca0000.0p-23, -0xb00000.0p-33,
+ 0xcb0000.0p-23, 0x840000.0p-34,
+ 0xcc0000.0p-23, -0xd00000.0p-34,
+ 0xcd0000.0p-23, 0x800000.0p-33,
+ 0xce0000.0p-23, -0xe00000.0p-35,
+ 0xcf0000.0p-23, 0xa60000.0p-33,
+ 0xd00000.0p-23, -0x800000.0p-35,
+ 0xd10000.0p-23, 0xb40000.0p-33,
+ 0xd20000.0p-23, -0x800000.0p-35,
+ 0xd30000.0p-23, 0xaa0000.0p-33,
+ 0xd40000.0p-23, -0xe00000.0p-35,
+ 0xd50000.0p-23, 0x880000.0p-33,
+ 0xd60000.0p-23, -0xd00000.0p-34,
+ 0xd70000.0p-23, 0x9c0000.0p-34,
+ 0xd80000.0p-23, -0xb00000.0p-33,
+ 0xd90000.0p-23, -0x800000.0p-38,
+ 0xda0000.0p-23, 0xa40000.0p-33,
+ 0xdb0000.0p-23, -0xdc0000.0p-34,
+ 0xdc0000.0p-23, 0xc00000.0p-35,
+ 0xdd0000.0p-23, 0xca0000.0p-33,
+ 0xde0000.0p-23, -0xb80000.0p-34,
+ 0xdf0000.0p-23, 0xd00000.0p-35,
+ 0xe00000.0p-23, 0xc00000.0p-33,
+ 0xe10000.0p-23, -0xf40000.0p-34,
+ 0xe20000.0p-23, 0x800000.0p-37,
+ 0xe30000.0p-23, 0x860000.0p-33,
+ 0xe40000.0p-23, -0xc80000.0p-33,
+ 0xe50000.0p-23, -0xa80000.0p-34,
+ 0xe60000.0p-23, 0xe00000.0p-36,
+ 0xe70000.0p-23, 0x880000.0p-33,
+ 0xe80000.0p-23, -0xe00000.0p-33,
+ 0xe90000.0p-23, -0xfc0000.0p-34,
+ 0xea0000.0p-23, -0x800000.0p-35,
+ 0xeb0000.0p-23, 0xe80000.0p-35,
+ 0xec0000.0p-23, 0x900000.0p-33,
+ 0xed0000.0p-23, 0xe20000.0p-33,
+ 0xee0000.0p-23, -0xac0000.0p-33,
+ 0xef0000.0p-23, -0xc80000.0p-34,
+ 0xf00000.0p-23, -0x800000.0p-35,
+ 0xf10000.0p-23, 0x800000.0p-35,
+ 0xf20000.0p-23, 0xb80000.0p-34,
+ 0xf30000.0p-23, 0x940000.0p-33,
+ 0xf40000.0p-23, 0xc80000.0p-33,
+ 0xf50000.0p-23, -0xf20000.0p-33,
+ 0xf60000.0p-23, -0xc80000.0p-33,
+ 0xf70000.0p-23, -0xa20000.0p-33,
+ 0xf80000.0p-23, -0x800000.0p-33,
+ 0xf90000.0p-23, -0xc40000.0p-34,
+ 0xfa0000.0p-23, -0x900000.0p-34,
+ 0xfb0000.0p-23, -0xc80000.0p-35,
+ 0xfc0000.0p-23, -0x800000.0p-35,
+ 0xfd0000.0p-23, -0x900000.0p-36,
+ 0xfe0000.0p-23, -0x800000.0p-37,
+ 0xff0000.0p-23, -0x800000.0p-39,
+ 0x800000.0p-22, 0,
+};
+#endif /* USE_UTAB */
+
+#ifdef STRUCT_RETURN
+#define RETURN1(rp, v) do { \
+ (rp)->hi = (v); \
+ (rp)->lo_set = 0; \
+ return; \
+} while (0)
+
+#define RETURN2(rp, h, l) do { \
+ (rp)->hi = (h); \
+ (rp)->lo = (l); \
+ (rp)->lo_set = 1; \
+ return; \
+} while (0)
+
+struct ld {
+ long double hi;
+ long double lo;
+ int lo_set;
+};
+#else
+#define RETURN1(rp, v) RETURNF(v)
+#define RETURN2(rp, h, l) RETURNI((h) + (l))
+#endif
+
+#ifdef STRUCT_RETURN
+static inline __always_inline void
+k_logl(long double x, struct ld *rp)
+#else
+long double
+logl(long double x)
+#endif
+{
+ long double d, val_hi, val_lo;
+ double dd, dk;
+ uint64_t lx, llx;
+ int i, k;
+ uint16_t hx;
+
+ EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+ k = -16383;
+#if 0 /* Hard to do efficiently. Don't do it until we support all modes. */
+ if (x == 1)
+ RETURN1(rp, 0); /* log(1) = +0 in all rounding modes */
+#endif
+ if (hx == 0 || hx >= 0x8000) { /* zero, negative or subnormal? */
+ if (((hx & 0x7fff) | lx | llx) == 0)
+ RETURN1(rp, -1 / zero); /* log(+-0) = -Inf */
+ if (hx != 0)
+ /* log(neg or NaN) = qNaN: */
+ RETURN1(rp, (x - x) / zero);
+ x *= 0x1.0p113; /* subnormal; scale up x */
+ EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+ k = -16383 - 113;
+ } else if (hx >= 0x7fff)
+ RETURN1(rp, x + x); /* log(Inf or NaN) = Inf or qNaN */
+#ifndef STRUCT_RETURN
+ ENTERI();
+#endif
+ k += hx;
+ dk = k;
+
+ /* Scale x to be in [1, 2). */
+ SET_LDBL_EXPSIGN(x, 0x3fff);
+
+ /* 0 <= i <= INTERVALS: */
+#define L2I (49 - LOG2_INTERVALS)
+ i = (lx + (1LL << (L2I - 2))) >> (L2I - 1);
+
+ /*
+ * -0.005280 < d < 0.004838. In particular, the infinite-
+ * precision |d| is <= 2**-7. Rounding of G(i) to 8 bits
+ * ensures that d is representable without extra precision for
+ * this bound on |d| (since when this calculation is expressed
+ * as x*G(i)-1, the multiplication needs as many extra bits as
+ * G(i) has and the subtraction cancels 8 bits). But for
+ * most i (107 cases out of 129), the infinite-precision |d|
+ * is <= 2**-8. G(i) is rounded to 9 bits for such i to give
+ * better accuracy (this works by improving the bound on |d|,
+ * which in turn allows rounding to 9 bits in more cases).
+ * This is only important when the original x is near 1 -- it
+ * lets us avoid using a special method to give the desired
+ * accuracy for such x.
+ */
+ if (0)
+ d = x * G(i) - 1;
+ else {
+#ifdef USE_UTAB
+ d = (x - H(i)) * G(i) + E(i);
+#else
+ long double x_hi;
+ double x_lo;
+
+ /*
+ * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly.
+ * G(i) has at most 9 bits, so the splitting point is not
+ * critical.
+ */
+ INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
+ llx & 0xffffffffff000000ULL);
+ x_lo = x - x_hi;
+ d = x_hi * G(i) - 1 + x_lo * G(i);
+#endif
+ }
+
+ /*
+ * Our algorithm depends on exact cancellation of F_lo(i) and
+ * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is
+ * at the end of the table. This and other technical complications
+ * make it difficult to avoid the double scaling in (dk*ln2) *
+ * log(base) for base != e without losing more accuracy and/or
+ * efficiency than is gained.
+ */
+ /*
+ * Use double precision operations wherever possible, since long
+ * double operations are emulated and are very slow on the only
+ * known machines that support ld128 (sparc64). Also, don't try
+ * to improve parallelism by increasing the number of operations,
+ * since any parallelism on such machines is needed for the
+ * emulation. Horner's method is good for this, and is also good
+ * for accuracy. Horner's method doesn't handle the `lo' term
+ * well, either for efficiency or accuracy. However, for accuracy
+ * we evaluate d * d * P2 separately to take advantage of
+ * by P2 being exact, and this gives a good place to sum the 'lo'
+ * term too.
+ */
+ dd = (double)d;
+ val_lo = d * d * d * (P3 +
+ d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
+ dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
+ dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo) + d * d * P2;
+ val_hi = d;
+#ifdef DEBUG
+ if (fetestexcept(FE_UNDERFLOW))
+ breakpoint();
+#endif
+
+ _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
+ RETURN2(rp, val_hi, val_lo);
+}
+
+long double
+log1pl(long double x)
+{
+ long double d, d_hi, f_lo, val_hi, val_lo;
+ long double f_hi, twopminusk;
+ double d_lo, dd, dk;
+ uint64_t lx, llx;
+ int i, k;
+ int16_t ax, hx;
+
+ DOPRINT_START(&x);
+ EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+ if (hx < 0x3fff) { /* x < 1, or x neg NaN */
+ ax = hx & 0x7fff;
+ if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */
+ if (ax == 0x3fff && (lx | llx) == 0)
+ RETURNP(-1 / zero); /* log1p(-1) = -Inf */
+ /* log1p(x < 1, or x NaN) = qNaN: */
+ RETURNP((x - x) / (x - x));
+ }
+ if (ax <= 0x3f8d) { /* |x| < 2**-113 */
+ if ((int)x == 0)
+ RETURNP(x); /* x with inexact if x != 0 */
+ }
+ f_hi = 1;
+ f_lo = x;
+ } else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */
+ RETURNP(x + x); /* log1p(Inf or NaN) = Inf or qNaN */
+ } else if (hx < 0x40e1) { /* 1 <= x < 2**226 */
+ f_hi = x;
+ f_lo = 1;
+ } else { /* 2**226 <= x < +Inf */
+ f_hi = x;
+ f_lo = 0; /* avoid underflow of the P3 term */
+ }
+ ENTERI();
+ x = f_hi + f_lo;
+ f_lo = (f_hi - x) + f_lo;
+
+ EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+ k = -16383;
+
+ k += hx;
+ dk = k;
+
+ SET_LDBL_EXPSIGN(x, 0x3fff);
+ twopminusk = 1;
+ SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
+ f_lo *= twopminusk;
+
+ i = (lx + (1LL << (L2I - 2))) >> (L2I - 1);
+
+ /*
+ * x*G(i)-1 (with a reduced x) can be represented exactly, as
+ * above, but now we need to evaluate the polynomial on d =
+ * (x+f_lo)*G(i)-1 and extra precision is needed for that.
+ * Since x+x_lo is a hi+lo decomposition and subtracting 1
+ * doesn't lose too many bits, an inexact calculation for
+ * f_lo*G(i) is good enough.
+ */
+ if (0)
+ d_hi = x * G(i) - 1;
+ else {
+#ifdef USE_UTAB
+ d_hi = (x - H(i)) * G(i) + E(i);
+#else
+ long double x_hi;
+ double x_lo;
+
+ INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
+ llx & 0xffffffffff000000ULL);
+ x_lo = x - x_hi;
+ d_hi = x_hi * G(i) - 1 + x_lo * G(i);
+#endif
+ }
+ d_lo = f_lo * G(i);
+
+ /*
+ * This is _2sumF(d_hi, d_lo) inlined. The condition
+ * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not
+ * always satisifed, so it is not clear that this works, but
+ * it works in practice. It works even if it gives a wrong
+ * normalized d_lo, since |d_lo| > |d_hi| implies that i is
+ * nonzero and d is tiny, so the F(i) term dominates d_lo.
+ * In float precision:
+ * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25.
+ * And if d is only a little tinier than that, we would have
+ * another underflow problem for the P3 term; this is also ruled
+ * out by exhaustive testing.)
+ */
+ d = d_hi + d_lo;
+ d_lo = d_hi - d + d_lo;
+ d_hi = d;
+
+ dd = (double)d;
+ val_lo = d * d * d * (P3 +
+ d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
+ dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
+ dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo + d_lo) + d * d * P2;
+ val_hi = d_hi;
+#ifdef DEBUG
+ if (fetestexcept(FE_UNDERFLOW))
+ breakpoint();
+#endif
+
+ _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
+ RETURN2PI(val_hi, val_lo);
+}
+
+#ifdef STRUCT_RETURN
+
+long double
+logl(long double x)
+{
+ struct ld r;
+
+ ENTERI();
+ DOPRINT_START(&x);
+ k_logl(x, &r);
+ RETURNSPI(&r);
+}
+
+/*
+ * 29+113 bit decompositions. The bits are distributed so that the products
+ * of the hi terms are exact in double precision. The types are chosen so
+ * that the products of the hi terms are done in at least double precision,
+ * without any explicit conversions. More natural choices would require a
+ * slow long double precision multiplication.
+ */
+static const double
+invln10_hi = 4.3429448176175356e-1, /* 0x1bcb7b15000000.0p-54 */
+invln2_hi = 1.4426950402557850e0; /* 0x17154765000000.0p-52 */
+static const long double
+invln10_lo = 1.41498268538580090791605082294397000e-10L, /* 0x137287195355baaafad33dc323ee3.0p-145L */
+invln2_lo = 6.33178418956604368501892137426645911e-10L; /* 0x15c17f0bbbe87fed0691d3e88eb57.0p-143L */
+
+long double
+log10l(long double x)
+{
+ struct ld r;
+ long double lo;
+ float hi;
+
+ ENTERI();
+ DOPRINT_START(&x);
+ k_logl(x, &r);
+ if (!r.lo_set)
+ RETURNPI(r.hi);
+ _2sumF(r.hi, r.lo);
+ hi = r.hi;
+ lo = r.lo + (r.hi - hi);
+ RETURN2PI(invln10_hi * hi,
+ (invln10_lo + invln10_hi) * lo + invln10_lo * hi);
+}
+
+long double
+log2l(long double x)
+{
+ struct ld r;
+ long double lo;
+ float hi;
+
+ ENTERI();
+ DOPRINT_START(&x);
+ k_logl(x, &r);
+ if (!r.lo_set)
+ RETURNPI(r.hi);
+ _2sumF(r.hi, r.lo);
+ hi = r.hi;
+ lo = r.lo + (r.hi - hi);
+ RETURN2PI(invln2_hi * hi,
+ (invln2_lo + invln2_hi) * lo + invln2_lo * hi);
+}
+
+#endif /* STRUCT_RETURN */
diff --git a/libm/upstream-freebsd/lib/msun/ld128/s_nanl.c b/libm/upstream-freebsd/lib/msun/ld128/s_nanl.c
new file mode 100644
index 0000000..0f74a13
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/ld128/s_nanl.c
@@ -0,0 +1,46 @@
+/*-
+ * Copyright (c) 2007 David Schultz
+ * 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 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 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 <math.h>
+
+#include "fpmath.h"
+#include "../src/math_private.h"
+
+long double
+nanl(const char *s)
+{
+ union {
+ union IEEEl2bits ieee;
+ uint32_t bits[4];
+ } u;
+
+ _scan_nan(u.bits, 4, s);
+ u.ieee.bits.exp = 0x7fff;
+ u.ieee.bits.manh |= 1ULL << 47; /* make it a quiet NaN */
+ return (u.ieee.e);
+}
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"
diff --git a/tests/math_test.cpp b/tests/math_test.cpp
index fc561bc..0bb072a 100644
--- a/tests/math_test.cpp
+++ b/tests/math_test.cpp
@@ -72,62 +72,97 @@
return u.d;
}
+long double ldouble_subnormal() {
+ union {
+ long double e;
+ unsigned char c[sizeof(long double)];
+ } u;
+
+ // Subnormals must have a zero exponent and non zero significand.
+ // On all supported representation the 17 bit (counting from either sides)
+ // is part of the significand so it should be enough to set that.
+ // It also applies for the case sizeof(double) = sizeof(long double)
+ for (unsigned int i = 0; i < sizeof(long double); i++) {
+ u.c[i] = 0x00;
+ }
+ u.c[sizeof(long double) - 3] = 0x80;
+ u.c[2] = 0x80;
+
+ return u.e;
+}
+
TEST(math, fpclassify) {
ASSERT_EQ(FP_INFINITE, fpclassify(INFINITY));
ASSERT_EQ(FP_INFINITE, fpclassify(HUGE_VALF));
ASSERT_EQ(FP_INFINITE, fpclassify(HUGE_VAL));
+ ASSERT_EQ(FP_INFINITE, fpclassify(HUGE_VALL));
ASSERT_EQ(FP_NAN, fpclassify(nanf("")));
ASSERT_EQ(FP_NAN, fpclassify(nan("")));
+ ASSERT_EQ(FP_NAN, fpclassify(nanl("")));
ASSERT_EQ(FP_NORMAL, fpclassify(1.0f));
ASSERT_EQ(FP_NORMAL, fpclassify(1.0));
+ ASSERT_EQ(FP_NORMAL, fpclassify(1.0l));
ASSERT_EQ(FP_SUBNORMAL, fpclassify(float_subnormal()));
ASSERT_EQ(FP_SUBNORMAL, fpclassify(double_subnormal()));
+ ASSERT_EQ(FP_SUBNORMAL, fpclassify(ldouble_subnormal()));
ASSERT_EQ(FP_ZERO, fpclassify(0.0f));
ASSERT_EQ(FP_ZERO, fpclassify(0.0));
+ ASSERT_EQ(FP_ZERO, fpclassify(0.0l));
}
TEST(math, isfinite) {
ASSERT_TRUE(test_capture_isfinite(123.0f));
ASSERT_TRUE(test_capture_isfinite(123.0));
+ ASSERT_TRUE(test_capture_isfinite(123.0l));
ASSERT_FALSE(test_capture_isfinite(HUGE_VALF));
ASSERT_FALSE(test_capture_isfinite(HUGE_VAL));
+ ASSERT_FALSE(test_capture_isfinite(HUGE_VALL));
}
TEST(math, isinf) {
ASSERT_FALSE(test_capture_isinf(123.0f));
ASSERT_FALSE(test_capture_isinf(123.0));
+ ASSERT_FALSE(test_capture_isinf(123.0l));
ASSERT_TRUE(test_capture_isinf(HUGE_VALF));
ASSERT_TRUE(test_capture_isinf(HUGE_VAL));
+ ASSERT_TRUE(test_capture_isinf(HUGE_VALL));
}
TEST(math, isnan) {
ASSERT_FALSE(test_capture_isnan(123.0f));
ASSERT_FALSE(test_capture_isnan(123.0));
+ ASSERT_FALSE(test_capture_isnan(123.0l));
ASSERT_TRUE(test_capture_isnan(nanf("")));
ASSERT_TRUE(test_capture_isnan(nan("")));
+ ASSERT_TRUE(test_capture_isnan(nanl("")));
}
TEST(math, isnormal) {
ASSERT_TRUE(isnormal(123.0f));
ASSERT_TRUE(isnormal(123.0));
+ ASSERT_TRUE(isnormal(123.0l));
ASSERT_FALSE(isnormal(float_subnormal()));
ASSERT_FALSE(isnormal(double_subnormal()));
+ ASSERT_FALSE(isnormal(ldouble_subnormal()));
}
// TODO: isgreater, isgreaterequals, isless, islessequal, islessgreater, isunordered
TEST(math, signbit) {
ASSERT_EQ(0, test_capture_signbit(0.0f));
ASSERT_EQ(0, test_capture_signbit(0.0));
+ ASSERT_EQ(0, test_capture_signbit(0.0l));
ASSERT_EQ(0, test_capture_signbit(1.0f));
ASSERT_EQ(0, test_capture_signbit(1.0));
+ ASSERT_EQ(0, test_capture_signbit(1.0l));
ASSERT_NE(0, test_capture_signbit(-1.0f));
ASSERT_NE(0, test_capture_signbit(-1.0));
+ ASSERT_NE(0, test_capture_signbit(-1.0l));
}
TEST(math, __fpclassifyd) {
@@ -158,9 +193,9 @@
#if defined(__BIONIC__)
EXPECT_EQ(FP_INFINITE, __fpclassifyl(HUGE_VALL));
EXPECT_EQ(FP_NAN, __fpclassifyl(nanl("")));
- EXPECT_EQ(FP_NORMAL, __fpclassifyl(1.0));
- EXPECT_EQ(FP_SUBNORMAL, __fpclassifyl(double_subnormal()));
- EXPECT_EQ(FP_ZERO, __fpclassifyl(0.0));
+ EXPECT_EQ(FP_NORMAL, __fpclassifyl(1.0l));
+ EXPECT_EQ(FP_SUBNORMAL, __fpclassifyl(ldouble_subnormal()));
+ EXPECT_EQ(FP_ZERO, __fpclassifyl(0.0l));
#else // __BIONIC__
GTEST_LOG_(INFO) << "This test does nothing.\n";
#endif // __BIONIC__
@@ -191,7 +226,7 @@
TEST(math, __isfinitel) {
#if defined(__BIONIC__)
- ASSERT_TRUE(__isfinitel(123.0f));
+ ASSERT_TRUE(__isfinitel(123.0l));
ASSERT_FALSE(__isfinitel(HUGE_VALL));
#else // __BIONIC__
GTEST_LOG_(INFO) << "This test does nothing.\n";
@@ -209,7 +244,7 @@
}
TEST(math, __isinfl) {
- ASSERT_FALSE(__isinfl(123.0));
+ ASSERT_FALSE(__isinfl(123.0l));
ASSERT_TRUE(__isinfl(HUGE_VALL));
}
@@ -219,7 +254,7 @@
}
TEST(math, __isnanl) {
- ASSERT_FALSE(__isnanl(123.0));
+ ASSERT_FALSE(__isnanl(123.0l));
ASSERT_TRUE(__isnanl(nanl("")));
}
@@ -248,8 +283,8 @@
TEST(math, __isnormall) {
#if defined(__BIONIC__)
- ASSERT_TRUE(__isnormall(123.0));
- ASSERT_FALSE(__isnormall(double_subnormal()));
+ ASSERT_TRUE(__isnormall(123.0l));
+ ASSERT_FALSE(__isnormall(ldouble_subnormal()));
#else // __BIONIC__
GTEST_LOG_(INFO) << "This test does nothing.\n";
#endif // __BIONIC__
@@ -268,13 +303,13 @@
}
TEST(math, __signbitl) {
- ASSERT_EQ(0, __signbitl(0.0));
- ASSERT_EQ(0, __signbitl(1.0));
- ASSERT_NE(0, __signbitl(-1.0));
+ ASSERT_EQ(0l, __signbitl(0.0l));
+ ASSERT_EQ(0l, __signbitl(1.0l));
+ ASSERT_NE(0l, __signbitl(-1.0l));
}
TEST(math, acos) {
- ASSERT_FLOAT_EQ(M_PI/2.0, acos(0.0));
+ ASSERT_DOUBLE_EQ(M_PI/2.0, acos(0.0));
}
TEST(math, acosf) {
@@ -282,11 +317,11 @@
}
TEST(math, acosl) {
- ASSERT_FLOAT_EQ(M_PI/2.0, acosl(0.0));
+ ASSERT_DOUBLE_EQ(M_PI/2.0l, acosl(0.0l));
}
TEST(math, asin) {
- ASSERT_FLOAT_EQ(0.0, asin(0.0));
+ ASSERT_DOUBLE_EQ(0.0, asin(0.0));
}
TEST(math, asinf) {
@@ -294,11 +329,11 @@
}
TEST(math, asinl) {
- ASSERT_FLOAT_EQ(0.0, asinl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, asinl(0.0l));
}
TEST(math, atan) {
- ASSERT_FLOAT_EQ(0.0, atan(0.0));
+ ASSERT_DOUBLE_EQ(0.0, atan(0.0));
}
TEST(math, atanf) {
@@ -306,11 +341,11 @@
}
TEST(math, atanl) {
- ASSERT_FLOAT_EQ(0.0, atanl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, atanl(0.0l));
}
TEST(math, atan2) {
- ASSERT_FLOAT_EQ(0.0, atan2(0.0, 0.0));
+ ASSERT_DOUBLE_EQ(0.0, atan2(0.0, 0.0));
}
TEST(math, atan2f) {
@@ -318,11 +353,11 @@
}
TEST(math, atan2l) {
- ASSERT_FLOAT_EQ(0.0, atan2l(0.0, 0.0));
+ ASSERT_DOUBLE_EQ(0.0l, atan2l(0.0l, 0.0l));
}
TEST(math, cos) {
- ASSERT_FLOAT_EQ(1.0, cos(0.0));
+ ASSERT_DOUBLE_EQ(1.0, cos(0.0));
}
TEST(math, cosf) {
@@ -330,7 +365,7 @@
}
TEST(math, cosl) {
- ASSERT_FLOAT_EQ(1.0, cosl(0.0));
+ ASSERT_DOUBLE_EQ(1.0l, cosl(0.0l));
}
TEST(math, sin) {
@@ -342,11 +377,11 @@
}
TEST(math, sinl) {
- ASSERT_FLOAT_EQ(0.0, sinl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, sinl(0.0l));
}
TEST(math, tan) {
- ASSERT_FLOAT_EQ(0.0, tan(0.0));
+ ASSERT_DOUBLE_EQ(0.0, tan(0.0));
}
TEST(math, tanf) {
@@ -354,11 +389,11 @@
}
TEST(math, tanl) {
- ASSERT_FLOAT_EQ(0.0, tanl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, tanl(0.0l));
}
TEST(math, acosh) {
- ASSERT_FLOAT_EQ(0.0, acosh(1.0));
+ ASSERT_DOUBLE_EQ(0.0, acosh(1.0));
}
TEST(math, acoshf) {
@@ -366,11 +401,11 @@
}
TEST(math, acoshl) {
- ASSERT_FLOAT_EQ(0.0, acoshl(1.0));
+ ASSERT_DOUBLE_EQ(0.0l, acoshl(1.0l));
}
TEST(math, asinh) {
- ASSERT_FLOAT_EQ(0.0, asinh(0.0));
+ ASSERT_DOUBLE_EQ(0.0, asinh(0.0));
}
TEST(math, asinhf) {
@@ -378,11 +413,11 @@
}
TEST(math, asinhl) {
- ASSERT_FLOAT_EQ(0.0, asinhl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, asinhl(0.0l));
}
TEST(math, atanh) {
- ASSERT_FLOAT_EQ(0.0, atanh(0.0));
+ ASSERT_DOUBLE_EQ(0.0, atanh(0.0));
}
TEST(math, atanhf) {
@@ -390,11 +425,11 @@
}
TEST(math, atanhl) {
- ASSERT_FLOAT_EQ(0.0, atanhl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, atanhl(0.0l));
}
TEST(math, cosh) {
- ASSERT_FLOAT_EQ(1.0, cosh(0.0));
+ ASSERT_DOUBLE_EQ(1.0, cosh(0.0));
}
TEST(math, coshf) {
@@ -402,11 +437,11 @@
}
TEST(math, coshl) {
- ASSERT_FLOAT_EQ(1.0, coshl(0.0));
+ ASSERT_DOUBLE_EQ(1.0l, coshl(0.0l));
}
TEST(math, sinh) {
- ASSERT_FLOAT_EQ(0.0, sinh(0.0));
+ ASSERT_DOUBLE_EQ(0.0, sinh(0.0));
}
TEST(math, sinhf) {
@@ -414,11 +449,11 @@
}
TEST(math, sinhl) {
- ASSERT_FLOAT_EQ(0.0, sinhl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, sinhl(0.0l));
}
TEST(math, tanh) {
- ASSERT_FLOAT_EQ(0.0, tanh(0.0));
+ ASSERT_DOUBLE_EQ(0.0, tanh(0.0));
}
TEST(math, tanhf) {
@@ -426,11 +461,11 @@
}
TEST(math, tanhl) {
- ASSERT_FLOAT_EQ(0.0, tanhl(0.0));
+ ASSERT_DOUBLE_EQ(0.0l, tanhl(0.0l));
}
TEST(math, log) {
- ASSERT_FLOAT_EQ(1.0, log(M_E));
+ ASSERT_DOUBLE_EQ(1.0, log(M_E));
}
TEST(math, logf) {
@@ -438,11 +473,11 @@
}
TEST(math, logl) {
- ASSERT_FLOAT_EQ(1.0, logl(M_E));
+ ASSERT_DOUBLE_EQ(1.0l, logl(M_E));
}
TEST(math, log2) {
- ASSERT_FLOAT_EQ(12.0, log2(4096.0));
+ ASSERT_DOUBLE_EQ(12.0, log2(4096.0));
}
TEST(math, log2f) {
@@ -450,11 +485,11 @@
}
TEST(math, log2l) {
- ASSERT_FLOAT_EQ(12.0, log2l(4096.0));
+ ASSERT_DOUBLE_EQ(12.0l, log2l(4096.0l));
}
TEST(math, log10) {
- ASSERT_FLOAT_EQ(3.0, log10(1000.0));
+ ASSERT_DOUBLE_EQ(3.0, log10(1000.0));
}
TEST(math, log10f) {
@@ -462,11 +497,11 @@
}
TEST(math, log10l) {
- ASSERT_FLOAT_EQ(3.0, log10l(1000.0));
+ ASSERT_DOUBLE_EQ(3.0l, log10l(1000.0l));
}
TEST(math, cbrt) {
- ASSERT_FLOAT_EQ(3.0, cbrt(27.0));
+ ASSERT_DOUBLE_EQ(3.0, cbrt(27.0));
}
TEST(math, cbrtf) {
@@ -474,11 +509,11 @@
}
TEST(math, cbrtl) {
- ASSERT_FLOAT_EQ(3.0, cbrtl(27.0));
+ ASSERT_DOUBLE_EQ(3.0l, cbrtl(27.0l));
}
TEST(math, sqrt) {
- ASSERT_FLOAT_EQ(2.0, sqrt(4.0));
+ ASSERT_DOUBLE_EQ(2.0, sqrt(4.0));
}
TEST(math, sqrtf) {
@@ -486,12 +521,12 @@
}
TEST(math, sqrtl) {
- ASSERT_FLOAT_EQ(2.0, sqrtl(4.0));
+ ASSERT_DOUBLE_EQ(2.0l, sqrtl(4.0l));
}
TEST(math, exp) {
- ASSERT_FLOAT_EQ(1.0, exp(0.0));
- ASSERT_FLOAT_EQ(M_E, exp(1.0));
+ ASSERT_DOUBLE_EQ(1.0, exp(0.0));
+ ASSERT_DOUBLE_EQ(M_E, exp(1.0));
}
TEST(math, expf) {
@@ -500,12 +535,12 @@
}
TEST(math, expl) {
- ASSERT_FLOAT_EQ(1.0, expl(0.0));
- ASSERT_FLOAT_EQ(M_E, expl(1.0));
+ ASSERT_DOUBLE_EQ(1.0l, expl(0.0l));
+ ASSERT_DOUBLE_EQ(M_E, expl(1.0l));
}
TEST(math, exp2) {
- ASSERT_FLOAT_EQ(8.0, exp2(3.0));
+ ASSERT_DOUBLE_EQ(8.0, exp2(3.0));
}
TEST(math, exp2f) {
@@ -513,11 +548,11 @@
}
TEST(math, exp2l) {
- ASSERT_FLOAT_EQ(8.0, exp2l(3.0));
+ ASSERT_DOUBLE_EQ(8.0l, exp2l(3.0l));
}
TEST(math, expm1) {
- ASSERT_FLOAT_EQ(M_E - 1.0, expm1(1.0));
+ ASSERT_DOUBLE_EQ(M_E - 1.0, expm1(1.0));
}
TEST(math, expm1f) {
@@ -525,14 +560,14 @@
}
TEST(math, expm1l) {
- ASSERT_FLOAT_EQ(M_E - 1.0, expm1l(1.0));
+ ASSERT_DOUBLE_EQ(M_E - 1.0l, expm1l(1.0l));
}
TEST(math, pow) {
ASSERT_TRUE(isnan(pow(nan(""), 3.0)));
- ASSERT_FLOAT_EQ(1.0, (pow(1.0, nan(""))));
+ ASSERT_DOUBLE_EQ(1.0, (pow(1.0, nan(""))));
ASSERT_TRUE(isnan(pow(2.0, nan(""))));
- ASSERT_FLOAT_EQ(8.0, pow(2.0, 3.0));
+ ASSERT_DOUBLE_EQ(8.0, pow(2.0, 3.0));
}
TEST(math, powf) {
@@ -543,14 +578,14 @@
}
TEST(math, powl) {
- ASSERT_TRUE(__isnanl(powl(nanl(""), 3.0)));
- ASSERT_FLOAT_EQ(1.0, (powl(1.0, nanl(""))));
- ASSERT_TRUE(__isnanl(powl(2.0, nanl(""))));
- ASSERT_FLOAT_EQ(8.0, powl(2.0, 3.0));
+ ASSERT_TRUE(__isnanl(powl(nanl(""), 3.0l)));
+ ASSERT_DOUBLE_EQ(1.0l, (powl(1.0l, nanl(""))));
+ ASSERT_TRUE(__isnanl(powl(2.0l, nanl(""))));
+ ASSERT_DOUBLE_EQ(8.0l, powl(2.0l, 3.0l));
}
TEST(math, ceil) {
- ASSERT_FLOAT_EQ(1.0, ceil(0.9));
+ ASSERT_DOUBLE_EQ(1.0, ceil(0.9));
}
TEST(math, ceilf) {
@@ -558,11 +593,11 @@
}
TEST(math, ceill) {
- ASSERT_FLOAT_EQ(1.0, ceill(0.9));
+ ASSERT_DOUBLE_EQ(1.0, ceill(0.9l));
}
TEST(math, floor) {
- ASSERT_FLOAT_EQ(1.0, floor(1.1));
+ ASSERT_DOUBLE_EQ(1.0, floor(1.1));
}
TEST(math, floorf) {
@@ -570,11 +605,11 @@
}
TEST(math, floorl) {
- ASSERT_FLOAT_EQ(1.0, floorl(1.1));
+ ASSERT_DOUBLE_EQ(1.0l, floorl(1.1l));
}
TEST(math, fabs) {
- ASSERT_FLOAT_EQ(1.0, fabs(-1.0));
+ ASSERT_DOUBLE_EQ(1.0, fabs(-1.0));
}
TEST(math, fabsf) {
@@ -582,11 +617,11 @@
}
TEST(math, fabsl) {
- ASSERT_FLOAT_EQ(1.0, fabsl(-1.0));
+ ASSERT_DOUBLE_EQ(1.0l, fabsl(-1.0l));
}
TEST(math, ldexp) {
- ASSERT_FLOAT_EQ(16.0, ldexp(2.0, 3.0));
+ ASSERT_DOUBLE_EQ(16.0, ldexp(2.0, 3.0));
}
TEST(math, ldexpf) {
@@ -594,11 +629,11 @@
}
TEST(math, ldexpl) {
- ASSERT_FLOAT_EQ(16.0, ldexpl(2.0, 3.0));
+ ASSERT_DOUBLE_EQ(16.0l, ldexpl(2.0l, 3.0));
}
TEST(math, fmod) {
- ASSERT_FLOAT_EQ(2.0, fmod(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(2.0, fmod(12.0, 10.0));
}
TEST(math, fmodf) {
@@ -606,11 +641,11 @@
}
TEST(math, fmodl) {
- ASSERT_FLOAT_EQ(2.0, fmodl(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(2.0l, fmodl(12.0l, 10.0l));
}
TEST(math, remainder) {
- ASSERT_FLOAT_EQ(2.0, remainder(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(2.0, remainder(12.0, 10.0));
}
TEST(math, remainderf) {
@@ -618,11 +653,11 @@
}
TEST(math, remainderl) {
- ASSERT_FLOAT_EQ(2.0, remainderl(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(2.0l, remainderl(12.0l, 10.0l));
}
TEST(math, drem) {
- ASSERT_FLOAT_EQ(2.0, drem(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(2.0, drem(12.0, 10.0));
}
TEST(math, dremf) {
@@ -630,9 +665,9 @@
}
TEST(math, fmax) {
- ASSERT_FLOAT_EQ(12.0, fmax(12.0, 10.0));
- ASSERT_FLOAT_EQ(12.0, fmax(12.0, nan("")));
- ASSERT_FLOAT_EQ(12.0, fmax(nan(""), 12.0));
+ ASSERT_DOUBLE_EQ(12.0, fmax(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(12.0, fmax(12.0, nan("")));
+ ASSERT_DOUBLE_EQ(12.0, fmax(nan(""), 12.0));
}
TEST(math, fmaxf) {
@@ -642,15 +677,15 @@
}
TEST(math, fmaxl) {
- ASSERT_FLOAT_EQ(12.0, fmaxl(12.0, 10.0));
- ASSERT_FLOAT_EQ(12.0, fmaxl(12.0, nanl("")));
- ASSERT_FLOAT_EQ(12.0, fmaxl(nanl(""), 12.0));
+ ASSERT_DOUBLE_EQ(12.0l, fmaxl(12.0l, 10.0l));
+ ASSERT_DOUBLE_EQ(12.0l, fmaxl(12.0l, nanl("")));
+ ASSERT_DOUBLE_EQ(12.0l, fmaxl(nanl(""), 12.0l));
}
TEST(math, fmin) {
- ASSERT_FLOAT_EQ(10.0, fmin(12.0, 10.0));
- ASSERT_FLOAT_EQ(12.0, fmin(12.0, nan("")));
- ASSERT_FLOAT_EQ(12.0, fmin(nan(""), 12.0));
+ ASSERT_DOUBLE_EQ(10.0, fmin(12.0, 10.0));
+ ASSERT_DOUBLE_EQ(12.0, fmin(12.0, nan("")));
+ ASSERT_DOUBLE_EQ(12.0, fmin(nan(""), 12.0));
}
TEST(math, fminf) {
@@ -660,13 +695,13 @@
}
TEST(math, fminl) {
- ASSERT_FLOAT_EQ(10.0, fminl(12.0, 10.0));
- ASSERT_FLOAT_EQ(12.0, fminl(12.0, nan("")));
- ASSERT_FLOAT_EQ(12.0, fminl(nan(""), 12.0));
+ ASSERT_DOUBLE_EQ(10.0l, fminl(12.0l, 10.0l));
+ ASSERT_DOUBLE_EQ(12.0l, fminl(12.0l, nanl("")));
+ ASSERT_DOUBLE_EQ(12.0l, fminl(nanl(""), 12.0l));
}
TEST(math, fma) {
- ASSERT_FLOAT_EQ(10.0, fma(2.0, 3.0, 4.0));
+ ASSERT_DOUBLE_EQ(10.0, fma(2.0, 3.0, 4.0));
}
TEST(math, fmaf) {
@@ -674,11 +709,11 @@
}
TEST(math, fmal) {
- ASSERT_FLOAT_EQ(10.0, fmal(2.0, 3.0, 4.0));
+ ASSERT_DOUBLE_EQ(10.0l, fmal(2.0l, 3.0l, 4.0l));
}
TEST(math, hypot) {
- ASSERT_FLOAT_EQ(5.0, hypot(3.0, 4.0));
+ ASSERT_DOUBLE_EQ(5.0, hypot(3.0, 4.0));
}
TEST(math, hypotf) {
@@ -686,11 +721,11 @@
}
TEST(math, hypotl) {
- ASSERT_FLOAT_EQ(5.0, hypotl(3.0, 4.0));
+ ASSERT_DOUBLE_EQ(5.0l, hypotl(3.0l, 4.0l));
}
TEST(math, erf) {
- ASSERT_FLOAT_EQ(0.84270078, erf(1.0));
+ ASSERT_DOUBLE_EQ(0.84270079294971489, erf(1.0));
}
TEST(math, erff) {
@@ -698,11 +733,11 @@
}
TEST(math, erfl) {
- ASSERT_FLOAT_EQ(0.84270078, erfl(1.0));
+ ASSERT_DOUBLE_EQ(0.84270079294971489l, erfl(1.0l));
}
TEST(math, erfc) {
- ASSERT_FLOAT_EQ(0.15729921, erfc(1.0));
+ ASSERT_DOUBLE_EQ(0.15729920705028513, erfc(1.0));
}
TEST(math, erfcf) {
@@ -710,27 +745,27 @@
}
TEST(math, erfcl) {
- ASSERT_FLOAT_EQ(0.15729921, erfcl(1.0));
+ ASSERT_DOUBLE_EQ(0.15729920705028513l, erfcl(1.0l));
}
TEST(math, lrint) {
fesetround(FE_UPWARD); // lrint/lrintf/lrintl obey the rounding mode.
ASSERT_EQ(1235, lrint(1234.01));
ASSERT_EQ(1235, lrintf(1234.01f));
- ASSERT_EQ(1235, lrintl(1234.01));
+ ASSERT_EQ(1235, lrintl(1234.01l));
fesetround(FE_TOWARDZERO); // lrint/lrintf/lrintl obey the rounding mode.
ASSERT_EQ(1234, lrint(1234.01));
ASSERT_EQ(1234, lrintf(1234.01f));
- ASSERT_EQ(1234, lrintl(1234.01));
+ ASSERT_EQ(1234, lrintl(1234.01l));
fesetround(FE_UPWARD); // llrint/llrintf/llrintl obey the rounding mode.
ASSERT_EQ(1235L, llrint(1234.01));
ASSERT_EQ(1235L, llrintf(1234.01f));
- ASSERT_EQ(1235L, llrintl(1234.01));
+ ASSERT_EQ(1235L, llrintl(1234.01l));
fesetround(FE_TOWARDZERO); // llrint/llrintf/llrintl obey the rounding mode.
ASSERT_EQ(1234L, llrint(1234.01));
ASSERT_EQ(1234L, llrintf(1234.01f));
- ASSERT_EQ(1234L, llrintl(1234.01));
+ ASSERT_EQ(1234L, llrintl(1234.01l));
}
TEST(math, rint) {
@@ -748,15 +783,15 @@
ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) != 0);
feclearexcept(FE_ALL_EXCEPT); // rint/rintf/rintl do set the FE_INEXACT flag.
- ASSERT_EQ(1234.0, rintl(1234.0));
+ ASSERT_EQ(1234.0, rintl(1234.0l));
ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
- ASSERT_EQ(1235.0, rintl(1234.01));
+ ASSERT_EQ(1235.0, rintl(1234.01l));
ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) != 0);
fesetround(FE_TOWARDZERO); // rint/rintf obey the rounding mode.
ASSERT_EQ(1234.0, rint(1234.01));
ASSERT_EQ(1234.0f, rintf(1234.01f));
- ASSERT_EQ(1234.0, rintl(1234.01));
+ ASSERT_EQ(1234.0, rintl(1234.01l));
}
TEST(math, nearbyint) {
@@ -774,29 +809,29 @@
ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
feclearexcept(FE_ALL_EXCEPT); // nearbyint/nearbyintf/nearbyintl don't set the FE_INEXACT flag.
- ASSERT_EQ(1234.0, nearbyintl(1234.0));
+ ASSERT_EQ(1234.0, nearbyintl(1234.0l));
ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
- ASSERT_EQ(1235.0, nearbyintl(1234.01));
+ ASSERT_EQ(1235.0, nearbyintl(1234.01l));
ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
fesetround(FE_TOWARDZERO); // nearbyint/nearbyintf/nearbyintl obey the rounding mode.
ASSERT_EQ(1234.0, nearbyint(1234.01));
ASSERT_EQ(1234.0f, nearbyintf(1234.01f));
- ASSERT_EQ(1234.0, nearbyintl(1234.01));
+ ASSERT_EQ(1234.0, nearbyintl(1234.01l));
}
TEST(math, lround) {
fesetround(FE_UPWARD); // lround ignores the rounding mode.
ASSERT_EQ(1234, lround(1234.01));
ASSERT_EQ(1234, lroundf(1234.01f));
- ASSERT_EQ(1234, lroundl(1234.01));
+ ASSERT_EQ(1234, lroundl(1234.01l));
}
TEST(math, llround) {
fesetround(FE_UPWARD); // llround ignores the rounding mode.
ASSERT_EQ(1234L, llround(1234.01));
ASSERT_EQ(1234L, llroundf(1234.01f));
- ASSERT_EQ(1234L, llroundl(1234.01));
+ ASSERT_EQ(1234L, llroundl(1234.01l));
}
TEST(math, ilogb) {
@@ -816,11 +851,11 @@
}
TEST(math, ilogbl) {
- ASSERT_EQ(FP_ILOGB0, ilogbl(0.0));
+ ASSERT_EQ(FP_ILOGB0, ilogbl(0.0l));
ASSERT_EQ(FP_ILOGBNAN, ilogbl(nanl("")));
ASSERT_EQ(INT_MAX, ilogbl(HUGE_VALL));
- ASSERT_EQ(0, ilogbl(1.0));
- ASSERT_EQ(3, ilogbl(10.0));
+ ASSERT_EQ(0l, ilogbl(1.0l));
+ ASSERT_EQ(3l, ilogbl(10.0l));
}
TEST(math, logb) {
@@ -840,18 +875,18 @@
}
TEST(math, logbl) {
- ASSERT_EQ(-HUGE_VAL, logbl(0.0));
+ ASSERT_EQ(-HUGE_VAL, logbl(0.0l));
ASSERT_TRUE(isnan(logbl(nanl(""))));
ASSERT_TRUE(isinf(logbl(HUGE_VALL)));
- ASSERT_EQ(0.0, logbl(1.0));
- ASSERT_EQ(3.0, logbl(10.0));
+ ASSERT_EQ(0.0l, logbl(1.0l));
+ ASSERT_EQ(3.0l, logbl(10.0l));
}
TEST(math, log1p) {
ASSERT_EQ(-HUGE_VAL, log1p(-1.0));
ASSERT_TRUE(isnan(log1p(nan(""))));
ASSERT_TRUE(isinf(log1p(HUGE_VAL)));
- ASSERT_FLOAT_EQ(1.0, log1p(M_E - 1.0));
+ ASSERT_DOUBLE_EQ(1.0, log1p(M_E - 1.0));
}
TEST(math, log1pf) {
@@ -862,16 +897,16 @@
}
TEST(math, log1pl) {
- ASSERT_EQ(-HUGE_VALL, log1pl(-1.0));
+ ASSERT_EQ(-HUGE_VALL, log1pl(-1.0l));
ASSERT_TRUE(isnan(log1pl(nanl(""))));
ASSERT_TRUE(isinf(log1pl(HUGE_VALL)));
- ASSERT_FLOAT_EQ(1.0, log1pl(M_E - 1.0));
+ ASSERT_DOUBLE_EQ(1.0l, log1pl(M_E - 1.0l));
}
TEST(math, fdim) {
- ASSERT_FLOAT_EQ(0.0, fdim(1.0, 1.0));
- ASSERT_FLOAT_EQ(1.0, fdim(2.0, 1.0));
- ASSERT_FLOAT_EQ(0.0, fdim(1.0, 2.0));
+ ASSERT_DOUBLE_EQ(0.0, fdim(1.0, 1.0));
+ ASSERT_DOUBLE_EQ(1.0, fdim(2.0, 1.0));
+ ASSERT_DOUBLE_EQ(0.0, fdim(1.0, 2.0));
}
TEST(math, fdimf) {
@@ -881,19 +916,19 @@
}
TEST(math, fdiml) {
- ASSERT_FLOAT_EQ(0.0, fdiml(1.0, 1.0));
- ASSERT_FLOAT_EQ(1.0, fdiml(2.0, 1.0));
- ASSERT_FLOAT_EQ(0.0, fdiml(1.0, 2.0));
+ ASSERT_DOUBLE_EQ(0.0l, fdiml(1.0l, 1.0l));
+ ASSERT_DOUBLE_EQ(1.0l, fdiml(2.0l, 1.0l));
+ ASSERT_DOUBLE_EQ(0.0l, fdiml(1.0l, 2.0l));
}
TEST(math, round) {
fesetround(FE_TOWARDZERO); // round ignores the rounding mode and always rounds away from zero.
- ASSERT_FLOAT_EQ(1.0, round(0.5));
- ASSERT_FLOAT_EQ(-1.0, round(-0.5));
- ASSERT_FLOAT_EQ(0.0, round(0.0));
- ASSERT_FLOAT_EQ(-0.0, round(-0.0));
+ ASSERT_DOUBLE_EQ(1.0, round(0.5));
+ ASSERT_DOUBLE_EQ(-1.0, round(-0.5));
+ ASSERT_DOUBLE_EQ(0.0, round(0.0));
+ ASSERT_DOUBLE_EQ(-0.0, round(-0.0));
ASSERT_TRUE(isnan(round(nan(""))));
- ASSERT_FLOAT_EQ(HUGE_VAL, round(HUGE_VAL));
+ ASSERT_DOUBLE_EQ(HUGE_VAL, round(HUGE_VAL));
}
TEST(math, roundf) {
@@ -908,22 +943,22 @@
TEST(math, roundl) {
fesetround(FE_TOWARDZERO); // roundl ignores the rounding mode and always rounds away from zero.
- ASSERT_FLOAT_EQ(1.0, roundl(0.5));
- ASSERT_FLOAT_EQ(-1.0, roundl(-0.5));
- ASSERT_FLOAT_EQ(0.0, roundl(0.0));
- ASSERT_FLOAT_EQ(-0.0, roundl(-0.0));
+ ASSERT_DOUBLE_EQ(1.0l, roundl(0.5l));
+ ASSERT_DOUBLE_EQ(-1.0l, roundl(-0.5l));
+ ASSERT_DOUBLE_EQ(0.0l, roundl(0.0l));
+ ASSERT_DOUBLE_EQ(-0.0l, roundl(-0.0l));
ASSERT_TRUE(isnan(roundl(nanl(""))));
- ASSERT_FLOAT_EQ(HUGE_VALL, roundl(HUGE_VALL));
+ ASSERT_DOUBLE_EQ(HUGE_VALL, roundl(HUGE_VALL));
}
TEST(math, trunc) {
fesetround(FE_UPWARD); // trunc ignores the rounding mode and always rounds toward zero.
- ASSERT_FLOAT_EQ(1.0, trunc(1.5));
- ASSERT_FLOAT_EQ(-1.0, trunc(-1.5));
- ASSERT_FLOAT_EQ(0.0, trunc(0.0));
- ASSERT_FLOAT_EQ(-0.0, trunc(-0.0));
+ ASSERT_DOUBLE_EQ(1.0, trunc(1.5));
+ ASSERT_DOUBLE_EQ(-1.0, trunc(-1.5));
+ ASSERT_DOUBLE_EQ(0.0, trunc(0.0));
+ ASSERT_DOUBLE_EQ(-0.0, trunc(-0.0));
ASSERT_TRUE(isnan(trunc(nan(""))));
- ASSERT_FLOAT_EQ(HUGE_VAL, trunc(HUGE_VAL));
+ ASSERT_DOUBLE_EQ(HUGE_VAL, trunc(HUGE_VAL));
}
TEST(math, truncf) {
@@ -938,18 +973,18 @@
TEST(math, truncl) {
fesetround(FE_UPWARD); // truncl ignores the rounding mode and always rounds toward zero.
- ASSERT_FLOAT_EQ(1.0, truncl(1.5));
- ASSERT_FLOAT_EQ(-1.0, truncl(-1.5));
- ASSERT_FLOAT_EQ(0.0, truncl(0.0));
- ASSERT_FLOAT_EQ(-0.0, truncl(-0.0));
+ ASSERT_DOUBLE_EQ(1.0l, truncl(1.5l));
+ ASSERT_DOUBLE_EQ(-1.0l, truncl(-1.5l));
+ ASSERT_DOUBLE_EQ(0.0l, truncl(0.0l));
+ ASSERT_DOUBLE_EQ(-0.0l, truncl(-0.0l));
ASSERT_TRUE(isnan(truncl(nan(""))));
- ASSERT_FLOAT_EQ(HUGE_VALL, truncl(HUGE_VALL));
+ ASSERT_DOUBLE_EQ(HUGE_VALL, truncl(HUGE_VALL));
}
TEST(math, nextafter) {
- ASSERT_FLOAT_EQ(0.0, nextafter(0.0, 0.0));
- ASSERT_FLOAT_EQ(1.4012985e-45, nextafter(0.0, 1.0));
- ASSERT_FLOAT_EQ(0.0, nextafter(0.0, -1.0));
+ ASSERT_DOUBLE_EQ(0.0, nextafter(0.0, 0.0));
+ ASSERT_DOUBLE_EQ(4.9406564584124654e-324, nextafter(0.0, 1.0));
+ ASSERT_DOUBLE_EQ(0.0, nextafter(0.0, -1.0));
}
TEST(math, nextafterf) {
@@ -959,9 +994,12 @@
}
TEST(math, nextafterl) {
- ASSERT_FLOAT_EQ(0.0, nextafterl(0.0, 0.0));
- ASSERT_FLOAT_EQ(1.4012985e-45, nextafterl(0.0, 1.0));
- ASSERT_FLOAT_EQ(0.0, nextafterl(0.0, -1.0));
+ ASSERT_DOUBLE_EQ(0.0l, nextafterl(0.0l, 0.0l));
+ // Use a runtime value to accomodate the case when
+ // sizeof(double) == sizeof(long double)
+ long double smallest_positive = ldexpl(1.0l, LDBL_MIN_EXP - LDBL_MANT_DIG);
+ ASSERT_DOUBLE_EQ(smallest_positive, nextafterl(0.0l, 1.0l));
+ ASSERT_DOUBLE_EQ(0.0l, nextafterl(0.0l, -1.0l));
}
// TODO: nexttoward
@@ -969,10 +1007,10 @@
// TODO: nexttowardl
TEST(math, copysign) {
- ASSERT_FLOAT_EQ(0.0, copysign(0.0, 1.0));
- ASSERT_FLOAT_EQ(-0.0, copysign(0.0, -1.0));
- ASSERT_FLOAT_EQ(2.0, copysign(2.0, 1.0));
- ASSERT_FLOAT_EQ(-2.0, copysign(2.0, -1.0));
+ ASSERT_DOUBLE_EQ(0.0, copysign(0.0, 1.0));
+ ASSERT_DOUBLE_EQ(-0.0, copysign(0.0, -1.0));
+ ASSERT_DOUBLE_EQ(2.0, copysign(2.0, 1.0));
+ ASSERT_DOUBLE_EQ(-2.0, copysign(2.0, -1.0));
}
TEST(math, copysignf) {
@@ -983,16 +1021,16 @@
}
TEST(math, copysignl) {
- ASSERT_FLOAT_EQ(0.0f, copysignl(0.0, 1.0));
- ASSERT_FLOAT_EQ(-0.0f, copysignl(0.0, -1.0));
- ASSERT_FLOAT_EQ(2.0f, copysignl(2.0, 1.0));
- ASSERT_FLOAT_EQ(-2.0f, copysignl(2.0, -1.0));
+ ASSERT_DOUBLE_EQ(0.0l, copysignl(0.0l, 1.0l));
+ ASSERT_DOUBLE_EQ(-0.0l, copysignl(0.0l, -1.0l));
+ ASSERT_DOUBLE_EQ(2.0l, copysignl(2.0l, 1.0l));
+ ASSERT_DOUBLE_EQ(-2.0l, copysignl(2.0l, -1.0l));
}
TEST(math, significand) {
- ASSERT_FLOAT_EQ(0.0, significand(0.0));
- ASSERT_FLOAT_EQ(1.2, significand(1.2));
- ASSERT_FLOAT_EQ(1.5375, significand(12.3));
+ ASSERT_DOUBLE_EQ(0.0, significand(0.0));
+ ASSERT_DOUBLE_EQ(1.2, significand(1.2));
+ ASSERT_DOUBLE_EQ(1.5375, significand(12.3));
}
TEST(math, significandf) {
@@ -1004,13 +1042,13 @@
extern "C" long double significandl(long double); // BSD's <math.h> doesn't declare this.
TEST(math, significandl) {
- ASSERT_FLOAT_EQ(0.0, significandl(0.0));
- ASSERT_FLOAT_EQ(1.2, significandl(1.2));
- ASSERT_FLOAT_EQ(1.5375, significandl(12.3));
+ ASSERT_DOUBLE_EQ(0.0l, significandl(0.0l));
+ ASSERT_DOUBLE_EQ(1.2l, significandl(1.2l));
+ ASSERT_DOUBLE_EQ(1.5375l, significandl(12.3l));
}
TEST(math, scalb) {
- ASSERT_FLOAT_EQ(12.0, scalb(3.0, 2.0));
+ ASSERT_DOUBLE_EQ(12.0, scalb(3.0, 2.0));
}
TEST(math, scalbf) {
@@ -1018,7 +1056,7 @@
}
TEST(math, scalbln) {
- ASSERT_FLOAT_EQ(12.0, scalbln(3.0, 2L));
+ ASSERT_DOUBLE_EQ(12.0, scalbln(3.0, 2L));
}
TEST(math, scalblnf) {
@@ -1026,11 +1064,11 @@
}
TEST(math, scalblnl) {
- ASSERT_FLOAT_EQ(12.0, scalblnl(3.0, 2L));
+ ASSERT_DOUBLE_EQ(12.0l, scalblnl(3.0l, 2L));
}
TEST(math, scalbn) {
- ASSERT_FLOAT_EQ(12.0, scalbn(3.0, 2));
+ ASSERT_DOUBLE_EQ(12.0, scalbn(3.0, 2));
}
TEST(math, scalbnf) {
@@ -1038,11 +1076,11 @@
}
TEST(math, scalbnl) {
- ASSERT_FLOAT_EQ(12.0, scalbnl(3.0, 2));
+ ASSERT_DOUBLE_EQ(12.0l, scalbnl(3.0l, 2));
}
TEST(math, gamma) {
- ASSERT_FLOAT_EQ(log(24.0), gamma(5.0));
+ ASSERT_DOUBLE_EQ(log(24.0), gamma(5.0));
}
TEST(math, gammaf) {
@@ -1052,7 +1090,7 @@
TEST(math, gamma_r) {
#if defined(__BIONIC__)
int sign;
- ASSERT_FLOAT_EQ(log(24.0), gamma_r(5.0, &sign));
+ ASSERT_DOUBLE_EQ(log(24.0), gamma_r(5.0, &sign));
ASSERT_EQ(1, sign);
#else // __BIONIC__
GTEST_LOG_(INFO) << "This test does nothing.\n";
@@ -1070,7 +1108,7 @@
}
TEST(math, lgamma) {
- ASSERT_FLOAT_EQ(log(24.0), lgamma(5.0));
+ ASSERT_DOUBLE_EQ(log(24.0), lgamma(5.0));
}
TEST(math, lgammaf) {
@@ -1078,12 +1116,12 @@
}
TEST(math, lgammal) {
- ASSERT_FLOAT_EQ(logl(24.0), lgammal(5.0));
+ ASSERT_DOUBLE_EQ(logl(24.0l), lgammal(5.0l));
}
TEST(math, lgamma_r) {
int sign;
- ASSERT_FLOAT_EQ(log(24.0), lgamma_r(5.0, &sign));
+ ASSERT_DOUBLE_EQ(log(24.0), lgamma_r(5.0, &sign));
ASSERT_EQ(1, sign);
}
@@ -1094,7 +1132,7 @@
}
TEST(math, tgamma) {
- ASSERT_FLOAT_EQ(24.0, tgamma(5.0));
+ ASSERT_DOUBLE_EQ(24.0, tgamma(5.0));
}
TEST(math, tgammaf) {
@@ -1102,12 +1140,12 @@
}
TEST(math, tgammal) {
- ASSERT_FLOAT_EQ(24.0, tgammal(5.0));
+ ASSERT_DOUBLE_EQ(24.0l, tgammal(5.0l));
}
TEST(math, j0) {
- ASSERT_FLOAT_EQ(1.0, j0(0.0));
- ASSERT_FLOAT_EQ(0.76519769, j0(1.0));
+ ASSERT_DOUBLE_EQ(1.0, j0(0.0));
+ ASSERT_DOUBLE_EQ(0.76519768655796661, j0(1.0));
}
TEST(math, j0f) {
@@ -1116,8 +1154,8 @@
}
TEST(math, j1) {
- ASSERT_FLOAT_EQ(0.0, j1(0.0));
- ASSERT_FLOAT_EQ(0.44005057, j1(1.0));
+ ASSERT_DOUBLE_EQ(0.0, j1(0.0));
+ ASSERT_DOUBLE_EQ(0.44005058574493355, j1(1.0));
}
TEST(math, j1f) {
@@ -1126,8 +1164,8 @@
}
TEST(math, jn) {
- ASSERT_FLOAT_EQ(0.0, jn(4, 0.0));
- ASSERT_FLOAT_EQ(0.0024766389, jn(4, 1.0));
+ ASSERT_DOUBLE_EQ(0.0, jn(4, 0.0));
+ ASSERT_DOUBLE_EQ(0.0024766389641099553, jn(4, 1.0));
}
TEST(math, jnf) {
@@ -1136,8 +1174,8 @@
}
TEST(math, y0) {
- ASSERT_FLOAT_EQ(-HUGE_VAL, y0(0.0));
- ASSERT_FLOAT_EQ(0.088256963, y0(1.0));
+ ASSERT_DOUBLE_EQ(-HUGE_VAL, y0(0.0));
+ ASSERT_DOUBLE_EQ(0.08825696421567697, y0(1.0));
}
TEST(math, y0f) {
@@ -1146,8 +1184,8 @@
}
TEST(math, y1) {
- ASSERT_FLOAT_EQ(-HUGE_VAL, y1(0.0));
- ASSERT_FLOAT_EQ(-0.78121281, y1(1.0));
+ ASSERT_DOUBLE_EQ(-HUGE_VAL, y1(0.0));
+ ASSERT_DOUBLE_EQ(-0.78121282130028868, y1(1.0));
}
TEST(math, y1f) {
@@ -1156,8 +1194,8 @@
}
TEST(math, yn) {
- ASSERT_FLOAT_EQ(-HUGE_VAL, yn(4, 0.0));
- ASSERT_FLOAT_EQ(-33.278423, yn(4, 1.0));
+ ASSERT_DOUBLE_EQ(-HUGE_VAL, yn(4, 0.0));
+ ASSERT_DOUBLE_EQ(-33.278423028972114, yn(4, 1.0));
}
TEST(math, ynf) {
@@ -1168,7 +1206,7 @@
TEST(math, frexp) {
int exp;
double dr = frexp(1024.0, &exp);
- ASSERT_FLOAT_EQ(1024.0, scalbn(dr, exp));
+ ASSERT_DOUBLE_EQ(1024.0, scalbn(dr, exp));
}
TEST(math, frexpf) {
@@ -1179,15 +1217,15 @@
TEST(math, frexpl) {
int exp;
- long double ldr = frexpl(1024.0, &exp);
- ASSERT_FLOAT_EQ(1024.0, scalbnl(ldr, exp));
+ long double ldr = frexpl(1024.0l, &exp);
+ ASSERT_DOUBLE_EQ(1024.0l, scalbnl(ldr, exp));
}
TEST(math, modf) {
double di;
double df = modf(123.456, &di);
- ASSERT_FLOAT_EQ(123.0, di);
- ASSERT_FLOAT_EQ(0.456, df);
+ ASSERT_DOUBLE_EQ(123.0, di);
+ ASSERT_DOUBLE_EQ(0.45600000000000307, df);
}
TEST(math, modff) {
@@ -1199,16 +1237,16 @@
TEST(math, modfl) {
long double ldi;
- long double ldf = modfl(123.456, &ldi);
- ASSERT_FLOAT_EQ(123.0, ldi);
- ASSERT_FLOAT_EQ(0.456, ldf);
+ long double ldf = modfl(123.456l, &ldi);
+ ASSERT_DOUBLE_EQ(123.0l, ldi);
+ ASSERT_DOUBLE_EQ(0.45600000000000002l, ldf);
}
TEST(math, remquo) {
int q;
double d = remquo(13.0, 4.0, &q);
ASSERT_EQ(3, q);
- ASSERT_FLOAT_EQ(1.0, d);
+ ASSERT_DOUBLE_EQ(1.0, d);
}
TEST(math, remquof) {
@@ -1220,9 +1258,9 @@
TEST(math, remquol) {
int q;
- long double ld = remquol(13.0, 4.0, &q);
- ASSERT_EQ(3, q);
- ASSERT_FLOAT_EQ(1.0, ld);
+ long double ld = remquol(13.0l, 4.0l, &q);
+ ASSERT_DOUBLE_EQ(3l, q);
+ ASSERT_DOUBLE_EQ(1.0l, ld);
}
// https://code.google.com/p/android/issues/detail?id=6697
