generic-library/bionic
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge "Fix LP64 libm for 128-bit long doubles"

Browse Source
This commit is contained in:
Calin Juravle 2014-03-24 12:13:24 +00:00 committed by Gerrit Code Review
commit 989a80b750
22 changed files with 2982 additions and 270 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

97
libm/Android.mk
View File

@ -179,42 +179,55 @@ libm_common_src_files += \
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 @@ libm_common_cflags := \
-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_SYSTEM_SHARED_LIBRARIES := libc
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)

31
libm/amd64/_fpmath.h
View File

@ -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)

37
libm/arm64/_fpmath.h
View File

@ -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)

26
libm/fake_long_double.c
View File

@ -22,6 +22,18 @@ int (isinf)(double a1) { return __isinf(a1); }
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 atanl(long double a1) { return atan(a1); }
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 frexpl(long double a1, int* a2) { return frexp(a1, a2); }
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 double logl(long double a1) { return log(a1); }
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__

9
libm/fpmath.h
View File

@ -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 @@ union IEEEd2bits {
* 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_

140
libm/upstream-freebsd/lib/msun/ld128/e_rem_pio2l.h Normal file
View File

@ -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;
}

100
libm/upstream-freebsd/lib/msun/ld128/invtrig.c Normal file
View File

@ -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;

113
libm/upstream-freebsd/lib/msun/ld128/invtrig.h Normal file
View File

@ -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])))))))))));
}

61
libm/upstream-freebsd/lib/msun/ld128/k_cosl.c Normal file
View File

@ -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));
}

59
libm/upstream-freebsd/lib/msun/ld128/k_sinl.c Normal file
View File

@ -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);
}

119
libm/upstream-freebsd/lib/msun/ld128/k_tanl.c Normal file
View File

@ -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);
}
}

427
libm/upstream-freebsd/lib/msun/ld128/s_exp2l.c Normal file
View File

@ -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);
}
}

494
libm/upstream-freebsd/lib/msun/ld128/s_expl.c Normal file
View File

@ -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);
}

737
libm/upstream-freebsd/lib/msun/ld128/s_logl.c Normal file
View File

@ -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 */

46
libm/upstream-freebsd/lib/msun/ld128/s_nanl.c Normal file
View File

@ -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);
}

89
libm/upstream-freebsd/lib/msun/src/e_acoshl.c Normal file
View File

@ -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)));
}
}

74
libm/upstream-freebsd/lib/msun/src/e_atanhl.c Normal file
View File

@ -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);
}

69
libm/upstream-freebsd/lib/msun/src/imprecise.c Normal file
View File

@ -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);

91
libm/upstream-freebsd/lib/msun/src/s_asinhl.c Normal file
View File

@ -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);
}

2
libm/upstream-freebsd/lib/msun/src/s_remquol.c
View File

@ -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.
* ====================================================
*/

1
libm/upstream-freebsd/lib/msun/src/s_rintl.c
View File

@ -29,6 +29,7 @@ __FBSDID("$FreeBSD$");
#include <float.h>
#include <math.h>
#include <stdint.h>
#include "fpmath.h"

430
tests/math_test.cpp
View File

File diff suppressed because it is too large Load Diff