1 /*
2 * Copyright (c) 2005, 2010, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation.
8 *
9 * This code is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
12 * version 2 for more details (a copy is included in the LICENSE file that
13 * accompanied this code).
14 *
15 * You should have received a copy of the GNU General Public License version
16 * 2 along with this work; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
18 *
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
20 * or visit www.oracle.com if you need additional information or have any
21 * questions.
22 *
39 // Enabling optimizations in this file causes incorrect code to be
40 // generated; can not figure out how to turn down optimization for one
41 // file in the IDE on Windows
42 #ifdef WIN32
43 # pragma optimize ( "", off )
44 #endif
45
46 #include <math.h>
47
48 // VM_LITTLE_ENDIAN is #defined appropriately in the Makefiles
49 // [jk] this is not 100% correct because the float word order may different
50 // from the byte order (e.g. on ARM)
51 #ifdef VM_LITTLE_ENDIAN
52 # define __HI(x) *(1+(int*)&x)
53 # define __LO(x) *(int*)&x
54 #else
55 # define __HI(x) *(int*)&x
56 # define __LO(x) *(1+(int*)&x)
57 #endif
58
59 #if !defined(AIX)
60 double copysign(double x, double y) {
61 __HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000);
62 return x;
63 }
64 #endif
65
66 /*
67 * ====================================================
68 * Copyright (c) 1998 Oracle and/or its affiliates. All rights reserved.
69 *
70 * Developed at SunSoft, a Sun Microsystems, Inc. business.
71 * Permission to use, copy, modify, and distribute this
72 * software is freely granted, provided that this notice
73 * is preserved.
74 * ====================================================
75 */
76
77 /*
78 * scalbn (double x, int n)
79 * scalbn(x,n) returns x* 2**n computed by exponent
80 * manipulation rather than by actually performing an
81 * exponentiation or a multiplication.
82 */
83
84 static const double
85 two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
86 twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
87 hugeX = 1.0e+300,
88 tiny = 1.0e-300;
89
90 #if !defined(AIX)
91 double scalbn (double x, int n) {
92 int k,hx,lx;
93 hx = __HI(x);
94 lx = __LO(x);
95 k = (hx&0x7ff00000)>>20; /* extract exponent */
96 if (k==0) { /* 0 or subnormal x */
97 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
98 x *= two54;
99 hx = __HI(x);
100 k = ((hx&0x7ff00000)>>20) - 54;
101 if (n< -50000) return tiny*x; /*underflow*/
102 }
103 if (k==0x7ff) return x+x; /* NaN or Inf */
104 k = k+n;
105 if (k > 0x7fe) return hugeX*copysign(hugeX,x); /* overflow */
106 if (k > 0) /* normal result */
107 {__HI(x) = (hx&0x800fffff)|(k<<20); return x;}
108 if (k <= -54) {
109 if (n > 50000) /* in case integer overflow in n+k */
110 return hugeX*copysign(hugeX,x); /*overflow*/
111 else return tiny*copysign(tiny,x); /*underflow*/
112 }
113 k += 54; /* subnormal result */
114 __HI(x) = (hx&0x800fffff)|(k<<20);
115 return x*twom54;
116 }
117 #endif
118
119 /* __ieee754_log(x)
120 * Return the logarithm of x
121 *
122 * Method :
123 * 1. Argument Reduction: find k and f such that
124 * x = 2^k * (1+f),
125 * where sqrt(2)/2 < 1+f < sqrt(2) .
126 *
127 * 2. Approximation of log(1+f).
128 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
129 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
130 * = 2s + s*R
131 * We use a special Reme algorithm on [0,0.1716] to generate
132 * a polynomial of degree 14 to approximate R The maximum error
133 * of this polynomial approximation is bounded by 2**-58.45. In
134 * other words,
135 * 2 4 6 8 10 12 14
136 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
137 * (the values of Lg1 to Lg7 are listed in the program)
702 n = j+(0x00100000>>(k+1));
703 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
704 t = zeroX;
705 __HI(t) = (n&~(0x000fffff>>k));
706 n = ((n&0x000fffff)|0x00100000)>>(20-k);
707 if(j<0) n = -n;
708 p_h -= t;
709 }
710 t = p_l+p_h;
711 __LO(t) = 0;
712 u = t*lg2_h;
713 v = (p_l-(t-p_h))*lg2+t*lg2_l;
714 z = u+v;
715 w = v-(z-u);
716 t = z*z;
717 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
718 r = (z*t1)/(t1-two)-(w+z*w);
719 z = one-(r-z);
720 j = __HI(z);
721 j += (n<<20);
722 if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
723 else __HI(z) += (n<<20);
724 return s*z;
725 }
726
727
728 JRT_LEAF(jdouble, SharedRuntime::dpow(jdouble x, jdouble y))
729 return __ieee754_pow(x, y);
730 JRT_END
731
732 #ifdef WIN32
733 # pragma optimize ( "", on )
734 #endif
|
1 /*
2 * Copyright (c) 2005, 2014, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation.
8 *
9 * This code is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
12 * version 2 for more details (a copy is included in the LICENSE file that
13 * accompanied this code).
14 *
15 * You should have received a copy of the GNU General Public License version
16 * 2 along with this work; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
18 *
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
20 * or visit www.oracle.com if you need additional information or have any
21 * questions.
22 *
39 // Enabling optimizations in this file causes incorrect code to be
40 // generated; can not figure out how to turn down optimization for one
41 // file in the IDE on Windows
42 #ifdef WIN32
43 # pragma optimize ( "", off )
44 #endif
45
46 #include <math.h>
47
48 // VM_LITTLE_ENDIAN is #defined appropriately in the Makefiles
49 // [jk] this is not 100% correct because the float word order may different
50 // from the byte order (e.g. on ARM)
51 #ifdef VM_LITTLE_ENDIAN
52 # define __HI(x) *(1+(int*)&x)
53 # define __LO(x) *(int*)&x
54 #else
55 # define __HI(x) *(int*)&x
56 # define __LO(x) *(1+(int*)&x)
57 #endif
58
59 static double copysignA(double x, double y) {
60 __HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000);
61 return x;
62 }
63
64 /*
65 * ====================================================
66 * Copyright (c) 1998 Oracle and/or its affiliates. All rights reserved.
67 *
68 * Developed at SunSoft, a Sun Microsystems, Inc. business.
69 * Permission to use, copy, modify, and distribute this
70 * software is freely granted, provided that this notice
71 * is preserved.
72 * ====================================================
73 */
74
75 /*
76 * scalbn (double x, int n)
77 * scalbn(x,n) returns x* 2**n computed by exponent
78 * manipulation rather than by actually performing an
79 * exponentiation or a multiplication.
80 */
81
82 static const double
83 two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
84 twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
85 hugeX = 1.0e+300,
86 tiny = 1.0e-300;
87
88 static double scalbnA (double x, int n) {
89 int k,hx,lx;
90 hx = __HI(x);
91 lx = __LO(x);
92 k = (hx&0x7ff00000)>>20; /* extract exponent */
93 if (k==0) { /* 0 or subnormal x */
94 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
95 x *= two54;
96 hx = __HI(x);
97 k = ((hx&0x7ff00000)>>20) - 54;
98 if (n< -50000) return tiny*x; /*underflow*/
99 }
100 if (k==0x7ff) return x+x; /* NaN or Inf */
101 k = k+n;
102 if (k > 0x7fe) return hugeX*copysignA(hugeX,x); /* overflow */
103 if (k > 0) /* normal result */
104 {__HI(x) = (hx&0x800fffff)|(k<<20); return x;}
105 if (k <= -54) {
106 if (n > 50000) /* in case integer overflow in n+k */
107 return hugeX*copysignA(hugeX,x); /*overflow*/
108 else return tiny*copysignA(tiny,x); /*underflow*/
109 }
110 k += 54; /* subnormal result */
111 __HI(x) = (hx&0x800fffff)|(k<<20);
112 return x*twom54;
113 }
114
115 /* __ieee754_log(x)
116 * Return the logarithm of x
117 *
118 * Method :
119 * 1. Argument Reduction: find k and f such that
120 * x = 2^k * (1+f),
121 * where sqrt(2)/2 < 1+f < sqrt(2) .
122 *
123 * 2. Approximation of log(1+f).
124 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
125 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
126 * = 2s + s*R
127 * We use a special Reme algorithm on [0,0.1716] to generate
128 * a polynomial of degree 14 to approximate R The maximum error
129 * of this polynomial approximation is bounded by 2**-58.45. In
130 * other words,
131 * 2 4 6 8 10 12 14
132 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
133 * (the values of Lg1 to Lg7 are listed in the program)
698 n = j+(0x00100000>>(k+1));
699 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
700 t = zeroX;
701 __HI(t) = (n&~(0x000fffff>>k));
702 n = ((n&0x000fffff)|0x00100000)>>(20-k);
703 if(j<0) n = -n;
704 p_h -= t;
705 }
706 t = p_l+p_h;
707 __LO(t) = 0;
708 u = t*lg2_h;
709 v = (p_l-(t-p_h))*lg2+t*lg2_l;
710 z = u+v;
711 w = v-(z-u);
712 t = z*z;
713 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
714 r = (z*t1)/(t1-two)-(w+z*w);
715 z = one-(r-z);
716 j = __HI(z);
717 j += (n<<20);
718 if((j>>20)<=0) z = scalbnA(z,n); /* subnormal output */
719 else __HI(z) += (n<<20);
720 return s*z;
721 }
722
723
724 JRT_LEAF(jdouble, SharedRuntime::dpow(jdouble x, jdouble y))
725 return __ieee754_pow(x, y);
726 JRT_END
727
728 #ifdef WIN32
729 # pragma optimize ( "", on )
730 #endif
|