1 /*
2 * Copyright (c) 1998, 2015, 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. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
62 private static final double INFINITY = Double.POSITIVE_INFINITY;
63
64 private FdLibm() {
65 throw new UnsupportedOperationException("No FdLibm instances for you.");
66 }
67
68 /**
69 * Return the low-order 32 bits of the double argument as an int.
70 */
71 private static int __LO(double x) {
72 long transducer = Double.doubleToRawLongBits(x);
73 return (int)transducer;
74 }
75
76 /**
77 * Return a double with its low-order bits of the second argument
78 * and the high-order bits of the first argument..
79 */
80 private static double __LO(double x, int low) {
81 long transX = Double.doubleToRawLongBits(x);
82 return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L)|low );
83 }
84
85 /**
86 * Return the high-order 32 bits of the double argument as an int.
87 */
88 private static int __HI(double x) {
89 long transducer = Double.doubleToRawLongBits(x);
90 return (int)(transducer >> 32);
91 }
92
93 /**
94 * Return a double with its high-order bits of the second argument
95 * and the low-order bits of the first argument..
96 */
97 private static double __HI(double x, int high) {
98 long transX = Double.doubleToRawLongBits(x);
99 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 );
100 }
101
102 /**
103 * cbrt(x)
104 * Return cube root of x
105 */
106 public static class Cbrt {
107 // unsigned
108 private static final int B1 = 715094163; /* B1 = (682-0.03306235651)*2**20 */
109 private static final int B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
110
111 private static final double C = 0x1.15f15f15f15f1p-1; // 19/35 ~= 5.42857142857142815906e-01
112 private static final double D = -0x1.691de2532c834p-1; // -864/1225 ~= 7.05306122448979611050e-01
113 private static final double E = 0x1.6a0ea0ea0ea0fp0; // 99/70 ~= 1.41428571428571436819e+00
114 private static final double F = 0x1.9b6db6db6db6ep0; // 45/28 ~= 1.60714285714285720630e+00
115 private static final double G = 0x1.6db6db6db6db7p-2; // 5/14 ~= 3.57142857142857150787e-01
116
117 public static strictfp double compute(double x) {
118 double t = 0.0;
119 double sign;
561 t = p_l + p_h;
562 t = __LO(t, 0);
563 u = t * LG2_H;
564 v = (p_l - (t - p_h)) * LG2 + t * LG2_L;
565 z = u + v;
566 w = v - (z - u);
567 t = z * z;
568 t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
569 r = (z * t1)/(t1 - 2.0) - (w + z * w);
570 z = 1.0 - (r - z);
571 j = __HI(z);
572 j += (n << 20);
573 if ((j >> 20) <= 0)
574 z = Math.scalb(z, n); // subnormal output
575 else {
576 int z_hi = __HI(z);
577 z_hi += (n << 20);
578 z = __HI(z, z_hi);
579 }
580 return s * z;
581 }
582 }
583 }
|
1 /*
2 * Copyright (c) 1998, 2016, 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. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
62 private static final double INFINITY = Double.POSITIVE_INFINITY;
63
64 private FdLibm() {
65 throw new UnsupportedOperationException("No FdLibm instances for you.");
66 }
67
68 /**
69 * Return the low-order 32 bits of the double argument as an int.
70 */
71 private static int __LO(double x) {
72 long transducer = Double.doubleToRawLongBits(x);
73 return (int)transducer;
74 }
75
76 /**
77 * Return a double with its low-order bits of the second argument
78 * and the high-order bits of the first argument..
79 */
80 private static double __LO(double x, int low) {
81 long transX = Double.doubleToRawLongBits(x);
82 return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L) |
83 (low & 0x0000_0000_FFFF_FFFFL));
84 }
85
86 /**
87 * Return the high-order 32 bits of the double argument as an int.
88 */
89 private static int __HI(double x) {
90 long transducer = Double.doubleToRawLongBits(x);
91 return (int)(transducer >> 32);
92 }
93
94 /**
95 * Return a double with its high-order bits of the second argument
96 * and the low-order bits of the first argument..
97 */
98 private static double __HI(double x, int high) {
99 long transX = Double.doubleToRawLongBits(x);
100 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL) |
101 ( ((long)high)) << 32 );
102 }
103
104 /**
105 * cbrt(x)
106 * Return cube root of x
107 */
108 public static class Cbrt {
109 // unsigned
110 private static final int B1 = 715094163; /* B1 = (682-0.03306235651)*2**20 */
111 private static final int B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
112
113 private static final double C = 0x1.15f15f15f15f1p-1; // 19/35 ~= 5.42857142857142815906e-01
114 private static final double D = -0x1.691de2532c834p-1; // -864/1225 ~= 7.05306122448979611050e-01
115 private static final double E = 0x1.6a0ea0ea0ea0fp0; // 99/70 ~= 1.41428571428571436819e+00
116 private static final double F = 0x1.9b6db6db6db6ep0; // 45/28 ~= 1.60714285714285720630e+00
117 private static final double G = 0x1.6db6db6db6db7p-2; // 5/14 ~= 3.57142857142857150787e-01
118
119 public static strictfp double compute(double x) {
120 double t = 0.0;
121 double sign;
563 t = p_l + p_h;
564 t = __LO(t, 0);
565 u = t * LG2_H;
566 v = (p_l - (t - p_h)) * LG2 + t * LG2_L;
567 z = u + v;
568 w = v - (z - u);
569 t = z * z;
570 t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
571 r = (z * t1)/(t1 - 2.0) - (w + z * w);
572 z = 1.0 - (r - z);
573 j = __HI(z);
574 j += (n << 20);
575 if ((j >> 20) <= 0)
576 z = Math.scalb(z, n); // subnormal output
577 else {
578 int z_hi = __HI(z);
579 z_hi += (n << 20);
580 z = __HI(z, z_hi);
581 }
582 return s * z;
583 }
584 }
585
586 /**
587 * Returns the exponential of x.
588 *
589 * Method
590 * 1. Argument reduction:
591 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
592 * Given x, find r and integer k such that
593 *
594 * x = k*ln2 + r, |r| <= 0.5*ln2.
595 *
596 * Here r will be represented as r = hi-lo for better
597 * accuracy.
598 *
599 * 2. Approximation of exp(r) by a special rational function on
600 * the interval [0,0.34658]:
601 * Write
602 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
603 * We use a special Reme algorithm on [0,0.34658] to generate
604 * a polynomial of degree 5 to approximate R. The maximum error
605 * of this polynomial approximation is bounded by 2**-59. In
606 * other words,
607 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
608 * (where z=r*r, and the values of P1 to P5 are listed below)
609 * and
610 * | 5 | -59
611 * | 2.0+P1*z+...+P5*z - R(z) | <= 2
612 * | |
613 * The computation of exp(r) thus becomes
614 * 2*r
615 * exp(r) = 1 + -------
616 * R - r
617 * r*R1(r)
618 * = 1 + r + ----------- (for better accuracy)
619 * 2 - R1(r)
620 * where
621 * 2 4 10
622 * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
623 *
624 * 3. Scale back to obtain exp(x):
625 * From step 1, we have
626 * exp(x) = 2^k * exp(r)
627 *
628 * Special cases:
629 * exp(INF) is INF, exp(NaN) is NaN;
630 * exp(-INF) is 0, and
631 * for finite argument, only exp(0)=1 is exact.
632 *
633 * Accuracy:
634 * according to an error analysis, the error is always less than
635 * 1 ulp (unit in the last place).
636 *
637 * Misc. info.
638 * For IEEE double
639 * if x > 7.09782712893383973096e+02 then exp(x) overflow
640 * if x < -7.45133219101941108420e+02 then exp(x) underflow
641 *
642 * Constants:
643 * The hexadecimal values are the intended ones for the following
644 * constants. The decimal values may be used, provided that the
645 * compiler will convert from decimal to binary accurately enough
646 * to produce the hexadecimal values shown.
647 */
648 static class Exp {
649 private static final double one = 1.0;
650 private static final double[] half = {0.5, -0.5,};
651 private static final double huge = 1.0e+300;
652 private static final double twom1000= 0x1.0p-1000; // 9.33263618503218878990e-302 = 2^-1000
653 private static final double o_threshold= 0x1.62e42fefa39efp9; // 7.09782712893383973096e+02
654 private static final double u_threshold= -0x1.74910d52d3051p9; // -7.45133219101941108420e+02;
655 private static final double[] ln2HI ={ 0x1.62e42feep-1, // 6.93147180369123816490e-01
656 -0x1.62e42feep-1}; // -6.93147180369123816490e-01
657 private static final double[] ln2LO ={ 0x1.a39ef35793c76p-33, // 1.90821492927058770002e-10
658 -0x1.a39ef35793c76p-33}; // -1.90821492927058770002e-10
659 private static final double invln2 = 0x1.71547652b82fep0; // 1.44269504088896338700e+00
660
661 private static final double P1 = 0x1.555555555553ep-3; // 1.66666666666666019037e-01
662 private static final double P2 = -0x1.6c16c16bebd93p-9; // -2.77777777770155933842e-03
663 private static final double P3 = 0x1.1566aaf25de2cp-14; // 6.61375632143793436117e-05
664 private static final double P4 = -0x1.bbd41c5d26bf1p-20; // -1.65339022054652515390e-06
665 private static final double P5 = 0x1.6376972bea4d0p-25; // 4.13813679705723846039e-08
666
667 // should be able to forgo strictfp due to controlled over/underflow
668 public static strictfp double compute(double x) {
669 double y;
670 double hi = 0.0;
671 double lo = 0.0;
672 double c;
673 double t;
674 int k = 0;
675 int xsb;
676 /*unsigned*/ int hx;
677
678 hx = __HI(x); /* high word of x */
679 xsb = (hx >> 31) & 1; /* sign bit of x */
680 hx &= 0x7fffffff; /* high word of |x| */
681
682 /* filter out non-finite argument */
683 if (hx >= 0x40862E42) { /* if |x| >= 709.78... */
684 if (hx >= 0x7ff00000) {
685 if (((hx & 0xfffff) | __LO(x)) != 0)
686 return x + x; /* NaN */
687 else
688 return (xsb == 0) ? x : 0.0; /* exp(+-inf) = {inf, 0} */
689 }
690 if (x > o_threshold)
691 return huge * huge; /* overflow */
692 if (x < u_threshold) // unsigned compare needed here?
693 return twom1000 * twom1000; /* underflow */
694 }
695
696 /* argument reduction */
697 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
698 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
699 hi = x - ln2HI[xsb];
700 lo=ln2LO[xsb];
701 k = 1 - xsb - xsb;
702 } else {
703 k = (int)(invln2 * x + half[xsb]);
704 t = k;
705 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
706 lo = t*ln2LO[0];
707 }
708 x = hi - lo;
709 } else if (hx < 0x3e300000) { /* when |x|<2**-28 */
710 if (huge + x > one)
711 return one + x; /* trigger inexact */
712 } else {
713 k = 0;
714 }
715
716 /* x is now in primary range */
717 t = x * x;
718 c = x - t*(P1 + t*(P2 + t*(P3 + t*(P4 + t*P5))));
719 if (k == 0)
720 return one - ((x*c)/(c - 2.0) - x);
721 else
722 y = one - ((lo - (x*c)/(2.0 - c)) - hi);
723
724 if(k >= -1021) {
725 y = __HI(y, __HI(y) + (k << 20)); /* add k to y's exponent */
726 return y;
727 } else {
728 y = __HI(y, __HI(y) + ((k + 1000) << 20)); /* add k to y's exponent */
729 return y * twom1000;
730 }
731 }
732 }
733 }
|