diff options
author | Jeremy Andrews <athenian200@outlook.com> | 2021-10-10 08:53:14 -0500 |
---|---|---|
committer | Jeremy Andrews <athenian200@outlook.com> | 2021-10-10 08:53:14 -0500 |
commit | 0ea4019528592afbe0ae918350dd259ef6e5719a (patch) | |
tree | 061bb883d883c3b2c045e68084b11f1d92980307 /modules | |
parent | 31dfa8ed97209ed4004f071c552efe499378ded6 (diff) | |
download | aura-central-0ea4019528592afbe0ae918350dd259ef6e5719a.tar.gz |
Issue %3003 - Move fdlibm to libs/
Diffstat (limited to 'modules')
59 files changed, 0 insertions, 7688 deletions
diff --git a/modules/fdlibm/README.mozilla b/modules/fdlibm/README.mozilla deleted file mode 100644 index d25ef8b8e..000000000 --- a/modules/fdlibm/README.mozilla +++ /dev/null @@ -1,19 +0,0 @@ -This is the fdlibm library imported from -https://github.com/freebsd/freebsd - -Upstream code can be viewed at - https://github.com/freebsd/freebsd/tree/master/lib/msun/src - -Each file is downloaded separately, as cloning whole repository takes so much -resources. - -The in-tree copy is updated by running - sh update.sh -or - sh update.sh <sha-commit> -from within the modules/fdlibm directory. - -Current version: [commit cf4707bb2f78ecf56ba350bdc24e3135b4339622 (2019-09-25T18:50:57Z)]. - -patches 01-18 fixes files to be usable within mozilla-central tree. -See https://bugzilla.mozilla.org/show_bug.cgi?id=933257 diff --git a/modules/fdlibm/import.sh b/modules/fdlibm/import.sh deleted file mode 100644 index 11acb064b..000000000 --- a/modules/fdlibm/import.sh +++ /dev/null @@ -1,111 +0,0 @@ -#!/bin/sh - -set -e - -BASE_URL=https://raw.githubusercontent.com/freebsd/freebsd/"${1}"/lib/msun/src - -download_source() { - REMOTE_FILENAME=$1 - LOCAL_FILENAME=$2 - while true; do - curl -o "src/${LOCAL_FILENAME}" "${BASE_URL}/${REMOTE_FILENAME}" && break - done -} - -mkdir -p src - -# headers -download_source math.h fdlibm.h -download_source math_private.h math_private.h - -# Math.acos -download_source e_acos.c e_acos.cpp - -# Math.acosh -download_source e_acosh.c e_acosh.cpp - -# Math.asin -download_source e_asin.c e_asin.cpp - -# Math.asinh -download_source s_asinh.c s_asinh.cpp - -# Math.atan -download_source s_atan.c s_atan.cpp - -# Math.atanh -download_source e_atanh.c e_atanh.cpp - -# Math.atan2 -download_source e_atan2.c e_atan2.cpp - -# Math.cbrt -download_source s_cbrt.c s_cbrt.cpp - -# Math.ceil -download_source s_ceil.c s_ceil.cpp -download_source s_ceilf.c s_ceilf.cpp - -# Math.cos (not used due to poor performance) - -# Math.cosh -download_source e_cosh.c e_cosh.cpp - -# Math.exp -download_source e_exp.c e_exp.cpp - -# Math.expm1 -download_source s_expm1.c s_expm1.cpp - -# Math.floor and Math.round -download_source s_floor.c s_floor.cpp - -# Math.fround -download_source s_floorf.c s_floorf.cpp - -# Math.hypot -download_source e_hypot.c e_hypot.cpp - -# Math.log -download_source e_log.c e_log.cpp - -# Math.log1p -download_source s_log1p.c s_log1p.cpp - -# Math.log10 -download_source e_log10.c e_log10.cpp -download_source k_log.h k_log.h - -# Math.log2 -download_source e_log2.c e_log2.cpp - -# Math.pow (not used due to poor performance) - -# Math.sin (not used due to poor performance) - -# Math.sinh -download_source e_sinh.c e_sinh.cpp - -# Math.sqrt (not used due to poor performance) - -# Math.tan (not used due to poor performance) - -# Math.tanh -download_source s_tanh.c s_tanh.cpp - -# Math.trunc -download_source s_trunc.c s_trunc.cpp -download_source s_truncf.c s_truncf.cpp - -# dependencies -download_source k_exp.c k_exp.cpp -download_source s_copysign.c s_copysign.cpp -download_source s_fabs.c s_fabs.cpp -download_source s_scalbn.c s_scalbn.cpp - -# These are not not used in Math.* functions, but used internally. -download_source e_pow.c e_pow.cpp - -download_source s_nearbyint.c s_nearbyint.cpp -download_source s_rint.c s_rint.cpp -download_source s_rintf.c s_rintf.cpp diff --git a/modules/fdlibm/moz.build b/modules/fdlibm/moz.build deleted file mode 100644 index f8d8f81d0..000000000 --- a/modules/fdlibm/moz.build +++ /dev/null @@ -1,6 +0,0 @@ -# -*- Mode: python; indent-tabs-mode: nil; tab-width: 40 -*- -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. - -DIRS += ['src'] diff --git a/modules/fdlibm/patches/01_remove_unused_declarations_from_fdlibm_h.patch b/modules/fdlibm/patches/01_remove_unused_declarations_from_fdlibm_h.patch deleted file mode 100644 index 4e5774ff0..000000000 --- a/modules/fdlibm/patches/01_remove_unused_declarations_from_fdlibm_h.patch +++ /dev/null @@ -1,507 +0,0 @@ -diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h ---- a/modules/fdlibm/src/fdlibm.h -+++ b/modules/fdlibm/src/fdlibm.h -@@ -12,499 +12,49 @@ - /* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - - #ifndef _MATH_H_ - #define _MATH_H_ - --#include <sys/cdefs.h> --#include <sys/_types.h> --#include <machine/_limits.h> -- --/* -- * ANSI/POSIX -- */ --extern const union __infinity_un { -- unsigned char __uc[8]; -- double __ud; --} __infinity; -- --extern const union __nan_un { -- unsigned char __uc[sizeof(float)]; -- float __uf; --} __nan; -- --#if __GNUC_PREREQ__(3, 3) || (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 800) --#define __MATH_BUILTIN_CONSTANTS --#endif -- --#if __GNUC_PREREQ__(3, 0) && !defined(__INTEL_COMPILER) --#define __MATH_BUILTIN_RELOPS --#endif -- --#ifdef __MATH_BUILTIN_CONSTANTS --#define HUGE_VAL __builtin_huge_val() --#else --#define HUGE_VAL (__infinity.__ud) --#endif -- --#if __ISO_C_VISIBLE >= 1999 --#define FP_ILOGB0 (-__INT_MAX) --#define FP_ILOGBNAN __INT_MAX -- --#ifdef __MATH_BUILTIN_CONSTANTS --#define HUGE_VALF __builtin_huge_valf() --#define HUGE_VALL __builtin_huge_vall() --#define INFINITY __builtin_inff() --#define NAN __builtin_nanf("") --#else --#define HUGE_VALF (float)HUGE_VAL --#define HUGE_VALL (long double)HUGE_VAL --#define INFINITY HUGE_VALF --#define NAN (__nan.__uf) --#endif /* __MATH_BUILTIN_CONSTANTS */ -- --#define MATH_ERRNO 1 --#define MATH_ERREXCEPT 2 --#define math_errhandling MATH_ERREXCEPT -- --#define FP_FAST_FMAF 1 -- --/* Symbolic constants to classify floating point numbers. */ --#define FP_INFINITE 0x01 --#define FP_NAN 0x02 --#define FP_NORMAL 0x04 --#define FP_SUBNORMAL 0x08 --#define FP_ZERO 0x10 -- --#if (__STDC_VERSION__ >= 201112L && defined(__clang__)) || \ -- __has_extension(c_generic_selections) --#define __fp_type_select(x, f, d, ld) _Generic((x), \ -- float: f(x), \ -- double: d(x), \ -- long double: ld(x), \ -- volatile float: f(x), \ -- volatile double: d(x), \ -- volatile long double: ld(x), \ -- volatile const float: f(x), \ -- volatile const double: d(x), \ -- volatile const long double: ld(x), \ -- const float: f(x), \ -- const double: d(x), \ -- const long double: ld(x)) --#elif __GNUC_PREREQ__(3, 1) && !defined(__cplusplus) --#define __fp_type_select(x, f, d, ld) __builtin_choose_expr( \ -- __builtin_types_compatible_p(__typeof(x), long double), ld(x), \ -- __builtin_choose_expr( \ -- __builtin_types_compatible_p(__typeof(x), double), d(x), \ -- __builtin_choose_expr( \ -- __builtin_types_compatible_p(__typeof(x), float), f(x), (void)0))) --#else --#define __fp_type_select(x, f, d, ld) \ -- ((sizeof(x) == sizeof(float)) ? f(x) \ -- : (sizeof(x) == sizeof(double)) ? d(x) \ -- : ld(x)) --#endif -- --#define fpclassify(x) \ -- __fp_type_select(x, __fpclassifyf, __fpclassifyd, __fpclassifyl) --#define isfinite(x) __fp_type_select(x, __isfinitef, __isfinite, __isfinitel) --#define isinf(x) __fp_type_select(x, __isinff, __isinf, __isinfl) --#define isnan(x) \ -- __fp_type_select(x, __inline_isnanf, __inline_isnan, __inline_isnanl) --#define isnormal(x) __fp_type_select(x, __isnormalf, __isnormal, __isnormall) -- --#ifdef __MATH_BUILTIN_RELOPS --#define isgreater(x, y) __builtin_isgreater((x), (y)) --#define isgreaterequal(x, y) __builtin_isgreaterequal((x), (y)) --#define isless(x, y) __builtin_isless((x), (y)) --#define islessequal(x, y) __builtin_islessequal((x), (y)) --#define islessgreater(x, y) __builtin_islessgreater((x), (y)) --#define isunordered(x, y) __builtin_isunordered((x), (y)) --#else --#define isgreater(x, y) (!isunordered((x), (y)) && (x) > (y)) --#define isgreaterequal(x, y) (!isunordered((x), (y)) && (x) >= (y)) --#define isless(x, y) (!isunordered((x), (y)) && (x) < (y)) --#define islessequal(x, y) (!isunordered((x), (y)) && (x) <= (y)) --#define islessgreater(x, y) (!isunordered((x), (y)) && \ -- ((x) > (y) || (y) > (x))) --#define isunordered(x, y) (isnan(x) || isnan(y)) --#endif /* __MATH_BUILTIN_RELOPS */ -- --#define signbit(x) __fp_type_select(x, __signbitf, __signbit, __signbitl) -- --typedef __double_t double_t; --typedef __float_t float_t; --#endif /* __ISO_C_VISIBLE >= 1999 */ -- --/* -- * XOPEN/SVID -- */ --#if __BSD_VISIBLE || __XSI_VISIBLE --#define M_E 2.7182818284590452354 /* e */ --#define M_LOG2E 1.4426950408889634074 /* log 2e */ --#define M_LOG10E 0.43429448190325182765 /* log 10e */ --#define M_LN2 0.69314718055994530942 /* log e2 */ --#define M_LN10 2.30258509299404568402 /* log e10 */ --#define M_PI 3.14159265358979323846 /* pi */ --#define M_PI_2 1.57079632679489661923 /* pi/2 */ --#define M_PI_4 0.78539816339744830962 /* pi/4 */ --#define M_1_PI 0.31830988618379067154 /* 1/pi */ --#define M_2_PI 0.63661977236758134308 /* 2/pi */ --#define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */ --#define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ --#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ -- --#define MAXFLOAT ((float)3.40282346638528860e+38) --extern int signgam; --#endif /* __BSD_VISIBLE || __XSI_VISIBLE */ -- --#if __BSD_VISIBLE --#if 0 --/* Old value from 4.4BSD-Lite math.h; this is probably better. */ --#define HUGE HUGE_VAL --#else --#define HUGE MAXFLOAT --#endif --#endif /* __BSD_VISIBLE */ -- --/* -- * Most of these functions depend on the rounding mode and have the side -- * effect of raising floating-point exceptions, so they are not declared -- * as __pure2. In C99, FENV_ACCESS affects the purity of these functions. -- */ --__BEGIN_DECLS --/* -- * ANSI/POSIX -- */ --int __fpclassifyd(double) __pure2; --int __fpclassifyf(float) __pure2; --int __fpclassifyl(long double) __pure2; --int __isfinitef(float) __pure2; --int __isfinite(double) __pure2; --int __isfinitel(long double) __pure2; --int __isinff(float) __pure2; --int __isinf(double) __pure2; --int __isinfl(long double) __pure2; --int __isnormalf(float) __pure2; --int __isnormal(double) __pure2; --int __isnormall(long double) __pure2; --int __signbit(double) __pure2; --int __signbitf(float) __pure2; --int __signbitl(long double) __pure2; -- --static __inline int --__inline_isnan(__const double __x) --{ -- -- return (__x != __x); --} -- --static __inline int --__inline_isnanf(__const float __x) --{ -- -- return (__x != __x); --} -- --static __inline int --__inline_isnanl(__const long double __x) --{ -- -- return (__x != __x); --} -- --/* -- * Version 2 of the Single UNIX Specification (UNIX98) defined isnan() and -- * isinf() as functions taking double. C99, and the subsequent POSIX revisions -- * (SUSv3, POSIX.1-2001, define it as a macro that accepts any real floating -- * point type. If we are targeting SUSv2 and C99 or C11 (or C++11) then we -- * expose the newer definition, assuming that the language spec takes -- * precedence over the operating system interface spec. -- */ --#if __XSI_VISIBLE > 0 && __XSI_VISIBLE < 600 && __ISO_C_VISIBLE < 1999 --#undef isinf --#undef isnan --int isinf(double); --int isnan(double); --#endif -- - double acos(double); - double asin(double); - double atan(double); - double atan2(double, double); --double cos(double); --double sin(double); --double tan(double); - - double cosh(double); - double sinh(double); - double tanh(double); - - double exp(double); --double frexp(double, int *); /* fundamentally !__pure2 */ --double ldexp(double, int); - double log(double); - double log10(double); --double modf(double, double *); /* fundamentally !__pure2 */ - - double pow(double, double); --double sqrt(double); -+double fabs(double); - --double ceil(double); --double fabs(double) __pure2; - double floor(double); --double fmod(double, double); -+double trunc(double); -+double ceil(double); - --/* -- * These functions are not in C90. -- */ --#if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE - double acosh(double); - double asinh(double); - double atanh(double); - double cbrt(double); --double erf(double); --double erfc(double); --double exp2(double); - double expm1(double); --double fma(double, double, double); - double hypot(double, double); --int ilogb(double) __pure2; --double lgamma(double); --long long llrint(double); --long long llround(double); - double log1p(double); - double log2(double); --double logb(double); --long lrint(double); --long lround(double); --double nan(const char *) __pure2; --double nextafter(double, double); --double remainder(double, double); --double remquo(double, double, int *); - double rint(double); --#endif /* __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE */ -- --#if __BSD_VISIBLE || __XSI_VISIBLE --double j0(double); --double j1(double); --double jn(int, double); --double y0(double); --double y1(double); --double yn(int, double); -- --#if __XSI_VISIBLE <= 500 || __BSD_VISIBLE --double gamma(double); --#endif -- --#if __XSI_VISIBLE <= 600 || __BSD_VISIBLE --double scalb(double, double); --#endif --#endif /* __BSD_VISIBLE || __XSI_VISIBLE */ -- --#if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 --double copysign(double, double) __pure2; --double fdim(double, double); --double fmax(double, double) __pure2; --double fmin(double, double) __pure2; -+double copysign(double, double); - double nearbyint(double); --double round(double); --double scalbln(double, long); - double scalbn(double, int); --double tgamma(double); --double trunc(double); --#endif -- --/* -- * BSD math library entry points -- */ --#if __BSD_VISIBLE --double drem(double, double); --int finite(double) __pure2; --int isnanf(float) __pure2; -- --/* -- * Reentrant version of gamma & lgamma; passes signgam back by reference -- * as the second argument; user must allocate space for signgam. -- */ --double gamma_r(double, int *); --double lgamma_r(double, int *); -- --/* -- * IEEE Test Vector -- */ --double significand(double); --#endif /* __BSD_VISIBLE */ -- --/* float versions of ANSI/POSIX functions */ --#if __ISO_C_VISIBLE >= 1999 --float acosf(float); --float asinf(float); --float atanf(float); --float atan2f(float, float); --float cosf(float); --float sinf(float); --float tanf(float); -- --float coshf(float); --float sinhf(float); --float tanhf(float); -- --float exp2f(float); --float expf(float); --float expm1f(float); --float frexpf(float, int *); /* fundamentally !__pure2 */ --int ilogbf(float) __pure2; --float ldexpf(float, int); --float log10f(float); --float log1pf(float); --float log2f(float); --float logf(float); --float modff(float, float *); /* fundamentally !__pure2 */ -- --float powf(float, float); --float sqrtf(float); - - float ceilf(float); --float fabsf(float) __pure2; - float floorf(float); --float fmodf(float, float); --float roundf(float); -- --float erff(float); --float erfcf(float); --float hypotf(float, float); --float lgammaf(float); --float tgammaf(float); - --float acoshf(float); --float asinhf(float); --float atanhf(float); --float cbrtf(float); --float logbf(float); --float copysignf(float, float) __pure2; --long long llrintf(float); --long long llroundf(float); --long lrintf(float); --long lroundf(float); --float nanf(const char *) __pure2; - float nearbyintf(float); --float nextafterf(float, float); --float remainderf(float, float); --float remquof(float, float, int *); - float rintf(float); --float scalblnf(float, long); --float scalbnf(float, int); - float truncf(float); - --float fdimf(float, float); --float fmaf(float, float, float); --float fmaxf(float, float) __pure2; --float fminf(float, float) __pure2; --#endif -- --/* -- * float versions of BSD math library entry points -- */ --#if __BSD_VISIBLE --float dremf(float, float); --int finitef(float) __pure2; --float gammaf(float); --float j0f(float); --float j1f(float); --float jnf(int, float); --float scalbf(float, float); --float y0f(float); --float y1f(float); --float ynf(int, float); -- --/* -- * Float versions of reentrant version of gamma & lgamma; passes -- * signgam back by reference as the second argument; user must -- * allocate space for signgam. -- */ --float gammaf_r(float, int *); --float lgammaf_r(float, int *); -- --/* -- * float version of IEEE Test Vector -- */ --float significandf(float); --#endif /* __BSD_VISIBLE */ -- --/* -- * long double versions of ISO/POSIX math functions -- */ --#if __ISO_C_VISIBLE >= 1999 --long double acoshl(long double); --long double acosl(long double); --long double asinhl(long double); --long double asinl(long double); --long double atan2l(long double, long double); --long double atanhl(long double); --long double atanl(long double); --long double cbrtl(long double); --long double ceill(long double); --long double copysignl(long double, long double) __pure2; --long double coshl(long double); --long double cosl(long double); --long double erfcl(long double); --long double erfl(long double); --long double exp2l(long double); --long double expl(long double); --long double expm1l(long double); --long double fabsl(long double) __pure2; --long double fdiml(long double, long double); --long double floorl(long double); --long double fmal(long double, long double, long double); --long double fmaxl(long double, long double) __pure2; --long double fminl(long double, long double) __pure2; --long double fmodl(long double, long double); --long double frexpl(long double, int *); /* fundamentally !__pure2 */ --long double hypotl(long double, long double); --int ilogbl(long double) __pure2; --long double ldexpl(long double, int); --long double lgammal(long double); --long long llrintl(long double); --long long llroundl(long double); --long double log10l(long double); --long double log1pl(long double); --long double log2l(long double); --long double logbl(long double); --long double logl(long double); --long lrintl(long double); --long lroundl(long double); --long double modfl(long double, long double *); /* fundamentally !__pure2 */ --long double nanl(const char *) __pure2; --long double nearbyintl(long double); --long double nextafterl(long double, long double); --double nexttoward(double, long double); --float nexttowardf(float, long double); --long double nexttowardl(long double, long double); --long double powl(long double, long double); --long double remainderl(long double, long double); --long double remquol(long double, long double, int *); --long double rintl(long double); --long double roundl(long double); --long double scalblnl(long double, long); --long double scalbnl(long double, int); --long double sinhl(long double); --long double sinl(long double); --long double sqrtl(long double); --long double tanhl(long double); --long double tanl(long double); --long double tgammal(long double); --long double truncl(long double); --#endif /* __ISO_C_VISIBLE >= 1999 */ -- --#if __BSD_VISIBLE --long double lgammal_r(long double, int *); --void sincos(double, double *, double *); --void sincosf(float, float *, float *); --void sincosl(long double, long double *, long double *); --#endif -- --__END_DECLS -- - #endif /* !_MATH_H_ */ diff --git a/modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch b/modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch deleted file mode 100644 index f103af128..000000000 --- a/modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch +++ /dev/null @@ -1,35 +0,0 @@ -diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h ---- a/modules/fdlibm/src/fdlibm.h -+++ b/modules/fdlibm/src/fdlibm.h -@@ -9,18 +9,18 @@ - * ==================================================== - */ - - /* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - --#ifndef _MATH_H_ --#define _MATH_H_ -+#ifndef mozilla_imported_fdlibm_h -+#define mozilla_imported_fdlibm_h - - double acos(double); - double asin(double); - double atan(double); - double atan2(double, double); - - double cosh(double); - double sinh(double); -@@ -52,9 +52,9 @@ double scalbn(double, int); - - float ceilf(float); - float floorf(float); - - float nearbyintf(float); - float rintf(float); - float truncf(float); - --#endif /* !_MATH_H_ */ -+#endif /* mozilla_imported_fdlibm_h */ diff --git a/modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch b/modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch deleted file mode 100644 index b6cd7ab7f..000000000 --- a/modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch +++ /dev/null @@ -1,34 +0,0 @@ -diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h ---- a/modules/fdlibm/src/fdlibm.h -+++ b/modules/fdlibm/src/fdlibm.h -@@ -12,16 +12,18 @@ - /* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - - #ifndef mozilla_imported_fdlibm_h - #define mozilla_imported_fdlibm_h - -+namespace fdlibm { -+ - double acos(double); - double asin(double); - double atan(double); - double atan2(double, double); - - double cosh(double); - double sinh(double); - double tanh(double); -@@ -52,9 +54,11 @@ double scalbn(double, int); - - float ceilf(float); - float floorf(float); - - float nearbyintf(float); - float rintf(float); - float truncf(float); - -+} /* namespace fdlibm */ -+ - #endif /* mozilla_imported_fdlibm_h */ diff --git a/modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch b/modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch deleted file mode 100644 index 153018902..000000000 --- a/modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch +++ /dev/null @@ -1,695 +0,0 @@ -diff --git a/modules/fdlibm/src/e_acos.cpp b/modules/fdlibm/src/e_acos.cpp ---- a/modules/fdlibm/src/e_acos.cpp -+++ b/modules/fdlibm/src/e_acos.cpp -@@ -35,17 +35,16 @@ - * if x is NaN, return x itself; - * if |x|>1, return NaN with invalid signal. - * - * Function needed: sqrt - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ - pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ - pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ - static volatile double - pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */ -diff --git a/modules/fdlibm/src/e_acosh.cpp b/modules/fdlibm/src/e_acosh.cpp ---- a/modules/fdlibm/src/e_acosh.cpp -+++ b/modules/fdlibm/src/e_acosh.cpp -@@ -26,17 +26,16 @@ - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - one = 1.0, - ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ - - double - __ieee754_acosh(double x) -diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp ---- a/modules/fdlibm/src/e_asin.cpp -+++ b/modules/fdlibm/src/e_asin.cpp -@@ -41,17 +41,16 @@ - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN with invalid signal. - * - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ - huge = 1.000e+300, - pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ - pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ - pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ -diff --git a/modules/fdlibm/src/e_atan2.cpp b/modules/fdlibm/src/e_atan2.cpp ---- a/modules/fdlibm/src/e_atan2.cpp -+++ b/modules/fdlibm/src/e_atan2.cpp -@@ -39,17 +39,16 @@ - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static volatile double - tiny = 1.0e-300; - static const double - zero = 0.0, - pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ - pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ -diff --git a/modules/fdlibm/src/e_atanh.cpp b/modules/fdlibm/src/e_atanh.cpp ---- a/modules/fdlibm/src/e_atanh.cpp -+++ b/modules/fdlibm/src/e_atanh.cpp -@@ -30,17 +30,16 @@ - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double one = 1.0, huge = 1e300; - static const double zero = 0.0; - - double - __ieee754_atanh(double x) - { -diff --git a/modules/fdlibm/src/e_cosh.cpp b/modules/fdlibm/src/e_cosh.cpp ---- a/modules/fdlibm/src/e_cosh.cpp -+++ b/modules/fdlibm/src/e_cosh.cpp -@@ -32,17 +32,16 @@ - * - * Special cases: - * cosh(x) is |x| if x is +INF, -INF, or NaN. - * only cosh(0)=1 is exact for finite x. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double one = 1.0, half=0.5, huge = 1.0e300; - - double - __ieee754_cosh(double x) - { - double t,w; -diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp ---- a/modules/fdlibm/src/e_exp.cpp -+++ b/modules/fdlibm/src/e_exp.cpp -@@ -73,17 +73,16 @@ - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - one = 1.0, - halF[2] = {0.5,-0.5,}, - o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ - u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ - ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ -diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp ---- a/modules/fdlibm/src/e_hypot.cpp -+++ b/modules/fdlibm/src/e_hypot.cpp -@@ -43,17 +43,16 @@ - * - * Accuracy: - * hypot(x,y) returns sqrt(x^2+y^2) with error less - * than 1 ulps (units in the last place) - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - double - __ieee754_hypot(double x, double y) - { - double a,b,t1,t2,y1,y2,w; - int32_t j,k,ha,hb; - -diff --git a/modules/fdlibm/src/e_log.cpp b/modules/fdlibm/src/e_log.cpp ---- a/modules/fdlibm/src/e_log.cpp -+++ b/modules/fdlibm/src/e_log.cpp -@@ -62,17 +62,16 @@ - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ - ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ - two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ - Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ - Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ -diff --git a/modules/fdlibm/src/e_log10.cpp b/modules/fdlibm/src/e_log10.cpp ---- a/modules/fdlibm/src/e_log10.cpp -+++ b/modules/fdlibm/src/e_log10.cpp -@@ -19,17 +19,16 @@ - * comments. - * - * log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) - * in not-quite-routine extra precision. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - #include "k_log.h" - - static const double - two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ - ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ - ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ - log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ -diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp ---- a/modules/fdlibm/src/e_log2.cpp -+++ b/modules/fdlibm/src/e_log2.cpp -@@ -21,17 +21,16 @@ - * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel, - * then does the combining and scaling steps - * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k - * in not-quite-routine extra precision. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - #include "k_log.h" - - static const double - two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ - ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ - ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ - -diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp ---- a/modules/fdlibm/src/e_pow.cpp -+++ b/modules/fdlibm/src/e_pow.cpp -@@ -53,17 +53,16 @@ - * Constants : - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - #include <float.h> --#include "math.h" - #include "math_private.h" - - static const double - bp[] = {1.0, 1.5,}, - dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ - dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ - zero = 0.0, - half = 0.5, -diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp ---- a/modules/fdlibm/src/e_sinh.cpp -+++ b/modules/fdlibm/src/e_sinh.cpp -@@ -29,17 +29,16 @@ - * - * Special cases: - * sinh(x) is |x| if x is +INF, -INF, or NaN. - * only sinh(0)=0 is exact for finite x. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double one = 1.0, shuge = 1.0e307; - - double - __ieee754_sinh(double x) - { - double t,h; -diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp ---- a/modules/fdlibm/src/k_exp.cpp -+++ b/modules/fdlibm/src/k_exp.cpp -@@ -26,17 +26,16 @@ - * SUCH DAMAGE. - */ - - #include <sys/cdefs.h> - __FBSDID("$FreeBSD$"); - - #include <complex.h> - --#include "math.h" - #include "math_private.h" - - static const uint32_t k = 1799; /* constant for reduction */ - static const double kln2 = 1246.97177782734161156; /* k * ln2 */ - - /* - * Compute exp(x), scaled to avoid spurious overflow. An exponent is - * returned separately in 'expt'. -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -15,16 +15,18 @@ - */ - - #ifndef _MATH_PRIVATE_H_ - #define _MATH_PRIVATE_H_ - - #include <sys/types.h> - #include <machine/endian.h> - -+#include "fdlibm.h" -+ - /* - * The original fdlibm code used statements like: - * n0 = ((*(int*)&one)>>29)^1; * index of high word * - * ix0 = *(n0+(int*)&x); * high word of x * - * ix1 = *((1-n0)+(int*)&x); * low word of x * - * to dig two 32 bit words out of the 64 bit IEEE floating point - * value. That is non-ANSI, and, moreover, the gcc instruction - * scheduler gets it wrong. We instead use the following macros. -diff --git a/modules/fdlibm/src/s_asinh.cpp b/modules/fdlibm/src/s_asinh.cpp ---- a/modules/fdlibm/src/s_asinh.cpp -+++ b/modules/fdlibm/src/s_asinh.cpp -@@ -21,17 +21,16 @@ - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log(x)+ln2)) for large |x|, else - * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else - * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ - ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ - huge= 1.00000000000000000000e+300; - - double -diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp ---- a/modules/fdlibm/src/s_atan.cpp -+++ b/modules/fdlibm/src/s_atan.cpp -@@ -30,17 +30,16 @@ - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double atanhi[] = { - 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ - 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ - 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ - 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ - }; -diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp ---- a/modules/fdlibm/src/s_cbrt.cpp -+++ b/modules/fdlibm/src/s_cbrt.cpp -@@ -11,17 +11,16 @@ - * - * Optimized by Bruce D. Evans. - */ - - #include <sys/cdefs.h> - __FBSDID("$FreeBSD$"); - - #include <float.h> --#include "math.h" - #include "math_private.h" - - /* cbrt(x) - * Return cube root of x - */ - static const u_int32_t - B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */ - B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */ -diff --git a/modules/fdlibm/src/s_ceil.cpp b/modules/fdlibm/src/s_ceil.cpp ---- a/modules/fdlibm/src/s_ceil.cpp -+++ b/modules/fdlibm/src/s_ceil.cpp -@@ -19,17 +19,16 @@ - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double huge = 1.0e300; - - double - ceil(double x) - { - int32_t i0,i1,j0; -diff --git a/modules/fdlibm/src/s_ceilf.cpp b/modules/fdlibm/src/s_ceilf.cpp ---- a/modules/fdlibm/src/s_ceilf.cpp -+++ b/modules/fdlibm/src/s_ceilf.cpp -@@ -11,17 +11,16 @@ - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - - #include <sys/cdefs.h> - __FBSDID("$FreeBSD$"); - --#include "math.h" - #include "math_private.h" - - static const float huge = 1.0e30; - - float - ceilf(float x) - { - int32_t i0,j0; -diff --git a/modules/fdlibm/src/s_copysign.cpp b/modules/fdlibm/src/s_copysign.cpp ---- a/modules/fdlibm/src/s_copysign.cpp -+++ b/modules/fdlibm/src/s_copysign.cpp -@@ -14,17 +14,16 @@ - __FBSDID("$FreeBSD$"); - - /* - * copysign(double x, double y) - * copysign(x,y) returns a value with the magnitude of x and - * with the sign bit of y. - */ - --#include "math.h" - #include "math_private.h" - - double - copysign(double x, double y) - { - u_int32_t hx,hy; - GET_HIGH_WORD(hx,x); - GET_HIGH_WORD(hy,y); -diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp ---- a/modules/fdlibm/src/s_expm1.cpp -+++ b/modules/fdlibm/src/s_expm1.cpp -@@ -105,17 +105,16 @@ - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - one = 1.0, - tiny = 1.0e-300, - o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */ - ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */ - ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */ -diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp ---- a/modules/fdlibm/src/s_fabs.cpp -+++ b/modules/fdlibm/src/s_fabs.cpp -@@ -12,17 +12,16 @@ - - #include <sys/cdefs.h> - __FBSDID("$FreeBSD$"); - - /* - * fabs(x) returns the absolute value of x. - */ - --#include "math.h" - #include "math_private.h" - - double - fabs(double x) - { - u_int32_t high; - GET_HIGH_WORD(high,x); - SET_HIGH_WORD(x,high&0x7fffffff); -diff --git a/modules/fdlibm/src/s_floor.cpp b/modules/fdlibm/src/s_floor.cpp ---- a/modules/fdlibm/src/s_floor.cpp -+++ b/modules/fdlibm/src/s_floor.cpp -@@ -19,17 +19,16 @@ - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double huge = 1.0e300; - - double - floor(double x) - { - int32_t i0,i1,j0; -diff --git a/modules/fdlibm/src/s_floorf.cpp b/modules/fdlibm/src/s_floorf.cpp ---- a/modules/fdlibm/src/s_floorf.cpp -+++ b/modules/fdlibm/src/s_floorf.cpp -@@ -20,17 +20,16 @@ - * floorf(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floorf(x). - */ - --#include "math.h" - #include "math_private.h" - - static const float huge = 1.0e30; - - float - floorf(float x) - { - int32_t i0,j0; -diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp ---- a/modules/fdlibm/src/s_log1p.cpp -+++ b/modules/fdlibm/src/s_log1p.cpp -@@ -75,17 +75,16 @@ - * if(u==1.0) return x ; else - * return log(u)*(x/(u-1.0)); - * - * See HP-15C Advanced Functions Handbook, p.193. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ - ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ - two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ - Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ - Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ -diff --git a/modules/fdlibm/src/s_nearbyint.cpp b/modules/fdlibm/src/s_nearbyint.cpp ---- a/modules/fdlibm/src/s_nearbyint.cpp -+++ b/modules/fdlibm/src/s_nearbyint.cpp -@@ -25,17 +25,17 @@ - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - - #include <sys/cdefs.h> - __FBSDID("$FreeBSD$"); - - #include <fenv.h> --#include <math.h> -+#include "math_private.h" - - /* - * We save and restore the floating-point environment to avoid raising - * an inexact exception. We can get away with using fesetenv() - * instead of feclearexcept()/feupdateenv() to restore the environment - * because the only exception defined for rint() is overflow, and - * rounding can't overflow as long as emax >= p. - * -diff --git a/modules/fdlibm/src/s_rint.cpp b/modules/fdlibm/src/s_rint.cpp ---- a/modules/fdlibm/src/s_rint.cpp -+++ b/modules/fdlibm/src/s_rint.cpp -@@ -20,17 +20,16 @@ - * Method: - * Using floating addition. - * Exception: - * Inexact flag raised if x not equal to rint(x). - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - TWO52[2]={ - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ - }; - -diff --git a/modules/fdlibm/src/s_rintf.cpp b/modules/fdlibm/src/s_rintf.cpp ---- a/modules/fdlibm/src/s_rintf.cpp -+++ b/modules/fdlibm/src/s_rintf.cpp -@@ -14,17 +14,16 @@ - */ - - #include <sys/cdefs.h> - __FBSDID("$FreeBSD$"); - - #include <float.h> - #include <stdint.h> - --#include "math.h" - #include "math_private.h" - - static const float - TWO23[2]={ - 8.3886080000e+06, /* 0x4b000000 */ - -8.3886080000e+06, /* 0xcb000000 */ - }; - -diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp ---- a/modules/fdlibm/src/s_scalbn.cpp -+++ b/modules/fdlibm/src/s_scalbn.cpp -@@ -17,17 +17,16 @@ - * scalbn (double x, int n) - * scalbn(x,n) returns x* 2**n computed by exponent - * manipulation rather than by actually performing an - * exponentiation or a multiplication. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double - two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ - twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ - huge = 1.0e+300, - tiny = 1.0e-300; - -diff --git a/modules/fdlibm/src/s_tanh.cpp b/modules/fdlibm/src/s_tanh.cpp ---- a/modules/fdlibm/src/s_tanh.cpp -+++ b/modules/fdlibm/src/s_tanh.cpp -@@ -34,17 +34,16 @@ - * - * Special cases: - * tanh(NaN) is NaN; - * only tanh(0)=0 is exact for finite argument. - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const volatile double tiny = 1.0e-300; - static const double one = 1.0, two = 2.0, huge = 1.0e300; - - double - tanh(double x) - { -diff --git a/modules/fdlibm/src/s_trunc.cpp b/modules/fdlibm/src/s_trunc.cpp ---- a/modules/fdlibm/src/s_trunc.cpp -+++ b/modules/fdlibm/src/s_trunc.cpp -@@ -19,17 +19,16 @@ - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to trunc(x). - */ - - #include <float.h> - --#include "math.h" - #include "math_private.h" - - static const double huge = 1.0e300; - - double - trunc(double x) - { - int32_t i0,i1,j0; -diff --git a/modules/fdlibm/src/s_truncf.cpp b/modules/fdlibm/src/s_truncf.cpp ---- a/modules/fdlibm/src/s_truncf.cpp -+++ b/modules/fdlibm/src/s_truncf.cpp -@@ -17,17 +17,16 @@ - * truncf(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to truncf(x). - */ - --#include "math.h" - #include "math_private.h" - - static const float huge = 1.0e30F; - - float - truncf(float x) - { - int32_t i0,j0; diff --git a/modules/fdlibm/patches/05_include_stdint_h_in_math_private_h.patch b/modules/fdlibm/patches/05_include_stdint_h_in_math_private_h.patch deleted file mode 100644 index 46355478f..000000000 --- a/modules/fdlibm/patches/05_include_stdint_h_in_math_private_h.patch +++ /dev/null @@ -1,21 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -12,16 +12,17 @@ - /* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - - #ifndef _MATH_PRIVATE_H_ - #define _MATH_PRIVATE_H_ - -+#include <stdint.h> - #include <sys/types.h> - #include <machine/endian.h> - - #include "fdlibm.h" - - /* - * The original fdlibm code used statements like: - * n0 = ((*(int*)&one)>>29)^1; * index of high word * diff --git a/modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch b/modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch deleted file mode 100644 index bd41b919a..000000000 --- a/modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch +++ /dev/null @@ -1,121 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -14,52 +14,38 @@ - * $FreeBSD$ - */ - - #ifndef _MATH_PRIVATE_H_ - #define _MATH_PRIVATE_H_ - - #include <stdint.h> - #include <sys/types.h> --#include <machine/endian.h> - - #include "fdlibm.h" - -+#include "mozilla/EndianUtils.h" -+ - /* - * The original fdlibm code used statements like: - * n0 = ((*(int*)&one)>>29)^1; * index of high word * - * ix0 = *(n0+(int*)&x); * high word of x * - * ix1 = *((1-n0)+(int*)&x); * low word of x * - * to dig two 32 bit words out of the 64 bit IEEE floating point - * value. That is non-ANSI, and, moreover, the gcc instruction - * scheduler gets it wrong. We instead use the following macros. - * Unlike the original code, we determine the endianness at compile - * time, not at run time; I don't see much benefit to selecting - * endianness at run time. - */ - --/* -- * A union which permits us to convert between a double and two 32 bit -- * ints. -- */ -- --#ifdef __arm__ --#if defined(__VFP_FP__) || defined(__ARM_EABI__) --#define IEEE_WORD_ORDER BYTE_ORDER --#else --#define IEEE_WORD_ORDER BIG_ENDIAN --#endif --#else /* __arm__ */ --#define IEEE_WORD_ORDER BYTE_ORDER --#endif -- - /* A union which permits us to convert between a long double and - four 32 bit ints. */ - --#if IEEE_WORD_ORDER == BIG_ENDIAN -+#if MOZ_BIG_ENDIAN - - typedef union - { - long double value; - struct { - u_int32_t mswhi; - u_int32_t mswlo; - u_int32_t lswhi; -@@ -68,17 +54,17 @@ typedef union - struct { - u_int64_t msw; - u_int64_t lsw; - } parts64; - } ieee_quad_shape_type; - - #endif - --#if IEEE_WORD_ORDER == LITTLE_ENDIAN -+#if MOZ_LITTLE_ENDIAN - - typedef union - { - long double value; - struct { - u_int32_t lswlo; - u_int32_t lswhi; - u_int32_t mswlo; -@@ -87,17 +73,22 @@ typedef union - struct { - u_int64_t lsw; - u_int64_t msw; - } parts64; - } ieee_quad_shape_type; - - #endif - --#if IEEE_WORD_ORDER == BIG_ENDIAN -+/* -+ * A union which permits us to convert between a double and two 32 bit -+ * ints. -+ */ -+ -+#if MOZ_BIG_ENDIAN - - typedef union - { - double value; - struct - { - u_int32_t msw; - u_int32_t lsw; -@@ -105,17 +96,17 @@ typedef union - struct - { - u_int64_t w; - } xparts; - } ieee_double_shape_type; - - #endif - --#if IEEE_WORD_ORDER == LITTLE_ENDIAN -+#if MOZ_LITTLE_ENDIAN - - typedef union - { - double value; - struct - { - u_int32_t lsw; - u_int32_t msw; diff --git a/modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch b/modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch deleted file mode 100644 index 07d3a500a..000000000 --- a/modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch +++ /dev/null @@ -1,54 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -872,16 +872,50 @@ irintl(long double x) - #define __ieee754_j1f j1f - #define __ieee754_y0f y0f - #define __ieee754_y1f y1f - #define __ieee754_jnf jnf - #define __ieee754_ynf ynf - #define __ieee754_remainderf remainderf - #define __ieee754_scalbf scalbf - -+#define acos fdlibm::acos -+#define asin fdlibm::asin -+#define atan fdlibm::atan -+#define atan2 fdlibm::atan2 -+#define cosh fdlibm::cosh -+#define sinh fdlibm::sinh -+#define tanh fdlibm::tanh -+#define exp fdlibm::exp -+#define log fdlibm::log -+#define log10 fdlibm::log10 -+#define pow fdlibm::pow -+#define ceil fdlibm::ceil -+#define ceilf fdlibm::ceilf -+#define fabs fdlibm::fabs -+#define floor fdlibm::floor -+#define acosh fdlibm::acosh -+#define asinh fdlibm::asinh -+#define atanh fdlibm::atanh -+#define cbrt fdlibm::cbrt -+#define expm1 fdlibm::expm1 -+#define hypot fdlibm::hypot -+#define log1p fdlibm::log1p -+#define log2 fdlibm::log2 -+#define scalb fdlibm::scalb -+#define copysign fdlibm::copysign -+#define scalbn fdlibm::scalbn -+#define trunc fdlibm::trunc -+#define truncf fdlibm::truncf -+#define floorf fdlibm::floorf -+#define nearbyint fdlibm::nearbyint -+#define nearbyintf fdlibm::nearbyintf -+#define rint fdlibm::rint -+#define rintf fdlibm::rintf -+ - /* fdlibm kernel function */ - int __kernel_rem_pio2(double*,double*,int,int,int); - - /* double precision kernel functions */ - #ifndef INLINE_REM_PIO2 - int __ieee754_rem_pio2(double,double*); - #endif - double __kernel_sin(double,double,int); diff --git a/modules/fdlibm/patches/08_remove_weak_reference_macro.patch b/modules/fdlibm/patches/08_remove_weak_reference_macro.patch deleted file mode 100644 index 0b55bbd84..000000000 --- a/modules/fdlibm/patches/08_remove_weak_reference_macro.patch +++ /dev/null @@ -1,385 +0,0 @@ -diff --git a/modules/fdlibm/src/e_acos.cpp b/modules/fdlibm/src/e_acos.cpp ---- a/modules/fdlibm/src/e_acos.cpp -+++ b/modules/fdlibm/src/e_acos.cpp -@@ -99,12 +99,8 @@ double - c = (z-df*df)/(s+df); - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - r = p/q; - w = r*s+c; - return 2.0*(df+w); - } - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(acos, acosl); --#endif -diff --git a/modules/fdlibm/src/e_acosh.cpp b/modules/fdlibm/src/e_acosh.cpp ---- a/modules/fdlibm/src/e_acosh.cpp -+++ b/modules/fdlibm/src/e_acosh.cpp -@@ -56,12 +56,8 @@ double - } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ - t=x*x; - return __ieee754_log(2.0*x-one/(x+sqrt(t-one))); - } else { /* 1<x<2 */ - t = x-one; - return log1p(t+sqrt(2.0*t+t*t)); - } - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(acosh, acoshl); --#endif -diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp ---- a/modules/fdlibm/src/e_asin.cpp -+++ b/modules/fdlibm/src/e_asin.cpp -@@ -105,12 +105,8 @@ double - c = (t-w*w)/(s+w); - r = p/q; - p = 2.0*s*r-(pio2_lo-2.0*c); - q = pio4_hi-2.0*w; - t = pio4_hi-(p-q); - } - if(hx>0) return t; else return -t; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(asin, asinl); --#endif -diff --git a/modules/fdlibm/src/e_atan2.cpp b/modules/fdlibm/src/e_atan2.cpp ---- a/modules/fdlibm/src/e_atan2.cpp -+++ b/modules/fdlibm/src/e_atan2.cpp -@@ -117,12 +117,8 @@ double - switch (m) { - case 0: return z ; /* atan(+,+) */ - case 1: return -z ; /* atan(-,+) */ - case 2: return pi-(z-pi_lo);/* atan(+,-) */ - default: /* case 3 */ - return (z-pi_lo)-pi;/* atan(-,-) */ - } - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(atan2, atan2l); --#endif -diff --git a/modules/fdlibm/src/e_atanh.cpp b/modules/fdlibm/src/e_atanh.cpp ---- a/modules/fdlibm/src/e_atanh.cpp -+++ b/modules/fdlibm/src/e_atanh.cpp -@@ -56,12 +56,8 @@ double - SET_HIGH_WORD(x,ix); - if(ix<0x3fe00000) { /* x < 0.5 */ - t = x+x; - t = 0.5*log1p(t+t*x/(one-x)); - } else - t = 0.5*log1p((x+x)/(one-x)); - if(hx>=0) return t; else return -t; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(atanh, atanhl); --#endif -diff --git a/modules/fdlibm/src/e_cosh.cpp b/modules/fdlibm/src/e_cosh.cpp ---- a/modules/fdlibm/src/e_cosh.cpp -+++ b/modules/fdlibm/src/e_cosh.cpp -@@ -73,12 +73,8 @@ double - - /* |x| in [log(maxdouble), overflowthresold] */ - if (ix<=0x408633CE) - return __ldexp_exp(fabs(x), -1); - - /* |x| > overflowthresold, cosh(x) overflow */ - return huge*huge; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(cosh, coshl); --#endif -diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp ---- a/modules/fdlibm/src/e_exp.cpp -+++ b/modules/fdlibm/src/e_exp.cpp -@@ -152,12 +152,8 @@ double - else y = one-((lo-(x*c)/(2.0-c))-hi); - if(k >= -1021) { - if (k==1024) return y*2.0*0x1p1023; - return y*twopk; - } else { - return y*twopk*twom1000; - } - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(exp, expl); --#endif -diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp ---- a/modules/fdlibm/src/e_hypot.cpp -+++ b/modules/fdlibm/src/e_hypot.cpp -@@ -119,12 +119,8 @@ double - if(k!=0) { - u_int32_t high; - t1 = 1.0; - GET_HIGH_WORD(high,t1); - SET_HIGH_WORD(t1,high+(k<<20)); - return t1*w; - } else return w; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(hypot, hypotl); --#endif -diff --git a/modules/fdlibm/src/e_log.cpp b/modules/fdlibm/src/e_log.cpp ---- a/modules/fdlibm/src/e_log.cpp -+++ b/modules/fdlibm/src/e_log.cpp -@@ -135,12 +135,8 @@ double - hfsq=0.5*f*f; - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); - } else { - if(k==0) return f-s*(f-R); else - return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); - } - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(log, logl); --#endif -diff --git a/modules/fdlibm/src/e_log10.cpp b/modules/fdlibm/src/e_log10.cpp ---- a/modules/fdlibm/src/e_log10.cpp -+++ b/modules/fdlibm/src/e_log10.cpp -@@ -82,12 +82,8 @@ double - * with some parallelism and it reduces the error for many args. - */ - w = y2 + val_hi; - val_lo += (y2 - w) + val_hi; - val_hi = w; - - return val_lo + val_hi; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(log10, log10l); --#endif -diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp ---- a/modules/fdlibm/src/e_log2.cpp -+++ b/modules/fdlibm/src/e_log2.cpp -@@ -105,12 +105,8 @@ double - - /* spadd(val_hi, val_lo, y), except for not using double_t: */ - w = y + val_hi; - val_lo += (y - w) + val_hi; - val_hi = w; - - return val_lo + val_hi; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(log2, log2l); --#endif -diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp ---- a/modules/fdlibm/src/e_pow.cpp -+++ b/modules/fdlibm/src/e_pow.cpp -@@ -302,12 +302,8 @@ double - r = (z*t1)/(t1-two)-(w+z*w); - z = one-(r-z); - GET_HIGH_WORD(j,z); - j += (n<<20); - if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */ - else SET_HIGH_WORD(z,j); - return s*z; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(pow, powl); --#endif -diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp ---- a/modules/fdlibm/src/e_sinh.cpp -+++ b/modules/fdlibm/src/e_sinh.cpp -@@ -67,12 +67,8 @@ double - - /* |x| in [log(maxdouble), overflowthresold] */ - if (ix<=0x408633CE) - return h*2.0*__ldexp_exp(fabs(x), -1); - - /* |x| > overflowthresold, sinh(x) overflow */ - return x*shuge; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(sinh, sinhl); --#endif -diff --git a/modules/fdlibm/src/s_asinh.cpp b/modules/fdlibm/src/s_asinh.cpp ---- a/modules/fdlibm/src/s_asinh.cpp -+++ b/modules/fdlibm/src/s_asinh.cpp -@@ -50,12 +50,8 @@ asinh(double x) - t = fabs(x); - w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; - w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t))); - } - if(hx>0) return w; else return -w; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(asinh, asinhl); --#endif -diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp ---- a/modules/fdlibm/src/s_atan.cpp -+++ b/modules/fdlibm/src/s_atan.cpp -@@ -112,12 +112,8 @@ atan(double x) - s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); - s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); - if (id<0) return x - x*(s1+s2); - else { - z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); - return (hx<0)? -z:z; - } - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(atan, atanl); --#endif -diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp ---- a/modules/fdlibm/src/s_cbrt.cpp -+++ b/modules/fdlibm/src/s_cbrt.cpp -@@ -106,12 +106,8 @@ cbrt(double x) - s=t*t; /* t*t is exact */ - r=x/s; /* error <= 0.5 ulps; |r| < |t| */ - w=t+t; /* t+t is exact */ - r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */ - t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ - - return(t); - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(cbrt, cbrtl); --#endif -diff --git a/modules/fdlibm/src/s_ceil.cpp b/modules/fdlibm/src/s_ceil.cpp ---- a/modules/fdlibm/src/s_ceil.cpp -+++ b/modules/fdlibm/src/s_ceil.cpp -@@ -65,12 +65,8 @@ ceil(double x) - } - } - i1 &= (~i); - } - } - INSERT_WORDS(x,i0,i1); - return x; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(ceil, ceill); --#endif -diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp ---- a/modules/fdlibm/src/s_expm1.cpp -+++ b/modules/fdlibm/src/s_expm1.cpp -@@ -210,12 +210,8 @@ expm1(double x) - SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */ - y = x-(e+t); - y += one; - y = y*twopk; - } - } - return y; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(expm1, expm1l); --#endif -diff --git a/modules/fdlibm/src/s_floor.cpp b/modules/fdlibm/src/s_floor.cpp ---- a/modules/fdlibm/src/s_floor.cpp -+++ b/modules/fdlibm/src/s_floor.cpp -@@ -66,12 +66,8 @@ floor(double x) - } - } - i1 &= (~i); - } - } - INSERT_WORDS(x,i0,i1); - return x; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(floor, floorl); --#endif -diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp ---- a/modules/fdlibm/src/s_log1p.cpp -+++ b/modules/fdlibm/src/s_log1p.cpp -@@ -168,12 +168,8 @@ log1p(double x) - return k*ln2_hi-((R-(k*ln2_lo+c))-f); - } - s = f/(2.0+f); - z = s*s; - R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f); - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(log1p, log1pl); --#endif -diff --git a/modules/fdlibm/src/s_rint.cpp b/modules/fdlibm/src/s_rint.cpp ---- a/modules/fdlibm/src/s_rint.cpp -+++ b/modules/fdlibm/src/s_rint.cpp -@@ -80,12 +80,8 @@ rint(double x) - if((i1&i)==0) return x; /* x is integral */ - i>>=1; - if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); - } - INSERT_WORDS(x,i0,i1); - STRICT_ASSIGN(double,w,TWO52[sx]+x); - return w-TWO52[sx]; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(rint, rintl); --#endif -diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp ---- a/modules/fdlibm/src/s_scalbn.cpp -+++ b/modules/fdlibm/src/s_scalbn.cpp -@@ -53,13 +53,8 @@ scalbn (double x, int n) - return huge*copysign(huge,x); /*overflow*/ - else - return tiny*copysign(tiny,x); /*underflow*/ - } - k += 54; /* subnormal result */ - SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); - return x*twom54; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(scalbn, ldexpl); --__weak_reference(scalbn, scalbnl); --#endif -diff --git a/modules/fdlibm/src/s_tanh.cpp b/modules/fdlibm/src/s_tanh.cpp ---- a/modules/fdlibm/src/s_tanh.cpp -+++ b/modules/fdlibm/src/s_tanh.cpp -@@ -72,12 +72,8 @@ tanh(double x) - z= -t/(t+two); - } - /* |x| >= 22, return +-1 */ - } else { - z = one - tiny; /* raise inexact flag */ - } - return (jx>=0)? z: -z; - } -- --#if (LDBL_MANT_DIG == 53) --__weak_reference(tanh, tanhl); --#endif -diff --git a/modules/fdlibm/src/s_trunc.cpp b/modules/fdlibm/src/s_trunc.cpp ---- a/modules/fdlibm/src/s_trunc.cpp -+++ b/modules/fdlibm/src/s_trunc.cpp -@@ -55,12 +55,8 @@ trunc(double x) - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) /* raise inexact flag */ - i1 &= (~i); - } - INSERT_WORDS(x,i0,i1); - return x; - } -- --#if LDBL_MANT_DIG == 53 --__weak_reference(trunc, truncl); --#endif diff --git a/modules/fdlibm/patches/09_comment_out_rcsid_variable.patch b/modules/fdlibm/patches/09_comment_out_rcsid_variable.patch deleted file mode 100644 index d520d9257..000000000 --- a/modules/fdlibm/patches/09_comment_out_rcsid_variable.patch +++ /dev/null @@ -1,812 +0,0 @@ -diff --git a/modules/fdlibm/src/e_acos.cpp b/modules/fdlibm/src/e_acos.cpp ---- a/modules/fdlibm/src/e_acos.cpp -+++ b/modules/fdlibm/src/e_acos.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_acos(x) - * Method : - * acos(x) = pi/2 - asin(x) - * acos(-x) = pi/2 + asin(x) - * For |x|<=0.5 - * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) - * For x>0.5 -diff --git a/modules/fdlibm/src/e_acosh.cpp b/modules/fdlibm/src/e_acosh.cpp ---- a/modules/fdlibm/src/e_acosh.cpp -+++ b/modules/fdlibm/src/e_acosh.cpp -@@ -7,18 +7,18 @@ - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_acosh(x) - * Method : - * Based on - * acosh(x) = log [ x + sqrt(x*x-1) ] - * we have - * acosh(x) := log(x)+ln2, if x is large; else - * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else -diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp ---- a/modules/fdlibm/src/e_asin.cpp -+++ b/modules/fdlibm/src/e_asin.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_asin(x) - * Method : - * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... - * we approximate asin(x) on [0,0.5] by - * asin(x) = x + x*x^2*R(x^2) - * where - * R(x^2) is a rational approximation of (asin(x)-x)/x^3 -diff --git a/modules/fdlibm/src/e_atan2.cpp b/modules/fdlibm/src/e_atan2.cpp ---- a/modules/fdlibm/src/e_atan2.cpp -+++ b/modules/fdlibm/src/e_atan2.cpp -@@ -7,18 +7,18 @@ - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_atan2(y,x) - * Method : - * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): - * ARG (x+iy) = arctan(y/x) ... if x > 0, - * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, - * -diff --git a/modules/fdlibm/src/e_atanh.cpp b/modules/fdlibm/src/e_atanh.cpp ---- a/modules/fdlibm/src/e_atanh.cpp -+++ b/modules/fdlibm/src/e_atanh.cpp -@@ -7,18 +7,18 @@ - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_atanh(x) - * Method : - * 1.Reduced x to positive by atanh(-x) = -atanh(x) - * 2.For x>=0.5 - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x -diff --git a/modules/fdlibm/src/e_cosh.cpp b/modules/fdlibm/src/e_cosh.cpp ---- a/modules/fdlibm/src/e_cosh.cpp -+++ b/modules/fdlibm/src/e_cosh.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_cosh(x) - * Method : - * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 - * 1. Replace x by |x| (cosh(x) = cosh(-x)). - * 2. - * [ exp(x) - 1 ]^2 - * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- -diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp ---- a/modules/fdlibm/src/e_exp.cpp -+++ b/modules/fdlibm/src/e_exp.cpp -@@ -5,18 +5,18 @@ - * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. - * - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - --#include <sys/cdefs.h> --__FBSDID("$FreeBSD$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_exp(x) - * Returns the exponential of x. - * - * Method - * 1. Argument reduction: - * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. - * Given x, find r and integer k such that -diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp ---- a/modules/fdlibm/src/e_hypot.cpp -+++ b/modules/fdlibm/src/e_hypot.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_hypot(x,y) - * - * Method : - * If (assume round-to-nearest) z=x*x+y*y - * has error less than sqrt(2)/2 ulp, than - * sqrt(z) has error less than 1 ulp (exercise). - * -diff --git a/modules/fdlibm/src/e_log.cpp b/modules/fdlibm/src/e_log.cpp ---- a/modules/fdlibm/src/e_log.cpp -+++ b/modules/fdlibm/src/e_log.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_log(x) - * Return the logrithm of x - * - * Method : - * 1. Argument Reduction: find k and f such that - * x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . -diff --git a/modules/fdlibm/src/e_log10.cpp b/modules/fdlibm/src/e_log10.cpp ---- a/modules/fdlibm/src/e_log10.cpp -+++ b/modules/fdlibm/src/e_log10.cpp -@@ -6,32 +6,32 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * Return the base 10 logarithm of x. See e_log.c and k_log.h for most - * comments. - * - * log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) - * in not-quite-routine extra precision. - */ - - #include <float.h> - - #include "math_private.h" - #include "k_log.h" -- -+ - static const double - two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ - ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ - ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ - log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ - log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ - - static const double zero = 0.0; -diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp ---- a/modules/fdlibm/src/e_log2.cpp -+++ b/modules/fdlibm/src/e_log2.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * Return the base 2 logarithm of x. See e_log.c and k_log.h for most - * comments. - * - * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel, - * then does the combining and scaling steps - * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k -diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp ---- a/modules/fdlibm/src/e_pow.cpp -+++ b/modules/fdlibm/src/e_pow.cpp -@@ -4,18 +4,18 @@ - * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. - * - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - --#include <sys/cdefs.h> --__FBSDID("$FreeBSD$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_pow(x,y) return x**y - * - * n - * Method: Let x = 2 * (1+f) - * 1. Compute and return log2(x) in two pieces: - * log2(x) = w1 + w2, - * where w1 has 53-24 = 29 bit trailing zeros. -diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp ---- a/modules/fdlibm/src/e_sinh.cpp -+++ b/modules/fdlibm/src/e_sinh.cpp -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* __ieee754_sinh(x) - * Method : - * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 - * 1. Replace x by |x| (sinh(-x) = -sinh(x)). - * 2. - * E + E/(E+1) - * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) -diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp ---- a/modules/fdlibm/src/k_exp.cpp -+++ b/modules/fdlibm/src/k_exp.cpp -@@ -21,22 +21,22 @@ - * 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 <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - #include <complex.h> - --#include "math_private.h" -+ #include "math_private.h" - - static const uint32_t k = 1799; /* constant for reduction */ - static const double kln2 = 1246.97177782734161156; /* k * ln2 */ - - /* - * Compute exp(x), scaled to avoid spurious overflow. An exponent is - * returned separately in 'expt'. - * -diff --git a/modules/fdlibm/src/k_log.h b/modules/fdlibm/src/k_log.h ---- a/modules/fdlibm/src/k_log.h -+++ b/modules/fdlibm/src/k_log.h -@@ -6,18 +6,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * k_log1p(f): - * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)]. - * - * The following describes the overall strategy for computing - * logarithms in base e. The argument reduction and adding the final - * term of the polynomial are done by the caller for increased accuracy -diff --git a/modules/fdlibm/src/s_asinh.cpp b/modules/fdlibm/src/s_asinh.cpp ---- a/modules/fdlibm/src/s_asinh.cpp -+++ b/modules/fdlibm/src/s_asinh.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* asinh(x) - * Method : - * Based on - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log(x)+ln2)) for large |x|, else -diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp ---- a/modules/fdlibm/src/s_atan.cpp -+++ b/modules/fdlibm/src/s_atan.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* atan(x) - * Method - * 1. Reduce x to positive by atan(x) = -atan(-x). - * 2. According to the integer k=4t+0.25 chopped, t=x, the argument - * is further reduced to one of the following intervals and the - * arctangent of t is evaluated by the corresponding formula: - * -diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp ---- a/modules/fdlibm/src/s_cbrt.cpp -+++ b/modules/fdlibm/src/s_cbrt.cpp -@@ -7,18 +7,18 @@ - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - #include <float.h> - #include "math_private.h" - - /* cbrt(x) - * Return cube root of x - */ - static const u_int32_t -diff --git a/modules/fdlibm/src/s_ceil.cpp b/modules/fdlibm/src/s_ceil.cpp ---- a/modules/fdlibm/src/s_ceil.cpp -+++ b/modules/fdlibm/src/s_ceil.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * ceil(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). -diff --git a/modules/fdlibm/src/s_ceilf.cpp b/modules/fdlibm/src/s_ceilf.cpp ---- a/modules/fdlibm/src/s_ceilf.cpp -+++ b/modules/fdlibm/src/s_ceilf.cpp -@@ -8,18 +8,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - #include "math_private.h" - - static const float huge = 1.0e30; - - float - ceilf(float x) - { -diff --git a/modules/fdlibm/src/s_copysign.cpp b/modules/fdlibm/src/s_copysign.cpp ---- a/modules/fdlibm/src/s_copysign.cpp -+++ b/modules/fdlibm/src/s_copysign.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * copysign(double x, double y) - * copysign(x,y) returns a value with the magnitude of x and - * with the sign bit of y. - */ - - #include "math_private.h" -diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp ---- a/modules/fdlibm/src/s_expm1.cpp -+++ b/modules/fdlibm/src/s_expm1.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* expm1(x) - * Returns exp(x)-1, the exponential of x minus 1. - * - * Method - * 1. Argument reduction: - * Given x, find r and integer k such that - * -diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp ---- a/modules/fdlibm/src/s_fabs.cpp -+++ b/modules/fdlibm/src/s_fabs.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * fabs(x) returns the absolute value of x. - */ - - #include "math_private.h" - - double -diff --git a/modules/fdlibm/src/s_floor.cpp b/modules/fdlibm/src/s_floor.cpp ---- a/modules/fdlibm/src/s_floor.cpp -+++ b/modules/fdlibm/src/s_floor.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * floor(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). -diff --git a/modules/fdlibm/src/s_floorf.cpp b/modules/fdlibm/src/s_floorf.cpp ---- a/modules/fdlibm/src/s_floorf.cpp -+++ b/modules/fdlibm/src/s_floorf.cpp -@@ -8,18 +8,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * floorf(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floorf(x). -diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp ---- a/modules/fdlibm/src/s_log1p.cpp -+++ b/modules/fdlibm/src/s_log1p.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* double log1p(double x) - * - * Method : - * 1. Argument Reduction: find k and f such that - * 1+x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * -diff --git a/modules/fdlibm/src/s_nearbyint.cpp b/modules/fdlibm/src/s_nearbyint.cpp ---- a/modules/fdlibm/src/s_nearbyint.cpp -+++ b/modules/fdlibm/src/s_nearbyint.cpp -@@ -21,18 +21,18 @@ - * 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 <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - #include <fenv.h> - #include "math_private.h" - - /* - * We save and restore the floating-point environment to avoid raising - * an inexact exception. We can get away with using fesetenv() - * instead of feclearexcept()/feupdateenv() to restore the environment -diff --git a/modules/fdlibm/src/s_rint.cpp b/modules/fdlibm/src/s_rint.cpp ---- a/modules/fdlibm/src/s_rint.cpp -+++ b/modules/fdlibm/src/s_rint.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * rint(x) - * Return x rounded to integral value according to the prevailing - * rounding mode. - * Method: - * Using floating addition. - * Exception: -diff --git a/modules/fdlibm/src/s_rintf.cpp b/modules/fdlibm/src/s_rintf.cpp ---- a/modules/fdlibm/src/s_rintf.cpp -+++ b/modules/fdlibm/src/s_rintf.cpp -@@ -8,18 +8,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - #include <float.h> - #include <stdint.h> - - #include "math_private.h" - - static const float - TWO23[2]={ -diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp ---- a/modules/fdlibm/src/s_scalbn.cpp -+++ b/modules/fdlibm/src/s_scalbn.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * scalbn (double x, int n) - * scalbn(x,n) returns x* 2**n computed by exponent - * manipulation rather than by actually performing an - * exponentiation or a multiplication. - */ - -diff --git a/modules/fdlibm/src/s_tanh.cpp b/modules/fdlibm/src/s_tanh.cpp ---- a/modules/fdlibm/src/s_tanh.cpp -+++ b/modules/fdlibm/src/s_tanh.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* Tanh(x) - * Return the Hyperbolic Tangent of x - * - * Method : - * x -x - * e - e - * 0. tanh(x) is defined to be ----------- -diff --git a/modules/fdlibm/src/s_trunc.cpp b/modules/fdlibm/src/s_trunc.cpp ---- a/modules/fdlibm/src/s_trunc.cpp -+++ b/modules/fdlibm/src/s_trunc.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * trunc(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to trunc(x). -diff --git a/modules/fdlibm/src/s_truncf.cpp b/modules/fdlibm/src/s_truncf.cpp ---- a/modules/fdlibm/src/s_truncf.cpp -+++ b/modules/fdlibm/src/s_truncf.cpp -@@ -5,18 +5,18 @@ - * - * 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$"); -+//#include <sys/cdefs.h> -+//__FBSDID("$FreeBSD$"); - - /* - * truncf(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to truncf(x). diff --git a/modules/fdlibm/patches/10_remove_unused_function_from_k_exp_cpp.patch b/modules/fdlibm/patches/10_remove_unused_function_from_k_exp_cpp.patch deleted file mode 100644 index 36aee9bb6..000000000 --- a/modules/fdlibm/patches/10_remove_unused_function_from_k_exp_cpp.patch +++ /dev/null @@ -1,55 +0,0 @@ -diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp ---- a/modules/fdlibm/src/k_exp.cpp -+++ b/modules/fdlibm/src/k_exp.cpp -@@ -24,18 +24,16 @@ - * 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 <complex.h> -- - #include "math_private.h" - - static const uint32_t k = 1799; /* constant for reduction */ - static const double kln2 = 1246.97177782734161156; /* k * ln2 */ - - /* - * Compute exp(x), scaled to avoid spurious overflow. An exponent is - * returned separately in 'expt'. -@@ -78,32 +76,8 @@ double - double exp_x, scale; - int ex_expt; - - exp_x = __frexp_exp(x, &ex_expt); - expt += ex_expt; - INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); - return (exp_x * scale); - } -- --double complex --__ldexp_cexp(double complex z, int expt) --{ -- double x, y, exp_x, scale1, scale2; -- int ex_expt, half_expt; -- -- x = creal(z); -- y = cimag(z); -- exp_x = __frexp_exp(x, &ex_expt); -- expt += ex_expt; -- -- /* -- * Arrange so that scale1 * scale2 == 2**expt. We use this to -- * compensate for scalbn being horrendously slow. -- */ -- half_expt = expt / 2; -- INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); -- half_expt = expt - half_expt; -- INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); -- -- return (CMPLX(cos(y) * exp_x * scale1 * scale2, -- sin(y) * exp_x * scale1 * scale2)); --} diff --git a/modules/fdlibm/patches/11_include_cfloat_to_use_flt_eval_method.patch b/modules/fdlibm/patches/11_include_cfloat_to_use_flt_eval_method.patch deleted file mode 100644 index a5ddcb1ca..000000000 --- a/modules/fdlibm/patches/11_include_cfloat_to_use_flt_eval_method.patch +++ /dev/null @@ -1,21 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -12,16 +12,17 @@ - /* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - - #ifndef _MATH_PRIVATE_H_ - #define _MATH_PRIVATE_H_ - -+#include <cfloat> - #include <stdint.h> - #include <sys/types.h> - - #include "fdlibm.h" - - #include "mozilla/EndianUtils.h" - - /* diff --git a/modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch b/modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch deleted file mode 100644 index 106d64dbc..000000000 --- a/modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch +++ /dev/null @@ -1,27 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -33,16 +33,23 @@ - * to dig two 32 bit words out of the 64 bit IEEE floating point - * value. That is non-ANSI, and, moreover, the gcc instruction - * scheduler gets it wrong. We instead use the following macros. - * Unlike the original code, we determine the endianness at compile - * time, not at run time; I don't see much benefit to selecting - * endianness at run time. - */ - -+#ifndef u_int32_t -+#define u_int32_t uint32_t -+#endif -+#ifndef u_int64_t -+#define u_int64_t uint64_t -+#endif -+ - /* A union which permits us to convert between a long double and - four 32 bit ints. */ - - #if MOZ_BIG_ENDIAN - - typedef union - { - long double value; diff --git a/modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch b/modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch deleted file mode 100644 index a6117e24b..000000000 --- a/modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch +++ /dev/null @@ -1,31 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -328,16 +328,27 @@ do { \ - if (sizeof(type) >= sizeof(long double)) \ - (lval) = (rval); \ - else { \ - __lval = (rval); \ - (lval) = __lval; \ - } \ - } while (0) - #endif -+#else -+#define STRICT_ASSIGN(type, lval, rval) do { \ -+ volatile type __lval; \ -+ \ -+ if (sizeof(type) >= sizeof(long double)) \ -+ (lval) = (rval); \ -+ else { \ -+ __lval = (rval); \ -+ (lval) = __lval; \ -+ } \ -+} while (0) - #endif /* FLT_EVAL_METHOD */ - - /* Support switching the mode to FP_PE if necessary. */ - #if defined(__i386__) && !defined(NO_FPSETPREC) - #define ENTERI() ENTERIT(long double) - #define ENTERIT(returntype) \ - returntype __retval; \ - fp_prec_t __oprec; \ diff --git a/modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch b/modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch deleted file mode 100644 index 85f511882..000000000 --- a/modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch +++ /dev/null @@ -1,47 +0,0 @@ -diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp ---- a/modules/fdlibm/src/e_exp.cpp -+++ b/modules/fdlibm/src/e_exp.cpp -@@ -146,14 +146,17 @@ double - if(k >= -1021) - INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0); - else - INSERT_WORDS(twopk,((u_int32_t)(0x3ff+(k+1000)))<<20, 0); - c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - if(k==0) return one-((x*c)/(c-2.0)-x); - else y = one-((lo-(x*c)/(2.0-c))-hi); - if(k >= -1021) { -- if (k==1024) return y*2.0*0x1p1023; -+ if (k==1024) { -+ double const_0x1p1023 = pow(2, 1023); -+ return y*2.0*const_0x1p1023; -+ } - return y*twopk; - } else { - return y*twopk*twom1000; - } - } -diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp ---- a/modules/fdlibm/src/s_expm1.cpp -+++ b/modules/fdlibm/src/s_expm1.cpp -@@ -192,17 +192,20 @@ expm1(double x) - e -= hxs; - if(k== -1) return 0.5*(x-e)-0.5; - if(k==1) { - if(x < -0.25) return -2.0*(e-(x+0.5)); - else return one+2.0*(x-e); - } - if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - y = one-(e-x); -- if (k == 1024) y = y*2.0*0x1p1023; -+ if (k == 1024) { -+ double const_0x1p1023 = pow(2, 1023); -+ y = y*2.0*const_0x1p1023; -+ } - else y = y*twopk; - return y-one; - } - t = one; - if(k<20) { - SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ - y = t-(e-x); - y = y*twopk; diff --git a/modules/fdlibm/patches/15_remove_unused_rintl_function_from_s_nearbyint_cpp.patch b/modules/fdlibm/patches/15_remove_unused_rintl_function_from_s_nearbyint_cpp.patch deleted file mode 100644 index b64e281aa..000000000 --- a/modules/fdlibm/patches/15_remove_unused_rintl_function_from_s_nearbyint_cpp.patch +++ /dev/null @@ -1,13 +0,0 @@ -diff --git a/modules/fdlibm/src/s_nearbyint.cpp b/modules/fdlibm/src/s_nearbyint.cpp ---- a/modules/fdlibm/src/s_nearbyint.cpp -+++ b/modules/fdlibm/src/s_nearbyint.cpp -@@ -53,9 +53,8 @@ fn(type x) \ - fegetenv(&env); \ - ret = rint(x); \ - fesetenv(&env); \ - return (ret); \ - } - - DECL(double, nearbyint, rint) - DECL(float, nearbyintf, rintf) --DECL(long double, nearbyintl, rintl) diff --git a/modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch b/modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch deleted file mode 100644 index 855f91a0d..000000000 --- a/modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch +++ /dev/null @@ -1,22 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -314,17 +314,17 @@ do { \ - /* The above works on non-i386 too, but we use this to check v. */ - #define LD80C(m, ex, v) { .e = (v), } - #endif - - #ifdef FLT_EVAL_METHOD - /* - * Attempt to get strict C99 semantics for assignment with non-C99 compilers. - */ --#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 -+#if !defined(_MSC_VER) && (FLT_EVAL_METHOD == 0 || __GNUC__ == 0) - #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) - #else - #define STRICT_ASSIGN(type, lval, rval) do { \ - volatile type __lval; \ - \ - if (sizeof(type) >= sizeof(long double)) \ - (lval) = (rval); \ - else { \ diff --git a/modules/fdlibm/patches/17_exp_exact_result_for_positive_one.patch b/modules/fdlibm/patches/17_exp_exact_result_for_positive_one.patch deleted file mode 100644 index 9cda2beef..000000000 --- a/modules/fdlibm/patches/17_exp_exact_result_for_positive_one.patch +++ /dev/null @@ -1,40 +0,0 @@ -diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp ---- a/modules/fdlibm/src/e_exp.cpp -+++ b/modules/fdlibm/src/e_exp.cpp -@@ -91,16 +91,18 @@ ln2LO[2] ={ 1.90821492927058770002e-10 - -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ - invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ - P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ - P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ - P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ - P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ - P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ - -+static const double E = 2.7182818284590452354; /* e */ -+ - static volatile double - huge = 1.0e+300, - twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/ - - double - __ieee754_exp(double x) /* default IEEE double exp */ - { - double y,hi=0.0,lo=0.0,c,t,twopk; -@@ -122,16 +124,17 @@ double - } - if(x > o_threshold) return huge*huge; /* overflow */ - if(x < u_threshold) return twom1000*twom1000; /* underflow */ - } - - /* argument reduction */ - if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ - if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ -+ if (x == 1.0) return E; - hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; - } else { - k = (int)(invln2*x+halF[xsb]); - t = k; - hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ - lo = t*ln2LO[0]; - } - STRICT_ASSIGN(double, x, hi - lo); diff --git a/modules/fdlibm/patches/18_use_stdlib_sqrt.patch b/modules/fdlibm/patches/18_use_stdlib_sqrt.patch deleted file mode 100644 index 1f87dd73b..000000000 --- a/modules/fdlibm/patches/18_use_stdlib_sqrt.patch +++ /dev/null @@ -1,255 +0,0 @@ -diff --git a/modules/fdlibm/src/e_acos.cpp b/modules/fdlibm/src/e_acos.cpp ---- a/modules/fdlibm/src/e_acos.cpp -+++ b/modules/fdlibm/src/e_acos.cpp -@@ -33,16 +33,17 @@ - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN with invalid signal. - * - * Function needed: sqrt - */ - -+#include <cmath> - #include <float.h> - - #include "math_private.h" - - static const double - one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ - pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ - pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ -@@ -82,23 +83,23 @@ double - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - r = p/q; - return pio2_hi - (x - (pio2_lo-x*r)); - } else if (hx<0) { /* x < -0.5 */ - z = (one+x)*0.5; - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); -- s = sqrt(z); -+ s = std::sqrt(z); - r = p/q; - w = r*s-pio2_lo; - return pi - 2.0*(s+w); - } else { /* x > 0.5 */ - z = (one-x)*0.5; -- s = sqrt(z); -+ s = std::sqrt(z); - df = s; - SET_LOW_WORD(df,0); - c = (z-df*df)/(s+df); - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - r = p/q; - w = r*s+c; - return 2.0*(df+w); -diff --git a/modules/fdlibm/src/e_acosh.cpp b/modules/fdlibm/src/e_acosh.cpp ---- a/modules/fdlibm/src/e_acosh.cpp -+++ b/modules/fdlibm/src/e_acosh.cpp -@@ -24,16 +24,17 @@ - * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else - * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - */ - -+#include <cmath> - #include <float.h> - - #include "math_private.h" - - static const double - one = 1.0, - ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ - -@@ -50,14 +51,14 @@ double - if(hx >=0x7ff00000) { /* x is inf of NaN */ - return x+x; - } else - return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */ - } else if(((hx-0x3ff00000)|lx)==0) { - return 0.0; /* acosh(1) = 0 */ - } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ - t=x*x; -- return __ieee754_log(2.0*x-one/(x+sqrt(t-one))); -+ return __ieee754_log(2.0*x-one/(x+std::sqrt(t-one))); - } else { /* 1<x<2 */ - t = x-one; -- return log1p(t+sqrt(2.0*t+t*t)); -+ return log1p(t+std::sqrt(2.0*t+t*t)); - } - } -diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp ---- a/modules/fdlibm/src/e_asin.cpp -+++ b/modules/fdlibm/src/e_asin.cpp -@@ -39,16 +39,17 @@ - * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN with invalid signal. - * - */ - -+#include <cmath> - #include <float.h> - - #include "math_private.h" - - static const double - one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ - huge = 1.000e+300, - pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ -@@ -90,17 +91,17 @@ double - w = p/q; - return x+x*w; - } - /* 1> |x|>= 0.5 */ - w = one-fabs(x); - t = w*0.5; - p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); -- s = sqrt(t); -+ s = std::sqrt(t); - if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ - w = p/q; - t = pio2_hi-(2.0*(s+s*w)-pio2_lo); - } else { - w = s; - SET_LOW_WORD(w,0); - c = (t-w*w)/(s+w); - r = p/q; -diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp ---- a/modules/fdlibm/src/e_hypot.cpp -+++ b/modules/fdlibm/src/e_hypot.cpp -@@ -41,16 +41,17 @@ - * hypot(x,y) is INF if x or y is +INF or -INF; else - * hypot(x,y) is NAN if x or y is NAN. - * - * Accuracy: - * hypot(x,y) returns sqrt(x^2+y^2) with error less - * than 1 ulps (units in the last place) - */ - -+#include <cmath> - #include <float.h> - - #include "math_private.h" - - double - __ieee754_hypot(double x, double y) - { - double a,b,t1,t2,y1,y2,w; -@@ -100,26 +101,26 @@ double - } - } - /* medium size a and b */ - w = a-b; - if (w>b) { - t1 = 0; - SET_HIGH_WORD(t1,ha); - t2 = a-t1; -- w = sqrt(t1*t1-(b*(-b)-t2*(a+t1))); -+ w = std::sqrt(t1*t1-(b*(-b)-t2*(a+t1))); - } else { - a = a+a; - y1 = 0; - SET_HIGH_WORD(y1,hb); - y2 = b - y1; - t1 = 0; - SET_HIGH_WORD(t1,ha+0x00100000); - t2 = a - t1; -- w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); -+ w = std::sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); - } - if(k!=0) { - u_int32_t high; - t1 = 1.0; - GET_HIGH_WORD(high,t1); - SET_HIGH_WORD(t1,high+(k<<20)); - return t1*w; - } else return w; -diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp ---- a/modules/fdlibm/src/e_pow.cpp -+++ b/modules/fdlibm/src/e_pow.cpp -@@ -52,16 +52,18 @@ - * - * Constants : - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -+#include <cmath> -+ - #include <float.h> - #include "math_private.h" - - static const double - bp[] = {1.0, 1.5,}, - dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ - dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ - zero = 0.0, -@@ -151,17 +153,17 @@ double - return (hy<0)?-y: zero; - } - if(iy==0x3ff00000) { /* y is +-1 */ - if(hy<0) return one/x; else return x; - } - if(hy==0x40000000) return x*x; /* y is 2 */ - if(hy==0x3fe00000) { /* y is 0.5 */ - if(hx>=0) /* x >= +0 */ -- return sqrt(x); -+ return std::sqrt(x); - } - } - - ax = fabs(x); - /* special value of x */ - if(lx==0) { - if(ix==0x7ff00000||ix==0||ix==0x3ff00000){ - z = ax; /*x is +-0,+-inf,+-1*/ -diff --git a/modules/fdlibm/src/s_asinh.cpp b/modules/fdlibm/src/s_asinh.cpp ---- a/modules/fdlibm/src/s_asinh.cpp -+++ b/modules/fdlibm/src/s_asinh.cpp -@@ -19,16 +19,17 @@ - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log(x)+ln2)) for large |x|, else - * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else - * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) - */ - -+#include <cmath> - #include <float.h> - - #include "math_private.h" - - static const double - one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ - ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ - huge= 1.00000000000000000000e+300; -@@ -43,15 +44,15 @@ asinh(double x) - if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ - if(ix< 0x3e300000) { /* |x|<2**-28 */ - if(huge+x>one) return x; /* return x inexact except 0 */ - } - if(ix>0x41b00000) { /* |x| > 2**28 */ - w = __ieee754_log(fabs(x))+ln2; - } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ - t = fabs(x); -- w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t)); -+ w = __ieee754_log(2.0*t+one/(std::sqrt(x*x+one)+t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; -- w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t))); -+ w =log1p(fabs(x)+t/(one+std::sqrt(one+t))); - } - if(hx>0) return w; else return -w; - } diff --git a/modules/fdlibm/patches/19_remove_unneeded_round_to_integer_helpers.patch b/modules/fdlibm/patches/19_remove_unneeded_round_to_integer_helpers.patch deleted file mode 100644 index 6d1baa23a..000000000 --- a/modules/fdlibm/patches/19_remove_unneeded_round_to_integer_helpers.patch +++ /dev/null @@ -1,130 +0,0 @@ -diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h ---- a/modules/fdlibm/src/math_private.h -+++ b/modules/fdlibm/src/math_private.h -@@ -586,126 +586,16 @@ CMPLXL(long double x, long double y) - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); - } - #endif - - #endif /* _COMPLEX_H */ - --/* -- * The rnint() family rounds to the nearest integer for a restricted range -- * range of args (up to about 2**MANT_DIG). We assume that the current -- * rounding mode is FE_TONEAREST so that this can be done efficiently. -- * Extra precision causes more problems in practice, and we only centralize -- * this here to reduce those problems, and have not solved the efficiency -- * problems. The exp2() family uses a more delicate version of this that -- * requires extracting bits from the intermediate value, so it is not -- * centralized here and should copy any solution of the efficiency problems. -- */ -- --static inline double --rnint(__double_t x) --{ -- /* -- * This casts to double to kill any extra precision. This depends -- * on the cast being applied to a double_t to avoid compiler bugs -- * (this is a cleaner version of STRICT_ASSIGN()). This is -- * inefficient if there actually is extra precision, but is hard -- * to improve on. We use double_t in the API to minimise conversions -- * for just calling here. Note that we cannot easily change the -- * magic number to the one that works directly with double_t, since -- * the rounding precision is variable at runtime on x86 so the -- * magic number would need to be variable. Assuming that the -- * rounding precision is always the default is too fragile. This -- * and many other complications will move when the default is -- * changed to FP_PE. -- */ -- return ((double)(x + 0x1.8p52) - 0x1.8p52); --} -- --static inline float --rnintf(__float_t x) --{ -- /* -- * As for rnint(), except we could just call that to handle the -- * extra precision case, usually without losing efficiency. -- */ -- return ((float)(x + 0x1.8p23F) - 0x1.8p23F); --} -- --#ifdef LDBL_MANT_DIG --/* -- * The complications for extra precision are smaller for rnintl() since it -- * can safely assume that the rounding precision has been increased from -- * its default to FP_PE on x86. We don't exploit that here to get small -- * optimizations from limiting the rangle to double. We just need it for -- * the magic number to work with long doubles. ld128 callers should use -- * rnint() instead of this if possible. ld80 callers should prefer -- * rnintl() since for amd64 this avoids swapping the register set, while -- * for i386 it makes no difference (assuming FP_PE), and for other arches -- * it makes little difference. -- */ --static inline long double --rnintl(long double x) --{ -- return (x + __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2 - -- __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2); --} --#endif /* LDBL_MANT_DIG */ -- --/* -- * irint() and i64rint() give the same result as casting to their integer -- * return type provided their arg is a floating point integer. They can -- * sometimes be more efficient because no rounding is required. -- */ --#if (defined(amd64) || defined(__i386__)) && defined(__GNUCLIKE_ASM) --#define irint(x) \ -- (sizeof(x) == sizeof(float) && \ -- sizeof(__float_t) == sizeof(long double) ? irintf(x) : \ -- sizeof(x) == sizeof(double) && \ -- sizeof(__double_t) == sizeof(long double) ? irintd(x) : \ -- sizeof(x) == sizeof(long double) ? irintl(x) : (int)(x)) --#else --#define irint(x) ((int)(x)) --#endif -- --#define i64rint(x) ((int64_t)(x)) /* only needed for ld128 so not opt. */ -- --#if defined(__i386__) && defined(__GNUCLIKE_ASM) --static __inline int --irintf(float x) --{ -- int n; -- -- __asm("fistl %0" : "=m" (n) : "t" (x)); -- return (n); --} -- --static __inline int --irintd(double x) --{ -- int n; -- -- __asm("fistl %0" : "=m" (n) : "t" (x)); -- return (n); --} --#endif -- --#if (defined(__amd64__) || defined(__i386__)) && defined(__GNUCLIKE_ASM) --static __inline int --irintl(long double x) --{ -- int n; -- -- __asm("fistl %0" : "=m" (n) : "t" (x)); -- return (n); --} --#endif -- - #ifdef DEBUG - #if defined(__amd64__) || defined(__i386__) - #define breakpoint() asm("int $3") - #else - #include <signal.h> - - #define breakpoint() raise(SIGTRAP) - #endif diff --git a/modules/fdlibm/src/e_acos.cpp b/modules/fdlibm/src/e_acos.cpp deleted file mode 100644 index 4f497b3b3..000000000 --- a/modules/fdlibm/src/e_acos.cpp +++ /dev/null @@ -1,107 +0,0 @@ - -/* @(#)e_acos.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$"); - -/* __ieee754_acos(x) - * Method : - * acos(x) = pi/2 - asin(x) - * acos(-x) = pi/2 + asin(x) - * For |x|<=0.5 - * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) - * For x>0.5 - * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) - * = 2asin(sqrt((1-x)/2)) - * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) - * = 2f + (2c + 2s*z*R(z)) - * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term - * for f so that f+c ~ sqrt(z). - * For x<-0.5 - * acos(x) = pi - 2asin(sqrt((1-|x|)/2)) - * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z) - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN with invalid signal. - * - * Function needed: sqrt - */ - -#include <cmath> -#include <float.h> - -#include "math_private.h" - -static const double -one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ -pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ -static volatile double -pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */ -static const double -pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ -pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ -pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ -pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ -pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ -pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ -qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ -qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ -qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ -qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ - -double -__ieee754_acos(double x) -{ - double z,p,q,r,w,s,c,df; - int32_t hx,ix; - GET_HIGH_WORD(hx,x); - ix = hx&0x7fffffff; - if(ix>=0x3ff00000) { /* |x| >= 1 */ - u_int32_t lx; - GET_LOW_WORD(lx,x); - if(((ix-0x3ff00000)|lx)==0) { /* |x|==1 */ - if(hx>0) return 0.0; /* acos(1) = 0 */ - else return pi+2.0*pio2_lo; /* acos(-1)= pi */ - } - return (x-x)/(x-x); /* acos(|x|>1) is NaN */ - } - if(ix<0x3fe00000) { /* |x| < 0.5 */ - if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/ - z = x*x; - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - r = p/q; - return pio2_hi - (x - (pio2_lo-x*r)); - } else if (hx<0) { /* x < -0.5 */ - z = (one+x)*0.5; - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - s = std::sqrt(z); - r = p/q; - w = r*s-pio2_lo; - return pi - 2.0*(s+w); - } else { /* x > 0.5 */ - z = (one-x)*0.5; - s = std::sqrt(z); - df = s; - SET_LOW_WORD(df,0); - c = (z-df*df)/(s+df); - p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); - q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - r = p/q; - w = r*s+c; - return 2.0*(df+w); - } -} diff --git a/modules/fdlibm/src/e_acosh.cpp b/modules/fdlibm/src/e_acosh.cpp deleted file mode 100644 index ce52d5aaa..000000000 --- a/modules/fdlibm/src/e_acosh.cpp +++ /dev/null @@ -1,64 +0,0 @@ - -/* @(#)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$"); - -/* __ieee754_acosh(x) - * Method : - * Based on - * acosh(x) = log [ x + sqrt(x*x-1) ] - * we have - * acosh(x) := log(x)+ln2, if x is large; else - * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else - * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - */ - -#include <cmath> -#include <float.h> - -#include "math_private.h" - -static const double -one = 1.0, -ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ - -double -__ieee754_acosh(double x) -{ - double t; - int32_t hx; - u_int32_t lx; - EXTRACT_WORDS(hx,lx,x); - if(hx<0x3ff00000) { /* x < 1 */ - return (x-x)/(x-x); - } else if(hx >=0x41b00000) { /* x > 2**28 */ - if(hx >=0x7ff00000) { /* x is inf of NaN */ - return x+x; - } else - return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */ - } else if(((hx-0x3ff00000)|lx)==0) { - return 0.0; /* acosh(1) = 0 */ - } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ - t=x*x; - return __ieee754_log(2.0*x-one/(x+std::sqrt(t-one))); - } else { /* 1<x<2 */ - t = x-one; - return log1p(t+std::sqrt(2.0*t+t*t)); - } -} diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp deleted file mode 100644 index e896bde9e..000000000 --- a/modules/fdlibm/src/e_asin.cpp +++ /dev/null @@ -1,113 +0,0 @@ - -/* @(#)e_asin.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$"); - -/* __ieee754_asin(x) - * Method : - * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... - * we approximate asin(x) on [0,0.5] by - * asin(x) = x + x*x^2*R(x^2) - * where - * R(x^2) is a rational approximation of (asin(x)-x)/x^3 - * and its remez error is bounded by - * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) - * - * For x in [0.5,1] - * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) - * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; - * then for x>0.98 - * asin(x) = pi/2 - 2*(s+s*z*R(z)) - * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) - * For x<=0.98, let pio4_hi = pio2_hi/2, then - * f = hi part of s; - * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) - * and - * asin(x) = pi/2 - 2*(s+s*z*R(z)) - * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) - * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN with invalid signal. - * - */ - -#include <cmath> -#include <float.h> - -#include "math_private.h" - -static const double -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -huge = 1.000e+300, -pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ -pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ -pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ - /* coefficient for R(x^2) */ -pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ -pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ -pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ -pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ -pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ -pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ -qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ -qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ -qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ -qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ - -double -__ieee754_asin(double x) -{ - double t=0.0,w,p,q,c,r,s; - int32_t hx,ix; - GET_HIGH_WORD(hx,x); - ix = hx&0x7fffffff; - if(ix>= 0x3ff00000) { /* |x|>= 1 */ - u_int32_t lx; - GET_LOW_WORD(lx,x); - if(((ix-0x3ff00000)|lx)==0) - /* asin(1)=+-pi/2 with inexact */ - return x*pio2_hi+x*pio2_lo; - return (x-x)/(x-x); /* asin(|x|>1) is NaN */ - } else if (ix<0x3fe00000) { /* |x|<0.5 */ - if(ix<0x3e500000) { /* if |x| < 2**-26 */ - if(huge+x>one) return x;/* return x with inexact if x!=0*/ - } - t = x*x; - p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); - w = p/q; - return x+x*w; - } - /* 1> |x|>= 0.5 */ - w = one-fabs(x); - t = w*0.5; - p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); - s = std::sqrt(t); - if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ - w = p/q; - t = pio2_hi-(2.0*(s+s*w)-pio2_lo); - } else { - w = s; - SET_LOW_WORD(w,0); - c = (t-w*w)/(s+w); - r = p/q; - p = 2.0*s*r-(pio2_lo-2.0*c); - q = pio4_hi-2.0*w; - t = pio4_hi-(p-q); - } - if(hx>0) return t; else return -t; -} diff --git a/modules/fdlibm/src/e_atan2.cpp b/modules/fdlibm/src/e_atan2.cpp deleted file mode 100644 index f45ad187f..000000000 --- a/modules/fdlibm/src/e_atan2.cpp +++ /dev/null @@ -1,124 +0,0 @@ - -/* @(#)e_atan2.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$"); - -/* __ieee754_atan2(y,x) - * Method : - * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): - * ARG (x+iy) = arctan(y/x) ... if x > 0, - * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, - * - * Special cases: - * - * ATAN2((anything), NaN ) is NaN; - * ATAN2(NAN , (anything) ) is NaN; - * ATAN2(+-0, +(anything but NaN)) is +-0 ; - * ATAN2(+-0, -(anything but NaN)) is +-pi ; - * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; - * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; - * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; - * ATAN2(+-INF,+INF ) is +-pi/4 ; - * ATAN2(+-INF,-INF ) is +-3pi/4; - * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <float.h> - -#include "math_private.h" - -static volatile double -tiny = 1.0e-300; -static const double -zero = 0.0, -pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ -pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ -pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */ -static volatile double -pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ - -double -__ieee754_atan2(double y, double x) -{ - double z; - int32_t k,m,hx,hy,ix,iy; - u_int32_t lx,ly; - - EXTRACT_WORDS(hx,lx,x); - ix = hx&0x7fffffff; - EXTRACT_WORDS(hy,ly,y); - iy = hy&0x7fffffff; - if(((ix|((lx|-lx)>>31))>0x7ff00000)|| - ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */ - return nan_mix(x, y); - if(hx==0x3ff00000&&lx==0) return atan(y); /* x=1.0 */ - m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ - - /* when y = 0 */ - if((iy|ly)==0) { - switch(m) { - case 0: - case 1: return y; /* atan(+-0,+anything)=+-0 */ - case 2: return pi+tiny;/* atan(+0,-anything) = pi */ - case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ - } - } - /* when x = 0 */ - if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; - - /* when x is INF */ - if(ix==0x7ff00000) { - if(iy==0x7ff00000) { - switch(m) { - case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ - case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ - case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ - case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ - } - } else { - switch(m) { - case 0: return zero ; /* atan(+...,+INF) */ - case 1: return -zero ; /* atan(-...,+INF) */ - case 2: return pi+tiny ; /* atan(+...,-INF) */ - case 3: return -pi-tiny ; /* atan(-...,-INF) */ - } - } - } - /* when y is INF */ - if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; - - /* compute y/x */ - k = (iy-ix)>>20; - if(k > 60) { /* |y/x| > 2**60 */ - z=pi_o_2+0.5*pi_lo; - m&=1; - } - else if(hx<0&&k<-60) z=0.0; /* 0 > |y|/x > -2**-60 */ - else z=atan(fabs(y/x)); /* safe to do y/x */ - switch (m) { - case 0: return z ; /* atan(+,+) */ - case 1: return -z ; /* atan(-,+) */ - case 2: return pi-(z-pi_lo);/* atan(+,-) */ - default: /* case 3 */ - return (z-pi_lo)-pi;/* atan(-,-) */ - } -} diff --git a/modules/fdlibm/src/e_atanh.cpp b/modules/fdlibm/src/e_atanh.cpp deleted file mode 100644 index a8f0f8deb..000000000 --- a/modules/fdlibm/src/e_atanh.cpp +++ /dev/null @@ -1,63 +0,0 @@ - -/* @(#)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$"); - -/* __ieee754_atanh(x) - * Method : - * 1.Reduced x to positive by atanh(-x) = -atanh(x) - * 2.For x>=0.5 - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * For x<0.5 - * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) - * - * Special cases: - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * - */ - -#include <float.h> - -#include "math_private.h" - -static const double one = 1.0, huge = 1e300; -static const double zero = 0.0; - -double -__ieee754_atanh(double x) -{ - double t; - int32_t hx,ix; - u_int32_t lx; - EXTRACT_WORDS(hx,lx,x); - ix = hx&0x7fffffff; - if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ - return (x-x)/(x-x); - if(ix==0x3ff00000) - return x/zero; - if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ - SET_HIGH_WORD(x,ix); - if(ix<0x3fe00000) { /* x < 0.5 */ - t = x+x; - t = 0.5*log1p(t+t*x/(one-x)); - } else - t = 0.5*log1p((x+x)/(one-x)); - if(hx>=0) return t; else return -t; -} diff --git a/modules/fdlibm/src/e_cosh.cpp b/modules/fdlibm/src/e_cosh.cpp deleted file mode 100644 index 42cb277d4..000000000 --- a/modules/fdlibm/src/e_cosh.cpp +++ /dev/null @@ -1,80 +0,0 @@ - -/* @(#)e_cosh.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$"); - -/* __ieee754_cosh(x) - * Method : - * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 - * 1. Replace x by |x| (cosh(x) = cosh(-x)). - * 2. - * [ exp(x) - 1 ]^2 - * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- - * 2*exp(x) - * - * exp(x) + 1/exp(x) - * ln2/2 <= x <= 22 : cosh(x) := ------------------- - * 2 - * 22 <= x <= lnovft : cosh(x) := exp(x)/2 - * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) - * ln2ovft < x : cosh(x) := huge*huge (overflow) - * - * Special cases: - * cosh(x) is |x| if x is +INF, -INF, or NaN. - * only cosh(0)=1 is exact for finite x. - */ - -#include <float.h> - -#include "math_private.h" - -static const double one = 1.0, half=0.5, huge = 1.0e300; - -double -__ieee754_cosh(double x) -{ - double t,w; - int32_t ix; - - /* High word of |x|. */ - GET_HIGH_WORD(ix,x); - ix &= 0x7fffffff; - - /* x is INF or NaN */ - if(ix>=0x7ff00000) return x*x; - - /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ - if(ix<0x3fd62e43) { - t = expm1(fabs(x)); - w = one+t; - if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ - return one+(t*t)/(w+w); - } - - /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ - if (ix < 0x40360000) { - t = __ieee754_exp(fabs(x)); - return half*t+half/t; - } - - /* |x| in [22, log(maxdouble)] return half*exp(|x|) */ - if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x)); - - /* |x| in [log(maxdouble), overflowthresold] */ - if (ix<=0x408633CE) - return __ldexp_exp(fabs(x), -1); - - /* |x| > overflowthresold, cosh(x) overflow */ - return huge*huge; -} diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp deleted file mode 100644 index 92af819ce..000000000 --- a/modules/fdlibm/src/e_exp.cpp +++ /dev/null @@ -1,165 +0,0 @@ - -/* @(#)e_exp.c 1.6 04/04/22 */ -/* - * ==================================================== - * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. - * - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -//#include <sys/cdefs.h> -//__FBSDID("$FreeBSD$"); - -/* __ieee754_exp(x) - * Returns the exponential of x. - * - * Method - * 1. Argument reduction: - * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. - * Given x, find r and integer k such that - * - * x = k*ln2 + r, |r| <= 0.5*ln2. - * - * Here r will be represented as r = hi-lo for better - * accuracy. - * - * 2. Approximation of exp(r) by a special rational function on - * the interval [0,0.34658]: - * Write - * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... - * We use a special Remes algorithm on [0,0.34658] to generate - * a polynomial of degree 5 to approximate R. The maximum error - * of this polynomial approximation is bounded by 2**-59. In - * other words, - * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 - * (where z=r*r, and the values of P1 to P5 are listed below) - * and - * | 5 | -59 - * | 2.0+P1*z+...+P5*z - R(z) | <= 2 - * | | - * The computation of exp(r) thus becomes - * 2*r - * exp(r) = 1 + ------- - * R - r - * r*R1(r) - * = 1 + r + ----------- (for better accuracy) - * 2 - R1(r) - * where - * 2 4 10 - * R1(r) = r - (P1*r + P2*r + ... + P5*r ). - * - * 3. Scale back to obtain exp(x): - * From step 1, we have - * exp(x) = 2^k * exp(r) - * - * Special cases: - * exp(INF) is INF, exp(NaN) is NaN; - * exp(-INF) is 0, and - * for finite argument, only exp(0)=1 is exact. - * - * Accuracy: - * according to an error analysis, the error is always less than - * 1 ulp (unit in the last place). - * - * Misc. info. - * For IEEE double - * if x > 7.09782712893383973096e+02 then exp(x) overflow - * if x < -7.45133219101941108420e+02 then exp(x) underflow - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <float.h> - -#include "math_private.h" - -static const double -one = 1.0, -halF[2] = {0.5,-0.5,}, -o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ -u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ -ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ - -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ -ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ - -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ -invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ -P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ -P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ -P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ -P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ -P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ - -static const double E = 2.7182818284590452354; /* e */ - -static volatile double -huge = 1.0e+300, -twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/ - -double -__ieee754_exp(double x) /* default IEEE double exp */ -{ - double y,hi=0.0,lo=0.0,c,t,twopk; - int32_t k=0,xsb; - u_int32_t hx; - - GET_HIGH_WORD(hx,x); - xsb = (hx>>31)&1; /* sign bit of x */ - hx &= 0x7fffffff; /* high word of |x| */ - - /* filter out non-finite argument */ - if(hx >= 0x40862E42) { /* if |x|>=709.78... */ - if(hx>=0x7ff00000) { - u_int32_t lx; - GET_LOW_WORD(lx,x); - if(((hx&0xfffff)|lx)!=0) - return x+x; /* NaN */ - else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ - } - if(x > o_threshold) return huge*huge; /* overflow */ - if(x < u_threshold) return twom1000*twom1000; /* underflow */ - } - - /* argument reduction */ - if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ - if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ - if (x == 1.0) return E; - hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; - } else { - k = (int)(invln2*x+halF[xsb]); - t = k; - hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ - lo = t*ln2LO[0]; - } - STRICT_ASSIGN(double, x, hi - lo); - } - else if(hx < 0x3e300000) { /* when |x|<2**-28 */ - if(huge+x>one) return one+x;/* trigger inexact */ - } - else k = 0; - - /* x is now in primary range */ - t = x*x; - if(k >= -1021) - INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0); - else - INSERT_WORDS(twopk,((u_int32_t)(0x3ff+(k+1000)))<<20, 0); - c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - if(k==0) return one-((x*c)/(c-2.0)-x); - else y = one-((lo-(x*c)/(2.0-c))-hi); - if(k >= -1021) { - if (k==1024) { - double const_0x1p1023 = pow(2, 1023); - return y*2.0*const_0x1p1023; - } - return y*twopk; - } else { - return y*twopk*twom1000; - } -} diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp deleted file mode 100644 index a23571150..000000000 --- a/modules/fdlibm/src/e_hypot.cpp +++ /dev/null @@ -1,127 +0,0 @@ - -/* @(#)e_hypot.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$"); - -/* __ieee754_hypot(x,y) - * - * Method : - * If (assume round-to-nearest) z=x*x+y*y - * has error less than sqrt(2)/2 ulp, than - * sqrt(z) has error less than 1 ulp (exercise). - * - * So, compute sqrt(x*x+y*y) with some care as - * follows to get the error below 1 ulp: - * - * Assume x>y>0; - * (if possible, set rounding to round-to-nearest) - * 1. if x > 2y use - * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y - * where x1 = x with lower 32 bits cleared, x2 = x-x1; else - * 2. if x <= 2y use - * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) - * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, - * y1= y with lower 32 bits chopped, y2 = y-y1. - * - * NOTE: scaling may be necessary if some argument is too - * large or too tiny - * - * Special cases: - * hypot(x,y) is INF if x or y is +INF or -INF; else - * hypot(x,y) is NAN if x or y is NAN. - * - * Accuracy: - * hypot(x,y) returns sqrt(x^2+y^2) with error less - * than 1 ulps (units in the last place) - */ - -#include <cmath> -#include <float.h> - -#include "math_private.h" - -double -__ieee754_hypot(double x, double y) -{ - double a,b,t1,t2,y1,y2,w; - int32_t j,k,ha,hb; - - GET_HIGH_WORD(ha,x); - ha &= 0x7fffffff; - GET_HIGH_WORD(hb,y); - hb &= 0x7fffffff; - if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} - a = fabs(a); - b = fabs(b); - if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */ - k=0; - if(ha > 0x5f300000) { /* a>2**500 */ - if(ha >= 0x7ff00000) { /* Inf or NaN */ - u_int32_t low; - /* Use original arg order iff result is NaN; quieten sNaNs. */ - w = fabsl(x+0.0L)-fabs(y+0); - GET_LOW_WORD(low,a); - if(((ha&0xfffff)|low)==0) w = a; - GET_LOW_WORD(low,b); - if(((hb^0x7ff00000)|low)==0) w = b; - return w; - } - /* scale a and b by 2**-600 */ - ha -= 0x25800000; hb -= 0x25800000; k += 600; - SET_HIGH_WORD(a,ha); - SET_HIGH_WORD(b,hb); - } - if(hb < 0x20b00000) { /* b < 2**-500 */ - if(hb <= 0x000fffff) { /* subnormal b or 0 */ - u_int32_t low; - GET_LOW_WORD(low,b); - if((hb|low)==0) return a; - t1=0; - SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */ - b *= t1; - a *= t1; - k -= 1022; - } else { /* scale a and b by 2^600 */ - ha += 0x25800000; /* a *= 2^600 */ - hb += 0x25800000; /* b *= 2^600 */ - k -= 600; - SET_HIGH_WORD(a,ha); - SET_HIGH_WORD(b,hb); - } - } - /* medium size a and b */ - w = a-b; - if (w>b) { - t1 = 0; - SET_HIGH_WORD(t1,ha); - t2 = a-t1; - w = std::sqrt(t1*t1-(b*(-b)-t2*(a+t1))); - } else { - a = a+a; - y1 = 0; - SET_HIGH_WORD(y1,hb); - y2 = b - y1; - t1 = 0; - SET_HIGH_WORD(t1,ha+0x00100000); - t2 = a - t1; - w = std::sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); - } - if(k!=0) { - u_int32_t high; - t1 = 1.0; - GET_HIGH_WORD(high,t1); - SET_HIGH_WORD(t1,high+(k<<20)); - return t1*w; - } else return w; -} diff --git a/modules/fdlibm/src/e_log.cpp b/modules/fdlibm/src/e_log.cpp deleted file mode 100644 index fa2da8fcb..000000000 --- a/modules/fdlibm/src/e_log.cpp +++ /dev/null @@ -1,142 +0,0 @@ - -/* @(#)e_log.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$"); - -/* __ieee754_log(x) - * Return the logrithm of x - * - * Method : - * 1. Argument Reduction: find k and f such that - * x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * 2. Approximation of log(1+f). - * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R - * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error - * of this polynomial approximation is bounded by 2**-58.45. In - * other words, - * 2 4 6 8 10 12 14 - * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s - * (the values of Lg1 to Lg7 are listed in the program) - * and - * | 2 14 | -58.45 - * | Lg1*s +...+Lg7*s - R(z) | <= 2 - * | | - * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. - * In order to guarantee error in log below 1ulp, we compute log - * by - * log(1+f) = f - s*(f - R) (if f is not too large) - * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) - * - * 3. Finally, log(x) = k*ln2 + log(1+f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) - * Here ln2 is split into two floating point number: - * ln2_hi + ln2_lo, - * where n*ln2_hi is always exact for |n| < 2000. - * - * Special cases: - * log(x) is NaN with signal if x < 0 (including -INF) ; - * log(+INF) is +INF; log(0) is -INF with signal; - * log(NaN) is that NaN with no signal. - * - * Accuracy: - * according to an error analysis, the error is always less than - * 1 ulp (unit in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <float.h> - -#include "math_private.h" - -static const double -ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ -ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ -two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ -Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ -Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ -Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ -Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ -Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ -Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ -Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ - -static const double zero = 0.0; -static volatile double vzero = 0.0; - -double -__ieee754_log(double x) -{ - double hfsq,f,s,z,R,w,t1,t2,dk; - int32_t k,hx,i,j; - u_int32_t lx; - - EXTRACT_WORDS(hx,lx,x); - - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/vzero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ - GET_HIGH_WORD(hx,x); - } - if (hx >= 0x7ff00000) return x+x; - k += (hx>>20)-1023; - hx &= 0x000fffff; - i = (hx+0x95f64)&0x100000; - SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ - k += (i>>20); - f = x-1.0; - if((0x000fffff&(2+hx))<3) { /* -2**-20 <= f < 2**-20 */ - if(f==zero) { - if(k==0) { - return zero; - } else { - dk=(double)k; - return dk*ln2_hi+dk*ln2_lo; - } - } - R = f*f*(0.5-0.33333333333333333*f); - if(k==0) return f-R; else {dk=(double)k; - return dk*ln2_hi-((R-dk*ln2_lo)-f);} - } - s = f/(2.0+f); - dk = (double)k; - z = s*s; - i = hx-0x6147a; - w = z*z; - j = 0x6b851-hx; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - i |= j; - R = t2+t1; - if(i>0) { - hfsq=0.5*f*f; - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); - } else { - if(k==0) return f-s*(f-R); else - return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); - } -} diff --git a/modules/fdlibm/src/e_log10.cpp b/modules/fdlibm/src/e_log10.cpp deleted file mode 100644 index ed6879885..000000000 --- a/modules/fdlibm/src/e_log10.cpp +++ /dev/null @@ -1,89 +0,0 @@ - -/* @(#)e_log10.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$"); - -/* - * Return the base 10 logarithm of x. See e_log.c and k_log.h for most - * comments. - * - * log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) - * in not-quite-routine extra precision. - */ - -#include <float.h> - -#include "math_private.h" -#include "k_log.h" - -static const double -two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ -ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ -log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ -log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ - -static const double zero = 0.0; -static volatile double vzero = 0.0; - -double -__ieee754_log10(double x) -{ - double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; - int32_t i,k,hx; - u_int32_t lx; - - EXTRACT_WORDS(hx,lx,x); - - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/vzero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ - GET_HIGH_WORD(hx,x); - } - if (hx >= 0x7ff00000) return x+x; - if (hx == 0x3ff00000 && lx == 0) - return zero; /* log(1) = +0 */ - k += (hx>>20)-1023; - hx &= 0x000fffff; - i = (hx+0x95f64)&0x100000; - SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ - k += (i>>20); - y = (double)k; - f = x - 1.0; - hfsq = 0.5*f*f; - r = k_log1p(f); - - /* See e_log2.c for most details. */ - hi = f - hfsq; - SET_LOW_WORD(hi,0); - lo = (f - hi) - hfsq + r; - val_hi = hi*ivln10hi; - y2 = y*log10_2hi; - val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi; - - /* - * Extra precision in for adding y*log10_2hi is not strictly needed - * since there is no very large cancellation near x = sqrt(2) or - * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs - * with some parallelism and it reduces the error for many args. - */ - w = y2 + val_hi; - val_lo += (y2 - w) + val_hi; - val_hi = w; - - return val_lo + val_hi; -} diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp deleted file mode 100644 index 5649fec44..000000000 --- a/modules/fdlibm/src/e_log2.cpp +++ /dev/null @@ -1,112 +0,0 @@ - -/* @(#)e_log10.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$"); - -/* - * Return the base 2 logarithm of x. See e_log.c and k_log.h for most - * comments. - * - * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel, - * then does the combining and scaling steps - * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k - * in not-quite-routine extra precision. - */ - -#include <float.h> - -#include "math_private.h" -#include "k_log.h" - -static const double -two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ -ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ - -static const double zero = 0.0; -static volatile double vzero = 0.0; - -double -__ieee754_log2(double x) -{ - double f,hfsq,hi,lo,r,val_hi,val_lo,w,y; - int32_t i,k,hx; - u_int32_t lx; - - EXTRACT_WORDS(hx,lx,x); - - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/vzero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ - GET_HIGH_WORD(hx,x); - } - if (hx >= 0x7ff00000) return x+x; - if (hx == 0x3ff00000 && lx == 0) - return zero; /* log(1) = +0 */ - k += (hx>>20)-1023; - hx &= 0x000fffff; - i = (hx+0x95f64)&0x100000; - SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ - k += (i>>20); - y = (double)k; - f = x - 1.0; - hfsq = 0.5*f*f; - r = k_log1p(f); - - /* - * f-hfsq must (for args near 1) be evaluated in extra precision - * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2). - * This is fairly efficient since f-hfsq only depends on f, so can - * be evaluated in parallel with R. Not combining hfsq with R also - * keeps R small (though not as small as a true `lo' term would be), - * so that extra precision is not needed for terms involving R. - * - * Compiler bugs involving extra precision used to break Dekker's - * theorem for spitting f-hfsq as hi+lo, unless double_t was used - * or the multi-precision calculations were avoided when double_t - * has extra precision. These problems are now automatically - * avoided as a side effect of the optimization of combining the - * Dekker splitting step with the clear-low-bits step. - * - * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra - * precision to avoid a very large cancellation when x is very near - * these values. Unlike the above cancellations, this problem is - * specific to base 2. It is strange that adding +-1 is so much - * harder than adding +-ln2 or +-log10_2. - * - * This uses Dekker's theorem to normalize y+val_hi, so the - * compiler bugs are back in some configurations, sigh. And I - * don't want to used double_t to avoid them, since that gives a - * pessimization and the support for avoiding the pessimization - * is not yet available. - * - * The multi-precision calculations for the multiplications are - * routine. - */ - hi = f - hfsq; - SET_LOW_WORD(hi,0); - lo = (f - hi) - hfsq + r; - val_hi = hi*ivln2hi; - val_lo = (lo+hi)*ivln2lo + lo*ivln2hi; - - /* spadd(val_hi, val_lo, y), except for not using double_t: */ - w = y + val_hi; - val_lo += (y - w) + val_hi; - val_hi = w; - - return val_lo + val_hi; -} diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp deleted file mode 100644 index c18226b8a..000000000 --- a/modules/fdlibm/src/e_pow.cpp +++ /dev/null @@ -1,311 +0,0 @@ -/* @(#)e_pow.c 1.5 04/04/22 SMI */ -/* - * ==================================================== - * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. - * - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -//#include <sys/cdefs.h> -//__FBSDID("$FreeBSD$"); - -/* __ieee754_pow(x,y) return x**y - * - * n - * Method: Let x = 2 * (1+f) - * 1. Compute and return log2(x) in two pieces: - * log2(x) = w1 + w2, - * where w1 has 53-24 = 29 bit trailing zeros. - * 2. Perform y*log2(x) = n+y' by simulating multi-precision - * arithmetic, where |y'|<=0.5. - * 3. Return x**y = 2**n*exp(y'*log2) - * - * Special cases: - * 1. (anything) ** 0 is 1 - * 2. (anything) ** 1 is itself - * 3. (anything) ** NAN is NAN except 1 ** NAN = 1 - * 4. NAN ** (anything except 0) is NAN - * 5. +-(|x| > 1) ** +INF is +INF - * 6. +-(|x| > 1) ** -INF is +0 - * 7. +-(|x| < 1) ** +INF is +0 - * 8. +-(|x| < 1) ** -INF is +INF - * 9. +-1 ** +-INF is 1 - * 10. +0 ** (+anything except 0, NAN) is +0 - * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 - * 12. +0 ** (-anything except 0, NAN) is +INF - * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF - * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) - * 15. +INF ** (+anything except 0,NAN) is +INF - * 16. +INF ** (-anything except 0,NAN) is +0 - * 17. -INF ** (anything) = -0 ** (-anything) - * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) - * 19. (-anything except 0 and inf) ** (non-integer) is NAN - * - * Accuracy: - * pow(x,y) returns x**y nearly rounded. In particular - * pow(integer,integer) - * always returns the correct integer provided it is - * representable. - * - * Constants : - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <cmath> - -#include <float.h> -#include "math_private.h" - -static const double -bp[] = {1.0, 1.5,}, -dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ -dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ -zero = 0.0, -half = 0.5, -qrtr = 0.25, -thrd = 3.3333333333333331e-01, /* 0x3fd55555, 0x55555555 */ -one = 1.0, -two = 2.0, -two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ -huge = 1.0e300, -tiny = 1.0e-300, - /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ -L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ -L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ -L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ -L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ -L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ -L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ -P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ -P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ -P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ -P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ -P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ -lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ -lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ -lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ -ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ -cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ -cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ -cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ -ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ -ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ -ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ - -double -__ieee754_pow(double x, double y) -{ - double z,ax,z_h,z_l,p_h,p_l; - double y1,t1,t2,r,s,t,u,v,w; - int32_t i,j,k,yisint,n; - int32_t hx,hy,ix,iy; - u_int32_t lx,ly; - - EXTRACT_WORDS(hx,lx,x); - EXTRACT_WORDS(hy,ly,y); - ix = hx&0x7fffffff; iy = hy&0x7fffffff; - - /* y==zero: x**0 = 1 */ - if((iy|ly)==0) return one; - - /* x==1: 1**y = 1, even if y is NaN */ - if (hx==0x3ff00000 && lx == 0) return one; - - /* y!=zero: result is NaN if either arg is NaN */ - if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || - iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) - return nan_mix(x, y); - - /* determine if y is an odd int when x < 0 - * yisint = 0 ... y is not an integer - * yisint = 1 ... y is an odd int - * yisint = 2 ... y is an even int - */ - yisint = 0; - if(hx<0) { - if(iy>=0x43400000) yisint = 2; /* even integer y */ - else if(iy>=0x3ff00000) { - k = (iy>>20)-0x3ff; /* exponent */ - if(k>20) { - j = ly>>(52-k); - if(((u_int32_t)j<<(52-k))==ly) yisint = 2-(j&1); - } else if(ly==0) { - j = iy>>(20-k); - if((j<<(20-k))==iy) yisint = 2-(j&1); - } - } - } - - /* special value of y */ - if(ly==0) { - if (iy==0x7ff00000) { /* y is +-inf */ - if(((ix-0x3ff00000)|lx)==0) - return one; /* (-1)**+-inf is 1 */ - else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ - return (hy>=0)? y: zero; - else /* (|x|<1)**-,+inf = inf,0 */ - return (hy<0)?-y: zero; - } - if(iy==0x3ff00000) { /* y is +-1 */ - if(hy<0) return one/x; else return x; - } - if(hy==0x40000000) return x*x; /* y is 2 */ - if(hy==0x3fe00000) { /* y is 0.5 */ - if(hx>=0) /* x >= +0 */ - return std::sqrt(x); - } - } - - ax = fabs(x); - /* special value of x */ - if(lx==0) { - if(ix==0x7ff00000||ix==0||ix==0x3ff00000){ - z = ax; /*x is +-0,+-inf,+-1*/ - if(hy<0) z = one/z; /* z = (1/|x|) */ - if(hx<0) { - if(((ix-0x3ff00000)|yisint)==0) { - z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) - z = -z; /* (x<0)**odd = -(|x|**odd) */ - } - return z; - } - } - - /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be - n = (hx>>31)+1; - but ANSI C says a right shift of a signed negative quantity is - implementation defined. */ - n = ((u_int32_t)hx>>31)-1; - - /* (x<0)**(non-int) is NaN */ - if((n|yisint)==0) return (x-x)/(x-x); - - s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ - if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */ - - /* |y| is huge */ - if(iy>0x41e00000) { /* if |y| > 2**31 */ - if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ - if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; - if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; - } - /* over/underflow if x is not close to one */ - if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny; - if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny; - /* now |1-x| is tiny <= 2**-20, suffice to compute - log(x) by x-x^2/2+x^3/3-x^4/4 */ - t = ax-one; /* t has 20 trailing zeros */ - w = (t*t)*(half-t*(thrd-t*qrtr)); - u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ - v = t*ivln2_l-w*ivln2; - t1 = u+v; - SET_LOW_WORD(t1,0); - t2 = v-(t1-u); - } else { - double ss,s2,s_h,s_l,t_h,t_l; - n = 0; - /* take care subnormal number */ - if(ix<0x00100000) - {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); } - n += ((ix)>>20)-0x3ff; - j = ix&0x000fffff; - /* determine interval */ - ix = j|0x3ff00000; /* normalize ix */ - if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ - else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ - else {k=0;n+=1;ix -= 0x00100000;} - SET_HIGH_WORD(ax,ix); - - /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ - u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ - v = one/(ax+bp[k]); - ss = u*v; - s_h = ss; - SET_LOW_WORD(s_h,0); - /* t_h=ax+bp[k] High */ - t_h = zero; - SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18)); - t_l = ax - (t_h-bp[k]); - s_l = v*((u-s_h*t_h)-s_h*t_l); - /* compute log(ax) */ - s2 = ss*ss; - r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); - r += s_l*(s_h+ss); - s2 = s_h*s_h; - t_h = 3+s2+r; - SET_LOW_WORD(t_h,0); - t_l = r-((t_h-3)-s2); - /* u+v = ss*(1+...) */ - u = s_h*t_h; - v = s_l*t_h+t_l*ss; - /* 2/(3log2)*(ss+...) */ - p_h = u+v; - SET_LOW_WORD(p_h,0); - p_l = v-(p_h-u); - z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ - z_l = cp_l*p_h+p_l*cp+dp_l[k]; - /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ - t = n; - t1 = (((z_h+z_l)+dp_h[k])+t); - SET_LOW_WORD(t1,0); - t2 = z_l-(((t1-t)-dp_h[k])-z_h); - } - - /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ - y1 = y; - SET_LOW_WORD(y1,0); - p_l = (y-y1)*t1+y*t2; - p_h = y1*t1; - z = p_l+p_h; - EXTRACT_WORDS(j,i,z); - if (j>=0x40900000) { /* z >= 1024 */ - if(((j-0x40900000)|i)!=0) /* if z > 1024 */ - return s*huge*huge; /* overflow */ - else { - if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ - } - } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */ - if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ - return s*tiny*tiny; /* underflow */ - else { - if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ - } - } - /* - * compute 2**(p_h+p_l) - */ - i = j&0x7fffffff; - k = (i>>20)-0x3ff; - n = 0; - if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ - n = j+(0x00100000>>(k+1)); - k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ - t = zero; - SET_HIGH_WORD(t,n&~(0x000fffff>>k)); - n = ((n&0x000fffff)|0x00100000)>>(20-k); - if(j<0) n = -n; - p_h -= t; - } - t = p_l+p_h; - SET_LOW_WORD(t,0); - u = t*lg2_h; - v = (p_l-(t-p_h))*lg2+t*lg2_l; - z = u+v; - w = v-(z-u); - t = z*z; - t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - r = (z*t1)/(t1-two)-(w+z*w); - z = one-(r-z); - GET_HIGH_WORD(j,z); - j += (n<<20); - if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */ - else SET_HIGH_WORD(z,j); - return s*z; -} diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp deleted file mode 100644 index c3418e687..000000000 --- a/modules/fdlibm/src/e_sinh.cpp +++ /dev/null @@ -1,74 +0,0 @@ - -/* @(#)e_sinh.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$"); - -/* __ieee754_sinh(x) - * Method : - * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 - * 1. Replace x by |x| (sinh(-x) = -sinh(x)). - * 2. - * E + E/(E+1) - * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) - * 2 - * - * 22 <= x <= lnovft : sinh(x) := exp(x)/2 - * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) - * ln2ovft < x : sinh(x) := x*shuge (overflow) - * - * Special cases: - * sinh(x) is |x| if x is +INF, -INF, or NaN. - * only sinh(0)=0 is exact for finite x. - */ - -#include <float.h> - -#include "math_private.h" - -static const double one = 1.0, shuge = 1.0e307; - -double -__ieee754_sinh(double x) -{ - double t,h; - int32_t ix,jx; - - /* High word of |x|. */ - GET_HIGH_WORD(jx,x); - ix = jx&0x7fffffff; - - /* x is INF or NaN */ - if(ix>=0x7ff00000) return x+x; - - h = 0.5; - if (jx<0) h = -h; - /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ - if (ix < 0x40360000) { /* |x|<22 */ - if (ix<0x3e300000) /* |x|<2**-28 */ - if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ - t = expm1(fabs(x)); - if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one)); - return h*(t+t/(t+one)); - } - - /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ - if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x)); - - /* |x| in [log(maxdouble), overflowthresold] */ - if (ix<=0x408633CE) - return h*2.0*__ldexp_exp(fabs(x), -1); - - /* |x| > overflowthresold, sinh(x) overflow */ - return x*shuge; -} diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h deleted file mode 100644 index 324e5d0b0..000000000 --- a/modules/fdlibm/src/fdlibm.h +++ /dev/null @@ -1,64 +0,0 @@ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - -#ifndef mozilla_imported_fdlibm_h -#define mozilla_imported_fdlibm_h - -namespace fdlibm { - -double acos(double); -double asin(double); -double atan(double); -double atan2(double, double); - -double cosh(double); -double sinh(double); -double tanh(double); - -double exp(double); -double log(double); -double log10(double); - -double pow(double, double); -double fabs(double); - -double floor(double); -double trunc(double); -double ceil(double); - -double acosh(double); -double asinh(double); -double atanh(double); -double cbrt(double); -double expm1(double); -double hypot(double, double); -double log1p(double); -double log2(double); -double rint(double); -double copysign(double, double); -double nearbyint(double); -double scalbn(double, int); - -float ceilf(float); -float floorf(float); - -float nearbyintf(float); -float rintf(float); -float truncf(float); - -} /* namespace fdlibm */ - -#endif /* mozilla_imported_fdlibm_h */ diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp deleted file mode 100644 index 9394c8fd8..000000000 --- a/modules/fdlibm/src/k_exp.cpp +++ /dev/null @@ -1,83 +0,0 @@ -/*- - * SPDX-License-Identifier: BSD-2-Clause-FreeBSD - * - * Copyright (c) 2011 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 "math_private.h" - -static const uint32_t k = 1799; /* constant for reduction */ -static const double kln2 = 1246.97177782734161156; /* k * ln2 */ - -/* - * Compute exp(x), scaled to avoid spurious overflow. An exponent is - * returned separately in 'expt'. - * - * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 - * Output: 2**1023 <= y < 2**1024 - */ -static double -__frexp_exp(double x, int *expt) -{ - double exp_x; - uint32_t hx; - - /* - * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to - * minimize |exp(kln2) - 2**k|. We also scale the exponent of - * exp_x to MAX_EXP so that the result can be multiplied by - * a tiny number without losing accuracy due to denormalization. - */ - exp_x = exp(x - kln2); - GET_HIGH_WORD(hx, exp_x); - *expt = (hx >> 20) - (0x3ff + 1023) + k; - SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); - return (exp_x); -} - -/* - * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt. - * They are intended for large arguments (real part >= ln(DBL_MAX)) - * where care is needed to avoid overflow. - * - * The present implementation is narrowly tailored for our hyperbolic and - * exponential functions. We assume expt is small (0 or -1), and the caller - * has filtered out very large x, for which overflow would be inevitable. - */ - -double -__ldexp_exp(double x, int expt) -{ - double exp_x, scale; - int ex_expt; - - exp_x = __frexp_exp(x, &ex_expt); - expt += ex_expt; - INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); - return (exp_x * scale); -} diff --git a/modules/fdlibm/src/k_log.h b/modules/fdlibm/src/k_log.h deleted file mode 100644 index 0efa020f6..000000000 --- a/modules/fdlibm/src/k_log.h +++ /dev/null @@ -1,100 +0,0 @@ - -/* @(#)e_log.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$"); - -/* - * k_log1p(f): - * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)]. - * - * The following describes the overall strategy for computing - * logarithms in base e. The argument reduction and adding the final - * term of the polynomial are done by the caller for increased accuracy - * when different bases are used. - * - * Method : - * 1. Argument Reduction: find k and f such that - * x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * 2. Approximation of log(1+f). - * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R - * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error - * of this polynomial approximation is bounded by 2**-58.45. In - * other words, - * 2 4 6 8 10 12 14 - * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s - * (the values of Lg1 to Lg7 are listed in the program) - * and - * | 2 14 | -58.45 - * | Lg1*s +...+Lg7*s - R(z) | <= 2 - * | | - * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. - * In order to guarantee error in log below 1ulp, we compute log - * by - * log(1+f) = f - s*(f - R) (if f is not too large) - * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) - * - * 3. Finally, log(x) = k*ln2 + log(1+f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) - * Here ln2 is split into two floating point number: - * ln2_hi + ln2_lo, - * where n*ln2_hi is always exact for |n| < 2000. - * - * Special cases: - * log(x) is NaN with signal if x < 0 (including -INF) ; - * log(+INF) is +INF; log(0) is -INF with signal; - * log(NaN) is that NaN with no signal. - * - * Accuracy: - * according to an error analysis, the error is always less than - * 1 ulp (unit in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -static const double -Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ -Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ -Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ -Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ -Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ -Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ -Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ - -/* - * We always inline k_log1p(), since doing so produces a - * substantial performance improvement (~40% on amd64). - */ -static inline double -k_log1p(double f) -{ - double hfsq,s,z,R,w,t1,t2; - - s = f/(2.0+f); - z = s*s; - w = z*z; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - R = t2+t1; - hfsq=0.5*f*f; - return s*(hfsq+R); -} diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h deleted file mode 100644 index d9ec44817..000000000 --- a/modules/fdlibm/src/math_private.h +++ /dev/null @@ -1,861 +0,0 @@ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - -#ifndef _MATH_PRIVATE_H_ -#define _MATH_PRIVATE_H_ - -#include <cfloat> -#include <stdint.h> -#include <sys/types.h> - -#include "fdlibm.h" - -#include "mozilla/EndianUtils.h" - -/* - * The original fdlibm code used statements like: - * n0 = ((*(int*)&one)>>29)^1; * index of high word * - * ix0 = *(n0+(int*)&x); * high word of x * - * ix1 = *((1-n0)+(int*)&x); * low word of x * - * to dig two 32 bit words out of the 64 bit IEEE floating point - * value. That is non-ANSI, and, moreover, the gcc instruction - * scheduler gets it wrong. We instead use the following macros. - * Unlike the original code, we determine the endianness at compile - * time, not at run time; I don't see much benefit to selecting - * endianness at run time. - */ - -#ifndef u_int32_t -#define u_int32_t uint32_t -#endif -#ifndef u_int64_t -#define u_int64_t uint64_t -#endif - -/* A union which permits us to convert between a long double and - four 32 bit ints. */ - -#if MOZ_BIG_ENDIAN - -typedef union -{ - long double value; - struct { - u_int32_t mswhi; - u_int32_t mswlo; - u_int32_t lswhi; - u_int32_t lswlo; - } parts32; - struct { - u_int64_t msw; - u_int64_t lsw; - } parts64; -} ieee_quad_shape_type; - -#endif - -#if MOZ_LITTLE_ENDIAN - -typedef union -{ - long double value; - struct { - u_int32_t lswlo; - u_int32_t lswhi; - u_int32_t mswlo; - u_int32_t mswhi; - } parts32; - struct { - u_int64_t lsw; - u_int64_t msw; - } parts64; -} ieee_quad_shape_type; - -#endif - -/* - * A union which permits us to convert between a double and two 32 bit - * ints. - */ - -#if MOZ_BIG_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t msw; - u_int32_t lsw; - } parts; - struct - { - u_int64_t w; - } xparts; -} ieee_double_shape_type; - -#endif - -#if MOZ_LITTLE_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t lsw; - u_int32_t msw; - } parts; - struct - { - u_int64_t w; - } xparts; -} ieee_double_shape_type; - -#endif - -/* Get two 32 bit ints from a double. */ - -#define EXTRACT_WORDS(ix0,ix1,d) \ -do { \ - ieee_double_shape_type ew_u; \ - ew_u.value = (d); \ - (ix0) = ew_u.parts.msw; \ - (ix1) = ew_u.parts.lsw; \ -} while (0) - -/* Get a 64-bit int from a double. */ -#define EXTRACT_WORD64(ix,d) \ -do { \ - ieee_double_shape_type ew_u; \ - ew_u.value = (d); \ - (ix) = ew_u.xparts.w; \ -} while (0) - -/* Get the more significant 32 bit int from a double. */ - -#define GET_HIGH_WORD(i,d) \ -do { \ - ieee_double_shape_type gh_u; \ - gh_u.value = (d); \ - (i) = gh_u.parts.msw; \ -} while (0) - -/* Get the less significant 32 bit int from a double. */ - -#define GET_LOW_WORD(i,d) \ -do { \ - ieee_double_shape_type gl_u; \ - gl_u.value = (d); \ - (i) = gl_u.parts.lsw; \ -} while (0) - -/* Set a double from two 32 bit ints. */ - -#define INSERT_WORDS(d,ix0,ix1) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.parts.msw = (ix0); \ - iw_u.parts.lsw = (ix1); \ - (d) = iw_u.value; \ -} while (0) - -/* Set a double from a 64-bit int. */ -#define INSERT_WORD64(d,ix) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.xparts.w = (ix); \ - (d) = iw_u.value; \ -} while (0) - -/* Set the more significant 32 bits of a double from an int. */ - -#define SET_HIGH_WORD(d,v) \ -do { \ - ieee_double_shape_type sh_u; \ - sh_u.value = (d); \ - sh_u.parts.msw = (v); \ - (d) = sh_u.value; \ -} while (0) - -/* Set the less significant 32 bits of a double from an int. */ - -#define SET_LOW_WORD(d,v) \ -do { \ - ieee_double_shape_type sl_u; \ - sl_u.value = (d); \ - sl_u.parts.lsw = (v); \ - (d) = sl_u.value; \ -} while (0) - -/* - * A union which permits us to convert between a float and a 32 bit - * int. - */ - -typedef union -{ - float value; - /* FIXME: Assumes 32 bit int. */ - unsigned int word; -} ieee_float_shape_type; - -/* Get a 32 bit int from a float. */ - -#define GET_FLOAT_WORD(i,d) \ -do { \ - ieee_float_shape_type gf_u; \ - gf_u.value = (d); \ - (i) = gf_u.word; \ -} while (0) - -/* Set a float from a 32 bit int. */ - -#define SET_FLOAT_WORD(d,i) \ -do { \ - ieee_float_shape_type sf_u; \ - sf_u.word = (i); \ - (d) = sf_u.value; \ -} while (0) - -/* - * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long - * double. - */ - -#define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \ -do { \ - union IEEEl2bits ew_u; \ - ew_u.e = (d); \ - (ix0) = ew_u.xbits.expsign; \ - (ix1) = ew_u.xbits.man; \ -} while (0) - -/* - * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit - * long double. - */ - -#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \ -do { \ - union IEEEl2bits ew_u; \ - ew_u.e = (d); \ - (ix0) = ew_u.xbits.expsign; \ - (ix1) = ew_u.xbits.manh; \ - (ix2) = ew_u.xbits.manl; \ -} while (0) - -/* Get expsign as a 16 bit int from a long double. */ - -#define GET_LDBL_EXPSIGN(i,d) \ -do { \ - union IEEEl2bits ge_u; \ - ge_u.e = (d); \ - (i) = ge_u.xbits.expsign; \ -} while (0) - -/* - * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int - * mantissa. - */ - -#define INSERT_LDBL80_WORDS(d,ix0,ix1) \ -do { \ - union IEEEl2bits iw_u; \ - iw_u.xbits.expsign = (ix0); \ - iw_u.xbits.man = (ix1); \ - (d) = iw_u.e; \ -} while (0) - -/* - * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints - * comprising the mantissa. - */ - -#define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \ -do { \ - union IEEEl2bits iw_u; \ - iw_u.xbits.expsign = (ix0); \ - iw_u.xbits.manh = (ix1); \ - iw_u.xbits.manl = (ix2); \ - (d) = iw_u.e; \ -} while (0) - -/* Set expsign of a long double from a 16 bit int. */ - -#define SET_LDBL_EXPSIGN(d,v) \ -do { \ - union IEEEl2bits se_u; \ - se_u.e = (d); \ - se_u.xbits.expsign = (v); \ - (d) = se_u.e; \ -} while (0) - -#ifdef __i386__ -/* Long double constants are broken on i386. */ -#define LD80C(m, ex, v) { \ - .xbits.man = __CONCAT(m, ULL), \ - .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \ -} -#else -/* The above works on non-i386 too, but we use this to check v. */ -#define LD80C(m, ex, v) { .e = (v), } -#endif - -#ifdef FLT_EVAL_METHOD -/* - * Attempt to get strict C99 semantics for assignment with non-C99 compilers. - */ -#if !defined(_MSC_VER) && (FLT_EVAL_METHOD == 0 || __GNUC__ == 0) -#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) -#else -#define STRICT_ASSIGN(type, lval, rval) do { \ - volatile type __lval; \ - \ - if (sizeof(type) >= sizeof(long double)) \ - (lval) = (rval); \ - else { \ - __lval = (rval); \ - (lval) = __lval; \ - } \ -} while (0) -#endif -#else -#define STRICT_ASSIGN(type, lval, rval) do { \ - volatile type __lval; \ - \ - if (sizeof(type) >= sizeof(long double)) \ - (lval) = (rval); \ - else { \ - __lval = (rval); \ - (lval) = __lval; \ - } \ -} while (0) -#endif /* FLT_EVAL_METHOD */ - -/* Support switching the mode to FP_PE if necessary. */ -#if defined(__i386__) && !defined(NO_FPSETPREC) -#define ENTERI() ENTERIT(long double) -#define ENTERIT(returntype) \ - returntype __retval; \ - fp_prec_t __oprec; \ - \ - if ((__oprec = fpgetprec()) != FP_PE) \ - fpsetprec(FP_PE) -#define RETURNI(x) do { \ - __retval = (x); \ - if (__oprec != FP_PE) \ - fpsetprec(__oprec); \ - RETURNF(__retval); \ -} while (0) -#define ENTERV() \ - fp_prec_t __oprec; \ - \ - if ((__oprec = fpgetprec()) != FP_PE) \ - fpsetprec(FP_PE) -#define RETURNV() do { \ - if (__oprec != FP_PE) \ - fpsetprec(__oprec); \ - return; \ -} while (0) -#else -#define ENTERI() -#define ENTERIT(x) -#define RETURNI(x) RETURNF(x) -#define ENTERV() -#define RETURNV() return -#endif - -/* Default return statement if hack*_t() is not used. */ -#define RETURNF(v) return (v) - -/* - * 2sum gives the same result as 2sumF without requiring |a| >= |b| or - * a == 0, but is slower. - */ -#define _2sum(a, b) do { \ - __typeof(a) __s, __w; \ - \ - __w = (a) + (b); \ - __s = __w - (a); \ - (b) = ((a) - (__w - __s)) + ((b) - __s); \ - (a) = __w; \ -} while (0) - -/* - * 2sumF algorithm. - * - * "Normalize" the terms in the infinite-precision expression a + b for - * the sum of 2 floating point values so that b is as small as possible - * relative to 'a'. (The resulting 'a' is the value of the expression in - * the same precision as 'a' and the resulting b is the rounding error.) - * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and - * exponent overflow or underflow must not occur. This uses a Theorem of - * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum" - * is apparently due to Skewchuk (1997). - * - * For this to always work, assignment of a + b to 'a' must not retain any - * extra precision in a + b. This is required by C standards but broken - * in many compilers. The brokenness cannot be worked around using - * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this - * algorithm would be destroyed by non-null strict assignments. (The - * compilers are correct to be broken -- the efficiency of all floating - * point code calculations would be destroyed similarly if they forced the - * conversions.) - * - * Fortunately, a case that works well can usually be arranged by building - * any extra precision into the type of 'a' -- 'a' should have type float_t, - * double_t or long double. b's type should be no larger than 'a's type. - * Callers should use these types with scopes as large as possible, to - * reduce their own extra-precision and efficiciency problems. In - * particular, they shouldn't convert back and forth just to call here. - */ -#ifdef DEBUG -#define _2sumF(a, b) do { \ - __typeof(a) __w; \ - volatile __typeof(a) __ia, __ib, __r, __vw; \ - \ - __ia = (a); \ - __ib = (b); \ - assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \ - \ - __w = (a) + (b); \ - (b) = ((a) - __w) + (b); \ - (a) = __w; \ - \ - /* The next 2 assertions are weak if (a) is already long double. */ \ - assert((long double)__ia + __ib == (long double)(a) + (b)); \ - __vw = __ia + __ib; \ - __r = __ia - __vw; \ - __r += __ib; \ - assert(__vw == (a) && __r == (b)); \ -} while (0) -#else /* !DEBUG */ -#define _2sumF(a, b) do { \ - __typeof(a) __w; \ - \ - __w = (a) + (b); \ - (b) = ((a) - __w) + (b); \ - (a) = __w; \ -} while (0) -#endif /* DEBUG */ - -/* - * Set x += c, where x is represented in extra precision as a + b. - * x must be sufficiently normalized and sufficiently larger than c, - * and the result is then sufficiently normalized. - * - * The details of ordering are that |a| must be >= |c| (so that (a, c) - * can be normalized without extra work to swap 'a' with c). The details of - * the normalization are that b must be small relative to the normalized 'a'. - * Normalization of (a, c) makes the normalized c tiny relative to the - * normalized a, so b remains small relative to 'a' in the result. However, - * b need not ever be tiny relative to 'a'. For example, b might be about - * 2**20 times smaller than 'a' to give about 20 extra bits of precision. - * That is usually enough, and adding c (which by normalization is about - * 2**53 times smaller than a) cannot change b significantly. However, - * cancellation of 'a' with c in normalization of (a, c) may reduce 'a' - * significantly relative to b. The caller must ensure that significant - * cancellation doesn't occur, either by having c of the same sign as 'a', - * or by having |c| a few percent smaller than |a|. Pre-normalization of - * (a, b) may help. - * - * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2 - * exercise 19). We gain considerable efficiency by requiring the terms to - * be sufficiently normalized and sufficiently increasing. - */ -#define _3sumF(a, b, c) do { \ - __typeof(a) __tmp; \ - \ - __tmp = (c); \ - _2sumF(__tmp, (a)); \ - (b) += (a); \ - (a) = __tmp; \ -} while (0) - -/* - * Common routine to process the arguments to nan(), nanf(), and nanl(). - */ -void _scan_nan(uint32_t *__words, int __num_words, const char *__s); - -/* - * Mix 0, 1 or 2 NaNs. First add 0 to each arg. This normally just turns - * signaling NaNs into quiet NaNs by setting a quiet bit. We do this - * because we want to never return a signaling NaN, and also because we - * don't want the quiet bit to affect the result. Then mix the converted - * args using the specified operation. - * - * When one arg is NaN, the result is typically that arg quieted. When both - * args are NaNs, the result is typically the quietening of the arg whose - * mantissa is largest after quietening. When neither arg is NaN, the - * result may be NaN because it is indeterminate, or finite for subsequent - * construction of a NaN as the indeterminate 0.0L/0.0L. - * - * Technical complications: the result in bits after rounding to the final - * precision might depend on the runtime precision and/or on compiler - * optimizations, especially when different register sets are used for - * different precisions. Try to make the result not depend on at least the - * runtime precision by always doing the main mixing step in long double - * precision. Try to reduce dependencies on optimizations by adding the - * the 0's in different precisions (unless everything is in long double - * precision). - */ -#define nan_mix(x, y) (nan_mix_op((x), (y), +)) -#define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0)) - -#ifdef _COMPLEX_H - -/* - * C99 specifies that complex numbers have the same representation as - * an array of two elements, where the first element is the real part - * and the second element is the imaginary part. - */ -typedef union { - float complex f; - float a[2]; -} float_complex; -typedef union { - double complex f; - double a[2]; -} double_complex; -typedef union { - long double complex f; - long double a[2]; -} long_double_complex; -#define REALPART(z) ((z).a[0]) -#define IMAGPART(z) ((z).a[1]) - -/* - * Inline functions that can be used to construct complex values. - * - * The C99 standard intends x+I*y to be used for this, but x+I*y is - * currently unusable in general since gcc introduces many overflow, - * underflow, sign and efficiency bugs by rewriting I*y as - * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. - * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted - * to -0.0+I*0.0. - * - * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL() - * to construct complex values. Compilers that conform to the C99 - * standard require the following functions to avoid the above issues. - */ - -#ifndef CMPLXF -static __inline float complex -CMPLXF(float x, float y) -{ - float_complex z; - - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); -} -#endif - -#ifndef CMPLX -static __inline double complex -CMPLX(double x, double y) -{ - double_complex z; - - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); -} -#endif - -#ifndef CMPLXL -static __inline long double complex -CMPLXL(long double x, long double y) -{ - long_double_complex z; - - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); -} -#endif - -#endif /* _COMPLEX_H */ - -#ifdef DEBUG -#if defined(__amd64__) || defined(__i386__) -#define breakpoint() asm("int $3") -#else -#include <signal.h> - -#define breakpoint() raise(SIGTRAP) -#endif -#endif - -/* Write a pari script to test things externally. */ -#ifdef DOPRINT -#include <stdio.h> - -#ifndef DOPRINT_SWIZZLE -#define DOPRINT_SWIZZLE 0 -#endif - -#ifdef DOPRINT_LD80 - -#define DOPRINT_START(xp) do { \ - uint64_t __lx; \ - uint16_t __hx; \ - \ - /* Hack to give more-problematic args. */ \ - EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \ - __lx ^= DOPRINT_SWIZZLE; \ - INSERT_LDBL80_WORDS(*xp, __hx, __lx); \ - printf("x = %.21Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#elif defined(DOPRINT_D64) - -#define DOPRINT_START(xp) do { \ - uint32_t __hx, __lx; \ - \ - EXTRACT_WORDS(__hx, __lx, *xp); \ - __lx ^= DOPRINT_SWIZZLE; \ - INSERT_WORDS(*xp, __hx, __lx); \ - printf("x = %.21Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#elif defined(DOPRINT_F32) - -#define DOPRINT_START(xp) do { \ - uint32_t __hx; \ - \ - GET_FLOAT_WORD(__hx, *xp); \ - __hx ^= DOPRINT_SWIZZLE; \ - SET_FLOAT_WORD(*xp, __hx); \ - printf("x = %.21Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */ - -#ifndef DOPRINT_SWIZZLE_HIGH -#define DOPRINT_SWIZZLE_HIGH 0 -#endif - -#define DOPRINT_START(xp) do { \ - uint64_t __lx, __llx; \ - uint16_t __hx; \ - \ - EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \ - __llx ^= DOPRINT_SWIZZLE; \ - __lx ^= DOPRINT_SWIZZLE_HIGH; \ - INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \ - printf("x = %.36Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#endif /* DOPRINT_LD80 */ - -#else /* !DOPRINT */ -#define DOPRINT_START(xp) -#define DOPRINT_END1(v) -#define DOPRINT_END2(hi, lo) -#endif /* DOPRINT */ - -#define RETURNP(x) do { \ - DOPRINT_END1(x); \ - RETURNF(x); \ -} while (0) -#define RETURNPI(x) do { \ - DOPRINT_END1(x); \ - RETURNI(x); \ -} while (0) -#define RETURN2P(x, y) do { \ - DOPRINT_END2((x), (y)); \ - RETURNF((x) + (y)); \ -} while (0) -#define RETURN2PI(x, y) do { \ - DOPRINT_END2((x), (y)); \ - RETURNI((x) + (y)); \ -} while (0) -#ifdef STRUCT_RETURN -#define RETURNSP(rp) do { \ - if (!(rp)->lo_set) \ - RETURNP((rp)->hi); \ - RETURN2P((rp)->hi, (rp)->lo); \ -} while (0) -#define RETURNSPI(rp) do { \ - if (!(rp)->lo_set) \ - RETURNPI((rp)->hi); \ - RETURN2PI((rp)->hi, (rp)->lo); \ -} while (0) -#endif -#define SUM2P(x, y) ({ \ - const __typeof (x) __x = (x); \ - const __typeof (y) __y = (y); \ - \ - DOPRINT_END2(__x, __y); \ - __x + __y; \ -}) - -/* - * ieee style elementary functions - * - * We rename functions here to improve other sources' diffability - * against fdlibm. - */ -#define __ieee754_sqrt sqrt -#define __ieee754_acos acos -#define __ieee754_acosh acosh -#define __ieee754_log log -#define __ieee754_log2 log2 -#define __ieee754_atanh atanh -#define __ieee754_asin asin -#define __ieee754_atan2 atan2 -#define __ieee754_exp exp -#define __ieee754_cosh cosh -#define __ieee754_fmod fmod -#define __ieee754_pow pow -#define __ieee754_lgamma lgamma -#define __ieee754_gamma gamma -#define __ieee754_lgamma_r lgamma_r -#define __ieee754_gamma_r gamma_r -#define __ieee754_log10 log10 -#define __ieee754_sinh sinh -#define __ieee754_hypot hypot -#define __ieee754_j0 j0 -#define __ieee754_j1 j1 -#define __ieee754_y0 y0 -#define __ieee754_y1 y1 -#define __ieee754_jn jn -#define __ieee754_yn yn -#define __ieee754_remainder remainder -#define __ieee754_scalb scalb -#define __ieee754_sqrtf sqrtf -#define __ieee754_acosf acosf -#define __ieee754_acoshf acoshf -#define __ieee754_logf logf -#define __ieee754_atanhf atanhf -#define __ieee754_asinf asinf -#define __ieee754_atan2f atan2f -#define __ieee754_expf expf -#define __ieee754_coshf coshf -#define __ieee754_fmodf fmodf -#define __ieee754_powf powf -#define __ieee754_lgammaf lgammaf -#define __ieee754_gammaf gammaf -#define __ieee754_lgammaf_r lgammaf_r -#define __ieee754_gammaf_r gammaf_r -#define __ieee754_log10f log10f -#define __ieee754_log2f log2f -#define __ieee754_sinhf sinhf -#define __ieee754_hypotf hypotf -#define __ieee754_j0f j0f -#define __ieee754_j1f j1f -#define __ieee754_y0f y0f -#define __ieee754_y1f y1f -#define __ieee754_jnf jnf -#define __ieee754_ynf ynf -#define __ieee754_remainderf remainderf -#define __ieee754_scalbf scalbf - -#define acos fdlibm::acos -#define asin fdlibm::asin -#define atan fdlibm::atan -#define atan2 fdlibm::atan2 -#define cosh fdlibm::cosh -#define sinh fdlibm::sinh -#define tanh fdlibm::tanh -#define exp fdlibm::exp -#define log fdlibm::log -#define log10 fdlibm::log10 -#define pow fdlibm::pow -#define ceil fdlibm::ceil -#define ceilf fdlibm::ceilf -#define fabs fdlibm::fabs -#define floor fdlibm::floor -#define acosh fdlibm::acosh -#define asinh fdlibm::asinh -#define atanh fdlibm::atanh -#define cbrt fdlibm::cbrt -#define expm1 fdlibm::expm1 -#define hypot fdlibm::hypot -#define log1p fdlibm::log1p -#define log2 fdlibm::log2 -#define scalb fdlibm::scalb -#define copysign fdlibm::copysign -#define scalbn fdlibm::scalbn -#define trunc fdlibm::trunc -#define truncf fdlibm::truncf -#define floorf fdlibm::floorf -#define nearbyint fdlibm::nearbyint -#define nearbyintf fdlibm::nearbyintf -#define rint fdlibm::rint -#define rintf fdlibm::rintf - -/* fdlibm kernel function */ -int __kernel_rem_pio2(double*,double*,int,int,int); - -/* double precision kernel functions */ -#ifndef INLINE_REM_PIO2 -int __ieee754_rem_pio2(double,double*); -#endif -double __kernel_sin(double,double,int); -double __kernel_cos(double,double); -double __kernel_tan(double,double,int); -double __ldexp_exp(double,int); -#ifdef _COMPLEX_H -double complex __ldexp_cexp(double complex,int); -#endif - -/* float precision kernel functions */ -#ifndef INLINE_REM_PIO2F -int __ieee754_rem_pio2f(float,double*); -#endif -#ifndef INLINE_KERNEL_SINDF -float __kernel_sindf(double); -#endif -#ifndef INLINE_KERNEL_COSDF -float __kernel_cosdf(double); -#endif -#ifndef INLINE_KERNEL_TANDF -float __kernel_tandf(double,int); -#endif -float __ldexp_expf(float,int); -#ifdef _COMPLEX_H -float complex __ldexp_cexpf(float complex,int); -#endif - -/* long double precision kernel functions */ -long double __kernel_sinl(long double, long double, int); -long double __kernel_cosl(long double, long double); -long double __kernel_tanl(long double, long double, int); - -#endif /* !_MATH_PRIVATE_H_ */ diff --git a/modules/fdlibm/src/moz.build b/modules/fdlibm/src/moz.build deleted file mode 100644 index 541f84039..000000000 --- a/modules/fdlibm/src/moz.build +++ /dev/null @@ -1,74 +0,0 @@ -# -*- Mode: python; indent-tabs-mode: nil; tab-width: 40 -*- -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. - -EXPORTS += [ - 'fdlibm.h', -] - -FINAL_LIBRARY = 'js' - -if CONFIG['CC_TYPE'] in ('clang', 'gcc'): - CXXFLAGS += [ - '-Wno-parentheses', - '-Wno-sign-compare', - ] - -if CONFIG['CC_TYPE'] == 'clang': - CXXFLAGS += [ - '-Wno-dangling-else', - ] - -if CONFIG['CC_TYPE'] in ('msvc', 'clang-cl'): - CXXFLAGS += [ - '-wd4146', # unary minus operator applied to unsigned type - '-wd4305', # truncation from 'double' to 'const float' - '-wd4723', # potential divide by 0 - '-wd4756', # overflow in constant arithmetic - ] - -if CONFIG['CC_TYPE'] == 'msvc': - CXXFLAGS += [ - '-wd4018', # signed/unsigned mismatch - ] - -if CONFIG['CC_TYPE'] == 'clang-cl': - CXXFLAGS += [ - '-Wno-sign-compare', # signed/unsigned mismatch - ] - -SOURCES += [ - 'e_acos.cpp', - 'e_acosh.cpp', - 'e_asin.cpp', - 'e_atan2.cpp', - 'e_atanh.cpp', - 'e_cosh.cpp', - 'e_exp.cpp', - 'e_hypot.cpp', - 'e_log.cpp', - 'e_log10.cpp', - 'e_log2.cpp', - 'e_pow.cpp', - 'e_sinh.cpp', - 'k_exp.cpp', - 's_asinh.cpp', - 's_atan.cpp', - 's_cbrt.cpp', - 's_ceil.cpp', - 's_ceilf.cpp', - 's_copysign.cpp', - 's_expm1.cpp', - 's_fabs.cpp', - 's_floor.cpp', - 's_floorf.cpp', - 's_log1p.cpp', - 's_nearbyint.cpp', - 's_rint.cpp', - 's_rintf.cpp', - 's_scalbn.cpp', - 's_tanh.cpp', - 's_trunc.cpp', - 's_truncf.cpp', -] diff --git a/modules/fdlibm/src/s_asinh.cpp b/modules/fdlibm/src/s_asinh.cpp deleted file mode 100644 index 7ecc396bb..000000000 --- a/modules/fdlibm/src/s_asinh.cpp +++ /dev/null @@ -1,58 +0,0 @@ -/* @(#)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$"); - -/* asinh(x) - * Method : - * Based on - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log(x)+ln2)) for large |x|, else - * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else - * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) - */ - -#include <cmath> -#include <float.h> - -#include "math_private.h" - -static const double -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ -huge= 1.00000000000000000000e+300; - -double -asinh(double x) -{ - double t,w; - int32_t hx,ix; - GET_HIGH_WORD(hx,x); - ix = hx&0x7fffffff; - if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ - if(ix< 0x3e300000) { /* |x|<2**-28 */ - if(huge+x>one) return x; /* return x inexact except 0 */ - } - if(ix>0x41b00000) { /* |x| > 2**28 */ - w = __ieee754_log(fabs(x))+ln2; - } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ - t = fabs(x); - w = __ieee754_log(2.0*t+one/(std::sqrt(x*x+one)+t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; - w =log1p(fabs(x)+t/(one+std::sqrt(one+t))); - } - if(hx>0) return w; else return -w; -} diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp deleted file mode 100644 index 21bc0d820..000000000 --- a/modules/fdlibm/src/s_atan.cpp +++ /dev/null @@ -1,119 +0,0 @@ -/* @(#)s_atan.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$"); - -/* atan(x) - * Method - * 1. Reduce x to positive by atan(x) = -atan(-x). - * 2. According to the integer k=4t+0.25 chopped, t=x, the argument - * is further reduced to one of the following intervals and the - * arctangent of t is evaluated by the corresponding formula: - * - * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) - * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) - * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) - * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) - * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <float.h> - -#include "math_private.h" - -static const double atanhi[] = { - 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ - 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ - 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ - 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ -}; - -static const double atanlo[] = { - 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ - 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ - 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ - 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ -}; - -static const double aT[] = { - 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ - -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ - 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ - -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ - 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ - -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ - 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ - -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ - 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ - -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ - 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ -}; - - static const double -one = 1.0, -huge = 1.0e300; - -double -atan(double x) -{ - double w,s1,s2,z; - int32_t ix,hx,id; - - GET_HIGH_WORD(hx,x); - ix = hx&0x7fffffff; - if(ix>=0x44100000) { /* if |x| >= 2^66 */ - u_int32_t low; - GET_LOW_WORD(low,x); - if(ix>0x7ff00000|| - (ix==0x7ff00000&&(low!=0))) - return x+x; /* NaN */ - if(hx>0) return atanhi[3]+*(volatile double *)&atanlo[3]; - else return -atanhi[3]-*(volatile double *)&atanlo[3]; - } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ - if (ix < 0x3e400000) { /* |x| < 2^-27 */ - if(huge+x>one) return x; /* raise inexact */ - } - id = -1; - } else { - x = fabs(x); - if (ix < 0x3ff30000) { /* |x| < 1.1875 */ - if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */ - id = 0; x = (2.0*x-one)/(2.0+x); - } else { /* 11/16<=|x|< 19/16 */ - id = 1; x = (x-one)/(x+one); - } - } else { - if (ix < 0x40038000) { /* |x| < 2.4375 */ - id = 2; x = (x-1.5)/(one+1.5*x); - } else { /* 2.4375 <= |x| < 2^66 */ - id = 3; x = -1.0/x; - } - }} - /* end of argument reduction */ - z = x*x; - w = z*z; - /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ - s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); - s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); - if (id<0) return x - x*(s1+s2); - else { - z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); - return (hx<0)? -z:z; - } -} diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp deleted file mode 100644 index fe3747e81..000000000 --- a/modules/fdlibm/src/s_cbrt.cpp +++ /dev/null @@ -1,113 +0,0 @@ -/* @(#)s_cbrt.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. - * ==================================================== - * - * Optimized by Bruce D. Evans. - */ - -//#include <sys/cdefs.h> -//__FBSDID("$FreeBSD$"); - -#include <float.h> -#include "math_private.h" - -/* cbrt(x) - * Return cube root of x - */ -static const u_int32_t - B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */ - B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */ - -/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */ -static const double -P0 = 1.87595182427177009643, /* 0x3ffe03e6, 0x0f61e692 */ -P1 = -1.88497979543377169875, /* 0xbffe28e0, 0x92f02420 */ -P2 = 1.621429720105354466140, /* 0x3ff9f160, 0x4a49d6c2 */ -P3 = -0.758397934778766047437, /* 0xbfe844cb, 0xbee751d9 */ -P4 = 0.145996192886612446982; /* 0x3fc2b000, 0xd4e4edd7 */ - -double -cbrt(double x) -{ - int32_t hx; - union { - double value; - uint64_t bits; - } u; - double r,s,t=0.0,w; - u_int32_t sign; - u_int32_t high,low; - - EXTRACT_WORDS(hx,low,x); - sign=hx&0x80000000; /* sign= sign(x) */ - hx ^=sign; - if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */ - - /* - * Rough cbrt to 5 bits: - * cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3) - * where e is integral and >= 0, m is real and in [0, 1), and "/" and - * "%" are integer division and modulus with rounding towards minus - * infinity. The RHS is always >= the LHS and has a maximum relative - * error of about 1 in 16. Adding a bias of -0.03306235651 to the - * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE - * floating point representation, for finite positive normal values, - * ordinary integer division of the value in bits magically gives - * almost exactly the RHS of the above provided we first subtract the - * exponent bias (1023 for doubles) and later add it back. We do the - * subtraction virtually to keep e >= 0 so that ordinary integer - * division rounds towards minus infinity; this is also efficient. - */ - if(hx<0x00100000) { /* zero or subnormal? */ - if((hx|low)==0) - return(x); /* cbrt(0) is itself */ - SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */ - t*=x; - GET_HIGH_WORD(high,t); - INSERT_WORDS(t,sign|((high&0x7fffffff)/3+B2),0); - } else - INSERT_WORDS(t,sign|(hx/3+B1),0); - - /* - * New cbrt to 23 bits: - * cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x) - * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r) - * to within 2**-23.5 when |r - 1| < 1/10. The rough approximation - * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this - * gives us bounds for r = t**3/x. - * - * Try to optimize for parallel evaluation as in k_tanf.c. - */ - r=(t*t)*(t/x); - t=t*((P0+r*(P1+r*P2))+((r*r)*r)*(P3+r*P4)); - - /* - * Round t away from zero to 23 bits (sloppily except for ensuring that - * the result is larger in magnitude than cbrt(x) but not much more than - * 2 23-bit ulps larger). With rounding towards zero, the error bound - * would be ~5/6 instead of ~4/6. With a maximum error of 2 23-bit ulps - * in the rounded t, the infinite-precision error in the Newton - * approximation barely affects third digit in the final error - * 0.667; the error in the rounded t can be up to about 3 23-bit ulps - * before the final error is larger than 0.667 ulps. - */ - u.value=t; - u.bits=(u.bits+0x80000000)&0xffffffffc0000000ULL; - t=u.value; - - /* one step Newton iteration to 53 bits with error < 0.667 ulps */ - s=t*t; /* t*t is exact */ - r=x/s; /* error <= 0.5 ulps; |r| < |t| */ - w=t+t; /* t+t is exact */ - r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */ - t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ - - return(t); -} diff --git a/modules/fdlibm/src/s_ceil.cpp b/modules/fdlibm/src/s_ceil.cpp deleted file mode 100644 index 67e9c1679..000000000 --- a/modules/fdlibm/src/s_ceil.cpp +++ /dev/null @@ -1,72 +0,0 @@ -/* @(#)s_ceil.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$"); - -/* - * ceil(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). - */ - -#include <float.h> - -#include "math_private.h" - -static const double huge = 1.0e300; - -double -ceil(double x) -{ - int32_t i0,i1,j0; - u_int32_t i,j; - EXTRACT_WORDS(i0,i1,x); - j0 = ((i0>>20)&0x7ff)-0x3ff; - if(j0<20) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;i1=0;} - else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;} - } - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x00100000)>>j0; - i0 &= (~i); i1=0; - } - } - } else if (j0>51) { - if(j0==0x400) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) { - if(j0==20) i0+=1; - else { - j = i1 + (1<<(52-j0)); - if(j<i1) i0+=1; /* got a carry */ - i1 = j; - } - } - i1 &= (~i); - } - } - INSERT_WORDS(x,i0,i1); - return x; -} diff --git a/modules/fdlibm/src/s_ceilf.cpp b/modules/fdlibm/src/s_ceilf.cpp deleted file mode 100644 index 7b52deeed..000000000 --- a/modules/fdlibm/src/s_ceilf.cpp +++ /dev/null @@ -1,51 +0,0 @@ -/* s_ceilf.c -- float version of s_ceil.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * 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$"); - -#include "math_private.h" - -static const float huge = 1.0e30; - -float -ceilf(float x) -{ - int32_t i0,j0; - u_int32_t i; - - GET_FLOAT_WORD(i0,x); - j0 = ((i0>>23)&0xff)-0x7f; - if(j0<23) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;} - else if(i0!=0) { i0=0x3f800000;} - } - } else { - i = (0x007fffff)>>j0; - if((i0&i)==0) return x; /* x is integral */ - if(huge+x>(float)0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x00800000)>>j0; - i0 &= (~i); - } - } - } else { - if(j0==0x80) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } - SET_FLOAT_WORD(x,i0); - return x; -} diff --git a/modules/fdlibm/src/s_copysign.cpp b/modules/fdlibm/src/s_copysign.cpp deleted file mode 100644 index b150106fb..000000000 --- a/modules/fdlibm/src/s_copysign.cpp +++ /dev/null @@ -1,32 +0,0 @@ -/* @(#)s_copysign.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$"); - -/* - * copysign(double x, double y) - * copysign(x,y) returns a value with the magnitude of x and - * with the sign bit of y. - */ - -#include "math_private.h" - -double -copysign(double x, double y) -{ - u_int32_t hx,hy; - GET_HIGH_WORD(hx,x); - GET_HIGH_WORD(hy,y); - SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000)); - return x; -} diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp deleted file mode 100644 index 90ebc1698..000000000 --- a/modules/fdlibm/src/s_expm1.cpp +++ /dev/null @@ -1,220 +0,0 @@ -/* @(#)s_expm1.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$"); - -/* expm1(x) - * Returns exp(x)-1, the exponential of x minus 1. - * - * Method - * 1. Argument reduction: - * Given x, find r and integer k such that - * - * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 - * - * Here a correction term c will be computed to compensate - * the error in r when rounded to a floating-point number. - * - * 2. Approximating expm1(r) by a special rational function on - * the interval [0,0.34658]: - * Since - * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ... - * we define R1(r*r) by - * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r) - * That is, - * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) - * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) - * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... - * We use a special Reme algorithm on [0,0.347] to generate - * a polynomial of degree 5 in r*r to approximate R1. The - * maximum error of this polynomial approximation is bounded - * by 2**-61. In other words, - * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 - * where Q1 = -1.6666666666666567384E-2, - * Q2 = 3.9682539681370365873E-4, - * Q3 = -9.9206344733435987357E-6, - * Q4 = 2.5051361420808517002E-7, - * Q5 = -6.2843505682382617102E-9; - * z = r*r, - * with error bounded by - * | 5 | -61 - * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 - * | | - * - * expm1(r) = exp(r)-1 is then computed by the following - * specific way which minimize the accumulation rounding error: - * 2 3 - * r r [ 3 - (R1 + R1*r/2) ] - * expm1(r) = r + --- + --- * [--------------------] - * 2 2 [ 6 - r*(3 - R1*r/2) ] - * - * To compensate the error in the argument reduction, we use - * expm1(r+c) = expm1(r) + c + expm1(r)*c - * ~ expm1(r) + c + r*c - * Thus c+r*c will be added in as the correction terms for - * expm1(r+c). Now rearrange the term to avoid optimization - * screw up: - * ( 2 2 ) - * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) - * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) - * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) - * ( ) - * - * = r - E - * 3. Scale back to obtain expm1(x): - * From step 1, we have - * expm1(x) = either 2^k*[expm1(r)+1] - 1 - * = or 2^k*[expm1(r) + (1-2^-k)] - * 4. Implementation notes: - * (A). To save one multiplication, we scale the coefficient Qi - * to Qi*2^i, and replace z by (x^2)/2. - * (B). To achieve maximum accuracy, we compute expm1(x) by - * (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) - * (ii) if k=0, return r-E - * (iii) if k=-1, return 0.5*(r-E)-0.5 - * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) - * else return 1.0+2.0*(r-E); - * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) - * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else - * (vii) return 2^k(1-((E+2^-k)-r)) - * - * Special cases: - * expm1(INF) is INF, expm1(NaN) is NaN; - * expm1(-INF) is -1, and - * for finite argument, only expm1(0)=0 is exact. - * - * Accuracy: - * according to an error analysis, the error is always less than - * 1 ulp (unit in the last place). - * - * Misc. info. - * For IEEE double - * if x > 7.09782712893383973096e+02 then expm1(x) overflow - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - -#include <float.h> - -#include "math_private.h" - -static const double -one = 1.0, -tiny = 1.0e-300, -o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */ -ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */ -ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */ -invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */ -/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */ -Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */ -Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ -Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ -Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ -Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ - -static volatile double huge = 1.0e+300; - -double -expm1(double x) -{ - double y,hi,lo,c,t,e,hxs,hfx,r1,twopk; - int32_t k,xsb; - u_int32_t hx; - - GET_HIGH_WORD(hx,x); - xsb = hx&0x80000000; /* sign bit of x */ - hx &= 0x7fffffff; /* high word of |x| */ - - /* filter out huge and non-finite argument */ - if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */ - if(hx >= 0x40862E42) { /* if |x|>=709.78... */ - if(hx>=0x7ff00000) { - u_int32_t low; - GET_LOW_WORD(low,x); - if(((hx&0xfffff)|low)!=0) - return x+x; /* NaN */ - else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - } - if(x > o_threshold) return huge*huge; /* overflow */ - } - if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */ - if(x+tiny<0.0) /* raise inexact */ - return tiny-one; /* return -1 */ - } - } - - /* argument reduction */ - if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ - if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ - if(xsb==0) - {hi = x - ln2_hi; lo = ln2_lo; k = 1;} - else - {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} - } else { - k = invln2*x+((xsb==0)?0.5:-0.5); - t = k; - hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ - lo = t*ln2_lo; - } - STRICT_ASSIGN(double, x, hi - lo); - c = (hi-x)-lo; - } - else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ - t = huge+x; /* return x with inexact flags when x!=0 */ - return x - (t-(huge+x)); - } - else k = 0; - - /* x is now in primary range */ - hfx = 0.5*x; - hxs = x*hfx; - r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); - t = 3.0-r1*hfx; - e = hxs*((r1-t)/(6.0 - x*t)); - if(k==0) return x - (x*e-hxs); /* c is 0 */ - else { - INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20,0); /* 2^k */ - e = (x*(e-c)-c); - e -= hxs; - if(k== -1) return 0.5*(x-e)-0.5; - if(k==1) { - if(x < -0.25) return -2.0*(e-(x+0.5)); - else return one+2.0*(x-e); - } - if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - y = one-(e-x); - if (k == 1024) { - double const_0x1p1023 = pow(2, 1023); - y = y*2.0*const_0x1p1023; - } - else y = y*twopk; - return y-one; - } - t = one; - if(k<20) { - SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ - y = t-(e-x); - y = y*twopk; - } else { - SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */ - y = x-(e+t); - y += one; - y = y*twopk; - } - } - return y; -} diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp deleted file mode 100644 index 6ca84d71b..000000000 --- a/modules/fdlibm/src/s_fabs.cpp +++ /dev/null @@ -1,29 +0,0 @@ -/* @(#)s_fabs.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$"); - -/* - * fabs(x) returns the absolute value of x. - */ - -#include "math_private.h" - -double -fabs(double x) -{ - u_int32_t high; - GET_HIGH_WORD(high,x); - SET_HIGH_WORD(x,high&0x7fffffff); - return x; -} diff --git a/modules/fdlibm/src/s_floor.cpp b/modules/fdlibm/src/s_floor.cpp deleted file mode 100644 index da57fc828..000000000 --- a/modules/fdlibm/src/s_floor.cpp +++ /dev/null @@ -1,73 +0,0 @@ -/* @(#)s_floor.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$"); - -/* - * floor(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). - */ - -#include <float.h> - -#include "math_private.h" - -static const double huge = 1.0e300; - -double -floor(double x) -{ - int32_t i0,i1,j0; - u_int32_t i,j; - EXTRACT_WORDS(i0,i1,x); - j0 = ((i0>>20)&0x7ff)-0x3ff; - if(j0<20) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=i1=0;} - else if(((i0&0x7fffffff)|i1)!=0) - { i0=0xbff00000;i1=0;} - } - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x00100000)>>j0; - i0 &= (~i); i1=0; - } - } - } else if (j0>51) { - if(j0==0x400) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) { - if(j0==20) i0+=1; - else { - j = i1+(1<<(52-j0)); - if(j<i1) i0 +=1 ; /* got a carry */ - i1=j; - } - } - i1 &= (~i); - } - } - INSERT_WORDS(x,i0,i1); - return x; -} diff --git a/modules/fdlibm/src/s_floorf.cpp b/modules/fdlibm/src/s_floorf.cpp deleted file mode 100644 index 88511f209..000000000 --- a/modules/fdlibm/src/s_floorf.cpp +++ /dev/null @@ -1,60 +0,0 @@ -/* s_floorf.c -- float version of s_floor.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * 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$"); - -/* - * floorf(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floorf(x). - */ - -#include "math_private.h" - -static const float huge = 1.0e30; - -float -floorf(float x) -{ - int32_t i0,j0; - u_int32_t i; - GET_FLOAT_WORD(i0,x); - j0 = ((i0>>23)&0xff)-0x7f; - if(j0<23) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=0;} - else if((i0&0x7fffffff)!=0) - { i0=0xbf800000;} - } - } else { - i = (0x007fffff)>>j0; - if((i0&i)==0) return x; /* x is integral */ - if(huge+x>(float)0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x00800000)>>j0; - i0 &= (~i); - } - } - } else { - if(j0==0x80) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } - SET_FLOAT_WORD(x,i0); - return x; -} diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp deleted file mode 100644 index afc6919c6..000000000 --- a/modules/fdlibm/src/s_log1p.cpp +++ /dev/null @@ -1,175 +0,0 @@ -/* @(#)s_log1p.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$"); - -/* double log1p(double x) - * - * Method : - * 1. Argument Reduction: find k and f such that - * 1+x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * Note. If k=0, then f=x is exact. However, if k!=0, then f - * may not be representable exactly. In that case, a correction - * term is need. Let u=1+x rounded. Let c = (1+x)-u, then - * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u), - * and add back the correction term c/u. - * (Note: when x > 2**53, one can simply return log(x)) - * - * 2. Approximation of log1p(f). - * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R - * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error - * of this polynomial approximation is bounded by 2**-58.45. In - * other words, - * 2 4 6 8 10 12 14 - * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s - * (the values of Lp1 to Lp7 are listed in the program) - * and - * | 2 14 | -58.45 - * | Lp1*s +...+Lp7*s - R(z) | <= 2 - * | | - * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. - * In order to guarantee error in log below 1ulp, we compute log - * by - * log1p(f) = f - (hfsq - s*(hfsq+R)). - * - * 3. Finally, log1p(x) = k*ln2 + log1p(f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) - * Here ln2 is split into two floating point number: - * ln2_hi + ln2_lo, - * where n*ln2_hi is always exact for |n| < 2000. - * - * Special cases: - * log1p(x) is NaN with signal if x < -1 (including -INF) ; - * log1p(+INF) is +INF; log1p(-1) is -INF with signal; - * log1p(NaN) is that NaN with no signal. - * - * Accuracy: - * according to an error analysis, the error is always less than - * 1 ulp (unit in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - * - * Note: Assuming log() return accurate answer, the following - * algorithm can be used to compute log1p(x) to within a few ULP: - * - * u = 1+x; - * if(u==1.0) return x ; else - * return log(u)*(x/(u-1.0)); - * - * See HP-15C Advanced Functions Handbook, p.193. - */ - -#include <float.h> - -#include "math_private.h" - -static const double -ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ -ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ -two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ -Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ -Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ -Lp3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ -Lp4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ -Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ -Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ -Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ - -static const double zero = 0.0; -static volatile double vzero = 0.0; - -double -log1p(double x) -{ - double hfsq,f,c,s,z,R,u; - int32_t k,hx,hu,ax; - - GET_HIGH_WORD(hx,x); - ax = hx&0x7fffffff; - - k = 1; - if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */ - if(ax>=0x3ff00000) { /* x <= -1.0 */ - if(x==-1.0) return -two54/vzero; /* log1p(-1)=+inf */ - else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ - } - if(ax<0x3e200000) { /* |x| < 2**-29 */ - if(two54+x>zero /* raise inexact */ - &&ax<0x3c900000) /* |x| < 2**-54 */ - return x; - else - return x - x*x*0.5; - } - if(hx>0||hx<=((int32_t)0xbfd2bec4)) { - k=0;f=x;hu=1;} /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ - } - if (hx >= 0x7ff00000) return x+x; - if(k!=0) { - if(hx<0x43400000) { - STRICT_ASSIGN(double,u,1.0+x); - GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; - c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ - c /= u; - } else { - u = x; - GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; - c = 0; - } - hu &= 0x000fffff; - /* - * The approximation to sqrt(2) used in thresholds is not - * critical. However, the ones used above must give less - * strict bounds than the one here so that the k==0 case is - * never reached from here, since here we have committed to - * using the correction term but don't use it if k==0. - */ - if(hu<0x6a09e) { /* u ~< sqrt(2) */ - SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ - } else { - k += 1; - SET_HIGH_WORD(u,hu|0x3fe00000); /* normalize u/2 */ - hu = (0x00100000-hu)>>2; - } - f = u-1.0; - } - hfsq=0.5*f*f; - if(hu==0) { /* |f| < 2**-20 */ - if(f==zero) { - if(k==0) { - return zero; - } else { - c += k*ln2_lo; - return k*ln2_hi+c; - } - } - R = hfsq*(1.0-0.66666666666666666*f); - if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); - } - s = f/(2.0+f); - z = s*s; - R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f); -} diff --git a/modules/fdlibm/src/s_nearbyint.cpp b/modules/fdlibm/src/s_nearbyint.cpp deleted file mode 100644 index 6c04212d3..000000000 --- a/modules/fdlibm/src/s_nearbyint.cpp +++ /dev/null @@ -1,60 +0,0 @@ -/*- - * SPDX-License-Identifier: BSD-2-Clause-FreeBSD - * - * Copyright (c) 2004 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 <fenv.h> -#include "math_private.h" - -/* - * We save and restore the floating-point environment to avoid raising - * an inexact exception. We can get away with using fesetenv() - * instead of feclearexcept()/feupdateenv() to restore the environment - * because the only exception defined for rint() is overflow, and - * rounding can't overflow as long as emax >= p. - * - * The volatile keyword is needed below because clang incorrectly assumes - * that rint won't raise any floating-point exceptions. Declaring ret volatile - * is sufficient to trick the compiler into doing the right thing. - */ -#define DECL(type, fn, rint) \ -type \ -fn(type x) \ -{ \ - volatile type ret; \ - fenv_t env; \ - \ - fegetenv(&env); \ - ret = rint(x); \ - fesetenv(&env); \ - return (ret); \ -} - -DECL(double, nearbyint, rint) -DECL(float, nearbyintf, rintf) diff --git a/modules/fdlibm/src/s_rint.cpp b/modules/fdlibm/src/s_rint.cpp deleted file mode 100644 index 19171f87f..000000000 --- a/modules/fdlibm/src/s_rint.cpp +++ /dev/null @@ -1,87 +0,0 @@ -/* @(#)s_rint.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$"); - -/* - * rint(x) - * Return x rounded to integral value according to the prevailing - * rounding mode. - * Method: - * Using floating addition. - * Exception: - * Inexact flag raised if x not equal to rint(x). - */ - -#include <float.h> - -#include "math_private.h" - -static const double -TWO52[2]={ - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ -}; - -double -rint(double x) -{ - int32_t i0,j0,sx; - u_int32_t i,i1; - double w,t; - EXTRACT_WORDS(i0,i1,x); - sx = (i0>>31)&1; - j0 = ((i0>>20)&0x7ff)-0x3ff; - if(j0<20) { - if(j0<0) { - if(((i0&0x7fffffff)|i1)==0) return x; - i1 |= (i0&0x0fffff); - i0 &= 0xfffe0000; - i0 |= ((i1|-i1)>>12)&0x80000; - SET_HIGH_WORD(x,i0); - STRICT_ASSIGN(double,w,TWO52[sx]+x); - t = w-TWO52[sx]; - GET_HIGH_WORD(i0,t); - SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); - return t; - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - i>>=1; - if(((i0&i)|i1)!=0) { - /* - * Some bit is set after the 0.5 bit. To avoid the - * possibility of errors from double rounding in - * w = TWO52[sx]+x, adjust the 0.25 bit to a lower - * guard bit. We do this for all j0<=51. The - * adjustment is trickiest for j0==18 and j0==19 - * since then it spans the word boundary. - */ - if(j0==19) i1 = 0x40000000; else - if(j0==18) i1 = 0x80000000; else - i0 = (i0&(~i))|((0x20000)>>j0); - } - } - } else if (j0>51) { - if(j0==0x400) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ - i>>=1; - if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); - } - INSERT_WORDS(x,i0,i1); - STRICT_ASSIGN(double,w,TWO52[sx]+x); - return w-TWO52[sx]; -} diff --git a/modules/fdlibm/src/s_rintf.cpp b/modules/fdlibm/src/s_rintf.cpp deleted file mode 100644 index 3a729b005..000000000 --- a/modules/fdlibm/src/s_rintf.cpp +++ /dev/null @@ -1,52 +0,0 @@ -/* s_rintf.c -- float version of s_rint.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * 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$"); - -#include <float.h> -#include <stdint.h> - -#include "math_private.h" - -static const float -TWO23[2]={ - 8.3886080000e+06, /* 0x4b000000 */ - -8.3886080000e+06, /* 0xcb000000 */ -}; - -float -rintf(float x) -{ - int32_t i0,j0,sx; - float w,t; - GET_FLOAT_WORD(i0,x); - sx = (i0>>31)&1; - j0 = ((i0>>23)&0xff)-0x7f; - if(j0<23) { - if(j0<0) { - if((i0&0x7fffffff)==0) return x; - STRICT_ASSIGN(float,w,TWO23[sx]+x); - t = w-TWO23[sx]; - GET_FLOAT_WORD(i0,t); - SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31)); - return t; - } - STRICT_ASSIGN(float,w,TWO23[sx]+x); - return w-TWO23[sx]; - } - if(j0==0x80) return x+x; /* inf or NaN */ - else return x; /* x is integral */ -} diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp deleted file mode 100644 index dfbcf5c57..000000000 --- a/modules/fdlibm/src/s_scalbn.cpp +++ /dev/null @@ -1,60 +0,0 @@ -/* @(#)s_scalbn.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$"); - -/* - * scalbn (double x, int n) - * scalbn(x,n) returns x* 2**n computed by exponent - * manipulation rather than by actually performing an - * exponentiation or a multiplication. - */ - -#include <float.h> - -#include "math_private.h" - -static const double -two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ -huge = 1.0e+300, -tiny = 1.0e-300; - -double -scalbn (double x, int n) -{ - int32_t k,hx,lx; - EXTRACT_WORDS(hx,lx,x); - k = (hx&0x7ff00000)>>20; /* extract exponent */ - if (k==0) { /* 0 or subnormal x */ - if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ - x *= two54; - GET_HIGH_WORD(hx,x); - k = ((hx&0x7ff00000)>>20) - 54; - if (n< -50000) return tiny*x; /*underflow*/ - } - if (k==0x7ff) return x+x; /* NaN or Inf */ - k = k+n; - if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ - if (k > 0) /* normal result */ - {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} - if (k <= -54) { - if (n > 50000) /* in case integer overflow in n+k */ - return huge*copysign(huge,x); /*overflow*/ - else - return tiny*copysign(tiny,x); /*underflow*/ - } - k += 54; /* subnormal result */ - SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); - return x*twom54; -} diff --git a/modules/fdlibm/src/s_tanh.cpp b/modules/fdlibm/src/s_tanh.cpp deleted file mode 100644 index 238973fce..000000000 --- a/modules/fdlibm/src/s_tanh.cpp +++ /dev/null @@ -1,79 +0,0 @@ -/* @(#)s_tanh.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$"); - -/* Tanh(x) - * Return the Hyperbolic Tangent of x - * - * Method : - * x -x - * e - e - * 0. tanh(x) is defined to be ----------- - * x -x - * e + e - * 1. reduce x to non-negative by tanh(-x) = -tanh(x). - * 2. 0 <= x < 2**-28 : tanh(x) := x with inexact if x != 0 - * -t - * 2**-28 <= x < 1 : tanh(x) := -----; t = expm1(-2x) - * t + 2 - * 2 - * 1 <= x < 22 : tanh(x) := 1 - -----; t = expm1(2x) - * t + 2 - * 22 <= x <= INF : tanh(x) := 1. - * - * Special cases: - * tanh(NaN) is NaN; - * only tanh(0)=0 is exact for finite argument. - */ - -#include <float.h> - -#include "math_private.h" - -static const volatile double tiny = 1.0e-300; -static const double one = 1.0, two = 2.0, huge = 1.0e300; - -double -tanh(double x) -{ - double t,z; - int32_t jx,ix; - - GET_HIGH_WORD(jx,x); - ix = jx&0x7fffffff; - - /* x is INF or NaN */ - if(ix>=0x7ff00000) { - if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */ - else return one/x-one; /* tanh(NaN) = NaN */ - } - - /* |x| < 22 */ - if (ix < 0x40360000) { /* |x|<22 */ - if (ix<0x3e300000) { /* |x|<2**-28 */ - if(huge+x>one) return x; /* tanh(tiny) = tiny with inexact */ - } - if (ix>=0x3ff00000) { /* |x|>=1 */ - t = expm1(two*fabs(x)); - z = one - two/(t+two); - } else { - t = expm1(-two*fabs(x)); - z= -t/(t+two); - } - /* |x| >= 22, return +-1 */ - } else { - z = one - tiny; /* raise inexact flag */ - } - return (jx>=0)? z: -z; -} diff --git a/modules/fdlibm/src/s_trunc.cpp b/modules/fdlibm/src/s_trunc.cpp deleted file mode 100644 index d2294a272..000000000 --- a/modules/fdlibm/src/s_trunc.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/* @(#)s_floor.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$"); - -/* - * trunc(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to trunc(x). - */ - -#include <float.h> - -#include "math_private.h" - -static const double huge = 1.0e300; - -double -trunc(double x) -{ - int32_t i0,i1,j0; - u_int32_t i; - EXTRACT_WORDS(i0,i1,x); - j0 = ((i0>>20)&0x7ff)-0x3ff; - if(j0<20) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* |x|<1, so return 0*sign(x) */ - i0 &= 0x80000000U; - i1 = 0; - } - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - i0 &= (~i); i1=0; - } - } - } else if (j0>51) { - if(j0==0x400) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) /* raise inexact flag */ - i1 &= (~i); - } - INSERT_WORDS(x,i0,i1); - return x; -} diff --git a/modules/fdlibm/src/s_truncf.cpp b/modules/fdlibm/src/s_truncf.cpp deleted file mode 100644 index 4853a4450..000000000 --- a/modules/fdlibm/src/s_truncf.cpp +++ /dev/null @@ -1,52 +0,0 @@ -/* @(#)s_floor.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$"); - -/* - * truncf(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to truncf(x). - */ - -#include "math_private.h" - -static const float huge = 1.0e30F; - -float -truncf(float x) -{ - int32_t i0,j0; - u_int32_t i; - GET_FLOAT_WORD(i0,x); - j0 = ((i0>>23)&0xff)-0x7f; - if(j0<23) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0F) /* |x|<1, so return 0*sign(x) */ - i0 &= 0x80000000; - } else { - i = (0x007fffff)>>j0; - if((i0&i)==0) return x; /* x is integral */ - if(huge+x>0.0F) /* raise inexact flag */ - i0 &= (~i); - } - } else { - if(j0==0x80) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } - SET_FLOAT_WORD(x,i0); - return x; -} diff --git a/modules/fdlibm/update.sh b/modules/fdlibm/update.sh deleted file mode 100644 index 90d417862..000000000 --- a/modules/fdlibm/update.sh +++ /dev/null @@ -1,40 +0,0 @@ -#!/bin/sh - -# Script to update the mozilla in-tree copy of the fdlibm library. -# Run this within the /modules/fdlibm directory of the source tree. - -set -e - -API_BASE_URL=https://api.github.com/repos/freebsd/freebsd - -get_commit() { - curl -s "${API_BASE_URL}/commits?path=lib/msun/src&per_page=1" \ - | python -c 'import json, sys; print(json.loads(sys.stdin.read())[0]["sha"])' -} -get_date() { - curl -s "${API_BASE_URL}/commits/${COMMIT}" \ - | python -c 'import json, sys; print(json.loads(sys.stdin.read())["commit"]["committer"]["date"])' -} - -mv ./src/moz.build ./src_moz.build -rm -rf src -if [ "$#" -eq 0 ]; then - COMMIT=$(get_commit) -else - COMMIT="$1" -fi -sh ./import.sh "${COMMIT}" -mv ./src_moz.build ./src/moz.build -COMMITDATE=$(get_date) -for FILE in $(ls patches/*.patch | sort); do - echo "Applying ${FILE} ..." - patch -p3 --no-backup-if-mismatch < ${FILE} -done -hg add src - -perl -p -i -e "s/\[commit [0-9a-f]{40} \(.{1,100}\)\]/[commit ${COMMIT} (${COMMITDATE})]/" README.mozilla - -echo "###" -echo "### Updated fdlibm/src to ${COMMIT}." -echo "### Remember to verify and commit the changes to source control!" -echo "###" |