RandomHuman
/
Mypal
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Update fdlibm to September 2019 version.

This commit is contained in:
Fedor 2019-12-25 15:46:14 +03:00
parent c5b70ecd3d
commit d9300b2bb0
38 changed files with 745 additions and 716 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
modules/fdlibm/README.mozilla
View File

@ -9,9 +9,11 @@ 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 f2287da07ac7a26ac08745cac66eec82ab9ba384].
Current version: [commit cf4707bb2f78ecf56ba350bdc24e3135b4339622 (2019-09-25T18:50:57Z)].
patches 01-14 fixes files to be usable within mozilla-central tree.
patches 01-18 fixes files to be usable within mozilla-central tree.
See https://bugzilla.mozilla.org/show_bug.cgi?id=933257

3
modules/fdlibm/import.sh
View File

@ -2,7 +2,7 @@
set -e
BASE_URL=https://raw.githubusercontent.com/freebsd/freebsd/master/lib/msun/src
BASE_URL=https://raw.githubusercontent.com/freebsd/freebsd/"${1}"/lib/msun/src
download_source() {
REMOTE_FILENAME=$1
@ -105,7 +105,6 @@ 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 e_sqrt.c e_sqrt.cpp
download_source s_nearbyint.c s_nearbyint.cpp
download_source s_rint.c s_rint.cpp

19
modules/fdlibm/patches/01_remove_unused_declarations_from_fdlibm_h.patch
View File

@ -1,7 +1,7 @@
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,496 +12,50 @@
@@ -12,499 +12,49 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
* $FreeBSD$
@ -242,15 +242,15 @@ diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h
-double modf(double, double *); /* fundamentally !__pure2 */
double pow(double, double);
double sqrt(double);
-double sqrt(double);
+double fabs(double);
+double floor(double);
+double trunc(double);
double ceil(double);
-double ceil(double);
-double fabs(double) __pure2;
-double floor(double);
double floor(double);
-double fmod(double, double);
+double trunc(double);
+double ceil(double);
-/*
- * These functions are not in C90.
@ -368,13 +368,13 @@ diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h
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);
@ -497,6 +497,9 @@ diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h
-
-#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

2
modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch
View File

@ -22,7 +22,7 @@ diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h
double cosh(double);
double sinh(double);
@@ -53,9 +53,9 @@ double scalbn(double, int);
@@ -52,9 +52,9 @@ double scalbn(double, int);
float ceilf(float);
float floorf(float);

2
modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch
View File

@ -20,7 +20,7 @@ diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h
double cosh(double);
double sinh(double);
double tanh(double);
@@ -53,9 +55,11 @@ double scalbn(double, int);
@@ -52,9 +54,11 @@ double scalbn(double, int);
float ceilf(float);
float floorf(float);

53
modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch
View File

@ -232,15 +232,15 @@ diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp
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,17 +52,16 @@
*
@@ -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
* 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"
@ -249,7 +249,7 @@ diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
zero = 0.0,
one = 1.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
@ -271,31 +271,10 @@ diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp
__ieee754_sinh(double x)
{
double t,h;
diff --git a/modules/fdlibm/src/e_sqrt.cpp b/modules/fdlibm/src/e_sqrt.cpp
--- a/modules/fdlibm/src/e_sqrt.cpp
+++ b/modules/fdlibm/src/e_sqrt.cpp
@@ -81,17 +81,16 @@
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
*
* Other methods : see the appended file at the end of the program below.
*---------------
*/
#include <float.h>
-#include "math.h"
#include "math_private.h"
static const double one = 1.0, tiny=1.0e-300;
double
__ieee754_sqrt(double x)
{
double z;
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,17 +24,16 @@
@@ -26,17 +26,16 @@
* SUCH DAMAGE.
*/
@ -380,8 +359,7 @@ diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp
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
@@ -10,17 +10,16 @@
* ====================================================
@@ -11,17 +11,16 @@
*
* Optimized by Bruce D. Evans.
*/
@ -389,6 +367,7 @@ diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <float.h>
-#include "math.h"
#include "math_private.h"
@ -485,10 +464,10 @@ diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp
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
@@ -13,17 +13,16 @@
#ifndef lint
static char rcsid[] = "$FreeBSD$";
#endif
@@ -12,17 +12,16 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* fabs(x) returns the absolute value of x.
@ -569,7 +548,7 @@ diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp
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
@@ -23,17 +23,17 @@
@@ -25,17 +25,17 @@
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
@ -633,13 +612,13 @@ diff --git a/modules/fdlibm/src/s_rintf.cpp b/modules/fdlibm/src/s_rintf.cpp
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
@@ -19,17 +19,16 @@ static char rcsid[] = "$FreeBSD$";
@@ -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 <sys/cdefs.h>
#include <float.h>
-#include "math.h"

63
modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch
View File

@ -1,7 +1,7 @@
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,20 +14,21 @@
@@ -14,52 +14,38 @@
* $FreeBSD$
*/
@ -24,15 +24,16 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
* 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.
@@ -36,27 +37,17 @@
* 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.
*/
-/*
- * 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
@ -43,7 +44,53 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
-#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
@ -53,7 +100,7 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
{
u_int32_t msw;
u_int32_t lsw;
@@ -64,17 +55,17 @@ typedef union
@@ -105,17 +96,17 @@ typedef union
struct
{
u_int64_t w;

3
modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch
View File

@ -1,7 +1,7 @@
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
@@ -724,16 +724,51 @@ irintl(long double x)
@@ -872,16 +872,50 @@ irintl(long double x)
#define __ieee754_j1f j1f
#define __ieee754_y0f y0f
#define __ieee754_y1f y1f
@ -21,7 +21,6 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
+#define log fdlibm::log
+#define log10 fdlibm::log10
+#define pow fdlibm::pow
+#define sqrt fdlibm::sqrt
+#define ceil fdlibm::ceil
+#define ceilf fdlibm::ceilf
+#define fabs fdlibm::fabs

48
modules/fdlibm/patches/08_remove_weak_reference_macro.patch
View File

@ -174,6 +174,22 @@ diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp
-#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
@ -190,30 +206,6 @@ diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp
-#if (LDBL_MANT_DIG == 53)
-__weak_reference(sinh, sinhl);
-#endif
diff --git a/modules/fdlibm/src/e_sqrt.cpp b/modules/fdlibm/src/e_sqrt.cpp
--- a/modules/fdlibm/src/e_sqrt.cpp
+++ b/modules/fdlibm/src/e_sqrt.cpp
@@ -182,20 +182,16 @@ double
ix0 = (q>>1)+0x3fe00000;
ix1 = q1>>1;
if ((q&1)==1) ix1 |= sign;
ix0 += (m <<20);
INSERT_WORDS(z,ix0,ix1);
return z;
}
-#if (LDBL_MANT_DIG == 53)
-__weak_reference(sqrt, sqrtl);
-#endif
-
/*
Other methods (use floating-point arithmetic)
-------------
(This is a copy of a drafted paper by Prof W. Kahan
and K.C. Ng, written in May, 1986)
Two algorithms are given here to implement sqrt(x)
(IEEE double precision arithmetic) in software.
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
@ -249,7 +241,7 @@ diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp
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
@@ -105,12 +105,8 @@ cbrt(double x)
@@ -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 */
@ -346,10 +338,10 @@ 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)
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*/
else
return tiny*copysign(tiny,x); /*underflow*/
}
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;

71
modules/fdlibm/patches/09_comment_out_rcsid_variable.patch
View File

@ -53,7 +53,7 @@ diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@ -284,7 +284,7 @@ diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp
* 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
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@ -325,34 +325,10 @@ diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp
* 2.
* E + E/(E+1)
* 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
diff --git a/modules/fdlibm/src/e_sqrt.cpp b/modules/fdlibm/src/e_sqrt.cpp
--- a/modules/fdlibm/src/e_sqrt.cpp
+++ b/modules/fdlibm/src/e_sqrt.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_sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
* | Use the hardware sqrt if you have one |
* ------------------------------------------
* Method:
* Bit by bit method using integer arithmetic. (Slow, but portable)
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
@@ -19,22 +19,22 @@
@@ -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
@ -467,13 +443,13 @@ diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp
+//#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 */
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
@ -573,12 +549,7 @@ diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp
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
@@ -1,22 +1,22 @@
-/* @(#)s_fabs.c 5.1 93/09/24 */
+ /* @(#)s_fabs.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -5,18 +5,18 @@
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
@ -587,10 +558,10 @@ diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp
* ====================================================
*/
#ifndef lint
-static char rcsid[] = "$FreeBSD$";
+ //static char rcsid[] = "$FreeBSD$";
#endif
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
/*
* fabs(x) returns the absolute value of x.
@ -598,6 +569,7 @@ diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp
#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
@ -673,7 +645,7 @@ diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp
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
@@ -19,18 +19,18 @@
@@ -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
@ -745,7 +717,8 @@ diff --git a/modules/fdlibm/src/s_rintf.cpp b/modules/fdlibm/src/s_rintf.cpp
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
@@ -6,27 +6,27 @@
@@ -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
@ -753,10 +726,10 @@ diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp
* ====================================================
*/
#ifndef lint
-static char rcsid[] = "$FreeBSD$";
+//static char rcsid[] = "$FreeBSD$";
#endif
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
/*
* scalbn (double x, int n)
@ -765,16 +738,6 @@ diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp
* exponentiation or a multiplication.
*/
-#include <sys/cdefs.h>
+//#include <sys/cdefs.h>
#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,
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

4
modules/fdlibm/patches/10_remove_unused_function_from_k_exp_cpp.patch
View File

@ -1,7 +1,7 @@
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
@@ -22,18 +22,16 @@
@@ -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.
@ -20,7 +20,7 @@ diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp
/*
* Compute exp(x), scaled to avoid spurious overflow. An exponent is
* returned separately in 'expt'.
@@ -76,32 +74,8 @@ double
@@ -78,32 +76,8 @@ double
double exp_x, scale;
int ex_expt;

10
modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch
View File

@ -1,7 +1,7 @@
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,21 @@
@@ -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.
@ -17,11 +17,11 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
+#define u_int64_t uint64_t
+#endif
+
/*
* A union which permits us to convert between a double and two 32 bit
* ints.
*/
/* 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;

8
modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch
View File

@ -1,7 +1,7 @@
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
@@ -285,16 +285,27 @@ do { \
@@ -328,16 +328,27 @@ do { \
if (sizeof(type) >= sizeof(long double)) \
(lval) = (rval); \
else { \
@ -25,7 +25,7 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
/* Support switching the mode to FP_PE if necessary. */
#if defined(__i386__) && !defined(NO_FPSETPREC)
#define ENTERI() \
long double __retval; \
#define ENTERI() ENTERIT(long double)
#define ENTERIT(returntype) \
returntype __retval; \
fp_prec_t __oprec; \
\

4
modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch
View File

@ -3,9 +3,9 @@ diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp
+++ b/modules/fdlibm/src/e_exp.cpp
@@ -146,14 +146,17 @@ double
if(k >= -1021)
INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0);
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0);
else
INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0);
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);

2
modules/fdlibm/patches/15_remove_unused_rintl_function_from_s_nearbyint_cpp.patch
View File

@ -1,7 +1,7 @@
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
@@ -51,9 +51,8 @@ fn(type x) \
@@ -53,9 +53,8 @@ fn(type x) \
fegetenv(&env); \
ret = rint(x); \
fesetenv(&env); \

3
modules/fdlibm/patches/15_use_safer_strict_assign_on_visual_studio.patch → modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch
View File

@ -1,7 +1,7 @@
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
@@ -271,17 +271,17 @@ do { \
@@ -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
@ -20,4 +20,3 @@ diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private
if (sizeof(type) >= sizeof(long double)) \
(lval) = (rval); \
else { \

40
modules/fdlibm/patches/17_exp_exact_result_for_positive_one.patch Normal file
View File

@ -0,0 +1,40 @@
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);

255
modules/fdlibm/patches/18_use_stdlib_sqrt.patch Normal file
View File

@ -0,0 +1,255 @@
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;
}

130
modules/fdlibm/patches/19_remove_unneeded_round_to_integer_helpers.patch Normal file
View File

@ -0,0 +1,130 @@
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

5
modules/fdlibm/src/e_acos.cpp
View File

@ -38,6 +38,7 @@
* Function needed: sqrt
*/
#include <cmath>
#include <float.h>
#include "math_private.h"
@ -87,13 +88,13 @@ __ieee754_acos(double x)
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);

5
modules/fdlibm/src/e_acosh.cpp
View File

@ -29,6 +29,7 @@
* acosh(NaN) is NaN without signal.
*/
#include <cmath>
#include <float.h>
#include "math_private.h"
@ -55,9 +56,9 @@ __ieee754_acosh(double x)
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));
}
}

5
modules/fdlibm/src/e_asin.cpp
View File

@ -6,7 +6,7 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@ -44,6 +44,7 @@
*
*/
#include <cmath>
#include <float.h>
#include "math_private.h"
@ -95,7 +96,7 @@ __ieee754_asin(double 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);

4
modules/fdlibm/src/e_atan2.cpp
View File

@ -69,8 +69,8 @@ __ieee754_atan2(double y, double x)
iy = hy&0x7fffffff;
if(((ix|((lx|-lx)>>31))>0x7ff00000)||
((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
return x+y;
if((hx-0x3ff00000|lx)==0) return atan(y); /* x=1.0 */
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 */

7
modules/fdlibm/src/e_exp.cpp
View File

@ -96,6 +96,8 @@ 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*/
@ -127,6 +129,7 @@ __ieee754_exp(double x) /* default IEEE double exp */
/* 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]);
@ -144,9 +147,9 @@ __ieee754_exp(double x) /* default IEEE double exp */
/* x is now in primary range */
t = x*x;
if(k >= -1021)
INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0);
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0);
else
INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0);
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);

7
modules/fdlibm/src/e_hypot.cpp
View File

@ -46,6 +46,7 @@
* than 1 ulps (units in the last place)
*/
#include <cmath>
#include <float.h>
#include "math_private.h"
@ -69,7 +70,7 @@ __ieee754_hypot(double x, double y)
if(ha >= 0x7ff00000) { /* Inf or NaN */
u_int32_t low;
/* Use original arg order iff result is NaN; quieten sNaNs. */
w = fabs(x+0.0)-fabs(y+0.0);
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);
@ -105,7 +106,7 @@ __ieee754_hypot(double x, double y)
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;
@ -114,7 +115,7 @@ __ieee754_hypot(double x, double y)
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;

54
modules/fdlibm/src/e_pow.cpp
View File

@ -4,7 +4,7 @@
* 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
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@ -19,7 +19,7 @@
* 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
* 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)
*
@ -47,16 +47,19 @@
* Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular
* pow(integer,integer)
* always returns the correct integer provided it is
* 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
* 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
@ -64,6 +67,9 @@ 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 */
@ -106,15 +112,15 @@ __ieee754_pow(double x, double y)
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
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 (x+0.0)+(y+0.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
@ -122,22 +128,22 @@ __ieee754_pow(double x, double y)
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<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((j<<(52-k))==ly) yisint = 2-(j&1);
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(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if(((ix-0x3ff00000)|lx)==0)
return one; /* (-1)**+-inf is 1 */
@ -145,14 +151,14 @@ __ieee754_pow(double x, double y)
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 sqrt(x);
return std::sqrt(x);
}
}
@ -165,13 +171,13 @@ __ieee754_pow(double x, double y)
if(hx<0) {
if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} else if(yisint==1)
} 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
@ -193,10 +199,10 @@ __ieee754_pow(double x, double y)
/* 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
/* 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)*(0.5-t*(0.3333333333333333333333-t*0.25));
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;
@ -233,9 +239,9 @@ __ieee754_pow(double x, double y)
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.0+s2+r;
t_h = 3+s2+r;
SET_LOW_WORD(t_h,0);
t_l = r-((t_h-3.0)-s2);
t_l = r-((t_h-3)-s2);
/* u+v = ss*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*ss;
@ -246,7 +252,7 @@ __ieee754_pow(double x, double y)
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 = (double)n;
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);
@ -286,7 +292,7 @@ __ieee754_pow(double x, double y)
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;

446
modules/fdlibm/src/e_sqrt.cpp
View File

@ -1,446 +0,0 @@
/* @(#)e_sqrt.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_sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
* | Use the hardware sqrt if you have one |
* ------------------------------------------
* Method:
* Bit by bit method using integer arithmetic. (Slow, but portable)
* 1. Normalization
* Scale x to y in [1,4) with even powers of 2:
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
* sqrt(x) = 2^k * sqrt(y)
* 2. Bit by bit computation
* Let q = sqrt(y) truncated to i bit after binary point (q = 1),
* i 0
* i+1 2
* s = 2*q , and y = 2 * ( y - q ). (1)
* i i i i
*
* To compute q from q , one checks whether
* i+1 i
*
* -(i+1) 2
* (q + 2 ) <= y. (2)
* i
* -(i+1)
* If (2) is false, then q = q ; otherwise q = q + 2 .
* i+1 i i+1 i
*
* With some algebric manipulation, it is not difficult to see
* that (2) is equivalent to
* -(i+1)
* s + 2 <= y (3)
* i i
*
* The advantage of (3) is that s and y can be computed by
* i i
* the following recurrence formula:
* if (3) is false
*
* s = s , y = y ; (4)
* i+1 i i+1 i
*
* otherwise,
* -i -(i+1)
* s = s + 2 , y = y - s - 2 (5)
* i+1 i i+1 i i
*
* One may easily use induction to prove (4) and (5).
* Note. Since the left hand side of (3) contain only i+2 bits,
* it does not necessary to do a full (53-bit) comparison
* in (3).
* 3. Final rounding
* After generating the 53 bits result, we compute one more bit.
* Together with the remainder, we can decide whether the
* result is exact, bigger than 1/2ulp, or less than 1/2ulp
* (it will never equal to 1/2ulp).
* The rounding mode can be detected by checking whether
* huge + tiny is equal to huge, and whether huge - tiny is
* equal to huge for some floating point number "huge" and "tiny".
*
* Special cases:
* sqrt(+-0) = +-0 ... exact
* sqrt(inf) = inf
* sqrt(-ve) = NaN ... with invalid signal
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
*
* Other methods : see the appended file at the end of the program below.
*---------------
*/
#include <float.h>
#include "math_private.h"
static const double one = 1.0, tiny=1.0e-300;
double
__ieee754_sqrt(double x)
{
double z;
int32_t sign = (int)0x80000000;
int32_t ix0,s0,q,m,t,i;
u_int32_t r,t1,s1,ix1,q1;
EXTRACT_WORDS(ix0,ix1,x);
/* take care of Inf and NaN */
if((ix0&0x7ff00000)==0x7ff00000) {
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */
}
/* take care of zero */
if(ix0<=0) {
if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
else if(ix0<0)
return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix0>>20);
if(m==0) { /* subnormal x */
while(ix0==0) {
m -= 21;
ix0 |= (ix1>>11); ix1 <<= 21;
}
for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
m -= i-1;
ix0 |= (ix1>>(32-i));
ix1 <<= i;
}
m -= 1023; /* unbias exponent */
ix0 = (ix0&0x000fffff)|0x00100000;
if(m&1){ /* odd m, double x to make it even */
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* r = moving bit from right to left */
while(r!=0) {
t = s0+r;
if(t<=ix0) {
s0 = t+r;
ix0 -= t;
q += r;
}
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
r>>=1;
}
r = sign;
while(r!=0) {
t1 = s1+r;
t = s0;
if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
s1 = t1+r;
if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
ix0 -= t;
if (ix1 < t1) ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
r>>=1;
}
/* use floating add to find out rounding direction */
if((ix0|ix1)!=0) {
z = one-tiny; /* trigger inexact flag */
if (z>=one) {
z = one+tiny;
if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
else if (z>one) {
if (q1==(u_int32_t)0xfffffffe) q+=1;
q1+=2;
} else
q1 += (q1&1);
}
}
ix0 = (q>>1)+0x3fe00000;
ix1 = q1>>1;
if ((q&1)==1) ix1 |= sign;
ix0 += (m <<20);
INSERT_WORDS(z,ix0,ix1);
return z;
}
/*
Other methods (use floating-point arithmetic)
-------------
(This is a copy of a drafted paper by Prof W. Kahan
and K.C. Ng, written in May, 1986)
Two algorithms are given here to implement sqrt(x)
(IEEE double precision arithmetic) in software.
Both supply sqrt(x) correctly rounded. The first algorithm (in
Section A) uses newton iterations and involves four divisions.
The second one uses reciproot iterations to avoid division, but
requires more multiplications. Both algorithms need the ability
to chop results of arithmetic operations instead of round them,
and the INEXACT flag to indicate when an arithmetic operation
is executed exactly with no roundoff error, all part of the
standard (IEEE 754-1985). The ability to perform shift, add,
subtract and logical AND operations upon 32-bit words is needed
too, though not part of the standard.
A. sqrt(x) by Newton Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
1 11 52 ...widths
------------------------------------------------------
x: |s| e | f |
------------------------------------------------------
msb lsb msb lsb ...order
------------------------ ------------------------
x0: |s| e | f1 | x1: | f2 |
------------------------ ------------------------
By performing shifts and subtracts on x0 and x1 (both regarded
as integers), we obtain an 8-bit approximation of sqrt(x) as
follows.
k := (x0>>1) + 0x1ff80000;
y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
Here k is a 32-bit integer and T1[] is an integer array containing
correction terms. Now magically the floating value of y (y's
leading 32-bit word is y0, the value of its trailing word is 0)
approximates sqrt(x) to almost 8-bit.
Value of T1:
static int T1[32]= {
0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
(2) Iterative refinement
Apply Heron's rule three times to y, we have y approximates
sqrt(x) to within 1 ulp (Unit in the Last Place):
y := (y+x/y)/2 ... almost 17 sig. bits
y := (y+x/y)/2 ... almost 35 sig. bits
y := y-(y-x/y)/2 ... within 1 ulp
Remark 1.
Another way to improve y to within 1 ulp is:
y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
2
(x-y )*y
y := y + 2* ---------- ...within 1 ulp
2
3y + x
This formula has one division fewer than the one above; however,
it requires more multiplications and additions. Also x must be
scaled in advance to avoid spurious overflow in evaluating the
expression 3y*y+x. Hence it is not recommended uless division
is slow. If division is very slow, then one should use the
reciproot algorithm given in section B.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
I := FALSE; ... reset INEXACT flag I
R := RZ; ... set rounding mode to round-toward-zero
z := x/y; ... chopped quotient, possibly inexact
If(not I) then { ... if the quotient is exact
if(z=y) {
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
} else {
z := z - ulp; ... special rounding
}
}
i := TRUE; ... sqrt(x) is inexact
If (r=RN) then z=z+ulp ... rounded-to-nearest
If (r=RP) then { ... round-toward-+inf
y = y+ulp; z=z+ulp;
}
y := y+z; ... chopped sum
y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
(4) Special cases
Square root of +inf, +-0, or NaN is itself;
Square root of a negative number is NaN with invalid signal.
B. sqrt(x) by Reciproot Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
(see section A). By performing shifs and subtracts on x0 and y0,
we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
k := 0x5fe80000 - (x0>>1);
y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
Here k is a 32-bit integer and T2[] is an integer array
containing correction terms. Now magically the floating
value of y (y's leading 32-bit word is y0, the value of
its trailing word y1 is set to zero) approximates 1/sqrt(x)
to almost 7.8-bit.
Value of T2:
static int T2[64]= {
0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
(2) Iterative refinement
Apply Reciproot iteration three times to y and multiply the
result by x to get an approximation z that matches sqrt(x)
to about 1 ulp. To be exact, we will have
-1ulp < sqrt(x)-z<1.0625ulp.
... set rounding mode to Round-to-nearest
y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
... special arrangement for better accuracy
z := x*y ... 29 bits to sqrt(x), with z*y<1
z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
(a) the term z*y in the final iteration is always less than 1;
(b) the error in the final result is biased upward so that
-1 ulp < sqrt(x) - z < 1.0625 ulp
instead of |sqrt(x)-z|<1.03125ulp.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
R := RZ; ... set rounding mode to round-toward-zero
switch(r) {
case RN: ... round-to-nearest
if(x<= z*(z-ulp)...chopped) z = z - ulp; else
if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
break;
case RZ:case RM: ... round-to-zero or round-to--inf
R:=RP; ... reset rounding mod to round-to-+inf
if(x<z*z ... rounded up) z = z - ulp; else
if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
break;
case RP: ... round-to-+inf
if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
if(x>z*z ...chopped) z = z+ulp;
break;
}
Remark 3. The above comparisons can be done in fixed point. For
example, to compare x and w=z*z chopped, it suffices to compare
x1 and w1 (the trailing parts of x and w), regarding them as
two's complement integers.
...Is z an exact square root?
To determine whether z is an exact square root of x, let z1 be the
trailing part of z, and also let x0 and x1 be the leading and
trailing parts of x.
If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
I := 1; ... Raise Inexact flag: z is not exact
else {
j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
k := z1 >> 26; ... get z's 25-th and 26-th
fraction bits
I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
}
R:= r ... restore rounded mode
return sqrt(x):=z.
If multiplication is cheaper then the foregoing red tape, the
Inexact flag can be evaluated by
I := i;
I := (z*z!=x) or I.
Note that z*z can overwrite I; this value must be sensed if it is
True.
Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
zero.
--------------------
z1: | f2 |
--------------------
bit 31 bit 0
Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
or even of logb(x) have the following relations:
-------------------------------------------------
bit 27,26 of z1 bit 1,0 of x1 logb(x)
-------------------------------------------------
00 00 odd and even
01 01 even
10 10 odd
10 00 even
11 01 even
-------------------------------------------------
(4) Special cases (see (4) of Section A).
*/

1
modules/fdlibm/src/fdlibm.h
View File

@ -33,7 +33,6 @@ double log(double);
double log10(double);
double pow(double, double);
double sqrt(double);
double fabs(double);
double floor(double);

2
modules/fdlibm/src/k_exp.cpp
View File

@ -1,4 +1,6 @@
/*-
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*

129
modules/fdlibm/src/math_private.h
View File

@ -45,6 +45,47 @@
#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.
@ -307,8 +348,9 @@ do { \
/* Support switching the mode to FP_PE if necessary. */
#if defined(__i386__) && !defined(NO_FPSETPREC)
#define ENTERI() \
long double __retval; \
#define ENTERI() ENTERIT(long double)
#define ENTERIT(returntype) \
returntype __retval; \
fp_prec_t __oprec; \
\
if ((__oprec = fpgetprec()) != FP_PE) \
@ -319,9 +361,22 @@ do { \
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(x)
#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. */
@ -436,6 +491,31 @@ do { \
*/
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
/*
@ -511,48 +591,6 @@ CMPLXL(long double x, long double y)
#endif /* _COMPLEX_H */
#ifdef __GNUCLIKE_ASM
/* Asm versions of some functions. */
#ifdef __amd64__
static __inline int
irint(double x)
{
int n;
asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
return (n);
}
#define HAVE_EFFICIENT_IRINT
#endif
#ifdef __i386__
static __inline int
irint(double x)
{
int n;
asm("fistl %0" : "=m" (n) : "t" (x));
return (n);
}
#define HAVE_EFFICIENT_IRINT
#endif
#if defined(__amd64__) || defined(__i386__)
static __inline int
irintl(long double x)
{
int n;
asm("fistl %0" : "=m" (n) : "t" (x));
return (n);
}
#define HAVE_EFFICIENT_IRINTL
#endif
#endif /* __GNUCLIKE_ASM */
#ifdef DEBUG
#if defined(__amd64__) || defined(__i386__)
#define breakpoint() asm("int $3")
@ -759,7 +797,6 @@ irintl(long double x)
#define log fdlibm::log
#define log10 fdlibm::log10
#define pow fdlibm::pow
#define sqrt fdlibm::sqrt
#define ceil fdlibm::ceil
#define ceilf fdlibm::ceilf
#define fabs fdlibm::fabs

18
modules/fdlibm/src/moz.build
View File

@ -10,26 +10,35 @@ EXPORTS += [
FINAL_LIBRARY = 'js'
if CONFIG['GNU_CXX']:
if CONFIG['CC_TYPE'] in ('clang', 'gcc'):
CXXFLAGS += [
'-Wno-parentheses',
'-Wno-sign-compare',
]
if CONFIG['CLANG_CXX']:
if CONFIG['CC_TYPE'] == 'clang':
CXXFLAGS += [
'-Wno-dangling-else',
]
if CONFIG['_MSC_VER']:
if CONFIG['CC_TYPE'] in ('msvc', 'clang-cl'):
CXXFLAGS += [
'-wd4018', # signed/unsigned mismatch
'-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',
@ -44,7 +53,6 @@ SOURCES += [
'e_log2.cpp',
'e_pow.cpp',
'e_sinh.cpp',
'e_sqrt.cpp',
'k_exp.cpp',
's_asinh.cpp',
's_atan.cpp',

5
modules/fdlibm/src/s_asinh.cpp
View File

@ -24,6 +24,7 @@
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
#include <cmath>
#include <float.h>
#include "math_private.h"
@ -48,10 +49,10 @@ asinh(double x)
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;
}

1
modules/fdlibm/src/s_cbrt.cpp
View File

@ -15,6 +15,7 @@
//#include <sys/cdefs.h>
//__FBSDID("$FreeBSD$");
#include <float.h>
#include "math_private.h"
/* cbrt(x)

2
modules/fdlibm/src/s_expm1.cpp
View File

@ -187,7 +187,7 @@ expm1(double x)
e = hxs*((r1-t)/(6.0 - x*t));
if(k==0) return x - (x*e-hxs); /* c is 0 */
else {
INSERT_WORDS(twopk,0x3ff00000+(k<<20),0); /* 2^k */
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;

7
modules/fdlibm/src/s_fabs.cpp
View File

@ -1,4 +1,4 @@
/* @(#)s_fabs.c 5.1 93/09/24 */
/* @(#)s_fabs.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,9 +10,8 @@
* ====================================================
*/
#ifndef lint
//static char rcsid[] = "$FreeBSD$";
#endif
//#include <sys/cdefs.h>
//__FBSDID("$FreeBSD$");
/*
* fabs(x) returns the absolute value of x.

2
modules/fdlibm/src/s_nearbyint.cpp
View File

@ -1,4 +1,6 @@
/*-
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
* Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*

12
modules/fdlibm/src/s_scalbn.cpp
View File

@ -10,9 +10,8 @@
* ====================================================
*/
#ifndef lint
//static char rcsid[] = "$FreeBSD$";
#endif
//#include <sys/cdefs.h>
//__FBSDID("$FreeBSD$");
/*
* scalbn (double x, int n)
@ -21,7 +20,6 @@
* exponentiation or a multiplication.
*/
//#include <sys/cdefs.h>
#include <float.h>
#include "math_private.h"
@ -50,10 +48,12 @@ scalbn (double x, int 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 (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*/
else
return tiny*copysign(tiny,x); /*underflow*/
}
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;

23
modules/fdlibm/update.sh
View File

@ -11,23 +11,28 @@ 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
BEFORE_COMMIT=$(get_commit)
sh ./import.sh
mv ./src_moz.build ./src/moz.build
COMMIT=$(get_commit)
if [ ${BEFORE_COMMIT} != ${COMMIT} ]; then
echo "Latest commit is changed during import. Please run again."
exit 1
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
patch -p3 < ${FILE}
echo "Applying ${FILE} ..."
patch -p3 --no-backup-if-mismatch < ${FILE}
done
hg add src
perl -p -i -e "s/\[commit [0-9a-f]{40}\]/[commit ${COMMIT}]/" README.mozilla
perl -p -i -e "s/\[commit [0-9a-f]{40} \(.{1,100}\)\]/[commit ${COMMIT} (${COMMITDATE})]/" README.mozilla
echo "###"
echo "### Updated fdlibm/src to ${COMMIT}."