Attachment #243675 for bug #272783

View | Details | Raw Unified | Return to bug 272783
Collapse All | Expand All

Lines 14-20 Link Here

(-)b/lib/msun/src/e_acos.c (-2 / +2 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_acos(x)	17	/* acos(x)
18	* Method :	18	* Method :
19	* acos(x) = pi/2 - asin(x)	19	* acos(x) = pi/2 - asin(x)
20	* acos(-x) = pi/2 + asin(x)	20	* acos(-x) = pi/2 + asin(x)
Lines 62-68 qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ Link Here
62	qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */	62	qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
63		63
64	double	64	double
65	__ieee754_acos(double x)	65	acos(double x)
66	{	66	{
67	double z,p,q,r,w,s,c,df;	67	double z,p,q,r,w,s,c,df;
68	int32_t hx,ix;	68	int32_t hx,ix;

Lines 32-38 pS2 = -8.6563630030e-03, Link Here

(-)b/lib/msun/src/e_acosf.c (-1 / +1 lines)
32	qS1 = -7.0662963390e-01;	32	qS1 = -7.0662963390e-01;
33		33
34	float	34	float
35	__ieee754_acosf(float x)	35	acosf(float x)
36	{	36	{
37	float z,p,q,r,w,s,c,df;	37	float z,p,q,r,w,s,c,df;
38	int32_t hx,ix;	38	int32_t hx,ix;




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* acosh(x)
 * Method :
 *	Based on 
 *		acosh(x) = log [ x + sqrt(x*x-1) ]

ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */

double
acosh(double x)
{
	double t;
	int32_t hx;

	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
	        return x+x;
	    } else 
		return 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 log(2.0*x-one/(x+sqrt(t-one)));
	} else {			/* 1<x<2 */
	    t = x-one;
	    return log1p(t+sqrt(2.0*t+t*t));




ln2	= 6.9314718246e-01;  /* 0x3f317218 */

float
acoshf(float x)
{
	float t;
	int32_t hx;

	    if(hx >=0x7f800000) {	/* x is inf of NaN */
	        return x+x;
	    } else
		return logf(x)+ln2;	/* acosh(huge)=log(2x) */
	} else if (hx==0x3f800000) {
	    return 0.0;			/* acosh(1) = 0 */
	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
	    t=x*x;
	    return logf((float)2.0*x-one/(x+sqrtf(t-one)));
	} else {			/* 1<x<2 */
	    t = x-one;
	    return log1pf(t+sqrtf((float)2.0*t+t*t));
	}
}

Lines 14-20 Link Here

(-)b/lib/msun/src/e_asin.c (-2 / +2 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_asin(x)	17	/* asin(x)
18	* Method :	18	* Method :
19	* Since asin(x) = x + x^3/6 + x^53/40 + x^715/336 + ...	19	* Since asin(x) = x + x^3/6 + x^53/40 + x^715/336 + ...
20	* we approximate asin(x) on [0,0.5] by	20	* we approximate asin(x) on [0,0.5] by
Lines 68-74 qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ Link Here
68	qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */	68	qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
69		69
70	double	70	double
71	__ieee754_asin(double x)	71	asin(double x)
72	{	72	{
73	double t=0.0,w,p,q,c,r,s;	73	double t=0.0,w,p,q,c,r,s;
74	int32_t hx,ix;	74	int32_t hx,ix;

Lines 32-38 static const double Link Here

(-)b/lib/msun/src/e_asinf.c (-1 / +1 lines)
32	pio2 = 1.570796326794896558e+00;	32	pio2 = 1.570796326794896558e+00;
33		33
34	float	34	float
35	__ieee754_asinf(float x)	35	asinf(float x)
36	{	36	{
37	double s;	37	double s;
38	float t,w,p,q;	38	float t,w,p,q;

Lines 15-21 Link Here

(-)b/lib/msun/src/e_atan2.c (-2 / +2 lines)
15	#include <sys/cdefs.h>	15	#include <sys/cdefs.h>
16	__FBSDID("$FreeBSD$");	16	__FBSDID("$FreeBSD$");
17		17
18	/* __ieee754_atan2(y,x)	18	/* atan2(y,x)
19	* Method :	19	* Method :
20	* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).	20	* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
21	* 2. Reduce x to positive by (if x and y are unexceptional):	21	* 2. Reduce x to positive by (if x and y are unexceptional):
Lines 58-64 static volatile double Link Here
58	pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */	58	pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
59		59
60	double	60	double
61	__ieee754_atan2(double y, double x)	61	atan2(double y, double x)
62	{	62	{
63	double z;	63	double z;
64	int32_t k,m,hx,hy,ix,iy;	64	int32_t k,m,hx,hy,ix,iy;

Lines 30-36 static volatile float Link Here

(-)b/lib/msun/src/e_atan2f.c (-1 / +1 lines)
30	pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */	30	pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
31		31
32	float	32	float
33	__ieee754_atan2f(float y, float x)	33	atan2f(float y, float x)
34	{	34	{
35	float z;	35	float z;
36	int32_t k,m,hx,hy,ix,iy;	36	int32_t k,m,hx,hy,ix,iy;

Lines 15-21 Link Here

(-)b/lib/msun/src/e_atanh.c (-2 / +2 lines)
15	#include <sys/cdefs.h>	15	#include <sys/cdefs.h>
16	__FBSDID("$FreeBSD$");	16	__FBSDID("$FreeBSD$");
17		17
18	/* __ieee754_atanh(x)	18	/* atanh(x)
19	* Method :	19	* Method :
20	* 1.Reduced x to positive by atanh(-x) = -atanh(x)	20	* 1.Reduced x to positive by atanh(-x) = -atanh(x)
21	* 2.For x>=0.5	21	* 2.For x>=0.5
Lines 42-48 static const double one = 1.0, huge = 1e300; Link Here
42	static const double zero = 0.0;	42	static const double zero = 0.0;
43		43
44	double	44	double
45	__ieee754_atanh(double x)	45	atanh(double x)
46	{	46	{
47	double t;	47	double t;
48	int32_t hx,ix;	48	int32_t hx,ix;

Lines 24-30 static const float one = 1.0, huge = 1e30; Link Here

(-)b/lib/msun/src/e_atanhf.c (-1 / +1 lines)
24	static const float zero = 0.0;	24	static const float zero = 0.0;
25		25
26	float	26	float
27	__ieee754_atanhf(float x)	27	atanhf(float x)
28	{	28	{
29	float t;	29	float t;
30	int32_t hx,ix;	30	int32_t hx,ix;




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* cosh(x)
 * Method : 
 * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
 *	1. Replace x by |x| (cosh(x) = cosh(-x)). 

static const double one = 1.0, half=0.5, huge = 1.0e300;

double
cosh(double x)
{
	double t,w;
	int32_t ix;


    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
	if (ix < 0x40360000) {
		t = exp(fabs(x));
		return half*t+half/t;
	}

    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
	if (ix < 0x40862E42)  return half*exp(fabs(x));

    /* |x| in [log(maxdouble), overflowthresold] */
	if (ix<=0x408633CE)




static const float one = 1.0, half=0.5, huge = 1.0e30;

float
coshf(float x)
{
	float t,w;
	int32_t ix;


    /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
	if (ix < 0x41100000) {
		t = expf(fabsf(x));
		return half*t+half/t;
	}

    /* |x| in [9, log(maxfloat)] return half*exp(|x|) */
	if (ix < 0x42b17217)  return half*expf(fabsf(x));

    /* |x| in [log(maxfloat), overflowthresold] */
	if (ix<=0x42b2d4fc)

Lines 13-19 Link Here

(-)b/lib/msun/src/e_exp.c (-2 / +2 lines)
13	#include <sys/cdefs.h>	13	#include <sys/cdefs.h>
14	__FBSDID("$FreeBSD$");	14	__FBSDID("$FreeBSD$");
15		15
16	/* __ieee754_exp(x)	16	/* exp(x)
17	* Returns the exponential of x.	17	* Returns the exponential of x.
18	*	18	*
19	* Method	19	* Method
Lines 102-108 huge = 1.0e+300, Link Here
102	twom1000= 9.33263618503218878990e-302; /* 2*-1000=0x01700000,0/	102	twom1000= 9.33263618503218878990e-302; /* 2*-1000=0x01700000,0/
103		103
104	double	104	double
105	__ieee754_exp(double x) /* default IEEE double exp */	105	exp(double x) /* default IEEE double exp */
106	{	106	{
107	double y,hi=0.0,lo=0.0,c,t,twopk;	107	double y,hi=0.0,lo=0.0,c,t,twopk;
108	int32_t k=0,xsb;	108	int32_t k=0,xsb;

Lines 43-49 huge = 1.0e+30, Link Here

(-)b/lib/msun/src/e_expf.c (-1 / +1 lines)
43	twom100 = 7.8886090522e-31; /* 2*-100=0x0d800000 /	43	twom100 = 7.8886090522e-31; /* 2*-100=0x0d800000 /
44		44
45	float	45	float
46	__ieee754_expf(float x)	46	expf(float x)
47	{	47	{
48	float y,hi=0.0,lo=0.0,c,t,twopk;	48	float y,hi=0.0,lo=0.0,c,t,twopk;
49	int32_t k=0,xsb;	49	int32_t k=0,xsb;

Lines 15-21 Link Here

(-)b/lib/msun/src/e_fmod.c (-2 / +2 lines)
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/*	17	/*
18	* __ieee754_fmod(x,y)	18	* fmod(x,y)
19	* Return x mod y in exact arithmetic	19	* Return x mod y in exact arithmetic
20	* Method: shift and subtract	20	* Method: shift and subtract
21	*/	21	*/
Lines 28-34 __FBSDID("$FreeBSD$"); Link Here
28	static const double one = 1.0, Zero[] = {0.0, -0.0,};	28	static const double one = 1.0, Zero[] = {0.0, -0.0,};
29		29
30	double	30	double
31	__ieee754_fmod(double x, double y)	31	fmod(double x, double y)
32	{	32	{
33	int32_t n,hx,hy,hz,ix,iy,sx,i;	33	int32_t n,hx,hy,hz,ix,iy,sx,i;
34	u_int32_t lx,ly,lz;	34	u_int32_t lx,ly,lz;

Lines 17-23 Link Here

(-)b/lib/msun/src/e_fmodf.c (-2 / +2 lines)
17	__FBSDID("$FreeBSD$");	17	__FBSDID("$FreeBSD$");
18		18
19	/*	19	/*
20	* __ieee754_fmodf(x,y)	20	* fmodf(x,y)
21	* Return x mod y in exact arithmetic	21	* Return x mod y in exact arithmetic
22	* Method: shift and subtract	22	* Method: shift and subtract
23	*/	23	*/
Lines 28-34 __FBSDID("$FreeBSD$"); Link Here
28	static const float one = 1.0, Zero[] = {0.0, -0.0,};	28	static const float one = 1.0, Zero[] = {0.0, -0.0,};
29		29
30	float	30	float
31	__ieee754_fmodf(float x, float y)	31	fmodf(float x, float y)
32	{	32	{
33	int32_t n,hx,hy,hz,ix,iy,sx,i;	33	int32_t n,hx,hy,hz,ix,iy,sx,i;
34		34




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* gamma(x)
 * Return the logarithm of the Gamma function of x.
 *
 * Method: call gamma_r
 */

#include "math.h"

extern int signgam;

double
gamma(double x)
{
	return gamma_r(x,&signgam);
}




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* gamma_r(x, signgamp)
 * Reentrant version of the logarithm of the Gamma function 
 * with user provide pointer for the sign of Gamma(x). 
 *
 * Method: See lgamma_r
 */

#include "math.h"
#include "math_private.h"

double
gamma_r(double x, int *signgamp)
{
	return lgamma_r(x,signgamp);
}




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* gammaf(x)
 * Return the logarithm of the Gamma function of x.
 *
 * Method: call gammaf_r
 */

#include "math.h"

extern int signgam;

float
gammaf(float x)
{
	return gammaf_r(x,&signgam);
}




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* gammaf_r(x, signgamp)
 * Reentrant version of the logarithm of the Gamma function
 * with user provide pointer for the sign of Gamma(x).
 *
 * Method: See lgammaf_r
 */

#include "math.h"
#include "math_private.h"

float
gammaf_r(float x, int *signgamp)
{
	return lgammaf_r(x,signgamp);
}

Lines 14-20 Link Here

(-)b/lib/msun/src/e_hypot.c (-2 / +2 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_hypot(x,y)	17	/* hypot(x,y)
18	*	18	*
19	* Method :	19	* Method :
20	* If (assume round-to-nearest) z=xx+yy	20	* If (assume round-to-nearest) z=xx+yy
Lines 52-58 __FBSDID("$FreeBSD$"); Link Here
52	#include "math_private.h"	52	#include "math_private.h"
53		53
54	double	54	double
55	__ieee754_hypot(double x, double y)	55	hypot(double x, double y)
56	{	56	{
57	double a,b,t1,t2,y1,y2,w;	57	double a,b,t1,t2,y1,y2,w;
58	int32_t j,k,ha,hb;	58	int32_t j,k,ha,hb;

Lines 20-26 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/e_hypotf.c (-3 / +3 lines)
20	#include "math_private.h"	20	#include "math_private.h"
21		21
22	float	22	float
23	__ieee754_hypotf(float x, float y)	23	hypotf(float x, float y)
24	{	24	{
25	float a,b,t1,t2,y1,y2,w;	25	float a,b,t1,t2,y1,y2,w;
26	int32_t j,k,ha,hb;	26	int32_t j,k,ha,hb;
Lines 67-80 __ieee754_hypotf(float x, float y) Link Here
67	if (w>b) {	67	if (w>b) {
68	SET_FLOAT_WORD(t1,ha&0xfffff000);	68	SET_FLOAT_WORD(t1,ha&0xfffff000);
69	t2 = a-t1;	69	t2 = a-t1;
70	w = __ieee754_sqrtf(t1t1-(b(-b)-t2*(a+t1)));	70	w = sqrtf(t1t1-(b(-b)-t2*(a+t1)));
71	} else {	71	} else {
72	a = a+a;	72	a = a+a;
73	SET_FLOAT_WORD(y1,hb&0xfffff000);	73	SET_FLOAT_WORD(y1,hb&0xfffff000);
74	y2 = b - y1;	74	y2 = b - y1;
75	SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000);	75	SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000);
76	t2 = a - t1;	76	t2 = a - t1;
77	w = __ieee754_sqrtf(t1y1-(w(-w)-(t1y2+t2b)));	77	w = sqrtf(t1y1-(w(-w)-(t1y2+t2b)));
78	}	78	}
79	if(k!=0) {	79	if(k!=0) {
80	SET_FLOAT_WORD(t1,(127+k)<<23);	80	SET_FLOAT_WORD(t1,(127+k)<<23);

Lines 13-19 Link Here

(-)b/lib/msun/src/e_j0.c (-5 / +5 lines)
13	#include <sys/cdefs.h>	13	#include <sys/cdefs.h>
14	__FBSDID("$FreeBSD$");	14	__FBSDID("$FreeBSD$");
15		15
16	/* __ieee754_j0(x), __ieee754_y0(x)	16	/* j0(x), y0(x)
17	* Bessel function of the first and second kinds of order zero.	17	* Bessel function of the first and second kinds of order zero.
18	* Method -- j0(x):	18	* Method -- j0(x):
19	* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...	19	* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
Lines 83-89 S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */ Link Here
83	static const double zero = 0, qrtr = 0.25;	83	static const double zero = 0, qrtr = 0.25;
84		84
85	double	85	double
86	__ieee754_j0(double x)	86	j0(double x)
87	{	87	{
88	double z, s,c,ss,cc,r,u,v;	88	double z, s,c,ss,cc,r,u,v;
89	int32_t hx,ix;	89	int32_t hx,ix;
Lines 143-149 v03 = 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */ Link Here
143	v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */	143	v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
144		144
145	double	145	double
146	__ieee754_y0(double x)	146	y0(double x)
147	{	147	{
148	double z, s,c,ss,cc,u,v;	148	double z, s,c,ss,cc,u,v;
149	int32_t hx,ix,lx;	149	int32_t hx,ix,lx;
Lines 192-203 __ieee754_y0(double x) Link Here
192	return z;	192	return z;
193	}	193	}
194	if(ix<=0x3e400000) { /* x < 2*-27 /	194	if(ix<=0x3e400000) { /* x < 2*-27 /
195	return(u00 + tpi*__ieee754_log(x));	195	return(u00 + tpi*log(x));
196	}	196	}
197	z = x*x;	197	z = x*x;
198	u = u00+z(u01+z(u02+z(u03+z(u04+z(u05+zu06)))));	198	u = u00+z(u01+z(u02+z(u03+z(u04+z(u05+zu06)))));
199	v = one+z(v01+z(v02+z(v03+zv04)));	199	v = one+z(v01+z(v02+z(v03+zv04)));
200	return(u/v + tpi(__ieee754_j0(x)__ieee754_log(x)));	200	return(u/v + tpi(j0(x)log(x)));
201	}	201	}
202		202
203	/* The asymptotic expansions of pzero is	203	/* The asymptotic expansions of pzero is

Lines 45-51 S04 = 1.1661400734e-09; /* 0x30a045e8 */ Link Here

(-)b/lib/msun/src/e_j0f.c (-4 / +4 lines)
45	static const float zero = 0, qrtr = 0.25;	45	static const float zero = 0, qrtr = 0.25;
46		46
47	float	47	float
48	__ieee754_j0f(float x)	48	j0f(float x)
49	{	49	{
50	float z, s,c,ss,cc,r,u,v;	50	float z, s,c,ss,cc,r,u,v;
51	int32_t hx,ix;	51	int32_t hx,ix;
Lines 105-111 v03 = 2.5915085189e-07, /* 0x348b216c */ Link Here
105	v04 = 4.4111031494e-10; /* 0x2ff280c2 */	105	v04 = 4.4111031494e-10; /* 0x2ff280c2 */
106		106
107	float	107	float
108	__ieee754_y0f(float x)	108	y0f(float x)
109	{	109	{
110	float z, s,c,ss,cc,u,v;	110	float z, s,c,ss,cc,u,v;
111	int32_t hx,ix;	111	int32_t hx,ix;
Lines 147-158 __ieee754_y0f(float x) Link Here
147	return z;	147	return z;
148	}	148	}
149	if(ix<=0x39000000) { /* x < 2*-13 /	149	if(ix<=0x39000000) { /* x < 2*-13 /
150	return(u00 + tpi*__ieee754_logf(x));	150	return(u00 + tpi*logf(x));
151	}	151	}
152	z = x*x;	152	z = x*x;
153	u = u00+z(u01+z(u02+z(u03+z(u04+z(u05+zu06)))));	153	u = u00+z(u01+z(u02+z(u03+z(u04+z(u05+zu06)))));
154	v = one+z(v01+z(v02+z(v03+zv04)));	154	v = one+z(v01+z(v02+z(v03+zv04)));
155	return(u/v + tpi(__ieee754_j0f(x)__ieee754_logf(x)));	155	return(u/v + tpi(j0f(x)logf(x)));
156	}	156	}
157		157
158	/* The asymptotic expansions of pzero is	158	/* The asymptotic expansions of pzero is

Lines 13-19 Link Here

(-)b/lib/msun/src/e_j1.c (-4 / +4 lines)
13	#include <sys/cdefs.h>	13	#include <sys/cdefs.h>
14	__FBSDID("$FreeBSD$");	14	__FBSDID("$FreeBSD$");
15		15
16	/* __ieee754_j1(x), __ieee754_y1(x)	16	/* j1(x), y1(x)
17	* Bessel function of the first and second kinds of order zero.	17	* Bessel function of the first and second kinds of order zero.
18	* Method -- j1(x):	18	* Method -- j1(x):
19	* 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...	19	* 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
Lines 84-90 s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */ Link Here
84	static const double zero = 0.0;	84	static const double zero = 0.0;
85		85
86	double	86	double
87	__ieee754_j1(double x)	87	j1(double x)
88	{	88	{
89	double z, s,c,ss,cc,r,u,v,y;	89	double z, s,c,ss,cc,r,u,v,y;
90	int32_t hx,ix;	90	int32_t hx,ix;
Lines 140-146 static const double V0[5] = { Link Here
140	};	140	};
141		141
142	double	142	double
143	__ieee754_y1(double x)	143	y1(double x)
144	{	144	{
145	double z, s,c,ss,cc,u,v;	145	double z, s,c,ss,cc,u,v;
146	int32_t hx,ix,lx;	146	int32_t hx,ix,lx;
Lines 190-196 __ieee754_y1(double x) Link Here
190	z = x*x;	190	z = x*x;
191	u = U0[0]+z(U0[1]+z(U0[2]+z(U0[3]+zU0[4])));	191	u = U0[0]+z(U0[1]+z(U0[2]+z(U0[3]+zU0[4])));
192	v = one+z(V0[0]+z(V0[1]+z(V0[2]+z(V0[3]+z*V0[4]))));	192	v = one+z(V0[0]+z(V0[1]+z(V0[2]+z(V0[3]+z*V0[4]))));
193	return(x(u/v) + tpi(__ieee754_j1(x)*__ieee754_log(x)-one/x));	193	return(x(u/v) + tpi(j1(x)*log(x)-one/x));
194	}	194	}
195		195
196	/* For x >= 8, the asymptotic expansions of pone is	196	/* For x >= 8, the asymptotic expansions of pone is

Lines 46-52 s05 = 1.2354227016e-11; /* 0x2d59567e */ Link Here

(-)b/lib/msun/src/e_j1f.c (-3 / +3 lines)
46	static const float zero = 0.0;	46	static const float zero = 0.0;
47		47
48	float	48	float
49	__ieee754_j1f(float x)	49	j1f(float x)
50	{	50	{
51	float z, s,c,ss,cc,r,u,v,y;	51	float z, s,c,ss,cc,r,u,v,y;
52	int32_t hx,ix;	52	int32_t hx,ix;
Lines 102-108 static const float V0[5] = { Link Here
102	};	102	};
103		103
104	float	104	float
105	__ieee754_y1f(float x)	105	y1f(float x)
106	{	106	{
107	float z, s,c,ss,cc,u,v;	107	float z, s,c,ss,cc,u,v;
108	int32_t hx,ix;	108	int32_t hx,ix;
Lines 145-151 __ieee754_y1f(float x) Link Here
145	z = x*x;	145	z = x*x;
146	u = U0[0]+z(U0[1]+z(U0[2]+z(U0[3]+zU0[4])));	146	u = U0[0]+z(U0[1]+z(U0[2]+z(U0[3]+zU0[4])));
147	v = one+z(V0[0]+z(V0[1]+z(V0[2]+z(V0[3]+z*V0[4]))));	147	v = one+z(V0[0]+z(V0[1]+z(V0[2]+z(V0[3]+z*V0[4]))));
148	return(x(u/v) + tpi(__ieee754_j1f(x)*__ieee754_logf(x)-one/x));	148	return(x(u/v) + tpi(j1f(x)*logf(x)-one/x));
149	}	149	}
150		150
151	/* For x >= 8, the asymptotic expansions of pone is	151	/* For x >= 8, the asymptotic expansions of pone is

Lines 14-20 Link Here

(-)b/lib/msun/src/e_jn.c (-14 / +14 lines)
14	__FBSDID("$FreeBSD$");	14	__FBSDID("$FreeBSD$");
15		15
16	/*	16	/*
17	* __ieee754_jn(n, x), __ieee754_yn(n, x)	17	* jn(n, x), yn(n, x)
18	* floating point Bessel's function of the 1st and 2nd kind	18	* floating point Bessel's function of the 1st and 2nd kind
19	* of order n	19	* of order n
20	*	20	*
Lines 51-57 one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */ Link Here
51	static const double zero = 0.00000000000000000000e+00;	51	static const double zero = 0.00000000000000000000e+00;
52		52
53	double	53	double
54	__ieee754_jn(int n, double x)	54	jn(int n, double x)
55	{	55	{
56	int32_t i,hx,ix,lx, sgn;	56	int32_t i,hx,ix,lx, sgn;
57	double a, b, c, s, temp, di;	57	double a, b, c, s, temp, di;
Lines 69-76 __ieee754_jn(int n, double x) Link Here
69	x = -x;	69	x = -x;
70	hx ^= 0x80000000;	70	hx ^= 0x80000000;
71	}	71	}
72	if(n==0) return(__ieee754_j0(x));	72	if(n==0) return(j0(x));
73	if(n==1) return(__ieee754_j1(x));	73	if(n==1) return(j1(x));
74	sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */	74	sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
75	x = fabs(x);	75	x = fabs(x);
76	if((ix\|lx)==0\|\|ix>=0x7ff00000) /* if x is 0 or inf */	76	if((ix\|lx)==0\|\|ix>=0x7ff00000) /* if x is 0 or inf */
Lines 100-107 __ieee754_jn(int n, double x) Link Here
100	}	100	}
101	b = invsqrtpi*temp/sqrt(x);	101	b = invsqrtpi*temp/sqrt(x);
102	} else {	102	} else {
103	a = __ieee754_j0(x);	103	a = j0(x);
104	b = __ieee754_j1(x);	104	b = j1(x);
105	for(i=1;i<n;i++){	105	for(i=1;i<n;i++){
106	temp = b;	106	temp = b;
107	b = b((double)(i+i)/x) - a; / avoid underflow */	107	b = b((double)(i+i)/x) - a; / avoid underflow */
Lines 177-183 __ieee754_jn(int n, double x) Link Here
177	*/	177	*/
178	tmp = n;	178	tmp = n;
179	v = two/x;	179	v = two/x;
180	tmp = tmp__ieee754_log(fabs(vtmp));	180	tmp = tmplog(fabs(vtmp));
181	if(tmp<7.09782712893383973096e+02) {	181	if(tmp<7.09782712893383973096e+02) {
182	for(i=n-1,di=(double)(i+i);i>0;i--){	182	for(i=n-1,di=(double)(i+i);i>0;i--){
183	temp = b;	183	temp = b;
Lines 201-208 __ieee754_jn(int n, double x) Link Here
201	}	201	}
202	}	202	}
203	}	203	}
204	z = __ieee754_j0(x);	204	z = j0(x);
205	w = __ieee754_j1(x);	205	w = j1(x);
206	if (fabs(z) >= fabs(w))	206	if (fabs(z) >= fabs(w))
207	b = (t*z/b);	207	b = (t*z/b);
208	else	208	else
Lines 213-219 __ieee754_jn(int n, double x) Link Here
213	}	213	}
214		214
215	double	215	double
216	__ieee754_yn(int n, double x)	216	yn(int n, double x)
217	{	217	{
218	int32_t i,hx,ix,lx;	218	int32_t i,hx,ix,lx;
219	int32_t sign;	219	int32_t sign;
Lines 232-239 __ieee754_yn(int n, double x) Link Here
232	n = -n;	232	n = -n;
233	sign = 1 - ((n&1)<<1);	233	sign = 1 - ((n&1)<<1);
234	}	234	}
235	if(n==0) return(__ieee754_y0(x));	235	if(n==0) return(y0(x));
236	if(n==1) return(sign*__ieee754_y1(x));	236	if(n==1) return(sign*y1(x));
237	if(ix==0x7ff00000) return zero;	237	if(ix==0x7ff00000) return zero;
238	if(ix>=0x52D00000) { /* x > 2*302 /	238	if(ix>=0x52D00000) { /* x > 2*302 /
239	/* (x >> n**2)	239	/* (x >> n**2)
Lines 259-266 __ieee754_yn(int n, double x) Link Here
259	b = invsqrtpi*temp/sqrt(x);	259	b = invsqrtpi*temp/sqrt(x);
260	} else {	260	} else {
261	u_int32_t high;	261	u_int32_t high;
262	a = __ieee754_y0(x);	262	a = y0(x);
263	b = __ieee754_y1(x);	263	b = y1(x);
264	/* quit if b is -inf */	264	/* quit if b is -inf */
265	GET_HIGH_WORD(high,b);	265	GET_HIGH_WORD(high,b);
266	for(i=1;i<n&&high!=0xfff00000;i++){	266	for(i=1;i<n&&high!=0xfff00000;i++){

Lines 32-38 one = 1.0000000000e+00; /* 0x3F800000 */ Link Here

(-)b/lib/msun/src/e_jnf.c (-13 / +13 lines)
32	static const float zero = 0.0000000000e+00;	32	static const float zero = 0.0000000000e+00;
33		33
34	float	34	float
35	__ieee754_jnf(int n, float x)	35	jnf(int n, float x)
36	{	36	{
37	int32_t i,hx,ix, sgn;	37	int32_t i,hx,ix, sgn;
38	float a, b, temp, di;	38	float a, b, temp, di;
Lines 50-65 __ieee754_jnf(int n, float x) Link Here
50	x = -x;	50	x = -x;
51	hx ^= 0x80000000;	51	hx ^= 0x80000000;
52	}	52	}
53	if(n==0) return(__ieee754_j0f(x));	53	if(n==0) return(j0f(x));
54	if(n==1) return(__ieee754_j1f(x));	54	if(n==1) return(j1f(x));
55	sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */	55	sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
56	x = fabsf(x);	56	x = fabsf(x);
57	if(ix==0\|\|ix>=0x7f800000) /* if x is 0 or inf */	57	if(ix==0\|\|ix>=0x7f800000) /* if x is 0 or inf */
58	b = zero;	58	b = zero;
59	else if((float)n<=x) {	59	else if((float)n<=x) {
60	/* Safe to use J(n+1,x)=2n/x J(n,x)-J(n-1,x) /	60	/* Safe to use J(n+1,x)=2n/x J(n,x)-J(n-1,x) /
61	a = __ieee754_j0f(x);	61	a = j0f(x);
62	b = __ieee754_j1f(x);	62	b = j1f(x);
63	for(i=1;i<n;i++){	63	for(i=1;i<n;i++){
64	temp = b;	64	temp = b;
65	b = b((float)(i+i)/x) - a; / avoid underflow */	65	b = b((float)(i+i)/x) - a; / avoid underflow */
Lines 134-140 __ieee754_jnf(int n, float x) Link Here
134	*/	134	*/
135	tmp = n;	135	tmp = n;
136	v = two/x;	136	v = two/x;
137	tmp = tmp__ieee754_logf(fabsf(vtmp));	137	tmp = tmplogf(fabsf(vtmp));
138	if(tmp<(float)8.8721679688e+01) {	138	if(tmp<(float)8.8721679688e+01) {
139	for(i=n-1,di=(float)(i+i);i>0;i--){	139	for(i=n-1,di=(float)(i+i);i>0;i--){
140	temp = b;	140	temp = b;
Lines 158-165 __ieee754_jnf(int n, float x) Link Here
158	}	158	}
159	}	159	}
160	}	160	}
161	z = __ieee754_j0f(x);	161	z = j0f(x);
162	w = __ieee754_j1f(x);	162	w = j1f(x);
163	if (fabsf(z) >= fabsf(w))	163	if (fabsf(z) >= fabsf(w))
164	b = (t*z/b);	164	b = (t*z/b);
165	else	165	else
Lines 170-176 __ieee754_jnf(int n, float x) Link Here
170	}	170	}
171		171
172	float	172	float
173	__ieee754_ynf(int n, float x)	173	ynf(int n, float x)
174	{	174	{
175	int32_t i,hx,ix,ib;	175	int32_t i,hx,ix,ib;
176	int32_t sign;	176	int32_t sign;
Lines 186-197 __ieee754_ynf(int n, float x) Link Here
186	n = -n;	186	n = -n;
187	sign = 1 - ((n&1)<<1);	187	sign = 1 - ((n&1)<<1);
188	}	188	}
189	if(n==0) return(__ieee754_y0f(x));	189	if(n==0) return(y0f(x));
190	if(n==1) return(sign*__ieee754_y1f(x));	190	if(n==1) return(sign*y1f(x));
191	if(ix==0x7f800000) return zero;	191	if(ix==0x7f800000) return zero;
192		192
193	a = __ieee754_y0f(x);	193	a = y0f(x);
194	b = __ieee754_y1f(x);	194	b = y1f(x);
195	/* quit if b is -inf */	195	/* quit if b is -inf */
196	GET_FLOAT_WORD(ib,b);	196	GET_FLOAT_WORD(ib,b);
197	for(i=1;i<n&&ib!=0xff800000;i++){	197	for(i=1;i<n&&ib!=0xff800000;i++){




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* lgamma(x)
 * Return the logarithm of the Gamma function of x.
 *
 * Method: call lgamma_r
 */

#include <float.h>

extern int signgam;

double
lgamma(double x)
{
	return lgamma_r(x,&signgam);
}

#if (LDBL_MANT_DIG == 53)




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* lgamma_r(x, signgamp)
 * Reentrant version of the logarithm of the Gamma function
 * with user provide pointer for the sign of Gamma(x).
 *



double
lgamma_r(double x, int *signgamp)
{
	double nadj,p,p1,p2,p3,q,r,t,w,y,z;
	int32_t hx;

	if(ix<0x3c700000) {	/* |x|<2**-56, return -log(|x|) */
	    if((ix|lx)==0)
	        return one/vzero;
	    return -log(fabs(x));
	}

    /* purge negative integers and start evaluation for other x < 0 */

		return one/vzero;
	    t = sin_pi(x);
	    if(t==zero) return one/vzero; /* -integer */
	    nadj = log(pi/fabs(t*x));
	    if(t<zero) *signgamp = -1;
	    x = -x;
	}

    /* for x < 2.0 */
	else if(ix<0x40000000) {
	    if(ix<=0x3feccccc) { 	/* lgamma(x) = lgamma(x+1)-log(x) */
		r = -log(x);
		if(ix>=0x3FE76944) {y = one-x; i= 0;}
		else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
	  	else {y = x; i=2;}

	    case 5: z *= (y+4);		/* FALLTHRU */
	    case 4: z *= (y+3);		/* FALLTHRU */
	    case 3: z *= (y+2);		/* FALLTHRU */
		    r += log(z); break;
	    }
    /* 8.0 <= x < 2**56 */
	} else if (ix < 0x43700000) {
	    t = log(x);
	    z = one/x;
	    y = z*z;
	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
	    r = (x-half)*(t-one)+w;
	} else
    /* 2**56 <= x <= inf */
	    r =  x*(log(x)-one);
	if(hx<0) r = nadj - r;
	return r;
}




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* lgammaf(x)
 * Return the logarithm of the Gamma function of x.
 *
 * Method: call lgammaf_r
 */

#include "math.h"

extern int signgam;

float
lgammaf(float x)
{
	return lgammaf_r(x,&signgam);
}

Lines 120-126 sin_pif(float x) Link Here

(-)b/lib/msun/src/e_lgammaf_r.c (-7 / +7 lines)
120		120
121		121
122	float	122	float
123	__ieee754_lgammaf_r(float x, int *signgamp)	123	lgammaf_r(float x, int *signgamp)
124	{	124	{
125	float nadj,p,p1,p2,q,r,t,w,y,z;	125	float nadj,p,p1,p2,q,r,t,w,y,z;
126	int32_t hx;	126	int32_t hx;
Lines 138-144 __ieee754_lgammaf_r(float x, int *signgamp) Link Here
138	if(ix<0x32000000) { /* \|x\|<2*-27, return -log(\|x\|) /	138	if(ix<0x32000000) { /* \|x\|<2*-27, return -log(\|x\|) /
139	if(ix==0)	139	if(ix==0)
140	return one/vzero;	140	return one/vzero;
141	return -__ieee754_logf(fabsf(x));	141	return -logf(fabsf(x));
142	}	142	}
143		143
144	/* purge negative integers and start evaluation for other x < 0 */	144	/* purge negative integers and start evaluation for other x < 0 */
Lines 148-154 __ieee754_lgammaf_r(float x, int *signgamp) Link Here
148	return one/vzero;	148	return one/vzero;
149	t = sin_pif(x);	149	t = sin_pif(x);
150	if(t==zero) return one/vzero; /* -integer */	150	if(t==zero) return one/vzero; /* -integer */
151	nadj = __ieee754_logf(pi/fabsf(t*x));	151	nadj = logf(pi/fabsf(t*x));
152	if(t<zero) *signgamp = -1;	152	if(t<zero) *signgamp = -1;
153	x = -x;	153	x = -x;
154	}	154	}
Lines 158-164 __ieee754_lgammaf_r(float x, int *signgamp) Link Here
158	/* for x < 2.0 */	158	/* for x < 2.0 */
159	else if(ix<0x40000000) {	159	else if(ix<0x40000000) {
160	if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */	160	if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
161	r = -__ieee754_logf(x);	161	r = -logf(x);
162	if(ix>=0x3f3b4a20) {y = one-x; i= 0;}	162	if(ix>=0x3f3b4a20) {y = one-x; i= 0;}
163	else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;}	163	else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;}
164	else {y = x; i=2;}	164	else {y = x; i=2;}
Lines 198-215 __ieee754_lgammaf_r(float x, int *signgamp) Link Here
198	case 5: z = (y+4); / FALLTHRU */	198	case 5: z = (y+4); / FALLTHRU */
199	case 4: z = (y+3); / FALLTHRU */	199	case 4: z = (y+3); / FALLTHRU */
200	case 3: z = (y+2); / FALLTHRU */	200	case 3: z = (y+2); / FALLTHRU */
201	r += __ieee754_logf(z); break;	201	r += logf(z); break;
202	}	202	}
203	/* 8.0 <= x < 2*27 /	203	/* 8.0 <= x < 2*27 /
204	} else if (ix < 0x4d000000) {	204	} else if (ix < 0x4d000000) {
205	t = __ieee754_logf(x);	205	t = logf(x);
206	z = one/x;	206	z = one/x;
207	y = z*z;	207	y = z*z;
208	w = w0+z(w1+yw2);	208	w = w0+z(w1+yw2);
209	r = (x-half)*(t-one)+w;	209	r = (x-half)*(t-one)+w;
210	} else	210	} else
211	/* 2*27 <= x <= inf /	211	/* 2*27 <= x <= inf /
212	r = x*(__ieee754_logf(x)-one);	212	r = x*(logf(x)-one);
213	if(hx<0) r = nadj - r;	213	if(hx<0) r = nadj - r;
214	return r;	214	return r;
215	}	215	}

Lines 14-20 Link Here

(-)b/lib/msun/src/e_log.c (-2 / +2 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_log(x)	17	/* log(x)
18	* Return the logrithm of x	18	* Return the logrithm of x
19	*	19	*
20	* Method :	20	* Method :
Lines 86-92 static const double zero = 0.0; Link Here
86	static volatile double vzero = 0.0;	86	static volatile double vzero = 0.0;
87		87
88	double	88	double
89	__ieee754_log(double x)	89	log(double x)
90	{	90	{
91	double hfsq,f,s,z,R,w,t1,t2,dk;	91	double hfsq,f,s,z,R,w,t1,t2,dk;
92	int32_t k,hx,i,j;	92	int32_t k,hx,i,j;

Lines 39-45 static const double zero = 0.0; Link Here

(-)b/lib/msun/src/e_log10.c (-1 / +1 lines)
39	static volatile double vzero = 0.0;	39	static volatile double vzero = 0.0;
40		40
41	double	41	double
42	__ieee754_log10(double x)	42	log10(double x)
43	{	43	{
44	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;	44	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
45	int32_t i,k,hx;	45	int32_t i,k,hx;

Lines 31-37 static const float zero = 0.0; Link Here

(-)b/lib/msun/src/e_log10f.c (-1 / +1 lines)
31	static volatile float vzero = 0.0;	31	static volatile float vzero = 0.0;
32		32
33	float	33	float
34	__ieee754_log10f(float x)	34	log10f(float x)
35	{	35	{
36	float f,hfsq,hi,lo,r,y;	36	float f,hfsq,hi,lo,r,y;
37	int32_t i,k,hx;	37	int32_t i,k,hx;

Lines 39-45 static const double zero = 0.0; Link Here

(-)b/lib/msun/src/e_log2.c (-1 / +1 lines)
39	static volatile double vzero = 0.0;	39	static volatile double vzero = 0.0;
40		40
41	double	41	double
42	__ieee754_log2(double x)	42	log2(double x)
43	{	43	{
44	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;	44	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
45	int32_t i,k,hx;	45	int32_t i,k,hx;

Lines 29-35 static const float zero = 0.0; Link Here

(-)b/lib/msun/src/e_log2f.c (-1 / +1 lines)
29	static volatile float vzero = 0.0;	29	static volatile float vzero = 0.0;
30		30
31	float	31	float
32	__ieee754_log2f(float x)	32	log2f(float x)
33	{	33	{
34	float f,hfsq,hi,lo,r,y;	34	float f,hfsq,hi,lo,r,y;
35	int32_t i,k,hx;	35	int32_t i,k,hx;

Lines 33-39 static const float zero = 0.0; Link Here

(-)b/lib/msun/src/e_logf.c (-1 / +1 lines)
33	static volatile float vzero = 0.0;	33	static volatile float vzero = 0.0;
34		34
35	float	35	float
36	__ieee754_logf(float x)	36	logf(float x)
37	{	37	{
38	float hfsq,f,s,z,R,w,t1,t2,dk;	38	float hfsq,f,s,z,R,w,t1,t2,dk;
39	int32_t k,ix,i,j;	39	int32_t k,ix,i,j;

Lines 12-18 Link Here

(-)b/lib/msun/src/e_pow.c (-2 / +2 lines)
12	#include <sys/cdefs.h>	12	#include <sys/cdefs.h>
13	__FBSDID("$FreeBSD$");	13	__FBSDID("$FreeBSD$");
14		14
15	/* __ieee754_pow(x,y) return x**y	15	/* pow(x,y) return x**y
16	*	16	*
17	* n	17	* n
18	* Method: Let x = 2 * (1+f)	18	* Method: Let x = 2 * (1+f)
Lines 98-104 ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ Link Here
98	ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/	98	ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
99		99
100	double	100	double
101	__ieee754_pow(double x, double y)	101	pow(double x, double y)
102	{	102	{
103	double z,ax,z_h,z_l,p_h,p_l;	103	double z,ax,z_h,z_l,p_h,p_l;
104	double y1,t1,t2,r,s,t,u,v,w;	104	double y1,t1,t2,r,s,t,u,v,w;

Lines 56-62 ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ Link Here

(-)b/lib/msun/src/e_powf.c (-2 / +2 lines)
56	ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/	56	ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
57		57
58	float	58	float
59	__ieee754_powf(float x, float y)	59	powf(float x, float y)
60	{	60	{
61	float z,ax,z_h,z_l,p_h,p_l;	61	float z,ax,z_h,z_l,p_h,p_l;
62	float y1,t1,t2,r,s,sn,t,u,v,w;	62	float y1,t1,t2,r,s,sn,t,u,v,w;
Lines 108-114 __ieee754_powf(float x, float y) Link Here
108	if(hy==0x40000000) return xx; / y is 2 */	108	if(hy==0x40000000) return xx; / y is 2 */
109	if(hy==0x3f000000) { /* y is 0.5 */	109	if(hy==0x3f000000) { /* y is 0.5 */
110	if(hx>=0) /* x >= +0 */	110	if(hx>=0) /* x >= +0 */
111	return __ieee754_sqrtf(x);	111	return sqrtf(x);
112	}	112	}
113		113
114	ax = fabsf(x);	114	ax = fabsf(x);




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* remainder(x,p)
 * Return :                  
 * 	returns  x REM p  =  x - [x/p]*p as if in infinite 
 * 	precise arithmetic, where [x/p] is the (infinite bit) 



double
remainder(double x, double p)
{
	int32_t hx,hp;
	u_int32_t sx,lx,lp;

	    return nan_mix_op(x, p, *)/nan_mix_op(x, p, *);


	if (hp<=0x7fdfffff) x = fmod(x,p+p);	/* now x < 2p */
	if (((hx-hp)|(lx-lp))==0) return zero*x;
	x  = fabs(x);
	p  = fabs(p);

Lines 23-29 static const float zero = 0.0; Link Here

(-)b/lib/msun/src/e_remainderf.c (-2 / +2 lines)
23		23
24		24
25	float	25	float
26	__ieee754_remainderf(float x, float p)	26	remainderf(float x, float p)
27	{	27	{
28	int32_t hx,hp;	28	int32_t hx,hp;
29	u_int32_t sx;	29	u_int32_t sx;
Lines 42-48 __ieee754_remainderf(float x, float p) Link Here
42	return nan_mix_op(x, p, )/nan_mix_op(x, p, );	42	return nan_mix_op(x, p, )/nan_mix_op(x, p, );
43		43
44		44
45	if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */	45	if (hp<=0x7effffff) x = fmodf(x,p+p); /* now x < 2p */
46	if ((hx-hp)==0) return zero*x;	46	if ((hx-hp)==0) return zero*x;
47	x = fabsf(x);	47	x = fabsf(x);
48	p = fabsf(p);	48	p = fabsf(p);




__FBSDID("$FreeBSD$");

/*
 * scalb(x, fn) is provide for
 * passing various standard test suite. One 
 * should use scalbn() instead.
 */


#ifdef _SCALB_INT
double
scalb(double x, int fn)
#else
double
scalb(double x, double fn)
#endif
{
#ifdef _SCALB_INT

Lines 21-30 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/e_scalbf.c (-2 / +2 lines)
21		21
22	#ifdef _SCALB_INT	22	#ifdef _SCALB_INT
23	float	23	float
24	__ieee754_scalbf(float x, int fn)	24	scalbf(float x, int fn)
25	#else	25	#else
26	float	26	float
27	__ieee754_scalbf(float x, float fn)	27	scalbf(float x, float fn)
28	#endif	28	#endif
29	{	29	{
30	#ifdef _SCALB_INT	30	#ifdef _SCALB_INT




#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

/* sinh(x)
 * Method : 
 * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
 *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 

static const double one = 1.0, shuge = 1.0e307;

double
sinh(double x)
{
	double t,h;
	int32_t ix,jx;

	}

    /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
	if (ix < 0x40862E42)  return h*exp(fabs(x));

    /* |x| in [log(maxdouble), overflowthresold] */
	if (ix<=0x408633CE)

Lines 22-28 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/e_sinhf.c (-2 / +2 lines)
22	static const float one = 1.0, shuge = 1.0e37;	22	static const float one = 1.0, shuge = 1.0e37;
23		23
24	float	24	float
25	__ieee754_sinhf(float x)	25	sinhf(float x)
26	{	26	{
27	float t,h;	27	float t,h;
28	int32_t ix,jx;	28	int32_t ix,jx;
Lines 45-51 __ieee754_sinhf(float x) Link Here
45	}	45	}
46		46
47	/* \|x\| in [9, logf(maxfloat)] return 0.5exp(\|x\|) /	47	/* \|x\| in [9, logf(maxfloat)] return 0.5exp(\|x\|) /
48	if (ix < 0x42b17217) return h*__ieee754_expf(fabsf(x));	48	if (ix < 0x42b17217) return h*expf(fabsf(x));
49		49
50	/* \|x\| in [logf(maxfloat), overflowthresold] */	50	/* \|x\| in [logf(maxfloat), overflowthresold] */
51	if (ix<=0x42b2d4fc)	51	if (ix<=0x42b2d4fc)





#ifdef USE_BUILTIN_SQRT
double
sqrt(double x)
{
	return (__builtin_sqrt(x));
}
#else
/* sqrt(x)
 * Return correctly rounded sqrt.
 *           ------------------------------------------
 *	     |  Use the hardware sqrt if you have one |

static	const double	one	= 1.0, tiny=1.0e-300;

double
sqrt(double x)
{
	double z;
	int32_t sign = (int)0x80000000;

Lines 22-28 static char rcsid[] = "$FreeBSD$"; Link Here

(-)b/lib/msun/src/e_sqrtf.c (-2 / +2 lines)
22		22
23	#ifdef USE_BUILTIN_SQRTF	23	#ifdef USE_BUILTIN_SQRTF
24	float	24	float
25	__ieee754_sqrtf(float x)	25	sqrtf(float x)
26	{	26	{
27	return (__builtin_sqrtf(x));	27	return (__builtin_sqrtf(x));
28	}	28	}
Lines 30-36 __ieee754_sqrtf(float x) Link Here
30	static const float one = 1.0, tiny=1.0e-30;	30	static const float one = 1.0, tiny=1.0e-30;
31		31
32	float	32	float
33	__ieee754_sqrtf(float x)	33	sqrtf(float x)
34	{	34	{
35	float z;	35	float z;
36	int32_t sign = (int)0x80000000;	36	int32_t sign = (int)0x80000000;




	__x + __y;			\
})

/*
 * ieee style elementary functions
 *
 * We rename functions here to improve other sources' diffability
 * against fdlibm.
 */
#define	__ieee754_sqrt	sqrt
#define	__ieee754_acos	acos
#define	__ieee754_acosh	acosh
#define	__ieee754_log	log
#define	__ieee754_log2	log2
#define	__ieee754_atanh	atanh
#define	__ieee754_asin	asin
#define	__ieee754_atan2	atan2
#define	__ieee754_exp	exp
#define	__ieee754_cosh	cosh
#define	__ieee754_fmod	fmod
#define	__ieee754_pow	pow
#define	__ieee754_lgamma lgamma
#define	__ieee754_gamma	gamma
#define	__ieee754_lgamma_r lgamma_r
#define	__ieee754_gamma_r gamma_r
#define	__ieee754_log10	log10
#define	__ieee754_sinh	sinh
#define	__ieee754_hypot	hypot
#define	__ieee754_j0	j0
#define	__ieee754_j1	j1
#define	__ieee754_y0	y0
#define	__ieee754_y1	y1
#define	__ieee754_jn	jn
#define	__ieee754_yn	yn
#define	__ieee754_remainder remainder
#define	__ieee754_scalb	scalb
#define	__ieee754_sqrtf	sqrtf
#define	__ieee754_acosf	acosf
#define	__ieee754_acoshf acoshf
#define	__ieee754_logf	logf
#define	__ieee754_atanhf atanhf
#define	__ieee754_asinf	asinf
#define	__ieee754_atan2f atan2f
#define	__ieee754_expf	expf
#define	__ieee754_coshf	coshf
#define	__ieee754_fmodf	fmodf
#define	__ieee754_powf	powf
#define	__ieee754_lgammaf lgammaf
#define	__ieee754_gammaf gammaf
#define	__ieee754_lgammaf_r lgammaf_r
#define	__ieee754_gammaf_r gammaf_r
#define	__ieee754_log10f log10f
#define	__ieee754_log2f log2f
#define	__ieee754_sinhf	sinhf
#define	__ieee754_hypotf hypotf
#define	__ieee754_j0f	j0f
#define	__ieee754_j1f	j1f
#define	__ieee754_y0f	y0f
#define	__ieee754_y1f	y1f
#define	__ieee754_jnf	jnf
#define	__ieee754_ynf	ynf
#define	__ieee754_remainderf remainderf
#define	__ieee754_scalbf scalbf

/* fdlibm kernel function */
int	__kernel_rem_pio2(double*,double*,int,int,int);


Lines 46-58 asinh(double x) Link Here

(-)b/lib/msun/src/s_asinh.c (-3 / +3 lines)
46	if(huge+x>one) return x; /* return x inexact except 0 */	46	if(huge+x>one) return x; /* return x inexact except 0 */
47	}	47	}
48	if(ix>0x41b00000) { /* \|x\| > 2*28 /	48	if(ix>0x41b00000) { /* \|x\| > 2*28 /
49	w = __ieee754_log(fabs(x))+ln2;	49	w = log(fabs(x))+ln2;
50	} else if (ix>0x40000000) { /* 2*28 > \|x\| > 2.0 /	50	} else if (ix>0x40000000) { /* 2*28 > \|x\| > 2.0 /
51	t = fabs(x);	51	t = fabs(x);
52	w = __ieee754_log(2.0t+one/(__ieee754_sqrt(xx+one)+t));	52	w = log(2.0t+one/(sqrt(xx+one)+t));
53	} else { /* 2.0 > \|x\| > 2*-28 /	53	} else { /* 2.0 > \|x\| > 2*-28 /
54	t = x*x;	54	t = x*x;
55	w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));	55	w =log1p(fabs(x)+t/(one+sqrt(one+t)));
56	}	56	}
57	if(hx>0) return w; else return -w;	57	if(hx>0) return w; else return -w;
58	}	58	}

Lines 36-48 asinhf(float x) Link Here

(-)b/lib/msun/src/s_asinhf.c (-3 / +3 lines)
36	if(huge+x>one) return x; /* return x inexact except 0 */	36	if(huge+x>one) return x; /* return x inexact except 0 */
37	}	37	}
38	if(ix>0x4d800000) { /* \|x\| > 2*28 /	38	if(ix>0x4d800000) { /* \|x\| > 2*28 /
39	w = __ieee754_logf(fabsf(x))+ln2;	39	w = logf(fabsf(x))+ln2;
40	} else if (ix>0x40000000) { /* 2*28 > \|x\| > 2.0 /	40	} else if (ix>0x40000000) { /* 2*28 > \|x\| > 2.0 /
41	t = fabsf(x);	41	t = fabsf(x);
42	w = __ieee754_logf((float)2.0t+one/(__ieee754_sqrtf(xx+one)+t));	42	w = logf((float)2.0t+one/(sqrtf(xx+one)+t));
43	} else { /* 2.0 > \|x\| > 2*-28 /	43	} else { /* 2.0 > \|x\| > 2*-28 /
44	t = x*x;	44	t = x*x;
45	w =log1pf(fabsf(x)+t/(one+__ieee754_sqrtf(one+t)));	45	w =log1pf(fabsf(x)+t/(one+sqrtf(one+t)));
46	}	46	}
47	if(hx>0) return w; else return -w;	47	if(hx>0) return w; else return -w;
48	}	48	}

Lines 238-244 erf(double x) Link Here

(-)b/lib/msun/src/s_erf.c (-2 / +2 lines)
238	}	238	}
239	z = x;	239	z = x;
240	SET_LOW_WORD(z,0);	240	SET_LOW_WORD(z,0);
241	r = __ieee754_exp(-zz-0.5625)__ieee754_exp((z-x)*(z+x)+R/S);	241	r = exp(-zz-0.5625)exp((z-x)*(z+x)+R/S);
242	if(hx>=0) return one-r/x; else return r/x-one;	242	if(hx>=0) return one-r/x; else return r/x-one;
243	}	243	}
244		244
Lines 297-303 erfc(double x) Link Here
297	}	297	}
298	z = x;	298	z = x;
299	SET_LOW_WORD(z,0);	299	SET_LOW_WORD(z,0);
300	r = __ieee754_exp(-zz-0.5625)__ieee754_exp((z-x)*(z+x)+R/S);	300	r = exp(-zz-0.5625)exp((z-x)*(z+x)+R/S);
301	if(hx>0) return r/x; else return two-r/x;	301	if(hx>0) return r/x; else return two-r/x;
302	} else {	302	} else {
303	if(hx>0) return tiny*tiny; else return two-tiny;	303	if(hx>0) return tiny*tiny; else return two-tiny;

Lines 25-29 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/s_significand.c (-1 / +1 lines)
25	double	25	double
26	significand(double x)	26	significand(double x)
27	{	27	{
28	return __ieee754_scalb(x,(double) -ilogb(x));	28	return scalb(x,(double) -ilogb(x));
29	}	29	}

Lines 22-26 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/s_significandf.c (-1 / +1 lines)
22	float	22	float
23	significandf(float x)	23	significandf(float x)
24	{	24	{
25	return __ieee754_scalbf(x,(float) -ilogbf(x));	25	return scalbf(x,(float) -ilogbf(x));
26	}	26	}

Return to bug 272783

Lines 15-21 Link Here

(-)b/lib/msun/src/e_acosh.c (-4 / +4 lines)
15	#include <sys/cdefs.h>	15	#include <sys/cdefs.h>
16	__FBSDID("$FreeBSD$");	16	__FBSDID("$FreeBSD$");
17		17
18	/* __ieee754_acosh(x)	18	/* acosh(x)
19	* Method :	19	* Method :
20	* Based on	20	* Based on
21	* acosh(x) = log [ x + sqrt(x*x-1) ]	21	* acosh(x) = log [ x + sqrt(x*x-1) ]
Lines 39-45 one = 1.0, Link Here
39	ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */	39	ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
40		40
41	double	41	double
42	__ieee754_acosh(double x)	42	acosh(double x)
43	{	43	{
44	double t;	44	double t;
45	int32_t hx;	45	int32_t hx;
Lines 51-62 __ieee754_acosh(double x) Link Here
51	if(hx >=0x7ff00000) { /* x is inf of NaN */	51	if(hx >=0x7ff00000) { /* x is inf of NaN */
52	return x+x;	52	return x+x;
53	} else	53	} else
54	return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */	54	return log(x)+ln2; /* acosh(huge)=log(2x) */
55	} else if(((hx-0x3ff00000)\|lx)==0) {	55	} else if(((hx-0x3ff00000)\|lx)==0) {
56	return 0.0; /* acosh(1) = 0 */	56	return 0.0; /* acosh(1) = 0 */
57	} else if (hx > 0x40000000) { /* 2*28 > x > 2 /	57	} else if (hx > 0x40000000) { /* 2*28 > x > 2 /
58	t=x*x;	58	t=x*x;
59	return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));	59	return log(2.0*x-one/(x+sqrt(t-one)));
60	} else { /* 1<x<2 */	60	} else { /* 1<x<2 */
61	t = x-one;	61	t = x-one;
62	return log1p(t+sqrt(2.0t+tt));	62	return log1p(t+sqrt(2.0t+tt));

Lines 24-30 one = 1.0, Link Here

(-)b/lib/msun/src/e_acoshf.c (-4 / +4 lines)
24	ln2 = 6.9314718246e-01; /* 0x3f317218 */	24	ln2 = 6.9314718246e-01; /* 0x3f317218 */
25		25
26	float	26	float
27	__ieee754_acoshf(float x)	27	acoshf(float x)
28	{	28	{
29	float t;	29	float t;
30	int32_t hx;	30	int32_t hx;
Lines 35-48 __ieee754_acoshf(float x) Link Here
35	if(hx >=0x7f800000) { /* x is inf of NaN */	35	if(hx >=0x7f800000) { /* x is inf of NaN */
36	return x+x;	36	return x+x;
37	} else	37	} else
38	return __ieee754_logf(x)+ln2; /* acosh(huge)=log(2x) */	38	return logf(x)+ln2; /* acosh(huge)=log(2x) */
39	} else if (hx==0x3f800000) {	39	} else if (hx==0x3f800000) {
40	return 0.0; /* acosh(1) = 0 */	40	return 0.0; /* acosh(1) = 0 */
41	} else if (hx > 0x40000000) { /* 2*28 > x > 2 /	41	} else if (hx > 0x40000000) { /* 2*28 > x > 2 /
42	t=x*x;	42	t=x*x;
43	return __ieee754_logf((float)2.0*x-one/(x+__ieee754_sqrtf(t-one)));	43	return logf((float)2.0*x-one/(x+sqrtf(t-one)));
44	} else { /* 1<x<2 */	44	} else { /* 1<x<2 */
45	t = x-one;	45	t = x-one;
46	return log1pf(t+__ieee754_sqrtf((float)2.0t+tt));	46	return log1pf(t+sqrtf((float)2.0t+tt));
47	}	47	}
48	}	48	}

Lines 14-20 Link Here

(-)b/lib/msun/src/e_cosh.c (-4 / +4 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_cosh(x)	17	/* cosh(x)
18	* Method :	18	* Method :
19	* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2	19	* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
20	* 1. Replace x by \|x\| (cosh(x) = cosh(-x)).	20	* 1. Replace x by \|x\| (cosh(x) = cosh(-x)).
Lines 43-49 __FBSDID("$FreeBSD$"); Link Here
43	static const double one = 1.0, half=0.5, huge = 1.0e300;	43	static const double one = 1.0, half=0.5, huge = 1.0e300;
44		44
45	double	45	double
46	__ieee754_cosh(double x)	46	cosh(double x)
47	{	47	{
48	double t,w;	48	double t,w;
49	int32_t ix;	49	int32_t ix;
Lines 65-76 __ieee754_cosh(double x) Link Here
65		65
66	/* \|x\| in [0.5ln2,22], return (exp(\|x\|)+1/exp(\|x\|)/2; /	66	/* \|x\| in [0.5ln2,22], return (exp(\|x\|)+1/exp(\|x\|)/2; /
67	if (ix < 0x40360000) {	67	if (ix < 0x40360000) {
68	t = __ieee754_exp(fabs(x));	68	t = exp(fabs(x));
69	return half*t+half/t;	69	return half*t+half/t;
70	}	70	}
71		71
72	/* \|x\| in [22, log(maxdouble)] return halfexp(\|x\|) /	72	/* \|x\| in [22, log(maxdouble)] return halfexp(\|x\|) /
73	if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x));	73	if (ix < 0x40862E42) return half*exp(fabs(x));
74		74
75	/* \|x\| in [log(maxdouble), overflowthresold] */	75	/* \|x\| in [log(maxdouble), overflowthresold] */
76	if (ix<=0x408633CE)	76	if (ix<=0x408633CE)

Lines 22-28 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/e_coshf.c (-3 / +3 lines)
22	static const float one = 1.0, half=0.5, huge = 1.0e30;	22	static const float one = 1.0, half=0.5, huge = 1.0e30;
23		23
24	float	24	float
25	__ieee754_coshf(float x)	25	coshf(float x)
26	{	26	{
27	float t,w;	27	float t,w;
28	int32_t ix;	28	int32_t ix;
Lines 43-54 __ieee754_coshf(float x) Link Here
43		43
44	/* \|x\| in [0.5ln2,9], return (exp(\|x\|)+1/exp(\|x\|))/2; /	44	/* \|x\| in [0.5ln2,9], return (exp(\|x\|)+1/exp(\|x\|))/2; /
45	if (ix < 0x41100000) {	45	if (ix < 0x41100000) {
46	t = __ieee754_expf(fabsf(x));	46	t = expf(fabsf(x));
47	return half*t+half/t;	47	return half*t+half/t;
48	}	48	}
49		49
50	/* \|x\| in [9, log(maxfloat)] return halfexp(\|x\|) /	50	/* \|x\| in [9, log(maxfloat)] return halfexp(\|x\|) /
51	if (ix < 0x42b17217) return half*__ieee754_expf(fabsf(x));	51	if (ix < 0x42b17217) return half*expf(fabsf(x));
52		52
53	/* \|x\| in [log(maxfloat), overflowthresold] */	53	/* \|x\| in [log(maxfloat), overflowthresold] */
54	if (ix<=0x42b2d4fc)	54	if (ix<=0x42b2d4fc)

Lines 15-24 Link Here

(-)b/lib/msun/src/e_gamma.c (-4 / +4 lines)
15	#include <sys/cdefs.h>	15	#include <sys/cdefs.h>
16	__FBSDID("$FreeBSD$");	16	__FBSDID("$FreeBSD$");
17		17
18	/* __ieee754_gamma(x)	18	/* gamma(x)
19	* Return the logarithm of the Gamma function of x.	19	* Return the logarithm of the Gamma function of x.
20	*	20	*
21	* Method: call __ieee754_gamma_r	21	* Method: call gamma_r
22	*/	22	*/
23		23
24	#include "math.h"	24	#include "math.h"
Lines 27-33 __FBSDID("$FreeBSD$"); Link Here
27	extern int signgam;	27	extern int signgam;
28		28
29	double	29	double
30	__ieee754_gamma(double x)	30	gamma(double x)
31	{	31	{
32	return __ieee754_gamma_r(x,&signgam);	32	return gamma_r(x,&signgam);
33	}	33	}

Lines 15-32 Link Here

(-)b/lib/msun/src/e_gamma_r.c (-4 / +4 lines)
15	#include <sys/cdefs.h>	15	#include <sys/cdefs.h>
16	__FBSDID("$FreeBSD$");	16	__FBSDID("$FreeBSD$");
17		17
18	/* __ieee754_gamma_r(x, signgamp)	18	/* gamma_r(x, signgamp)
19	* Reentrant version of the logarithm of the Gamma function	19	* Reentrant version of the logarithm of the Gamma function
20	* with user provide pointer for the sign of Gamma(x).	20	* with user provide pointer for the sign of Gamma(x).
21	*	21	*
22	* Method: See __ieee754_lgamma_r	22	* Method: See lgamma_r
23	*/	23	*/
24		24
25	#include "math.h"	25	#include "math.h"
26	#include "math_private.h"	26	#include "math_private.h"
27		27
28	double	28	double
29	__ieee754_gamma_r(double x, int *signgamp)	29	gamma_r(double x, int *signgamp)
30	{	30	{
31	return __ieee754_lgamma_r(x,signgamp);	31	return lgamma_r(x,signgamp);
32	}	32	}

Lines 16-25 Link Here

(-)b/lib/msun/src/e_gammaf.c (-4 / +4 lines)
16	#include <sys/cdefs.h>	16	#include <sys/cdefs.h>
17	__FBSDID("$FreeBSD$");	17	__FBSDID("$FreeBSD$");
18		18
19	/* __ieee754_gammaf(x)	19	/* gammaf(x)
20	* Return the logarithm of the Gamma function of x.	20	* Return the logarithm of the Gamma function of x.
21	*	21	*
22	* Method: call __ieee754_gammaf_r	22	* Method: call gammaf_r
23	*/	23	*/
24		24
25	#include "math.h"	25	#include "math.h"
Lines 28-34 __FBSDID("$FreeBSD$"); Link Here
28	extern int signgam;	28	extern int signgam;
29		29
30	float	30	float
31	__ieee754_gammaf(float x)	31	gammaf(float x)
32	{	32	{
33	return __ieee754_gammaf_r(x,&signgam);	33	return gammaf_r(x,&signgam);
34	}	34	}

Lines 16-33 Link Here

(-)b/lib/msun/src/e_gammaf_r.c (-4 / +4 lines)
16	#include <sys/cdefs.h>	16	#include <sys/cdefs.h>
17	__FBSDID("$FreeBSD$");	17	__FBSDID("$FreeBSD$");
18		18
19	/* __ieee754_gammaf_r(x, signgamp)	19	/* gammaf_r(x, signgamp)
20	* Reentrant version of the logarithm of the Gamma function	20	* Reentrant version of the logarithm of the Gamma function
21	* with user provide pointer for the sign of Gamma(x).	21	* with user provide pointer for the sign of Gamma(x).
22	*	22	*
23	* Method: See __ieee754_lgammaf_r	23	* Method: See lgammaf_r
24	*/	24	*/
25		25
26	#include "math.h"	26	#include "math.h"
27	#include "math_private.h"	27	#include "math_private.h"
28		28
29	float	29	float
30	__ieee754_gammaf_r(float x, int *signgamp)	30	gammaf_r(float x, int *signgamp)
31	{	31	{
32	return __ieee754_lgammaf_r(x,signgamp);	32	return lgammaf_r(x,signgamp);
33	}	33	}

Lines 15-24 Link Here

(-)b/lib/msun/src/e_lgamma.c (-4 / +4 lines)
15	#include <sys/cdefs.h>	15	#include <sys/cdefs.h>
16	__FBSDID("$FreeBSD$");	16	__FBSDID("$FreeBSD$");
17		17
18	/* __ieee754_lgamma(x)	18	/* lgamma(x)
19	* Return the logarithm of the Gamma function of x.	19	* Return the logarithm of the Gamma function of x.
20	*	20	*
21	* Method: call __ieee754_lgamma_r	21	* Method: call lgamma_r
22	*/	22	*/
23		23
24	#include <float.h>	24	#include <float.h>
Lines 29-37 __FBSDID("$FreeBSD$"); Link Here
29	extern int signgam;	29	extern int signgam;
30		30
31	double	31	double
32	__ieee754_lgamma(double x)	32	lgamma(double x)
33	{	33	{
34	return __ieee754_lgamma_r(x,&signgam);	34	return lgamma_r(x,&signgam);
35	}	35	}
36		36
37	#if (LDBL_MANT_DIG == 53)	37	#if (LDBL_MANT_DIG == 53)

Lines 13-19 Link Here

(-)b/lib/msun/src/e_lgamma_r.c (-8 / +8 lines)
13	#include <sys/cdefs.h>	13	#include <sys/cdefs.h>
14	__FBSDID("$FreeBSD$");	14	__FBSDID("$FreeBSD$");
15		15
16	/* __ieee754_lgamma_r(x, signgamp)	16	/* lgamma_r(x, signgamp)
17	* Reentrant version of the logarithm of the Gamma function	17	* Reentrant version of the logarithm of the Gamma function
18	* with user provide pointer for the sign of Gamma(x).	18	* with user provide pointer for the sign of Gamma(x).
19	*	19	*
Lines 199-205 sin_pi(double x) Link Here
199		199
200		200
201	double	201	double
202	__ieee754_lgamma_r(double x, int *signgamp)	202	lgamma_r(double x, int *signgamp)
203	{	203	{
204	double nadj,p,p1,p2,p3,q,r,t,w,y,z;	204	double nadj,p,p1,p2,p3,q,r,t,w,y,z;
205	int32_t hx;	205	int32_t hx;
Lines 217-223 __ieee754_lgamma_r(double x, int *signgamp) Link Here
217	if(ix<0x3c700000) { /* \|x\|<2*-56, return -log(\|x\|) /	217	if(ix<0x3c700000) { /* \|x\|<2*-56, return -log(\|x\|) /
218	if((ix\|lx)==0)	218	if((ix\|lx)==0)
219	return one/vzero;	219	return one/vzero;
220	return -__ieee754_log(fabs(x));	220	return -log(fabs(x));
221	}	221	}
222		222
223	/* purge negative integers and start evaluation for other x < 0 */	223	/* purge negative integers and start evaluation for other x < 0 */
Lines 227-233 __ieee754_lgamma_r(double x, int *signgamp) Link Here
227	return one/vzero;	227	return one/vzero;
228	t = sin_pi(x);	228	t = sin_pi(x);
229	if(t==zero) return one/vzero; /* -integer */	229	if(t==zero) return one/vzero; /* -integer */
230	nadj = __ieee754_log(pi/fabs(t*x));	230	nadj = log(pi/fabs(t*x));
231	if(t<zero) *signgamp = -1;	231	if(t<zero) *signgamp = -1;
232	x = -x;	232	x = -x;
233	}	233	}
Lines 237-243 __ieee754_lgamma_r(double x, int *signgamp) Link Here
237	/* for x < 2.0 */	237	/* for x < 2.0 */
238	else if(ix<0x40000000) {	238	else if(ix<0x40000000) {
239	if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */	239	if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
240	r = -__ieee754_log(x);	240	r = -log(x);
241	if(ix>=0x3FE76944) {y = one-x; i= 0;}	241	if(ix>=0x3FE76944) {y = one-x; i= 0;}
242	else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}	242	else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
243	else {y = x; i=2;}	243	else {y = x; i=2;}
Lines 282-299 __ieee754_lgamma_r(double x, int *signgamp) Link Here
282	case 5: z = (y+4); / FALLTHRU */	282	case 5: z = (y+4); / FALLTHRU */
283	case 4: z = (y+3); / FALLTHRU */	283	case 4: z = (y+3); / FALLTHRU */
284	case 3: z = (y+2); / FALLTHRU */	284	case 3: z = (y+2); / FALLTHRU */
285	r += __ieee754_log(z); break;	285	r += log(z); break;
286	}	286	}
287	/* 8.0 <= x < 2*56 /	287	/* 8.0 <= x < 2*56 /
288	} else if (ix < 0x43700000) {	288	} else if (ix < 0x43700000) {
289	t = __ieee754_log(x);	289	t = log(x);
290	z = one/x;	290	z = one/x;
291	y = z*z;	291	y = z*z;
292	w = w0+z(w1+y(w2+y(w3+y(w4+y(w5+yw6)))));	292	w = w0+z(w1+y(w2+y(w3+y(w4+y(w5+yw6)))));
293	r = (x-half)*(t-one)+w;	293	r = (x-half)*(t-one)+w;
294	} else	294	} else
295	/* 2*56 <= x <= inf /	295	/* 2*56 <= x <= inf /
296	r = x*(__ieee754_log(x)-one);	296	r = x*(log(x)-one);
297	if(hx<0) r = nadj - r;	297	if(hx<0) r = nadj - r;
298	return r;	298	return r;
299	}	299	}

Lines 16-25 Link Here

(-)b/lib/msun/src/e_lgammaf.c (-4 / +4 lines)
16	#include <sys/cdefs.h>	16	#include <sys/cdefs.h>
17	__FBSDID("$FreeBSD$");	17	__FBSDID("$FreeBSD$");
18		18
19	/* __ieee754_lgammaf(x)	19	/* lgammaf(x)
20	* Return the logarithm of the Gamma function of x.	20	* Return the logarithm of the Gamma function of x.
21	*	21	*
22	* Method: call __ieee754_lgammaf_r	22	* Method: call lgammaf_r
23	*/	23	*/
24		24
25	#include "math.h"	25	#include "math.h"
Lines 28-34 __FBSDID("$FreeBSD$"); Link Here
28	extern int signgam;	28	extern int signgam;
29		29
30	float	30	float
31	__ieee754_lgammaf(float x)	31	lgammaf(float x)
32	{	32	{
33	return __ieee754_lgammaf_r(x,&signgam);	33	return lgammaf_r(x,&signgam);
34	}	34	}

Lines 14-20 Link Here

(-)b/lib/msun/src/e_remainder.c (-3 / +3 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_remainder(x,p)	17	/* remainder(x,p)
18	* Return :	18	* Return :
19	* returns x REM p = x - [x/p]*p as if in infinite	19	* returns x REM p = x - [x/p]*p as if in infinite
20	* precise arithmetic, where [x/p] is the (infinite bit)	20	* precise arithmetic, where [x/p] is the (infinite bit)
Lines 32-38 static const double zero = 0.0; Link Here
32		32
33		33
34	double	34	double
35	__ieee754_remainder(double x, double p)	35	remainder(double x, double p)
36	{	36	{
37	int32_t hx,hp;	37	int32_t hx,hp;
38	u_int32_t sx,lx,lp;	38	u_int32_t sx,lx,lp;
Lines 52-58 __ieee754_remainder(double x, double p) Link Here
52	return nan_mix_op(x, p, )/nan_mix_op(x, p, );	52	return nan_mix_op(x, p, )/nan_mix_op(x, p, );
53		53
54		54
55	if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */	55	if (hp<=0x7fdfffff) x = fmod(x,p+p); /* now x < 2p */
56	if (((hx-hp)\|(lx-lp))==0) return zero*x;	56	if (((hx-hp)\|(lx-lp))==0) return zero*x;
57	x = fabs(x);	57	x = fabs(x);
58	p = fabs(p);	58	p = fabs(p);

Lines 15-21 Link Here

(-)b/lib/msun/src/e_scalb.c (-3 / +3 lines)
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/*	17	/*
18	* __ieee754_scalb(x, fn) is provide for	18	* scalb(x, fn) is provide for
19	* passing various standard test suite. One	19	* passing various standard test suite. One
20	* should use scalbn() instead.	20	* should use scalbn() instead.
21	*/	21	*/
Lines 25-34 __FBSDID("$FreeBSD$"); Link Here
25		25
26	#ifdef _SCALB_INT	26	#ifdef _SCALB_INT
27	double	27	double
28	__ieee754_scalb(double x, int fn)	28	scalb(double x, int fn)
29	#else	29	#else
30	double	30	double
31	__ieee754_scalb(double x, double fn)	31	scalb(double x, double fn)
32	#endif	32	#endif
33	{	33	{
34	#ifdef _SCALB_INT	34	#ifdef _SCALB_INT

Lines 14-20 Link Here

(-)b/lib/msun/src/e_sinh.c (-3 / +3 lines)
14	#include <sys/cdefs.h>	14	#include <sys/cdefs.h>
15	__FBSDID("$FreeBSD$");	15	__FBSDID("$FreeBSD$");
16		16
17	/* __ieee754_sinh(x)	17	/* sinh(x)
18	* Method :	18	* Method :
19	* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2	19	* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
20	* 1. Replace x by \|x\| (sinh(-x) = -sinh(x)).	20	* 1. Replace x by \|x\| (sinh(-x) = -sinh(x)).
Lines 40-46 __FBSDID("$FreeBSD$"); Link Here
40	static const double one = 1.0, shuge = 1.0e307;	40	static const double one = 1.0, shuge = 1.0e307;
41		41
42	double	42	double
43	__ieee754_sinh(double x)	43	sinh(double x)
44	{	44	{
45	double t,h;	45	double t,h;
46	int32_t ix,jx;	46	int32_t ix,jx;
Lines 64-70 __ieee754_sinh(double x) Link Here
64	}	64	}
65		65
66	/* \|x\| in [22, log(maxdouble)] return 0.5exp(\|x\|) /	66	/* \|x\| in [22, log(maxdouble)] return 0.5exp(\|x\|) /
67	if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x));	67	if (ix < 0x40862E42) return h*exp(fabs(x));
68		68
69	/* \|x\| in [log(maxdouble), overflowthresold] */	69	/* \|x\| in [log(maxdouble), overflowthresold] */
70	if (ix<=0x408633CE)	70	if (ix<=0x408633CE)

Lines 21-32 __FBSDID("$FreeBSD$"); Link Here

(-)b/lib/msun/src/e_sqrt.c (-3 / +3 lines)
21		21
22	#ifdef USE_BUILTIN_SQRT	22	#ifdef USE_BUILTIN_SQRT
23	double	23	double
24	__ieee754_sqrt(double x)	24	sqrt(double x)
25	{	25	{
26	return (__builtin_sqrt(x));	26	return (__builtin_sqrt(x));
27	}	27	}
28	#else	28	#else
29	/* __ieee754_sqrt(x)	29	/* sqrt(x)
30	* Return correctly rounded sqrt.	30	* Return correctly rounded sqrt.
31	* ------------------------------------------	31	* ------------------------------------------
32	* \| Use the hardware sqrt if you have one \|	32	* \| Use the hardware sqrt if you have one \|
Lines 99-105 __ieee754_sqrt(double x) Link Here
99	static const double one = 1.0, tiny=1.0e-300;	99	static const double one = 1.0, tiny=1.0e-300;
100		100
101	double	101	double
102	__ieee754_sqrt(double x)	102	sqrt(double x)
103	{	103	{
104	double z;	104	double z;
105	int32_t sign = (int)0x80000000;	105	int32_t sign = (int)0x80000000;

Lines 822-888 irintl(long double x) Link Here

(-)b/lib/msun/src/math_private.h (-61 lines)
822	__x + __y; \	822	__x + __y; \
823	})	823	})
824		824
825	/*
826	* ieee style elementary functions
827	*
828	* We rename functions here to improve other sources' diffability
829	* against fdlibm.
830	*/
831	#define __ieee754_sqrt sqrt
832	#define __ieee754_acos acos
833	#define __ieee754_acosh acosh
834	#define __ieee754_log log
835	#define __ieee754_log2 log2
836	#define __ieee754_atanh atanh
837	#define __ieee754_asin asin
838	#define __ieee754_atan2 atan2
839	#define __ieee754_exp exp
840	#define __ieee754_cosh cosh
841	#define __ieee754_fmod fmod
842	#define __ieee754_pow pow
843	#define __ieee754_lgamma lgamma
844	#define __ieee754_gamma gamma
845	#define __ieee754_lgamma_r lgamma_r
846	#define __ieee754_gamma_r gamma_r
847	#define __ieee754_log10 log10
848	#define __ieee754_sinh sinh
849	#define __ieee754_hypot hypot
850	#define __ieee754_j0 j0
851	#define __ieee754_j1 j1
852	#define __ieee754_y0 y0
853	#define __ieee754_y1 y1
854	#define __ieee754_jn jn
855	#define __ieee754_yn yn
856	#define __ieee754_remainder remainder
857	#define __ieee754_scalb scalb
858	#define __ieee754_sqrtf sqrtf
859	#define __ieee754_acosf acosf
860	#define __ieee754_acoshf acoshf
861	#define __ieee754_logf logf
862	#define __ieee754_atanhf atanhf
863	#define __ieee754_asinf asinf
864	#define __ieee754_atan2f atan2f
865	#define __ieee754_expf expf
866	#define __ieee754_coshf coshf
867	#define __ieee754_fmodf fmodf
868	#define __ieee754_powf powf
869	#define __ieee754_lgammaf lgammaf
870	#define __ieee754_gammaf gammaf
871	#define __ieee754_lgammaf_r lgammaf_r
872	#define __ieee754_gammaf_r gammaf_r
873	#define __ieee754_log10f log10f
874	#define __ieee754_log2f log2f
875	#define __ieee754_sinhf sinhf
876	#define __ieee754_hypotf hypotf
877	#define __ieee754_j0f j0f
878	#define __ieee754_j1f j1f
879	#define __ieee754_y0f y0f
880	#define __ieee754_y1f y1f
881	#define __ieee754_jnf jnf
882	#define __ieee754_ynf ynf
883	#define __ieee754_remainderf remainderf
884	#define __ieee754_scalbf scalbf
885
886	/* fdlibm kernel function */	825	/* fdlibm kernel function */
887	int __kernel_rem_pio2(double,double,int,int,int);	826	int __kernel_rem_pio2(double,double,int,int,int);
888		827