WebSVN - shark - Blame - Rev 1619 - /shark/branches/xen/libc/arch/x86/libm/msun/src/s

Rev	Author	Line No.	Line
2	pj	1	/* s_erff.c -- float version of s_erf.c.
		2	* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
		3	*/
		4
		5	/*
		6	* ====================================================
		7	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
		8	*
		9	* Developed at SunPro, a Sun Microsystems, Inc. business.
		10	* Permission to use, copy, modify, and distribute this
		11	* software is freely granted, provided that this notice
		12	* is preserved.
		13	* ====================================================
		14	*/
		15
		16	#ifndef lint
		17	static char rcsid[] = "$\Id: s_erff.c,v 1.2 1995/05/30 05:49:32 rgrimes Exp $";
		18	#endif
		19
		20	#include "math.h"
		21	#include "math_private.h"
		22
		23	#ifdef __STDC__
		24	static const float
		25	#else
		26	static float
		27	#endif
		28	tiny = 1e-30,
		29	half= 5.0000000000e-01, /* 0x3F000000 */
		30	one = 1.0000000000e+00, /* 0x3F800000 */
		31	two = 2.0000000000e+00, /* 0x40000000 */
		32	/* c = (subfloat)0.84506291151 */
		33	erx = 8.4506291151e-01, /* 0x3f58560b */
		34	/*
		35	* Coefficients for approximation to erf on [0,0.84375]
		36	*/
		37	efx = 1.2837916613e-01, /* 0x3e0375d4 */
		38	efx8= 1.0270333290e+00, /* 0x3f8375d4 */
		39	pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
		40	pp1 = -3.2504209876e-01, /* 0xbea66beb */
		41	pp2 = -2.8481749818e-02, /* 0xbce9528f */
		42	pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
		43	pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
		44	qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
		45	qq2 = 6.5022252500e-02, /* 0x3d852a63 */
		46	qq3 = 5.0813062117e-03, /* 0x3ba68116 */
		47	qq4 = 1.3249473704e-04, /* 0x390aee49 */
		48	qq5 = -3.9602282413e-06, /* 0xb684e21a */
		49	/*
		50	* Coefficients for approximation to erf in [0.84375,1.25]
		51	*/
		52	pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
		53	pa1 = 4.1485610604e-01, /* 0x3ed46805 */
		54	pa2 = -3.7220788002e-01, /* 0xbebe9208 */
		55	pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
		56	pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
		57	pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
		58	pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
		59	qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
		60	qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
		61	qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
		62	qa4 = 1.2617121637e-01, /* 0x3e013307 */
		63	qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
		64	qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
		65	/*
		66	* Coefficients for approximation to erfc in [1.25,1/0.35]
		67	*/
		68	ra0 = -9.8649440333e-03, /* 0xbc21a093 */
		69	ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
		70	ra2 = -1.0558626175e+01, /* 0xc128f022 */
		71	ra3 = -6.2375331879e+01, /* 0xc2798057 */
		72	ra4 = -1.6239666748e+02, /* 0xc322658c */
		73	ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
		74	ra6 = -8.1287437439e+01, /* 0xc2a2932b */
		75	ra7 = -9.8143291473e+00, /* 0xc11d077e */
		76	sa1 = 1.9651271820e+01, /* 0x419d35ce */
		77	sa2 = 1.3765776062e+02, /* 0x4309a863 */
		78	sa3 = 4.3456588745e+02, /* 0x43d9486f */
		79	sa4 = 6.4538726807e+02, /* 0x442158c9 */
		80	sa5 = 4.2900814819e+02, /* 0x43d6810b */
		81	sa6 = 1.0863500214e+02, /* 0x42d9451f */
		82	sa7 = 6.5702495575e+00, /* 0x40d23f7c */
		83	sa8 = -6.0424413532e-02, /* 0xbd777f97 */
		84	/*
		85	* Coefficients for approximation to erfc in [1/.35,28]
		86	*/
		87	rb0 = -9.8649431020e-03, /* 0xbc21a092 */
		88	rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
		89	rb2 = -1.7757955551e+01, /* 0xc18e104b */
		90	rb3 = -1.6063638306e+02, /* 0xc320a2ea */
		91	rb4 = -6.3756646729e+02, /* 0xc41f6441 */
		92	rb5 = -1.0250950928e+03, /* 0xc480230b */
		93	rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
		94	sb1 = 3.0338060379e+01, /* 0x41f2b459 */
		95	sb2 = 3.2579251099e+02, /* 0x43a2e571 */
		96	sb3 = 1.5367296143e+03, /* 0x44c01759 */
		97	sb4 = 3.1998581543e+03, /* 0x4547fdbb */
		98	sb5 = 2.5530502930e+03, /* 0x451f90ce */
		99	sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
		100	sb7 = -2.2440952301e+01; /* 0xc1b38712 */
		101
		102	#ifdef __STDC__
		103	float erff(float x)
		104	#else
		105	float erff(x)
		106	float x;
		107	#endif
		108	{
		109	int32_t hx,ix,i;
		110	float R,S,P,Q,s,y,z,r;
		111	GET_FLOAT_WORD(hx,x);
		112	ix = hx&0x7fffffff;
		113	if(ix>=0x7f800000) { /* erf(nan)=nan */
		114	i = ((u_int32_t)hx>>31)<<1;
		115	return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */
		116	}
		117
		118	if(ix < 0x3f580000) { /* \|x\|<0.84375 */
		119	if(ix < 0x31800000) { /* \|x\|<2*-28 /
		120	if (ix < 0x04000000)
		121	/avoid underflow /
		122	return (float)0.125((float)8.0x+efx8*x);
		123	return x + efx*x;
		124	}
		125	z = x*x;
		126	r = pp0+z(pp1+z(pp2+z(pp3+zpp4)));
		127	s = one+z(qq1+z(qq2+z(qq3+z(qq4+z*qq5))));
		128	y = r/s;
		129	return x + x*y;
		130	}
		131	if(ix < 0x3fa00000) { /* 0.84375 <= \|x\| < 1.25 */
		132	s = fabsf(x)-one;
		133	P = pa0+s(pa1+s(pa2+s(pa3+s(pa4+s(pa5+spa6)))));
		134	Q = one+s(qa1+s(qa2+s(qa3+s(qa4+s(qa5+sqa6)))));
		135	if(hx>=0) return erx + P/Q; else return -erx - P/Q;
		136	}
		137	if (ix >= 0x40c00000) { /* inf>\|x\|>=6 */
		138	if(hx>=0) return one-tiny; else return tiny-one;
		139	}
		140	x = fabsf(x);
		141	s = one/(x*x);
		142	if(ix< 0x4036DB6E) { /* \|x\| < 1/0.35 */
		143	R=ra0+s(ra1+s(ra2+s(ra3+s(ra4+s*(
		144	ra5+s(ra6+sra7))))));
		145	S=one+s(sa1+s(sa2+s(sa3+s(sa4+s*(
		146	sa5+s(sa6+s(sa7+s*sa8)))))));
		147	} else { /* \|x\| >= 1/0.35 */
		148	R=rb0+s(rb1+s(rb2+s(rb3+s(rb4+s*(
		149	rb5+s*rb6)))));
		150	S=one+s(sb1+s(sb2+s(sb3+s(sb4+s*(
		151	sb5+s(sb6+ssb7))))));
		152	}
		153	GET_FLOAT_WORD(ix,x);
		154	SET_FLOAT_WORD(z,ix&0xfffff000);
		155	r = __ieee754_expf(-zz-(float)0.5625)__ieee754_expf((z-x)*(z+x)+R/S);
		156	if(hx>=0) return one-r/x; else return r/x-one;
		157	}
		158
		159	#ifdef __STDC__
		160	float erfcf(float x)
		161	#else
		162	float erfcf(x)
		163	float x;
		164	#endif
		165	{
		166	int32_t hx,ix;
		167	float R,S,P,Q,s,y,z,r;
		168	GET_FLOAT_WORD(hx,x);
		169	ix = hx&0x7fffffff;
		170	if(ix>=0x7f800000) { /* erfc(nan)=nan */
		171	/* erfc(+-inf)=0,2 */
		172	return (float)(((u_int32_t)hx>>31)<<1)+one/x;
		173	}
		174
		175	if(ix < 0x3f580000) { /* \|x\|<0.84375 */
		176	if(ix < 0x23800000) /* \|x\|<2*-56 /
		177	return one-x;
		178	z = x*x;
		179	r = pp0+z(pp1+z(pp2+z(pp3+zpp4)));
		180	s = one+z(qq1+z(qq2+z(qq3+z(qq4+z*qq5))));
		181	y = r/s;
		182	if(hx < 0x3e800000) { /* x<1/4 */
		183	return one-(x+x*y);
		184	} else {
		185	r = x*y;
		186	r += (x-half);
		187	return half - r ;
		188	}
		189	}
		190	if(ix < 0x3fa00000) { /* 0.84375 <= \|x\| < 1.25 */
		191	s = fabsf(x)-one;
		192	P = pa0+s(pa1+s(pa2+s(pa3+s(pa4+s(pa5+spa6)))));
		193	Q = one+s(qa1+s(qa2+s(qa3+s(qa4+s(qa5+sqa6)))));
		194	if(hx>=0) {
		195	z = one-erx; return z - P/Q;
		196	} else {
		197	z = erx+P/Q; return one+z;
		198	}
		199	}
		200	if (ix < 0x41e00000) { /* \|x\|<28 */
		201	x = fabsf(x);
		202	s = one/(x*x);
		203	if(ix< 0x4036DB6D) { /* \|x\| < 1/.35 ~ 2.857143*/
		204	R=ra0+s(ra1+s(ra2+s(ra3+s(ra4+s*(
		205	ra5+s(ra6+sra7))))));
		206	S=one+s(sa1+s(sa2+s(sa3+s(sa4+s*(
		207	sa5+s(sa6+s(sa7+s*sa8)))))));
		208	} else { /* \|x\| >= 1/.35 ~ 2.857143 */
		209	if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
		210	R=rb0+s(rb1+s(rb2+s(rb3+s(rb4+s*(
		211	rb5+s*rb6)))));
		212	S=one+s(sb1+s(sb2+s(sb3+s(sb4+s*(
		213	sb5+s(sb6+ssb7))))));
		214	}
		215	GET_FLOAT_WORD(ix,x);
		216	SET_FLOAT_WORD(z,ix&0xfffff000);
		217	r = __ieee754_expf(-zz-(float)0.5625)
		218	__ieee754_expf((z-x)*(z+x)+R/S);
		219	if(hx>0) return r/x; else return two-r/x;
		220	} else {
		221	if(hx>0) return tiny*tiny; else return two-tiny;
		222	}
		223	}

Subversion Repositories shark

(root)/shark/branches/xen/libc/arch/x86/libm/msun/src/s_erff.c - Rev 1619