8301202: Port fdlibm log to Java

Reviewed-by: bpb

src/java.base/share/classes/java/lang/FdLibm.java

@@ -747,6 +747,146 @@ class FdLibm {
         }
     }
 
+    /**
+     * Return the (natural) logarithm of x
+     *
+     * Method :
+     *   1. Argument Reduction: find k and f such that
+     *                      x = 2^k * (1+f),
+     *         where  sqrt(2)/2 < 1+f < sqrt(2) .
+     *
+     *   2. Approximation of log(1+f).
+     *      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+     *               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+     *               = 2s + s*R
+     *      We use a special Reme algorithm on [0,0.1716] to generate
+     *      a polynomial of degree 14 to approximate R The maximum error
+     *      of this polynomial approximation is bounded by 2**-58.45. In
+     *      other words,
+     *                      2      4      6      8      10      12      14
+     *          R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+     *      (the values of Lg1 to Lg7 are listed in the program)
+     *      and
+     *          |      2          14          |     -58.45
+     *          | Lg1*s +...+Lg7*s    -  R(z) | <= 2
+     *          |                             |
+     *      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+     *      In order to guarantee error in log below 1ulp, we compute log
+     *      by
+     *              log(1+f) = f - s*(f - R)        (if f is not too large)
+     *              log(1+f) = f - (hfsq - s*(hfsq+R)).     (better accuracy)
+     *
+     *      3. Finally,  log(x) = k*ln2 + log(1+f).
+     *                          = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+     *         Here ln2 is split into two floating point number:
+     *                      ln2_hi + ln2_lo,
+     *         where n*ln2_hi is always exact for |n| < 2000.
+     *
+     * Special cases:
+     *      log(x) is NaN with signal if x < 0 (including -INF) ;
+     *      log(+INF) is +INF; log(0) is -INF with signal;
+     *      log(NaN) is that NaN with no signal.
+     *
+     * Accuracy:
+     *      according to an error analysis, the error is always less than
+     *      1 ulp (unit in the last place).
+     *
+     * Constants:
+     * 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.
+     */
+    static final class Log {
+        private Log() {throw new UnsupportedOperationException();}
+
+        private static final double
+            ln2_hi = 0x1.62e42feep-1,       // 6.93147180369123816490e-01
+            ln2_lo = 0x1.a39ef35793c76p-33, // 1.90821492927058770002e-10
+
+            Lg1    = 0x1.5555555555593p-1,  // 6.666666666666735130e-01
+            Lg2    = 0x1.999999997fa04p-2,  // 3.999999999940941908e-01
+            Lg3    = 0x1.2492494229359p-2,  // 2.857142874366239149e-01
+            Lg4    = 0x1.c71c51d8e78afp-3,  // 2.222219843214978396e-01
+            Lg5    = 0x1.7466496cb03dep-3,  // 1.818357216161805012e-01
+            Lg6    = 0x1.39a09d078c69fp-3,  // 1.531383769920937332e-01
+            Lg7    = 0x1.2f112df3e5244p-3;  // 1.479819860511658591e-01
+
+        private static final double zero = 0.0;
+
+        static double compute(double x) {
+            double hfsq, f, s, z, R, w, t1, t2, dk;
+            int k, hx, i, j;
+            /*unsigned*/ int lx;
+
+            hx = __HI(x);           // high word of x
+            lx = __LO(x);           // low  word of x
+
+            k=0;
+            if (hx < 0x0010_0000) {                  // x < 2**-1022
+                if (((hx & 0x7fff_ffff) | lx) == 0) { // log(+-0) = -inf
+                    return -TWO54/zero;
+                }
+                if (hx < 0) {                        // log(-#) = NaN
+                    return (x - x)/zero;
+                }
+                k -= 54;
+                x *= TWO54;    // subnormal number, scale up x
+                hx = __HI(x);  // high word of x
+            }
+            if (hx >= 0x7ff0_0000) {
+                return x + x;
+            }
+            k += (hx >> 20) - 1023;
+            hx &= 0x000f_ffff;
+            i = (hx + 0x9_5f64) & 0x10_0000;
+            x =__HI(x, hx | (i ^ 0x3ff0_0000));  // normalize x or x/2
+            k += (i >> 20);
+            f = x - 1.0;
+            if ((0x000f_ffff & (2 + hx)) < 3) {// |f| < 2**-20
+                if (f == zero) {
+                    if (k == 0) {
+                        return zero;
+                    } else {
+                        dk = (double)k;
+                        return dk*ln2_hi + dk*ln2_lo;
+                    }
+                }
+                R = f*f*(0.5 - 0.33333333333333333*f);
+                if (k == 0) {
+                    return f - R;
+                } else {
+                    dk = (double)k;
+                    return dk*ln2_hi - ((R - dk*ln2_lo) - f);
+                }
+            }
+            s = f/(2.0 + f);
+            dk = (double)k;
+            z = s*s;
+            i = hx - 0x6_147a;
+            w = z*z;
+            j = 0x6b851 - hx;
+            t1= w*(Lg2 + w*(Lg4 + w*Lg6));
+            t2= z*(Lg1 + w*(Lg3 + w*(Lg5 + w*Lg7)));
+            i |= j;
+            R = t2 + t1;
+            if (i > 0) {
+                hfsq = 0.5*f*f;
+                if (k == 0) {
+                    return f-(hfsq - s*(hfsq + R));
+                } else {
+                    return dk*ln2_hi - ((hfsq - (s*(hfsq + R) + dk*ln2_lo)) - f);
+                }
+            } else {
+                if (k == 0) {
+                    return f - s*(f - R);
+                } else {
+                    return dk*ln2_hi - ((s*(f - R) - dk*ln2_lo) - f);
+                }
+            }
+        }
+    }
+
     /**
      * Return the base 10 logarithm of x
      *