procyberian/gcc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[multiple changes]

Browse source 
Fri May 28 22:20:03 1999  Anthony Green  <green@cygnus.com>
	* java/lang/fdlibm.h: Don't use __uint32_t.  Include mprec.h.
	* java/lang/e_log.c: Don't use __uint32_t.
1999-05-27  Eric Christopher <echristo@cygnus.com>
	* configure: Rebuilt
	* configure.in: Fixed ISO C9X and namespace collision with __uint32_t
	* acconfig.h: Rebuilt
	* include/config.h.in: Rebuilt
	* java/lang/mprec.h, java/lang/e_acos.c, java/lang/e_asin.c,
 	java/lang/e_atan2.c, java/lang/e_exp.c, java/lang/e_fmod.c,
 	e_log.c, java/lang/e_pow.c, java/lang/e_rem_pio2.c,
 	java/lang/e_remainder.c, java/lang/e_sqrt.c, java/lang/fdlibm.h,
 	k_tan.c, java/lang/mprec.h, java/lang/s_atan.c,
 	java/lang/s_ceil.c, java/lang/s_copysign.c, java/lang/s_fabs.c,
 	s_floor.c, java/lang/s_rint.c, java/lang/sf_rint.c: Fixed ISO C9X
 	and namespace collision with __uint32_t

From-SVN: r27729
This commit is contained in:
Tom Tromey 1999-06-24 20:06:09 +00:00
parent fe574d5d92
commit 0d16618c58
31 changed files with 614 additions and 526 deletions
Download patch file Download diff file Expand all files Collapse all files

40
libjava/java/lang/e_exp.c
View file

@ -6,7 +6,7 @@
*
* Developed at SunPro, 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.
* ====================================================
*/
@ -19,36 +19,36 @@
* Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
* Given x, find r and integer k such that
*
* x = k*ln2 + r, |r| <= 0.5*ln2.
* x = k*ln2 + r, |r| <= 0.5*ln2.
*
* Here r will be represented as r = hi-lo for better
* Here r will be represented as r = hi-lo for better
* accuracy.
*
* 2. Approximation of exp(r) by a special rational function on
* the interval [0,0.34658]:
* Write
* R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
* We use a special Reme algorithm on [0,0.34658] to generate
* a polynomial of degree 5 to approximate R. The maximum error
* We use a special Reme algorithm on [0,0.34658] to generate
* a polynomial of degree 5 to approximate R. The maximum error
* of this polynomial approximation is bounded by 2**-59. In
* other words,
* R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
* (where z=r*r, and the values of P1 to P5 are listed below)
* and
* | 5 | -59
* | 2.0+P1*z+...+P5*z - R(z) | <= 2
* | 2.0+P1*z+...+P5*z - R(z) | <= 2
* | |
* The computation of exp(r) thus becomes
* 2*r
* exp(r) = 1 + -------
* R - r
* r*R1(r)
* r*R1(r)
* = 1 + r + ----------- (for better accuracy)
* 2 - R1(r)
* where
* 2 4 10
* R1(r) = r - (P1*r + P2*r + ... + P5*r ).
*
*
* 3. Scale back to obtain exp(x):
* From step 1, we have
* exp(x) = 2^k * exp(r)
@ -63,13 +63,13 @@
* 1 ulp (unit in the last place).
*
* Misc. info.
* For IEEE double
* For IEEE double
* if x > 7.09782712893383973096e+02 then exp(x) overflow
* if x < -7.45133219101941108420e+02 then exp(x) underflow
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
@ -109,8 +109,8 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
#endif
{
double y,hi,lo,c,t;
__int32_t k,xsb;
__uint32_t hx;
int32_t k,xsb;
uint32_t hx;
GET_HIGH_WORD(hx,x);
xsb = (hx>>31)&1; /* sign bit of x */
@ -119,9 +119,9 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
/* filter out non-finite argument */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
__uint32_t lx;
uint32_t lx;
GET_LOW_WORD(lx,x);
if(((hx&0xfffff)|lx)!=0)
if(((hx&0xfffff)|lx)!=0)
return x+x; /* NaN */
else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
}
@ -130,7 +130,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
}
/* argument reduction */
if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
} else {
@ -140,7 +140,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
lo = t*ln2LO[0];
}
x = hi - lo;
}
}
else if(hx < 0x3e300000) { /* when |x|<2**-28 */
if(huge+x>one) return one+x;/* trigger inexact */
}
@ -149,15 +149,15 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
/* x is now in primary range */
t = x*x;
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if(k==0) return one-((x*c)/(c-2.0)-x);
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);
if(k >= -1021) {
__uint32_t hy;
uint32_t hy;
GET_HIGH_WORD(hy,y);
SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */
return y;
} else {
__uint32_t hy;
uint32_t hy;
GET_HIGH_WORD(hy,y);
SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */
return y*twom1000;