Magic-Mirror/gecko-dev
1
0
Fork 0
mirror of https://github.com/mozilla/gecko-dev.git synced 2025-02-26 12:20:56 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Bug 1564117 - Part 1: Use std::sqrt in fdlibm. r=arai

Browse Source 
std::sqrt is guaranteed to have an error less than 0.5 ulp, so there's no need
to use fdlibm's slower implementation.

Differential Revision: https://phabricator.services.mozilla.com/D37256

--HG--
extra : moz-landing-system : lando
This commit is contained in:
André Bargull 2019-07-08 16:52:45 +00:00
parent 33a8904ee9
commit eebb17f171
18 changed files with 276 additions and 534 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
modules/fdlibm/README.mozilla
View File

@ -13,5 +13,5 @@ from within the modules/fdlibm directory.
Current version: [commit b21ccf63f28a3a4692d8a31419e0a725a1b1a800 (2018-02-14T07:59:30Z)].
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

1
modules/fdlibm/import.sh
View File

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

4
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,499 +12,50 @@
@@ -12,499 +12,49 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
* $FreeBSD$
@ -242,7 +242,7 @@ 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 ceil(double);

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

@ -271,27 +271,6 @@ 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

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
@@ -736,16 +736,51 @@ irintl(long double x)
@@ -742,16 +742,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

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

@ -190,30 +190,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

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

@ -325,30 +325,6 @@ 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

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 "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,
one = 1.0,
@@ -147,17 +149,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;
}

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

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

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

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

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

@ -57,6 +57,8 @@
* to produce the hexadecimal values shown.
*/
#include <cmath>
#include "math_private.h"
static const double
@ -152,7 +154,7 @@ __ieee754_pow(double x, double y)
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);
}
}

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

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

@ -771,7 +771,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

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

@ -46,7 +46,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;
}