112 static double __ieee754_log(double x) { |
113 static double __ieee754_log(double x) { |
113 double hfsq,f,s,z,R,w,t1,t2,dk; |
114 double hfsq,f,s,z,R,w,t1,t2,dk; |
114 int k,hx,i,j; |
115 int k,hx,i,j; |
115 unsigned lx; |
116 unsigned lx; |
116 |
117 |
117 hx = __HI(x); /* high word of x */ |
118 hx = high(x); /* high word of x */ |
118 lx = __LO(x); /* low word of x */ |
119 lx = low(x); /* low word of x */ |
119 |
120 |
120 k=0; |
121 k=0; |
121 if (hx < 0x00100000) { /* x < 2**-1022 */ |
122 if (hx < 0x00100000) { /* x < 2**-1022 */ |
122 if (((hx&0x7fffffff)|lx)==0) |
123 if (((hx&0x7fffffff)|lx)==0) |
123 return -two54/zero; /* log(+-0)=-inf */ |
124 return -two54/zero; /* log(+-0)=-inf */ |
124 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
125 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
125 k -= 54; x *= two54; /* subnormal number, scale up x */ |
126 k -= 54; x *= two54; /* subnormal number, scale up x */ |
126 hx = __HI(x); /* high word of x */ |
127 hx = high(x); /* high word of x */ |
127 } |
128 } |
128 if (hx >= 0x7ff00000) return x+x; |
129 if (hx >= 0x7ff00000) return x+x; |
129 k += (hx>>20)-1023; |
130 k += (hx>>20)-1023; |
130 hx &= 0x000fffff; |
131 hx &= 0x000fffff; |
131 i = (hx+0x95f64)&0x100000; |
132 i = (hx+0x95f64)&0x100000; |
132 __HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */ |
133 set_high(&x, hx|(i^0x3ff00000)); /* normalize x or x/2 */ |
133 k += (i>>20); |
134 k += (i>>20); |
134 f = x-1.0; |
135 f = x-1.0; |
135 if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ |
136 if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ |
136 if(f==zero) { |
137 if(f==zero) { |
137 if (k==0) return zero; |
138 if (k==0) return zero; |
206 static double __ieee754_log10(double x) { |
207 static double __ieee754_log10(double x) { |
207 double y,z; |
208 double y,z; |
208 int i,k,hx; |
209 int i,k,hx; |
209 unsigned lx; |
210 unsigned lx; |
210 |
211 |
211 hx = __HI(x); /* high word of x */ |
212 hx = high(x); /* high word of x */ |
212 lx = __LO(x); /* low word of x */ |
213 lx = low(x); /* low word of x */ |
213 |
214 |
214 k=0; |
215 k=0; |
215 if (hx < 0x00100000) { /* x < 2**-1022 */ |
216 if (hx < 0x00100000) { /* x < 2**-1022 */ |
216 if (((hx&0x7fffffff)|lx)==0) |
217 if (((hx&0x7fffffff)|lx)==0) |
217 return -two54/zero; /* log(+-0)=-inf */ |
218 return -two54/zero; /* log(+-0)=-inf */ |
218 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
219 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
219 k -= 54; x *= two54; /* subnormal number, scale up x */ |
220 k -= 54; x *= two54; /* subnormal number, scale up x */ |
220 hx = __HI(x); /* high word of x */ |
221 hx = high(x); /* high word of x */ |
221 } |
222 } |
222 if (hx >= 0x7ff00000) return x+x; |
223 if (hx >= 0x7ff00000) return x+x; |
223 k += (hx>>20)-1023; |
224 k += (hx>>20)-1023; |
224 i = ((unsigned)k&0x80000000)>>31; |
225 i = ((unsigned)k&0x80000000)>>31; |
225 hx = (hx&0x000fffff)|((0x3ff-i)<<20); |
226 hx = (hx&0x000fffff)|((0x3ff-i)<<20); |
226 y = (double)(k+i); |
227 y = (double)(k+i); |
227 __HI(x) = hx; |
228 set_high(&x, hx); |
228 z = y*log10_2lo + ivln10*__ieee754_log(x); |
229 z = y*log10_2lo + ivln10*__ieee754_log(x); |
229 return z+y*log10_2hi; |
230 return z+y*log10_2hi; |
230 } |
231 } |
231 |
232 |
232 JRT_LEAF(jdouble, SharedRuntime::dlog10(jdouble x)) |
233 JRT_LEAF(jdouble, SharedRuntime::dlog10(jdouble x)) |
317 static double __ieee754_exp(double x) { |
318 static double __ieee754_exp(double x) { |
318 double y,hi=0,lo=0,c,t; |
319 double y,hi=0,lo=0,c,t; |
319 int k=0,xsb; |
320 int k=0,xsb; |
320 unsigned hx; |
321 unsigned hx; |
321 |
322 |
322 hx = __HI(x); /* high word of x */ |
323 hx = high(x); /* high word of x */ |
323 xsb = (hx>>31)&1; /* sign bit of x */ |
324 xsb = (hx>>31)&1; /* sign bit of x */ |
324 hx &= 0x7fffffff; /* high word of |x| */ |
325 hx &= 0x7fffffff; /* high word of |x| */ |
325 |
326 |
326 /* filter out non-finite argument */ |
327 /* filter out non-finite argument */ |
327 if(hx >= 0x40862E42) { /* if |x|>=709.78... */ |
328 if(hx >= 0x40862E42) { /* if |x|>=709.78... */ |
328 if(hx>=0x7ff00000) { |
329 if(hx>=0x7ff00000) { |
329 if(((hx&0xfffff)|__LO(x))!=0) |
330 if(((hx&0xfffff)|low(x))!=0) |
330 return x+x; /* NaN */ |
331 return x+x; /* NaN */ |
331 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ |
332 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ |
332 } |
333 } |
333 if(x > o_threshold) return hugeX*hugeX; /* overflow */ |
334 if(x > o_threshold) return hugeX*hugeX; /* overflow */ |
334 if(x < u_threshold) return twom1000*twom1000; /* underflow */ |
335 if(x < u_threshold) return twom1000*twom1000; /* underflow */ |
445 int i0,i1,i,j,k,yisint,n; |
446 int i0,i1,i,j,k,yisint,n; |
446 int hx,hy,ix,iy; |
447 int hx,hy,ix,iy; |
447 unsigned lx,ly; |
448 unsigned lx,ly; |
448 |
449 |
449 i0 = ((*(int*)&one)>>29)^1; i1=1-i0; |
450 i0 = ((*(int*)&one)>>29)^1; i1=1-i0; |
450 hx = __HI(x); lx = __LO(x); |
451 hx = high(x); lx = low(x); |
451 hy = __HI(y); ly = __LO(y); |
452 hy = high(y); ly = low(y); |
452 ix = hx&0x7fffffff; iy = hy&0x7fffffff; |
453 ix = hx&0x7fffffff; iy = hy&0x7fffffff; |
453 |
454 |
454 /* y==zero: x**0 = 1 */ |
455 /* y==zero: x**0 = 1 */ |
455 if((iy|ly)==0) return one; |
456 if((iy|ly)==0) return one; |
456 |
457 |
546 t = ax-one; /* t has 20 trailing zeros */ |
547 t = ax-one; /* t has 20 trailing zeros */ |
547 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); |
548 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); |
548 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ |
549 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ |
549 v = t*ivln2_l-w*ivln2; |
550 v = t*ivln2_l-w*ivln2; |
550 t1 = u+v; |
551 t1 = u+v; |
551 __LO(t1) = 0; |
552 set_low(&t1, 0); |
552 t2 = v-(t1-u); |
553 t2 = v-(t1-u); |
553 } else { |
554 } else { |
554 double ss,s2,s_h,s_l,t_h,t_l; |
555 double ss,s2,s_h,s_l,t_h,t_l; |
555 n = 0; |
556 n = 0; |
556 /* take care subnormal number */ |
557 /* take care subnormal number */ |
557 if(ix<0x00100000) |
558 if(ix<0x00100000) |
558 {ax *= two53; n -= 53; ix = __HI(ax); } |
559 {ax *= two53; n -= 53; ix = high(ax); } |
559 n += ((ix)>>20)-0x3ff; |
560 n += ((ix)>>20)-0x3ff; |
560 j = ix&0x000fffff; |
561 j = ix&0x000fffff; |
561 /* determine interval */ |
562 /* determine interval */ |
562 ix = j|0x3ff00000; /* normalize ix */ |
563 ix = j|0x3ff00000; /* normalize ix */ |
563 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ |
564 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ |
564 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ |
565 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ |
565 else {k=0;n+=1;ix -= 0x00100000;} |
566 else {k=0;n+=1;ix -= 0x00100000;} |
566 __HI(ax) = ix; |
567 set_high(&ax, ix); |
567 |
568 |
568 /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ |
569 /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ |
569 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ |
570 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ |
570 v = one/(ax+bp[k]); |
571 v = one/(ax+bp[k]); |
571 ss = u*v; |
572 ss = u*v; |
572 s_h = ss; |
573 s_h = ss; |
573 __LO(s_h) = 0; |
574 set_low(&s_h, 0); |
574 /* t_h=ax+bp[k] High */ |
575 /* t_h=ax+bp[k] High */ |
575 t_h = zeroX; |
576 t_h = zeroX; |
576 __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); |
577 set_high(&t_h, ((ix>>1)|0x20000000)+0x00080000+(k<<18)); |
577 t_l = ax - (t_h-bp[k]); |
578 t_l = ax - (t_h-bp[k]); |
578 s_l = v*((u-s_h*t_h)-s_h*t_l); |
579 s_l = v*((u-s_h*t_h)-s_h*t_l); |
579 /* compute log(ax) */ |
580 /* compute log(ax) */ |
580 s2 = ss*ss; |
581 s2 = ss*ss; |
581 r = s2*s2*(L1X+s2*(L2X+s2*(L3X+s2*(L4X+s2*(L5X+s2*L6X))))); |
582 r = s2*s2*(L1X+s2*(L2X+s2*(L3X+s2*(L4X+s2*(L5X+s2*L6X))))); |
582 r += s_l*(s_h+ss); |
583 r += s_l*(s_h+ss); |
583 s2 = s_h*s_h; |
584 s2 = s_h*s_h; |
584 t_h = 3.0+s2+r; |
585 t_h = 3.0+s2+r; |
585 __LO(t_h) = 0; |
586 set_low(&t_h, 0); |
586 t_l = r-((t_h-3.0)-s2); |
587 t_l = r-((t_h-3.0)-s2); |
587 /* u+v = ss*(1+...) */ |
588 /* u+v = ss*(1+...) */ |
588 u = s_h*t_h; |
589 u = s_h*t_h; |
589 v = s_l*t_h+t_l*ss; |
590 v = s_l*t_h+t_l*ss; |
590 /* 2/(3log2)*(ss+...) */ |
591 /* 2/(3log2)*(ss+...) */ |
591 p_h = u+v; |
592 p_h = u+v; |
592 __LO(p_h) = 0; |
593 set_low(&p_h, 0); |
593 p_l = v-(p_h-u); |
594 p_l = v-(p_h-u); |
594 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ |
595 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ |
595 z_l = cp_l*p_h+p_l*cp+dp_l[k]; |
596 z_l = cp_l*p_h+p_l*cp+dp_l[k]; |
596 /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ |
597 /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ |
597 t = (double)n; |
598 t = (double)n; |
598 t1 = (((z_h+z_l)+dp_h[k])+t); |
599 t1 = (((z_h+z_l)+dp_h[k])+t); |
599 __LO(t1) = 0; |
600 set_low(&t1, 0); |
600 t2 = z_l-(((t1-t)-dp_h[k])-z_h); |
601 t2 = z_l-(((t1-t)-dp_h[k])-z_h); |
601 } |
602 } |
602 |
603 |
603 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ |
604 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ |
604 y1 = y; |
605 y1 = y; |
605 __LO(y1) = 0; |
606 set_low(&y1, 0); |
606 p_l = (y-y1)*t1+y*t2; |
607 p_l = (y-y1)*t1+y*t2; |
607 p_h = y1*t1; |
608 p_h = y1*t1; |
608 z = p_l+p_h; |
609 z = p_l+p_h; |
609 j = __HI(z); |
610 j = high(z); |
610 i = __LO(z); |
611 i = low(z); |
611 if (j>=0x40900000) { /* z >= 1024 */ |
612 if (j>=0x40900000) { /* z >= 1024 */ |
612 if(((j-0x40900000)|i)!=0) /* if z > 1024 */ |
613 if(((j-0x40900000)|i)!=0) /* if z > 1024 */ |
613 return s*hugeX*hugeX; /* overflow */ |
614 return s*hugeX*hugeX; /* overflow */ |
614 else { |
615 else { |
615 if(p_l+ovt>z-p_h) return s*hugeX*hugeX; /* overflow */ |
616 if(p_l+ovt>z-p_h) return s*hugeX*hugeX; /* overflow */ |
629 n = 0; |
630 n = 0; |
630 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ |
631 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ |
631 n = j+(0x00100000>>(k+1)); |
632 n = j+(0x00100000>>(k+1)); |
632 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ |
633 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ |
633 t = zeroX; |
634 t = zeroX; |
634 __HI(t) = (n&~(0x000fffff>>k)); |
635 set_high(&t, (n&~(0x000fffff>>k))); |
635 n = ((n&0x000fffff)|0x00100000)>>(20-k); |
636 n = ((n&0x000fffff)|0x00100000)>>(20-k); |
636 if(j<0) n = -n; |
637 if(j<0) n = -n; |
637 p_h -= t; |
638 p_h -= t; |
638 } |
639 } |
639 t = p_l+p_h; |
640 t = p_l+p_h; |
640 __LO(t) = 0; |
641 set_low(&t, 0); |
641 u = t*lg2_h; |
642 u = t*lg2_h; |
642 v = (p_l-(t-p_h))*lg2+t*lg2_l; |
643 v = (p_l-(t-p_h))*lg2+t*lg2_l; |
643 z = u+v; |
644 z = u+v; |
644 w = v-(z-u); |
645 w = v-(z-u); |
645 t = z*z; |
646 t = z*z; |
646 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); |
647 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); |
647 r = (z*t1)/(t1-two)-(w+z*w); |
648 r = (z*t1)/(t1-two)-(w+z*w); |
648 z = one-(r-z); |
649 z = one-(r-z); |
649 j = __HI(z); |
650 j = high(z); |
650 j += (n<<20); |
651 j += (n<<20); |
651 if((j>>20)<=0) z = scalbnA(z,n); /* subnormal output */ |
652 if((j>>20)<=0) z = scalbnA(z,n); /* subnormal output */ |
652 else __HI(z) += (n<<20); |
653 else set_high(&z, high(z) + (n<<20)); |
653 return s*z; |
654 return s*z; |
654 } |
655 } |
655 |
656 |
656 |
657 |
657 JRT_LEAF(jdouble, SharedRuntime::dpow(jdouble x, jdouble y)) |
658 JRT_LEAF(jdouble, SharedRuntime::dpow(jdouble x, jdouble y)) |