Merge "libm: sync with upstream FreeBSD."

This commit is contained in:
Elliott Hughes 2019-10-25 17:19:57 +00:00 committed by Gerrit Code Review
commit b21d968281
10 changed files with 39 additions and 45 deletions

View file

@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_exp.c 352710 2019-09-25 18:50:57Z dim $");
/* __ieee754_exp(x)
* Returns the exponential of x.
@ -145,9 +145,9 @@ __ieee754_exp(double x) /* default IEEE double exp */
/* x is now in primary range */
t = x*x;
if(k >= -1021)
INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0);
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0);
else
INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0);
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+(k+1000)))<<20, 0);
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);

View file

@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_expf.c 352710 2019-09-25 18:50:57Z dim $");
#include <float.h>
@ -83,9 +83,9 @@ __ieee754_expf(float x)
/* x is now in primary range */
t = x*x;
if(k >= -125)
SET_FLOAT_WORD(twopk,0x3f800000+(k<<23));
SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+k))<<23);
else
SET_FLOAT_WORD(twopk,0x3f800000+((k+100)<<23));
SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+(k+100)))<<23);
c = x - t*(P1+t*P2);
if(k==0) return one-((x*c)/(c-(float)2.0)-x);
else y = one-((lo-(x*c)/((float)2.0-c))-hi);

View file

@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 343023 2019-01-14 15:48:35Z pfg $");
__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 343953 2019-02-10 08:46:07Z peterj $");
/* __ieee754_j0(x), __ieee754_y0(x)
* Bessel function of the first and second kinds of order zero.
@ -93,8 +93,7 @@ __ieee754_j0(double x)
if(ix>=0x7ff00000) return one/(x*x);
x = fabs(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
sincos(x, &s, &c);
ss = s-c;
cc = s+c;
if(ix<0x7fe00000) { /* Make sure x+x does not overflow. */
@ -173,8 +172,7 @@ __ieee754_y0(double x)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
s = sin(x);
c = cos(x);
sincos(x, &s, &c);
ss = s-c;
cc = s+c;
/*

View file

@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD: head/lib/msun/src/e_j0f.c 343023 2019-01-14 15:48:35Z pfg $");
__FBSDID("$FreeBSD: head/lib/msun/src/e_j0f.c 343953 2019-02-10 08:46:07Z peterj $");
/*
* See e_j0.c for complete comments.
@ -55,8 +55,7 @@ __ieee754_j0f(float x)
if(ix>=0x7f800000) return one/(x*x);
x = fabsf(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
sincosf(x, &s, &c);
ss = s-c;
cc = s+c;
if(ix<0x7f000000) { /* Make sure x+x does not overflow. */
@ -128,8 +127,7 @@ __ieee754_y0f(float x)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
s = sinf(x);
c = cosf(x);
sincosf(x, &s, &c);
ss = s-c;
cc = s+c;
/*

View file

@ -11,7 +11,7 @@
*/
#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 343953 2019-02-10 08:46:07Z peterj $");
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
@ -94,8 +94,7 @@ __ieee754_j1(double x)
if(ix>=0x7ff00000) return one/x;
y = fabs(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(y);
c = cos(y);
sincos(y, &s, &c);
ss = -s-c;
cc = s-c;
if(ix<0x7fe00000) { /* make sure y+y not overflow */
@ -159,8 +158,7 @@ __ieee754_y1(double x)
/* y1(x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
sincos(x, &s, &c);
ss = -s-c;
cc = s-c;
if(ix<0x7fe00000) { /* make sure x+x not overflow */

View file

@ -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 343953 2019-02-10 08:46:07Z peterj $");
/*
* See e_j1.c for complete comments.
@ -56,8 +56,7 @@ __ieee754_j1f(float x)
if(ix>=0x7f800000) return one/x;
y = fabsf(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(y);
c = cosf(y);
sincosf(y, &s, &c);
ss = -s-c;
cc = s-c;
if(ix<0x7f000000) { /* make sure y+y not overflow */
@ -114,8 +113,7 @@ __ieee754_y1f(float x)
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
sincosf(x, &s, &c);
ss = -s-c;
cc = s-c;
if(ix<0x7f000000) { /* make sure x+x not overflow */

View file

@ -11,7 +11,7 @@
*/
#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 343953 2019-02-10 08:46:07Z peterj $");
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
@ -54,7 +54,7 @@ double
__ieee754_jn(int n, double x)
{
int32_t i,hx,ix,lx, sgn;
double a, b, temp, di;
double a, b, c, s, temp, di;
double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@ -91,11 +91,12 @@ __ieee754_jn(int n, double x)
* 2 -s+c -c-s
* 3 s+c c-s
*/
sincos(x, &s, &c);
switch(n&3) {
case 0: temp = cos(x)+sin(x); break;
case 1: temp = -cos(x)+sin(x); break;
case 2: temp = -cos(x)-sin(x); break;
case 3: temp = cos(x)-sin(x); break;
case 0: temp = c+s; break;
case 1: temp = -c+s; break;
case 2: temp = -c-s; break;
case 3: temp = c-s; break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
@ -216,7 +217,7 @@ __ieee754_yn(int n, double x)
{
int32_t i,hx,ix,lx;
int32_t sign;
double a, b, temp;
double a, b, c, s, temp;
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
@ -248,11 +249,12 @@ __ieee754_yn(int n, double x)
* 2 -s+c -c-s
* 3 s+c c-s
*/
sincos(x, &s, &c);
switch(n&3) {
case 0: temp = sin(x)-cos(x); break;
case 1: temp = -sin(x)-cos(x); break;
case 2: temp = -sin(x)+cos(x); break;
case 3: temp = sin(x)+cos(x); break;
case 0: temp = s-c; break;
case 1: temp = -s-c; break;
case 2: temp = -s+c; break;
case 3: temp = s+c; break;
}
b = invsqrtpi*temp/sqrt(x);
} else {

View file

@ -11,7 +11,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_expm1.c 352710 2019-09-25 18:50:57Z dim $");
/* expm1(x)
* Returns exp(x)-1, the exponential of x minus 1.
@ -188,7 +188,7 @@ expm1(double x)
e = hxs*((r1-t)/(6.0 - x*t));
if(k==0) return x - (x*e-hxs); /* c is 0 */
else {
INSERT_WORDS(twopk,0x3ff00000+(k<<20),0); /* 2^k */
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20,0); /* 2^k */
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return 0.5*(x-e)-0.5;

View file

@ -14,7 +14,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_expm1f.c 352710 2019-09-25 18:50:57Z dim $");
#include <float.h>
@ -94,7 +94,7 @@ expm1f(float x)
e = hxs*((r1-t)/((float)6.0 - x*t));
if(k==0) return x - (x*e-hxs); /* c is 0 */
else {
SET_FLOAT_WORD(twopk,0x3f800000+(k<<23)); /* 2^k */
SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+k))<<23); /* 2^k */
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return (float)0.5*(x-e)-(float)0.5;

View file

@ -32,7 +32,7 @@
#include <math.h>
#ifndef type
__FBSDID("$FreeBSD: head/lib/msun/src/s_lround.c 326219 2017-11-26 02:00:33Z pfg $");
__FBSDID("$FreeBSD: head/lib/msun/src/s_lround.c 353329 2019-10-08 21:39:51Z brooks $");
#define type double
#define roundit round
#define dtype long
@ -49,9 +49,9 @@ __FBSDID("$FreeBSD: head/lib/msun/src/s_lround.c 326219 2017-11-26 02:00:33Z pfg
* that everything is in range. At compile time, INRANGE(x) should reduce to
* two floating-point comparisons in the former case, or TRUE otherwise.
*/
static const type dtype_min = DTYPE_MIN - 0.5;
static const type dtype_max = DTYPE_MAX + 0.5;
#define INRANGE(x) (dtype_max - DTYPE_MAX != 0.5 || \
static const type dtype_min = (type)DTYPE_MIN - 0.5;
static const type dtype_max = (type)DTYPE_MAX + 0.5;
#define INRANGE(x) (dtype_max - (type)DTYPE_MAX != 0.5 || \
((x) > dtype_min && (x) < dtype_max))
dtype