Use a faster but less accurate log algorithm for computing Geotile Y coordinate (#87515)

This commit introduces a new algorithm to ESSloppyMath to compute logarithm (base e)

docs/changelog/87515.yaml

@@ -0,0 +1,5 @@
+pr: 87515
+summary: Use a faster but less accurate log algorithm for computing Geotile Y coordinate
+area: Geo
+type: enhancement
+issues: []

libs/core/src/main/java/org/elasticsearch/core/ESSloppyMath.java

@@ -22,4 +22,8 @@ public class ESSloppyMath {
     public static double atan(double value) {
         return FastMath.atan(value);
     }
+
+    public static double log(double value) {
+        return FastMath.log(value);
+    }
 }

libs/core/src/main/java/org/elasticsearch/core/FastMath.java

@@ -48,6 +48,14 @@ final class FastMath {
     private static final double TWO_POW_N28 = Double.longBitsToDouble(0x3E30000000000000L);
     private static final double TWO_POW_66 = Double.longBitsToDouble(0x4410000000000000L);
     private static final double LOG_DOUBLE_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
+    // First double value (from zero) such as (value+-1/value == value).
+    private static final double TWO_POW_27 = Double.longBitsToDouble(0x41A0000000000000L);
+    private static final double TWO_POW_52 = Double.longBitsToDouble(0x4330000000000000L);
+    // Smallest double normal value.
+    private static final double MIN_DOUBLE_NORMAL = Double.longBitsToDouble(0x0010000000000000L); // 2.2250738585072014E-308
+    private static final int MIN_DOUBLE_EXPONENT = -1074;
+    private static final int MAX_DOUBLE_EXPONENT = 1023;
+    private static final double LOG_2 = StrictMath.log(2.0);
 
     // --------------------------------------------------------------------------
     // CONSTANTS AND TABLES FOR ATAN
@@ -84,6 +92,22 @@ final class FastMath {
     private static final double ATAN_AT9 = Double.longBitsToDouble(0xbfa2b4442c6a6c2fL); // -3.65315727442169155270e-02
     private static final double ATAN_AT10 = Double.longBitsToDouble(0x3f90ad3ae322da11L); // 1.62858201153657823623e-02
 
+    // --------------------------------------------------------------------------
+    // CONSTANTS AND TABLES FOR LOG AND LOG1P
+    // --------------------------------------------------------------------------
+
+    private static final int LOG_BITS = 12;
+    private static final int LOG_TAB_SIZE = (1 << LOG_BITS);
+    private static final double[] logXLogTab = new double[LOG_TAB_SIZE];
+    private static final double[] logXTab = new double[LOG_TAB_SIZE];
+    private static final double[] logXInvTab = new double[LOG_TAB_SIZE];
+
+    // --------------------------------------------------------------------------
+    // TABLE FOR POWERS OF TWO
+    // --------------------------------------------------------------------------
+
+    private static final double[] twoPowTab = new double[(MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT) + 1];
+
     static {
         // atan
         for (int i = 0; i < ATAN_TABS_SIZE; i++) {
@@ -99,6 +123,23 @@ final class FastMath {
             atanDer3DivF3Tab[i] = ((-2 + 6 * x * x) * onePlusXSqInv3) * ONE_DIV_F3;
             atanDer4DivF4Tab[i] = ((24 * x * (1 - x * x)) * onePlusXSqInv4) * ONE_DIV_F4;
         }
+
+        // log
+
+        for (int i = 0; i < LOG_TAB_SIZE; i++) {
+            // Exact to use inverse of tab size, since it is a power of two.
+            double x = 1 + i * (1.0 / LOG_TAB_SIZE);
+            logXLogTab[i] = StrictMath.log(x);
+            logXTab[i] = x;
+            logXInvTab[i] = 1 / x;
+        }
+
+        // twoPow
+
+        for (int i = MIN_DOUBLE_EXPONENT; i <= MAX_DOUBLE_EXPONENT; i++) {
+            twoPowTab[i - MIN_DOUBLE_EXPONENT] = StrictMath.pow(2.0, i);
+        }
+
     }
 
     /**
@@ -175,4 +216,69 @@ final class FastMath {
         }
     }
 
+    /**
+     * @param value A double value.
+     * @return Value logarithm (base e).
+     */
+    public static double log(double value) {
+        if (value > 0.0) {
+            if (value == Double.POSITIVE_INFINITY) {
+                return Double.POSITIVE_INFINITY;
+            }
+
+            // For normal values not close to 1.0, we use the following formula:
+            // log(value)
+            // = log(2^exponent*1.mantissa)
+            // = log(2^exponent) + log(1.mantissa)
+            // = exponent * log(2) + log(1.mantissa)
+            // = exponent * log(2) + log(1.mantissaApprox) + log(1.mantissa/1.mantissaApprox)
+            // = exponent * log(2) + log(1.mantissaApprox) + log(1+epsilon)
+            // = exponent * log(2) + log(1.mantissaApprox) + epsilon-epsilon^2/2+epsilon^3/3-epsilon^4/4+...
+            // with:
+            // 1.mantissaApprox <= 1.mantissa,
+            // log(1.mantissaApprox) in table,
+            // epsilon = (1.mantissa/1.mantissaApprox)-1
+            //
+            // To avoid bad relative error for small results,
+            // values close to 1.0 are treated aside, with the formula:
+            // log(x) = z*(2+z^2*((2.0/3)+z^2*((2.0/5))+z^2*((2.0/7))+...)))
+            // with z=(x-1)/(x+1)
+
+            double h;
+            if (value > 0.95) {
+                if (value < 1.14) {
+                    double z = (value - 1.0) / (value + 1.0);
+                    double z2 = z * z;
+                    return z * (2 + z2 * ((2.0 / 3) + z2 * ((2.0 / 5) + z2 * ((2.0 / 7) + z2 * ((2.0 / 9) + z2 * ((2.0 / 11)))))));
+                }
+                h = 0.0;
+            } else if (value < MIN_DOUBLE_NORMAL) {
+                // Ensuring value is normal.
+                value *= TWO_POW_52;
+                // log(x*2^52)
+                // = log(x)-ln(2^52)
+                // = log(x)-52*ln(2)
+                h = -52 * LOG_2;
+            } else {
+                h = 0.0;
+            }
+
+            int valueBitsHi = (int) (Double.doubleToRawLongBits(value) >> 32);
+            int valueExp = (valueBitsHi >> 20) - MAX_DOUBLE_EXPONENT;
+            // Getting the first LOG_BITS bits of the mantissa.
+            int xIndex = ((valueBitsHi << 12) >>> (32 - LOG_BITS));
+
+            // 1.mantissa/1.mantissaApprox - 1
+            double z = (value * twoPowTab[-valueExp - MIN_DOUBLE_EXPONENT]) * logXInvTab[xIndex] - 1;
+
+            z *= (1 - z * ((1.0 / 2) - z * ((1.0 / 3))));
+
+            return h + valueExp * LOG_2 + (logXLogTab[xIndex] + z);
+
+        } else if (value == 0.0) {
+            return Double.NEGATIVE_INFINITY;
+        } else { // value < 0.0, or value is NaN
+            return Double.NaN;
+        }
+    }
 }

libs/core/src/test/java/org/elasticsearch/common/util/ESSloppyMathTests.java

@@ -11,6 +11,7 @@ package org.elasticsearch.common.util;
 import org.elasticsearch.test.ESTestCase;
 
 import static org.elasticsearch.core.ESSloppyMath.atan;
+import static org.elasticsearch.core.ESSloppyMath.log;
 import static org.elasticsearch.core.ESSloppyMath.sinh;
 
 public class ESSloppyMathTests extends ESTestCase {
@@ -19,6 +20,9 @@ public class ESSloppyMathTests extends ESTestCase {
     static double ATAN_DELTA = 1E-15;
     // accuracy for sinh(x)
     static double SINH_DELTA = 1E-15; // for small x
+    // accuracy for log(x)
+    static double LOG_DELTA = 1E-12;
+    static double LOG_DELTA_1 = 1E-14; // near 1.0 we can be more accurate
 
     public void testAtan() {
         assertTrue(Double.isNaN(atan(Double.NaN)));
@@ -42,4 +46,17 @@ public class ESSloppyMathTests extends ESTestCase {
             assertEquals(StrictMath.sinh(d), sinh(d), SINH_DELTA);
         }
     }
+
+    public void testLog() {
+        assertTrue(Double.isNaN(log(Double.NaN)));
+        assertEquals(Double.NaN, log(Double.NEGATIVE_INFINITY), LOG_DELTA);
+        assertEquals(Double.POSITIVE_INFINITY, log(Double.POSITIVE_INFINITY), LOG_DELTA);
+        for (int i = 0; i < 10000; i++) {
+            double d = randomDoubleBetween(-Double.MAX_VALUE, Double.MAX_VALUE, true);
+            assertEquals(StrictMath.log(d), log(d), LOG_DELTA);
+            // Do extra testing around 1.0 due to the special cases near 0, 0.95 and 1.14
+            d = randomDoubleBetween(0, 2, true);
+            assertEquals(StrictMath.log(d), log(d), LOG_DELTA_1);
+        }
+    }
 }

server/src/main/java/org/elasticsearch/search/aggregations/bucket/geogrid/GeoTileUtils.java

@@ -133,7 +133,7 @@ public final class GeoTileUtils {
      */
     public static int getYTile(double latitude, long tiles) {
         double latSin = SloppyMath.cos(PI_DIV_2 - Math.toRadians(normalizeLat(latitude)));
-        int yTile = (int) Math.floor((0.5 - (Math.log((1 + latSin) / (1 - latSin)) / (4 * Math.PI))) * tiles);
+        int yTile = (int) Math.floor((0.5 - (ESSloppyMath.log((1 + latSin) / (1 - latSin)) / (4 * Math.PI))) * tiles);
 
         if (yTile < 0) {
             yTile = 0;