Revert "Update to FreeBSD libm r336523."
This reverts commit f86ee10278116d211cfb08679416444e0e714cb5.
Incorrect result for fmodf(3.0f, 0f) = 1.0f breaks ART tests.
Bug: 111710419
Test: art/test/testrunner/testrunner.py -b -t 436-rem-float --target
Change-Id: I7eae68fb92740db33415d16418447bcbbd98ecba
diff --git a/libm/upstream-freebsd/lib/msun/src/catrig.c b/libm/upstream-freebsd/lib/msun/src/catrig.c
index 4a2d076..025076f 100644
--- a/libm/upstream-freebsd/lib/msun/src/catrig.c
+++ b/libm/upstream-freebsd/lib/msun/src/catrig.c
@@ -27,7 +27,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 327232 2017-12-27 03:23:41Z eadler $");
#include <complex.h>
#include <float.h>
@@ -300,7 +300,7 @@
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
- return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -384,7 +384,7 @@
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
- return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -601,7 +601,7 @@
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
- return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
diff --git a/libm/upstream-freebsd/lib/msun/src/catrigf.c b/libm/upstream-freebsd/lib/msun/src/catrigf.c
index e251871..344290a 100644
--- a/libm/upstream-freebsd/lib/msun/src/catrigf.c
+++ b/libm/upstream-freebsd/lib/msun/src/catrigf.c
@@ -41,7 +41,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <float.h>
@@ -163,7 +163,7 @@
return (CMPLXF(y, x + x));
if (y == 0)
return (CMPLXF(x + x, y));
- return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -221,7 +221,7 @@
return (CMPLXF(x + x, -y));
if (x == 0)
return (CMPLXF(pio2_hi + pio2_lo, y + y));
- return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -359,7 +359,7 @@
if (isinf(y))
return (CMPLXF(copysignf(0, x),
copysignf(pio2_hi + pio2_lo, y)));
- return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
diff --git a/libm/upstream-freebsd/lib/msun/src/catrigl.c b/libm/upstream-freebsd/lib/msun/src/catrigl.c
index 63b068a..960c1ca 100644
--- a/libm/upstream-freebsd/lib/msun/src/catrigl.c
+++ b/libm/upstream-freebsd/lib/msun/src/catrigl.c
@@ -40,7 +40,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/catrigl.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/catrigl.c 323003 2017-08-29 22:32:29Z rlibby $");
#include <complex.h>
#include <float.h>
@@ -182,7 +182,7 @@
return (CMPLXL(y, x + x));
if (y == 0)
return (CMPLXL(x + x, y));
- return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -241,7 +241,7 @@
return (CMPLXL(x + x, -y));
if (x == 0)
return (CMPLXL(pio2_hi + pio2_lo, y + y));
- return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -380,7 +380,7 @@
if (isinf(y))
return (CMPLXL(copysignl(0, x),
copysignl(pio2_hi + pio2_lo, y)));
- return (CMPLXL(nan_mix(x, y), nan_mix(x, y)));
+ return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atan2.c b/libm/upstream-freebsd/lib/msun/src/e_atan2.c
index bef9df8..ee5b15d 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atan2.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atan2.c
@@ -13,7 +13,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_atan2.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_atan2.c 329259 2018-02-14 07:59:30Z eadler $");
/* __ieee754_atan2(y,x)
* Method :
@@ -70,7 +70,7 @@
iy = hy&0x7fffffff;
if(((ix|((lx|-lx)>>31))>0x7ff00000)||
((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
- return nan_mix(x, y);
+ return 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) */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atan2f.c b/libm/upstream-freebsd/lib/msun/src/e_atan2f.c
index 0ce9bc0..fc77bff 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atan2f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atan2f.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_atan2f.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
@@ -41,7 +41,7 @@
iy = hy&0x7fffffff;
if((ix>0x7f800000)||
(iy>0x7f800000)) /* x or y is NaN */
- return nan_mix(x, y);
+ return x+y;
if(hx==0x3f800000) return atanf(y); /* x=1.0 */
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atan2l.c b/libm/upstream-freebsd/lib/msun/src/e_atan2l.c
index 8f80acc..0326482 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atan2l.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atan2l.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_atan2l.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
/*
* See comments in e_atan2.c.
@@ -62,7 +62,7 @@
((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */
(expty==BIAS+LDBL_MAX_EXP &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */
- return nan_mix(x, y);
+ return x+y;
if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
return atanl(y); /* x=1.0 */
m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_fmod.c b/libm/upstream-freebsd/lib/msun/src/e_fmod.c
index 36f98d1..bad066a 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_fmod.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_fmod.c
@@ -12,7 +12,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_fmod.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_fmod.c 305380 2016-09-04 12:01:32Z bde $");
/*
* __ieee754_fmod(x,y)
@@ -42,7 +42,7 @@
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
- return nan_mix(x, y)/nan_mix(x, y);
+ return (x*y)/(x*y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
if(lx==ly)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_fmodf.c b/libm/upstream-freebsd/lib/msun/src/e_fmodf.c
index 5a103f2..52ce373 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_fmodf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_fmodf.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_fmodf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
/*
* __ieee754_fmodf(x,y)
@@ -41,7 +41,7 @@
/* purge off exception values */
if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */
(hy>0x7f800000)) /* or y is NaN */
- return nan_mix(x, y)/nan_mix(x, y);
+ return (x*y)/(x*y);
if(hx<hy) return x; /* |x|<|y| return x */
if(hx==hy)
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
diff --git a/libm/upstream-freebsd/lib/msun/src/e_fmodl.c b/libm/upstream-freebsd/lib/msun/src/e_fmodl.c
index 6001df2..e315f76 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_fmodl.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_fmodl.c
@@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_fmodl.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include <float.h>
#include <stdint.h>
@@ -79,7 +79,7 @@
(ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(uy.bits.exp == BIAS + LDBL_MAX_EXP &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
- return nan_mix(x, y)/nan_mix(x, y);
+ return (x*y)/(x*y);
if(ux.bits.exp<=uy.bits.exp) {
if((ux.bits.exp<uy.bits.exp) ||
(ux.bits.manh<=uy.bits.manh &&
diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypot.c b/libm/upstream-freebsd/lib/msun/src/e_hypot.c
index 0c16767..2398e98 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypot.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypot.c
@@ -12,7 +12,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_hypot.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
/* __ieee754_hypot(x,y)
*
@@ -70,7 +70,7 @@
if(ha >= 0x7ff00000) { /* Inf or NaN */
u_int32_t low;
/* Use original arg order iff result is NaN; quieten sNaNs. */
- w = fabsl(x+0.0L)-fabs(y+0);
+ w = fabs(x+0.0)-fabs(y+0.0);
GET_LOW_WORD(low,a);
if(((ha&0xfffff)|low)==0) w = a;
GET_LOW_WORD(low,b);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypotf.c b/libm/upstream-freebsd/lib/msun/src/e_hypotf.c
index 79e4697..6d083e4 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypotf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypotf.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_hypotf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
@@ -37,7 +37,7 @@
if(ha > 0x58800000) { /* a>2**50 */
if(ha >= 0x7f800000) { /* Inf or NaN */
/* Use original arg order iff result is NaN; quieten sNaNs. */
- w = fabsl(x+0.0L)-fabsf(y+0);
+ w = fabsf(x+0.0F)-fabsf(y+0.0F);
if(ha == 0x7f800000) w = a;
if(hb == 0x7f800000) w = b;
return w;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypotl.c b/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
index 5603b6e..7b5ab89 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
@@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_hypotl.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
/* long double version of hypot(). See e_hypot.c for most comments. */
@@ -64,7 +64,7 @@
if(ha >= ESW(MAX_EXP)) { /* Inf or NaN */
man_t manh, manl;
/* Use original arg order iff result is NaN; quieten sNaNs. */
- w = fabsl(x+0.0L)-fabsl(y+0);
+ w = fabsl(x+0.0)-fabsl(y+0.0);
GET_LDBL_MAN(manh,manl,a);
if (manh == LDBL_NBIT && manl == 0) w = a;
GET_LDBL_MAN(manh,manl,b);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j0.c b/libm/upstream-freebsd/lib/msun/src/e_j0.c
index cfa71c2..36e72c2 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j0.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j0.c
@@ -1,3 +1,4 @@
+
/* @(#)e_j0.c 1.3 95/01/18 */
/*
* ====================================================
@@ -5,13 +6,13 @@
*
* 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.
* ====================================================
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 336089 2018-07-08 16:26:13Z markj $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 283032 2015-05-17 16:27:06Z kargl $");
/* __ieee754_j0(x), __ieee754_y0(x)
* Bessel function of the first and second kinds of order zero.
@@ -32,20 +33,20 @@
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
- *
+ *
* 3 Special cases
* j0(nan)= nan
* j0(0) = 1
* j0(inf) = 0
- *
+ *
* Method -- y0(x):
* 1. For x<2.
- * Since
+ * Since
* y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
* therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
* We use the following function to approximate y0,
* y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
- * where
+ * where
* U(z) = u00 + u01*z + ... + u06*z^6
* V(z) = 1 + v01*z + ... + v04*z^4
* with absolute approximation error bounded by 2**-72.
@@ -70,7 +71,7 @@
one = 1.0,
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
-/* R0/S0 on [0, 2.00] */
+ /* R0/S0 on [0, 2.00] */
R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
@@ -156,7 +157,7 @@
* y0(Inf) = 0.
* y0(-Inf) = NaN and raise invalid exception.
*/
- if(ix>=0x7ff00000) return vone/(x+x*x);
+ if(ix>=0x7ff00000) return vone/(x+x*x);
/* y0(+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* y0(x<0) = NaN and raise invalid exception. */
@@ -292,7 +293,7 @@
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return one+ r/s;
}
-
+
/* For x >= 8, the asymptotic expansions of qzero is
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j1.c b/libm/upstream-freebsd/lib/msun/src/e_j1.c
index 78bb329..b11ac2d 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j1.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j1.c
@@ -1,3 +1,4 @@
+
/* @(#)e_j1.c 1.3 95/01/18 */
/*
* ====================================================
@@ -5,13 +6,13 @@
*
* 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.
* ====================================================
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_j1.c 336089 2018-07-08 16:26:13Z markj $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_j1.c 283032 2015-05-17 16:27:06Z kargl $");
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
@@ -33,16 +34,16 @@
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
- *
+ *
* 3 Special cases
* j1(nan)= nan
* j1(0) = 0
* j1(inf) = 0
- *
+ *
* Method -- y1(x):
- * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
+ * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
* 2. For x<2.
- * Since
+ * Since
* y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
* therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
* We use the following function to approximate y1,
@@ -153,7 +154,7 @@
* y1(Inf) = 0.
* y1(-Inf) = NaN and raise invalid exception.
*/
- if(ix>=0x7ff00000) return vone/(x+x*x);
+ if(ix>=0x7ff00000) return vone/(x+x*x);
/* y1(+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* y1(x<0) = NaN and raise invalid exception. */
@@ -185,10 +186,10 @@
z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
}
return z;
- }
+ }
if(ix<=0x3c900000) { /* x < 2**-54 */
return(-tpi/x);
- }
+ }
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
@@ -286,7 +287,7 @@
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return one+ r/s;
}
-
+
/* For x >= 8, the asymptotic expansions of qone is
* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j1f.c b/libm/upstream-freebsd/lib/msun/src/e_j1f.c
index 3abe201..0cca823 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j1f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j1f.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_j1f.c 336089 2018-07-08 16:26:13Z markj $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_j1f.c 283032 2015-05-17 16:27:06Z kargl $");
/*
* See e_j1.c for complete comments.
@@ -32,7 +32,7 @@
one = 1.0,
invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
-/* R0/S0 on [0,2] */
+ /* R0/S0 on [0,2] */
r00 = -6.2500000000e-02, /* 0xbd800000 */
r01 = 1.4070566976e-03, /* 0x3ab86cfd */
r02 = -1.5995563444e-05, /* 0xb7862e36 */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_jn.c b/libm/upstream-freebsd/lib/msun/src/e_jn.c
index 58ec905..a1130c5 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_jn.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_jn.c
@@ -1,3 +1,4 @@
+
/* @(#)e_jn.c 1.4 95/01/18 */
/*
* ====================================================
@@ -5,19 +6,19 @@
*
* 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.
* ====================================================
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_jn.c 336089 2018-07-08 16:26:13Z markj $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_jn.c 279856 2015-03-10 17:10:54Z kargl $");
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
* floating point Bessel's function of the 1st and 2nd kind
* of order n
- *
+ *
* Special cases:
* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
@@ -36,6 +37,7 @@
* yn(n,x) is similar in all respects, except
* that forward recursion is used for all
* values of n>1.
+ *
*/
#include "math.h"
@@ -64,7 +66,7 @@
ix = 0x7fffffff&hx;
/* if J(n,NaN) is NaN */
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
- if(n<0){
+ if(n<0){
n = -n;
x = -x;
hx ^= 0x80000000;
@@ -75,14 +77,14 @@
x = fabs(x);
if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
b = zero;
- else if((double)n<=x) {
+ else if((double)n<=x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if(ix>=0x52D00000) { /* x > 2**302 */
- /* (x >> n**2)
+ /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
- * Let s=sin(x), c=cos(x),
- * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
+ * Let s=sin(x), c=cos(x),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
@@ -98,7 +100,7 @@
case 3: temp = cos(x)-sin(x); break;
}
b = invsqrtpi*temp/sqrt(x);
- } else {
+ } else {
a = __ieee754_j0(x);
b = __ieee754_j1(x);
for(i=1;i<n;i++){
@@ -109,7 +111,7 @@
}
} else {
if(ix<0x3e100000) { /* x < 2**-29 */
- /* x is tiny, return the first Taylor expansion of J(n,x)
+ /* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if(n>33) /* underflow */
@@ -124,14 +126,14 @@
}
} else {
/* use backward recurrence */
- /* x x^2 x^2
+ /* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
- * 1 1 1
+ * 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
- * -- - ------ - ------ -
+ * -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
@@ -147,9 +149,9 @@
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
- * When Q(k) > 1e4 good for single
- * When Q(k) > 1e9 good for double
- * When Q(k) > 1e17 good for quadruple
+ * When Q(k) > 1e4 good for single
+ * When Q(k) > 1e9 good for double
+ * When Q(k) > 1e17 good for quadruple
*/
/* determine k */
double t,v;
@@ -235,11 +237,11 @@
if(n==1) return(sign*__ieee754_y1(x));
if(ix==0x7ff00000) return zero;
if(ix>=0x52D00000) { /* x > 2**302 */
- /* (x >> n**2)
+ /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
- * Let s=sin(x), c=cos(x),
- * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
+ * Let s=sin(x), c=cos(x),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
diff --git a/libm/upstream-freebsd/lib/msun/src/e_pow.c b/libm/upstream-freebsd/lib/msun/src/e_pow.c
index 411dd52..d088b41 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_pow.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_pow.c
@@ -10,7 +10,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_pow.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/e_pow.c 326482 2017-12-03 01:56:03Z emaste $");
/* __ieee754_pow(x,y) return x**y
*
@@ -57,7 +57,6 @@
* to produce the hexadecimal values shown.
*/
-#include <float.h>
#include "math.h"
#include "math_private.h"
@@ -66,9 +65,6 @@
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 */
@@ -119,7 +115,7 @@
/* y!=zero: result is NaN if either arg is NaN */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
- return nan_mix(x, y);
+ return (x+0.0)+(y+0.0);
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
@@ -201,7 +197,7 @@
/* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax-one; /* t has 20 trailing zeros */
- w = (t*t)*(half-t*(thrd-t*qrtr));
+ w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
@@ -238,9 +234,9 @@
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+ss);
s2 = s_h*s_h;
- t_h = 3+s2+r;
+ t_h = 3.0+s2+r;
SET_LOW_WORD(t_h,0);
- t_l = r-((t_h-3)-s2);
+ t_l = r-((t_h-3.0)-s2);
/* u+v = ss*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*ss;
@@ -251,7 +247,7 @@
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
- t = n;
+ t = (double)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);
@@ -308,7 +304,3 @@
else SET_HIGH_WORD(z,j);
return s*z;
}
-
-#if (LDBL_MANT_DIG == 53)
-__weak_reference(pow, powl);
-#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/e_powf.c b/libm/upstream-freebsd/lib/msun/src/e_powf.c
index 13dd4eb..5c46478 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_powf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_powf.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_powf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
@@ -24,9 +24,6 @@
dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
zero = 0.0,
-half = 0.5,
-qrtr = 0.25,
-thrd = 3.33333343e-01, /* 0x3eaaaaab */
one = 1.0,
two = 2.0,
two24 = 16777216.0, /* 0x4b800000 */
@@ -76,7 +73,7 @@
/* y!=zero: result is NaN if either arg is NaN */
if(ix > 0x7f800000 ||
iy > 0x7f800000)
- return nan_mix(x, y);
+ return (x+0.0F)+(y+0.0F);
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
@@ -141,7 +138,7 @@
/* 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-1; /* t has 20 trailing zeros */
- w = (t*t)*(half-t*(thrd-t*qrtr));
+ w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
u = ivln2_h*t; /* ivln2_h has 16 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
@@ -180,10 +177,10 @@
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+s);
s2 = s_h*s_h;
- t_h = 3+s2+r;
+ t_h = (float)3.0+s2+r;
GET_FLOAT_WORD(is,t_h);
SET_FLOAT_WORD(t_h,is&0xfffff000);
- t_l = r-((t_h-3)-s2);
+ t_l = r-((t_h-(float)3.0)-s2);
/* u+v = s*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*s;
@@ -195,7 +192,7 @@
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) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
- t = n;
+ t = (float)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
GET_FLOAT_WORD(is,t1);
SET_FLOAT_WORD(t1,is&0xfffff000);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_remainder.c b/libm/upstream-freebsd/lib/msun/src/e_remainder.c
index a4d277f..9be513b 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_remainder.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_remainder.c
@@ -12,7 +12,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_remainder.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
/* __ieee754_remainder(x,p)
* Return :
@@ -45,11 +45,11 @@
hx &= 0x7fffffff;
/* purge off exception values */
- if(((hp|lp)==0)|| /* p = 0 */
- (hx>=0x7ff00000)|| /* x not finite */
+ if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
+ if((hx>=0x7ff00000)|| /* x not finite */
((hp>=0x7ff00000)&& /* p is NaN */
(((hp-0x7ff00000)|lp)!=0)))
- return nan_mix(x, p)/nan_mix(x, p);
+ return ((long double)x*p)/((long double)x*p);
if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_remainderf.c b/libm/upstream-freebsd/lib/msun/src/e_remainderf.c
index 89cd932..b0014ae 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_remainderf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_remainderf.c
@@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/e_remainderf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
@@ -36,10 +36,10 @@
hx &= 0x7fffffff;
/* purge off exception values */
- if((hp==0)|| /* p = 0 */
- (hx>=0x7f800000)|| /* x not finite */
+ if(hp==0) return (x*p)/(x*p); /* p = 0 */
+ if((hx>=0x7f800000)|| /* x not finite */
((hp>0x7f800000))) /* p is NaN */
- return nan_mix(x, p)/nan_mix(x, p);
+ return ((long double)x*p)/((long double)x*p);
if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */
diff --git a/libm/upstream-freebsd/lib/msun/src/math_private.h b/libm/upstream-freebsd/lib/msun/src/math_private.h
index 7832304..637a09a 100644
--- a/libm/upstream-freebsd/lib/msun/src/math_private.h
+++ b/libm/upstream-freebsd/lib/msun/src/math_private.h
@@ -11,7 +11,7 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math_private.h 336362 2018-07-17 07:42:14Z bde $
+ * $FreeBSD: head/lib/msun/src/math_private.h 319047 2017-05-28 06:13:38Z mmel $
*/
#ifndef _MATH_PRIVATE_H_
@@ -48,47 +48,6 @@
#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
-
-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 IEEE_WORD_ORDER == 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
-
#if IEEE_WORD_ORDER == BIG_ENDIAN
typedef union
@@ -335,9 +294,8 @@
/* Support switching the mode to FP_PE if necessary. */
#if defined(__i386__) && !defined(NO_FPSETPREC)
-#define ENTERI() ENTERIT(long double)
-#define ENTERIT(returntype) \
- returntype __retval; \
+#define ENTERI() \
+ long double __retval; \
fp_prec_t __oprec; \
\
if ((__oprec = fpgetprec()) != FP_PE) \
@@ -360,7 +318,6 @@
} while (0)
#else
#define ENTERI()
-#define ENTERIT(x)
#define RETURNI(x) RETURNF(x)
#define ENTERV()
#define RETURNV() return
@@ -478,24 +435,6 @@
*/
void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
-/*
- * Mix 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 addition. The result is typically the arg whose mantissa
- * bits (considered as in integer) are largest.
- *
- * Technical complications: the result in bits might depend on the 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 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) (((x) + 0.0L) + ((y) + 0))
-
#ifdef _COMPLEX_H
/*
@@ -846,7 +785,4 @@
long double __kernel_cosl(long double, long double);
long double __kernel_tanl(long double, long double, int);
-long double __p1evll(long double, void *, int);
-long double __polevll(long double, void *, int);
-
#endif /* !_MATH_PRIVATE_H_ */
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ccosh.c b/libm/upstream-freebsd/lib/msun/src/s_ccosh.c
index da568f4..534f4b6 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ccosh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ccosh.c
@@ -39,7 +39,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_ccosh.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_ccosh.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
@@ -146,8 +146,7 @@
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
- return (CMPLX(((long double)x * x) * (y - y),
- ((long double)x + x) * (y - y)));
+ return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
}
double complex
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ccoshf.c b/libm/upstream-freebsd/lib/msun/src/s_ccoshf.c
index 4aeea63..3f4423e 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ccoshf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ccoshf.c
@@ -31,7 +31,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_ccoshf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_ccoshf.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
@@ -92,8 +92,7 @@
return (CMPLXF(INFINITY * cosf(y), x * sinf(y)));
}
- return (CMPLXF(((long double)x * x) * (y - y),
- ((long double)x + x) * (y - y)));
+ return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
}
float complex
diff --git a/libm/upstream-freebsd/lib/msun/src/s_clog.c b/libm/upstream-freebsd/lib/msun/src/s_clog.c
deleted file mode 100644
index 2c31095..0000000
--- a/libm/upstream-freebsd/lib/msun/src/s_clog.c
+++ /dev/null
@@ -1,155 +0,0 @@
-/*-
- * Copyright (c) 2013 Bruce D. Evans
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice unmodified, this list of conditions, and the following
- * disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
- * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
- * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_clog.c 333577 2018-05-13 09:54:34Z kib $");
-
-#include <complex.h>
-#include <float.h>
-
-#include "fpmath.h"
-#include "math.h"
-#include "math_private.h"
-
-#define MANT_DIG DBL_MANT_DIG
-#define MAX_EXP DBL_MAX_EXP
-#define MIN_EXP DBL_MIN_EXP
-
-static const double
-ln2_hi = 6.9314718055829871e-1, /* 0x162e42fefa0000.0p-53 */
-ln2_lo = 1.6465949582897082e-12; /* 0x1cf79abc9e3b3a.0p-92 */
-
-double complex
-clog(double complex z)
-{
- double_t ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl, sh, sl, t;
- double x, y, v;
- uint32_t hax, hay;
- int kx, ky;
-
- x = creal(z);
- y = cimag(z);
- v = atan2(y, x);
-
- ax = fabs(x);
- ay = fabs(y);
- if (ax < ay) {
- t = ax;
- ax = ay;
- ay = t;
- }
-
- GET_HIGH_WORD(hax, ax);
- kx = (hax >> 20) - 1023;
- GET_HIGH_WORD(hay, ay);
- ky = (hay >> 20) - 1023;
-
- /* Handle NaNs and Infs using the general formula. */
- if (kx == MAX_EXP || ky == MAX_EXP)
- return (CMPLX(log(hypot(x, y)), v));
-
- /* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
- if (ax == 1) {
- if (ky < (MIN_EXP - 1) / 2)
- return (CMPLX((ay / 2) * ay, v));
- return (CMPLX(log1p(ay * ay) / 2, v));
- }
-
- /* Avoid underflow when ax is not small. Also handle zero args. */
- if (kx - ky > MANT_DIG || ay == 0)
- return (CMPLX(log(ax), v));
-
- /* Avoid overflow. */
- if (kx >= MAX_EXP - 1)
- return (CMPLX(log(hypot(x * 0x1p-1022, y * 0x1p-1022)) +
- (MAX_EXP - 2) * ln2_lo + (MAX_EXP - 2) * ln2_hi, v));
- if (kx >= (MAX_EXP - 1) / 2)
- return (CMPLX(log(hypot(x, y)), v));
-
- /* Reduce inaccuracies and avoid underflow when ax is denormal. */
- if (kx <= MIN_EXP - 2)
- return (CMPLX(log(hypot(x * 0x1p1023, y * 0x1p1023)) +
- (MIN_EXP - 2) * ln2_lo + (MIN_EXP - 2) * ln2_hi, v));
-
- /* Avoid remaining underflows (when ax is small but not denormal). */
- if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
- return (CMPLX(log(hypot(x, y)), v));
-
- /* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
- t = (double)(ax * (0x1p27 + 1));
- axh = (double)(ax - t) + t;
- axl = ax - axh;
- ax2h = ax * ax;
- ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
- t = (double)(ay * (0x1p27 + 1));
- ayh = (double)(ay - t) + t;
- ayl = ay - ayh;
- ay2h = ay * ay;
- ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
-
- /*
- * When log(|z|) is far from 1, accuracy in calculating the sum
- * of the squares is not very important since log() reduces
- * inaccuracies. We depended on this to use the general
- * formula when log(|z|) is very far from 1. When log(|z|) is
- * moderately far from 1, we go through the extra-precision
- * calculations to reduce branches and gain a little accuracy.
- *
- * When |z| is near 1, we subtract 1 and use log1p() and don't
- * leave it to log() to subtract 1, since we gain at least 1 bit
- * of accuracy in this way.
- *
- * When |z| is very near 1, subtracting 1 can cancel almost
- * 3*MANT_DIG bits. We arrange that subtracting 1 is exact in
- * doubled precision, and then do the rest of the calculation
- * in sloppy doubled precision. Although large cancellations
- * often lose lots of accuracy, here the final result is exact
- * in doubled precision if the large calculation occurs (because
- * then it is exact in tripled precision and the cancellation
- * removes enough bits to fit in doubled precision). Thus the
- * result is accurate in sloppy doubled precision, and the only
- * significant loss of accuracy is when it is summed and passed
- * to log1p().
- */
- sh = ax2h;
- sl = ay2h;
- _2sumF(sh, sl);
- if (sh < 0.5 || sh >= 3)
- return (CMPLX(log(ay2l + ax2l + sl + sh) / 2, v));
- sh -= 1;
- _2sum(sh, sl);
- _2sum(ax2l, ay2l);
- /* Briggs-Kahan algorithm (except we discard the final low term): */
- _2sum(sh, ax2l);
- _2sum(sl, ay2l);
- t = ax2l + sl;
- _2sumF(sh, t);
- return (CMPLX(log1p(ay2l + t + sh) / 2, v));
-}
-
-#if (LDBL_MANT_DIG == 53)
-__weak_reference(clog, clogl);
-#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_clogf.c b/libm/upstream-freebsd/lib/msun/src/s_clogf.c
deleted file mode 100644
index d94977d..0000000
--- a/libm/upstream-freebsd/lib/msun/src/s_clogf.c
+++ /dev/null
@@ -1,151 +0,0 @@
-/*-
- * Copyright (c) 2013 Bruce D. Evans
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice unmodified, this list of conditions, and the following
- * disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
- * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
- * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_clogf.c 333577 2018-05-13 09:54:34Z kib $");
-
-#include <complex.h>
-#include <float.h>
-
-#include "fpmath.h"
-#include "math.h"
-#include "math_private.h"
-
-#define MANT_DIG FLT_MANT_DIG
-#define MAX_EXP FLT_MAX_EXP
-#define MIN_EXP FLT_MIN_EXP
-
-static const float
-ln2f_hi = 6.9314575195e-1, /* 0xb17200.0p-24 */
-ln2f_lo = 1.4286067653e-6; /* 0xbfbe8e.0p-43 */
-
-float complex
-clogf(float complex z)
-{
- float_t ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl, sh, sl, t;
- float x, y, v;
- uint32_t hax, hay;
- int kx, ky;
-
- x = crealf(z);
- y = cimagf(z);
- v = atan2f(y, x);
-
- ax = fabsf(x);
- ay = fabsf(y);
- if (ax < ay) {
- t = ax;
- ax = ay;
- ay = t;
- }
-
- GET_FLOAT_WORD(hax, ax);
- kx = (hax >> 23) - 127;
- GET_FLOAT_WORD(hay, ay);
- ky = (hay >> 23) - 127;
-
- /* Handle NaNs and Infs using the general formula. */
- if (kx == MAX_EXP || ky == MAX_EXP)
- return (CMPLXF(logf(hypotf(x, y)), v));
-
- /* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
- if (hax == 0x3f800000) {
- if (ky < (MIN_EXP - 1) / 2)
- return (CMPLXF((ay / 2) * ay, v));
- return (CMPLXF(log1pf(ay * ay) / 2, v));
- }
-
- /* Avoid underflow when ax is not small. Also handle zero args. */
- if (kx - ky > MANT_DIG || hay == 0)
- return (CMPLXF(logf(ax), v));
-
- /* Avoid overflow. */
- if (kx >= MAX_EXP - 1)
- return (CMPLXF(logf(hypotf(x * 0x1p-126F, y * 0x1p-126F)) +
- (MAX_EXP - 2) * ln2f_lo + (MAX_EXP - 2) * ln2f_hi, v));
- if (kx >= (MAX_EXP - 1) / 2)
- return (CMPLXF(logf(hypotf(x, y)), v));
-
- /* Reduce inaccuracies and avoid underflow when ax is denormal. */
- if (kx <= MIN_EXP - 2)
- return (CMPLXF(logf(hypotf(x * 0x1p127F, y * 0x1p127F)) +
- (MIN_EXP - 2) * ln2f_lo + (MIN_EXP - 2) * ln2f_hi, v));
-
- /* Avoid remaining underflows (when ax is small but not denormal). */
- if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
- return (CMPLXF(logf(hypotf(x, y)), v));
-
- /* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
- t = (float)(ax * (0x1p12F + 1));
- axh = (float)(ax - t) + t;
- axl = ax - axh;
- ax2h = ax * ax;
- ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
- t = (float)(ay * (0x1p12F + 1));
- ayh = (float)(ay - t) + t;
- ayl = ay - ayh;
- ay2h = ay * ay;
- ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
-
- /*
- * When log(|z|) is far from 1, accuracy in calculating the sum
- * of the squares is not very important since log() reduces
- * inaccuracies. We depended on this to use the general
- * formula when log(|z|) is very far from 1. When log(|z|) is
- * moderately far from 1, we go through the extra-precision
- * calculations to reduce branches and gain a little accuracy.
- *
- * When |z| is near 1, we subtract 1 and use log1p() and don't
- * leave it to log() to subtract 1, since we gain at least 1 bit
- * of accuracy in this way.
- *
- * When |z| is very near 1, subtracting 1 can cancel almost
- * 3*MANT_DIG bits. We arrange that subtracting 1 is exact in
- * doubled precision, and then do the rest of the calculation
- * in sloppy doubled precision. Although large cancellations
- * often lose lots of accuracy, here the final result is exact
- * in doubled precision if the large calculation occurs (because
- * then it is exact in tripled precision and the cancellation
- * removes enough bits to fit in doubled precision). Thus the
- * result is accurate in sloppy doubled precision, and the only
- * significant loss of accuracy is when it is summed and passed
- * to log1p().
- */
- sh = ax2h;
- sl = ay2h;
- _2sumF(sh, sl);
- if (sh < 0.5F || sh >= 3)
- return (CMPLXF(logf(ay2l + ax2l + sl + sh) / 2, v));
- sh -= 1;
- _2sum(sh, sl);
- _2sum(ax2l, ay2l);
- /* Briggs-Kahan algorithm (except we discard the final low term): */
- _2sum(sh, ax2l);
- _2sum(sl, ay2l);
- t = ax2l + sl;
- _2sumF(sh, t);
- return (CMPLXF(log1pf(ay2l + t + sh) / 2, v));
-}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_clogl.c b/libm/upstream-freebsd/lib/msun/src/s_clogl.c
deleted file mode 100644
index 9f56ca3..0000000
--- a/libm/upstream-freebsd/lib/msun/src/s_clogl.c
+++ /dev/null
@@ -1,168 +0,0 @@
-/*-
- * Copyright (c) 2013 Bruce D. Evans
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice unmodified, this list of conditions, and the following
- * disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
- * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
- * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_clogl.c 333577 2018-05-13 09:54:34Z kib $");
-
-#include <complex.h>
-#include <float.h>
-#ifdef __i386__
-#include <ieeefp.h>
-#endif
-
-#include "fpmath.h"
-#include "math.h"
-#include "math_private.h"
-
-#define MANT_DIG LDBL_MANT_DIG
-#define MAX_EXP LDBL_MAX_EXP
-#define MIN_EXP LDBL_MIN_EXP
-
-static const double
-ln2_hi = 6.9314718055829871e-1; /* 0x162e42fefa0000.0p-53 */
-
-#if LDBL_MANT_DIG == 64
-#define MULT_REDUX 0x1p32 /* exponent MANT_DIG / 2 rounded up */
-static const double
-ln2l_lo = 1.6465949582897082e-12; /* 0x1cf79abc9e3b3a.0p-92 */
-#elif LDBL_MANT_DIG == 113
-#define MULT_REDUX 0x1p57
-static const long double
-ln2l_lo = 1.64659495828970812809844307550013433e-12L; /* 0x1cf79abc9e3b39803f2f6af40f343.0p-152L */
-#else
-#error "Unsupported long double format"
-#endif
-
-long double complex
-clogl(long double complex z)
-{
- long double ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl;
- long double sh, sl, t;
- long double x, y, v;
- uint16_t hax, hay;
- int kx, ky;
-
- ENTERIT(long double complex);
-
- x = creall(z);
- y = cimagl(z);
- v = atan2l(y, x);
-
- ax = fabsl(x);
- ay = fabsl(y);
- if (ax < ay) {
- t = ax;
- ax = ay;
- ay = t;
- }
-
- GET_LDBL_EXPSIGN(hax, ax);
- kx = hax - 16383;
- GET_LDBL_EXPSIGN(hay, ay);
- ky = hay - 16383;
-
- /* Handle NaNs and Infs using the general formula. */
- if (kx == MAX_EXP || ky == MAX_EXP)
- RETURNI(CMPLXL(logl(hypotl(x, y)), v));
-
- /* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
- if (ax == 1) {
- if (ky < (MIN_EXP - 1) / 2)
- RETURNI(CMPLXL((ay / 2) * ay, v));
- RETURNI(CMPLXL(log1pl(ay * ay) / 2, v));
- }
-
- /* Avoid underflow when ax is not small. Also handle zero args. */
- if (kx - ky > MANT_DIG || ay == 0)
- RETURNI(CMPLXL(logl(ax), v));
-
- /* Avoid overflow. */
- if (kx >= MAX_EXP - 1)
- RETURNI(CMPLXL(logl(hypotl(x * 0x1p-16382L, y * 0x1p-16382L)) +
- (MAX_EXP - 2) * ln2l_lo + (MAX_EXP - 2) * ln2_hi, v));
- if (kx >= (MAX_EXP - 1) / 2)
- RETURNI(CMPLXL(logl(hypotl(x, y)), v));
-
- /* Reduce inaccuracies and avoid underflow when ax is denormal. */
- if (kx <= MIN_EXP - 2)
- RETURNI(CMPLXL(logl(hypotl(x * 0x1p16383L, y * 0x1p16383L)) +
- (MIN_EXP - 2) * ln2l_lo + (MIN_EXP - 2) * ln2_hi, v));
-
- /* Avoid remaining underflows (when ax is small but not denormal). */
- if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
- RETURNI(CMPLXL(logl(hypotl(x, y)), v));
-
- /* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
- t = (long double)(ax * (MULT_REDUX + 1));
- axh = (long double)(ax - t) + t;
- axl = ax - axh;
- ax2h = ax * ax;
- ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
- t = (long double)(ay * (MULT_REDUX + 1));
- ayh = (long double)(ay - t) + t;
- ayl = ay - ayh;
- ay2h = ay * ay;
- ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
-
- /*
- * When log(|z|) is far from 1, accuracy in calculating the sum
- * of the squares is not very important since log() reduces
- * inaccuracies. We depended on this to use the general
- * formula when log(|z|) is very far from 1. When log(|z|) is
- * moderately far from 1, we go through the extra-precision
- * calculations to reduce branches and gain a little accuracy.
- *
- * When |z| is near 1, we subtract 1 and use log1p() and don't
- * leave it to log() to subtract 1, since we gain at least 1 bit
- * of accuracy in this way.
- *
- * When |z| is very near 1, subtracting 1 can cancel almost
- * 3*MANT_DIG bits. We arrange that subtracting 1 is exact in
- * doubled precision, and then do the rest of the calculation
- * in sloppy doubled precision. Although large cancellations
- * often lose lots of accuracy, here the final result is exact
- * in doubled precision if the large calculation occurs (because
- * then it is exact in tripled precision and the cancellation
- * removes enough bits to fit in doubled precision). Thus the
- * result is accurate in sloppy doubled precision, and the only
- * significant loss of accuracy is when it is summed and passed
- * to log1p().
- */
- sh = ax2h;
- sl = ay2h;
- _2sumF(sh, sl);
- if (sh < 0.5 || sh >= 3)
- RETURNI(CMPLXL(logl(ay2l + ax2l + sl + sh) / 2, v));
- sh -= 1;
- _2sum(sh, sl);
- _2sum(ax2l, ay2l);
- /* Briggs-Kahan algorithm (except we discard the final low term): */
- _2sum(sh, ax2l);
- _2sum(sl, ay2l);
- t = ax2l + sl;
- _2sumF(sh, t);
- RETURNI(CMPLXL(log1pl(ay2l + t + sh) / 2, v));
-}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cpow.c b/libm/upstream-freebsd/lib/msun/src/s_cpow.c
deleted file mode 100644
index 9be5c51..0000000
--- a/libm/upstream-freebsd/lib/msun/src/s_cpow.c
+++ /dev/null
@@ -1,74 +0,0 @@
-/*-
- * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- */
-
-/* cpow
- *
- * Complex power function
- *
- *
- *
- * SYNOPSIS:
- *
- * double complex cpow();
- * double complex a, z, w;
- *
- * w = cpow (a, z);
- *
- *
- *
- * DESCRIPTION:
- *
- * Raises complex A to the complex Zth power.
- * Definition is per AMS55 # 4.2.8,
- * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE -10,+10 30000 9.4e-15 1.5e-15
- *
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_cpow.c 336299 2018-07-15 00:23:10Z mmacy $");
-
-#include <complex.h>
-#include <float.h>
-#include <math.h>
-
-double complex
-cpow(double complex a, double complex z)
-{
- double complex w;
- double x, y, r, theta, absa, arga;
-
- x = creal (z);
- y = cimag (z);
- absa = cabs (a);
- if (absa == 0.0) {
- return (0.0 + 0.0 * I);
- }
- arga = carg (a);
- r = pow (absa, x);
- theta = x * arga;
- if (y != 0.0) {
- r = r * exp (-y * arga);
- theta = theta + y * log (absa);
- }
- w = r * cos (theta) + (r * sin (theta)) * I;
- return (w);
-}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cpowf.c b/libm/upstream-freebsd/lib/msun/src/s_cpowf.c
deleted file mode 100644
index 6e27f57..0000000
--- a/libm/upstream-freebsd/lib/msun/src/s_cpowf.c
+++ /dev/null
@@ -1,73 +0,0 @@
-/*-
- * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- */
-
-/* cpowf
- *
- * Complex power function
- *
- *
- *
- * SYNOPSIS:
- *
- * float complex cpowf();
- * float complex a, z, w;
- *
- * w = cpowf (a, z);
- *
- *
- *
- * DESCRIPTION:
- *
- * Raises complex A to the complex Zth power.
- * Definition is per AMS55 # 4.2.8,
- * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE -10,+10 30000 9.4e-15 1.5e-15
- *
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_cpowf.c 336299 2018-07-15 00:23:10Z mmacy $");
-
-#include <complex.h>
-#include <math.h>
-
-float complex
-cpowf(float complex a, float complex z)
-{
- float complex w;
- float x, y, r, theta, absa, arga;
-
- x = crealf(z);
- y = cimagf(z);
- absa = cabsf (a);
- if (absa == 0.0f) {
- return (0.0f + 0.0f * I);
- }
- arga = cargf (a);
- r = powf (absa, x);
- theta = x * arga;
- if (y != 0.0f) {
- r = r * expf (-y * arga);
- theta = theta + y * logf (absa);
- }
- w = r * cosf (theta) + (r * sinf (theta)) * I;
- return (w);
-}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cpowl.c b/libm/upstream-freebsd/lib/msun/src/s_cpowl.c
deleted file mode 100644
index fcb5c7f..0000000
--- a/libm/upstream-freebsd/lib/msun/src/s_cpowl.c
+++ /dev/null
@@ -1,73 +0,0 @@
-/*-
- * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- */
-
-/* cpowl
- *
- * Complex power function
- *
- *
- *
- * SYNOPSIS:
- *
- * long double complex cpowl();
- * long double complex a, z, w;
- *
- * w = cpowl (a, z);
- *
- *
- *
- * DESCRIPTION:
- *
- * Raises complex A to the complex Zth power.
- * Definition is per AMS55 # 4.2.8,
- * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE -10,+10 30000 9.4e-15 1.5e-15
- *
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_cpowl.c 336299 2018-07-15 00:23:10Z mmacy $");
-
-#include <complex.h>
-#include <math.h>
-
-long double complex
-cpowl(long double complex a, long double complex z)
-{
- long double complex w;
- long double x, y, r, theta, absa, arga;
-
- x = creall(z);
- y = cimagl(z);
- absa = cabsl(a);
- if (absa == 0.0L) {
- return (0.0L + 0.0L * I);
- }
- arga = cargl(a);
- r = powl(absa, x);
- theta = x * arga;
- if (y != 0.0L) {
- r = r * expl(-y * arga);
- theta = theta + y * logl(absa);
- }
- w = r * cosl(theta) + (r * sinl(theta)) * I;
- return (w);
-}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_csinh.c b/libm/upstream-freebsd/lib/msun/src/s_csinh.c
index 7408b13..5f8380e 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csinh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csinh.c
@@ -39,7 +39,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csinh.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csinh.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
@@ -145,8 +145,7 @@
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
- return (CMPLX(((long double)x + x) * (y - y),
- ((long double)x * x) * (y - y)));
+ return (CMPLX((x + x) * (y - y), (x * x) * (y - y)));
}
double complex
diff --git a/libm/upstream-freebsd/lib/msun/src/s_csinhf.c b/libm/upstream-freebsd/lib/msun/src/s_csinhf.c
index 6e29f19..7f58b55 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csinhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csinhf.c
@@ -31,7 +31,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csinhf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csinhf.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
@@ -92,8 +92,7 @@
return (CMPLXF(x * cosf(y), INFINITY * sinf(y)));
}
- return (CMPLXF(((long double)x + x) * (y - y),
- ((long double)x * x) * (y - y)));
+ return (CMPLXF((x + x) * (y - y), (x * x) * (y - y)));
}
float complex
diff --git a/libm/upstream-freebsd/lib/msun/src/s_csqrt.c b/libm/upstream-freebsd/lib/msun/src/s_csqrt.c
index 215b7d7..ad8ef4c 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csqrt.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csqrt.c
@@ -27,7 +27,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrt.c 336488 2018-07-19 15:04:10Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrt.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <float.h>
@@ -35,14 +35,25 @@
#include "math_private.h"
-/* For avoiding overflow for components >= DBL_MAX / (1 + sqrt(2)). */
+/*
+ * gcc doesn't implement complex multiplication or division correctly,
+ * so we need to handle infinities specially. We turn on this pragma to
+ * notify conforming c99 compilers that the fast-but-incorrect code that
+ * gcc generates is acceptable, since the special cases have already been
+ * handled.
+ */
+#pragma STDC CX_LIMITED_RANGE ON
+
+/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
#define THRESH 0x1.a827999fcef32p+1022
double complex
csqrt(double complex z)
{
double complex result;
- double a, b, rx, ry, scale, t;
+ double a, b;
+ double t;
+ int scale;
a = creal(z);
b = cimag(z);
@@ -54,7 +65,7 @@
return (CMPLX(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
- return (CMPLX(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
+ return (CMPLX(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@@ -68,46 +79,34 @@
else
return (CMPLX(a, copysign(b - b, b)));
}
- if (isnan(b)) {
- t = (a - a) / (a - a); /* raise invalid */
- return (CMPLX(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
- }
+ /*
+ * The remaining special case (b is NaN) is handled just fine by
+ * the normal code path below.
+ */
/* Scale to avoid overflow. */
if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
- /*
- * Don't scale a or b if this might give (spurious)
- * underflow. Then the unscaled value is an equivalent
- * infinitesmal (or 0).
- */
- if (fabs(a) >= 0x1p-1020)
- a *= 0.25;
- if (fabs(b) >= 0x1p-1020)
- b *= 0.25;
- scale = 2;
- } else {
+ a *= 0.25;
+ b *= 0.25;
scale = 1;
- }
-
- /* Scale to reduce inaccuracies when both components are denormal. */
- if (fabs(a) < 0x1p-1022 && fabs(b) < 0x1p-1022) {
- a *= 0x1p54;
- b *= 0x1p54;
- scale = 0x1p-27;
+ } else {
+ scale = 0;
}
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
- rx = scale * t;
- ry = scale * b / (2 * t);
+ result = CMPLX(t, b / (2 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
- rx = scale * fabs(b) / (2 * t);
- ry = copysign(scale * t, b);
+ result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
}
- return (CMPLX(rx, ry));
+ /* Rescale. */
+ if (scale)
+ return (result * 2);
+ else
+ return (result);
}
#if LDBL_MANT_DIG == 53
diff --git a/libm/upstream-freebsd/lib/msun/src/s_csqrtf.c b/libm/upstream-freebsd/lib/msun/src/s_csqrtf.c
index 16d3266..195e3b7 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csqrtf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csqrtf.c
@@ -27,21 +27,27 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtf.c 336412 2018-07-17 12:01:59Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtf.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
#include "math_private.h"
+/*
+ * gcc doesn't implement complex multiplication or division correctly,
+ * so we need to handle infinities specially. We turn on this pragma to
+ * notify conforming c99 compilers that the fast-but-incorrect code that
+ * gcc generates is acceptable, since the special cases have already been
+ * handled.
+ */
+#pragma STDC CX_LIMITED_RANGE ON
+
float complex
csqrtf(float complex z)
{
+ float a = crealf(z), b = cimagf(z);
double t;
- float a, b;
-
- a = creal(z);
- b = cimag(z);
/* Handle special cases. */
if (z == 0)
@@ -50,7 +56,7 @@
return (CMPLXF(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
- return (CMPLXF(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
+ return (CMPLXF(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@@ -64,10 +70,10 @@
else
return (CMPLXF(a, copysignf(b - b, b)));
}
- if (isnan(b)) {
- t = (a - a) / (a - a); /* raise invalid */
- return (CMPLXF(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
- }
+ /*
+ * The remaining special case (b is NaN) is handled just fine by
+ * the normal code path below.
+ */
/*
* We compute t in double precision to avoid overflow and to
@@ -76,9 +82,9 @@
*/
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
- return (CMPLXF(t, b / (2 * t)));
+ return (CMPLXF(t, b / (2.0 * t)));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
- return (CMPLXF(fabsf(b) / (2 * t), copysignf(t, b)));
+ return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)));
}
}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c b/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
index e58e4ae..a975833 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
@@ -27,7 +27,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtl.c 336488 2018-07-19 15:04:10Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtl.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <float.h>
@@ -36,28 +36,32 @@
#include "math_private.h"
/*
- * Several thresholds require a 15-bit exponent and also the usual bias.
- * s_logl.c and e_hypotl have less hard-coding but end up requiring the
- * same for the exponent and more for the mantissa.
+ * gcc doesn't implement complex multiplication or division correctly,
+ * so we need to handle infinities specially. We turn on this pragma to
+ * notify conforming c99 compilers that the fast-but-incorrect code that
+ * gcc generates is acceptable, since the special cases have already been
+ * handled.
+ */
+#pragma STDC CX_LIMITED_RANGE ON
+
+/*
+ * We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)).
+ * Rather than determining the fully precise value at which we might
+ * overflow, just use a threshold of approximately LDBL_MAX / 4.
*/
#if LDBL_MAX_EXP != 0x4000
#error "Unsupported long double format"
-#endif
-
-/*
- * Overflow must be avoided for components >= LDBL_MAX / (1 + sqrt(2)).
- * The precise threshold is nontrivial to determine and spell, so use a
- * lower threshold of approximaely LDBL_MAX / 4, and don't use LDBL_MAX
- * to spell this since LDBL_MAX is broken on i386 (it overflows in 53-bit
- * precision).
- */
+#else
#define THRESH 0x1p16382L
+#endif
long double complex
csqrtl(long double complex z)
{
long double complex result;
- long double a, b, rx, ry, scale, t;
+ long double a, b;
+ long double t;
+ int scale;
a = creall(z);
b = cimagl(z);
@@ -69,7 +73,7 @@
return (CMPLXL(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
- return (CMPLXL(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
+ return (CMPLXL(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@@ -83,44 +87,32 @@
else
return (CMPLXL(a, copysignl(b - b, b)));
}
- if (isnan(b)) {
- t = (a - a) / (a - a); /* raise invalid */
- return (CMPLXL(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
- }
+ /*
+ * The remaining special case (b is NaN) is handled just fine by
+ * the normal code path below.
+ */
/* Scale to avoid overflow. */
if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {
- /*
- * Don't scale a or b if this might give (spurious)
- * underflow. Then the unscaled value is an equivalent
- * infinitesmal (or 0).
- */
- if (fabsl(a) >= 0x1p-16380L)
- a *= 0.25;
- if (fabsl(b) >= 0x1p-16380L)
- b *= 0.25;
- scale = 2;
- } else {
+ a *= 0.25;
+ b *= 0.25;
scale = 1;
- }
-
- /* Scale to reduce inaccuracies when both components are denormal. */
- if (fabsl(a) < 0x1p-16382L && fabsl(b) < 0x1p-16382L) {
- a *= 0x1p64;
- b *= 0x1p64;
- scale = 0x1p-32;
+ } else {
+ scale = 0;
}
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrtl((a + hypotl(a, b)) * 0.5);
- rx = scale * t;
- ry = scale * b / (2 * t);
+ result = CMPLXL(t, b / (2 * t));
} else {
t = sqrtl((-a + hypotl(a, b)) * 0.5);
- rx = scale * fabsl(b) / (2 * t);
- ry = copysignl(scale * t, b);
+ result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b));
}
- return (CMPLXL(rx, ry));
+ /* Rescale. */
+ if (scale)
+ return (result * 2);
+ else
+ return (result);
}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ctanh.c b/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
index d005baf..bd16109 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
@@ -66,7 +66,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanh.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanh.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
@@ -104,8 +104,8 @@
*/
if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */
- return (CMPLX(nan_mix(x, y),
- y == 0 ? y : nan_mix(x, y)));
+ return (CMPLX((x + 0) * (y + 0),
+ y == 0 ? y : (x + 0) * (y + 0)));
SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c b/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
index 06825a3..1124fc2 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
@@ -31,7 +31,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanhf.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanhf.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
#include <math.h>
@@ -53,8 +53,8 @@
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
- return (CMPLXF(nan_mix(x, y),
- y == 0 ? y : nan_mix(x, y)));
+ return (CMPLXF((x + 0) * (y + 0),
+ y == 0 ? y : (x + 0) * (y + 0)));
SET_FLOAT_WORD(x, hx - 0x40000000);
return (CMPLXF(x,
copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))));
diff --git a/libm/upstream-freebsd/lib/msun/src/s_remquo.c b/libm/upstream-freebsd/lib/msun/src/s_remquo.c
index d007d1f..d811c69 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_remquo.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_remquo.c
@@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_remquo.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include <float.h>
@@ -44,7 +44,7 @@
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
- return nan_mix(x, y)/nan_mix(x, y);
+ return (x*y)/(x*y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) {
q = 0;
diff --git a/libm/upstream-freebsd/lib/msun/src/s_remquof.c b/libm/upstream-freebsd/lib/msun/src/s_remquof.c
index f376b13..f7b4c00 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_remquof.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_remquof.c
@@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_remquof.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
@@ -41,7 +41,7 @@
/* purge off exception values */
if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
- return nan_mix(x, y)/nan_mix(x, y);
+ return (x*y)/(x*y);
if(hx<hy) {
q = 0;
goto fixup; /* |x|<|y| return x or x-y */
diff --git a/libm/upstream-freebsd/lib/msun/src/s_remquol.c b/libm/upstream-freebsd/lib/msun/src/s_remquol.c
index 77ebe43..712651c 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_remquol.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_remquol.c
@@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_remquol.c 336362 2018-07-17 07:42:14Z bde $");
+__FBSDID("$FreeBSD$");
#include <float.h>
#include <stdint.h>
@@ -86,7 +86,7 @@
(ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(uy.bits.exp == BIAS + LDBL_MAX_EXP &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
- return nan_mix(x, y)/nan_mix(x, y);
+ return (x*y)/(x*y);
if(ux.bits.exp<=uy.bits.exp) {
if((ux.bits.exp<uy.bits.exp) ||
(ux.bits.manh<=uy.bits.manh &&