Quick benchmark for different ways to compute distance between points on the globe.

Distance.java

package com.mapr.bench;

import com.google.caliper.Benchmark;
import com.google.caliper.runner.Running;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;

public class Distance {
    @Benchmark
    double timeHaversinFast(int reps) {
        double sum = 0;
        for (int i = 0; i < reps; i++) {
            double x = 0.2 * i / reps;
            sum += haversinFast(32.1, 0.0, 32.1, x);
        }
        return sum;
    }

    @Benchmark
    double timeHaversin(int reps) {
        double sum = 0;
        for (int i = 0; i < reps; i++) {
            double x = 0.2 * i / reps;
            sum += haversin(32.1, 0, 32.1, x);
        }
        return sum;
    }

    @Benchmark
    double timeDot(int reps) {
        double sum = 0;
        for (int i = 0; i < reps; i++) {
            double x = 0.2 * i / reps;
            sum += dist(32.1, 0, 32.1, x);
        }
        return sum;
    }


    public void testTrig() {
        for (int i = 0; i < 10000; i++) {
            double x = 0.05 * i / 10000;
            Assert.assertTrue(Math.abs(polyAcos(x) - Math.acos(x)) < 1e-7);
        }
        long t0 = System.nanoTime();
        double sum1 = 0;
        for (int i = 0; i < 10000; i++) {
            double x = 0.2 * i / 10000;
            sum1 += dist(32.1, 0, 32.1, x);
        }
        long t1 = System.nanoTime();
        long t2 = System.nanoTime();
//        System.out.printf("%.5f\t%.5f\n%.1f\t%.1f\n", sum1, sum2, (t1 - t0) / 10000.0, (t2 - t1) / 10000.0);

    }

    double dist(double lat1, double long1, double lat2, double long2) {
        lat1 *= TO_RADIANS;
        long1 *= TO_RADIANS;
        lat2 *= TO_RADIANS;
        long2 *= TO_RADIANS;
        //        double x1 = 0;
        double y1 = FastMath.cos(lat1);
        double z1 = FastMath.sin(lat1);
        double q2 = FastMath.cos(lat2);

//        double x2 = FastMath.sin(long1 - long2) * q2;
        double y2 = FastMath.cos(long1 - long2) * q2;
        double z2 = FastMath.sin(lat2);

        return polyAcos(FastMath.sqrt(y1 * y2 + z1 * z2)) * TO_KILOMETERS;
    }

    double polyAsin(double x) {
        double z = x * x;
        return x *
                (z * (z * (z * (z * (z * (z * (1288287 * z +
                        1600830) + 2063880) + 2802800) + 4118400) + 6918912) + 15375360) + 92252160) / 92252160;
    }

    double polyAcos(double x) {
        double z = x * x;
//        return (x * (z * (z * (z * (z * (z * (-800415 * z - 1031940) - 1401400) - 2059200) - 3459456) - 7687680) - 46126080) + 23063040 * Math.PI / 46126080);

        return (x * (z * (z * (z * (z * (z * (z * (-1288287.0 * z -
                1600830.0) - 2063880.0) - 2802800.0) - 4118400.0) - 6918912.0) - 15375360.0) - 92252160.0) + 46126080.0 * Math.PI) / 92252160.0;
    }

    /**
     * Returns the distance between two points in decimal degrees.
     *
     * @param lat1 Latitude of the first point.
     * @param lon1 Longitude of the first point.
     * @param lat2 Latitude of the second point.
     * @param lon2 Longitude of the second point.
     * @return distance in kilometers.
     */
    public static double haversin(double lat1, double lon1, double lat2, double lon2) {
        double x1 = lat1 * TO_RADIANS;
        double x2 = lat2 * TO_RADIANS;
        double h1 = (1 - Math.cos(x1 - x2)) / 2;
        double h2 = (1 - Math.cos((lon1 - lon2) * TO_RADIANS)) / 2;
        double h = h1 + Math.cos(x1) * Math.cos(x2) * h2;
        return TO_KILOMETERS * 2 * Math.asin(Math.min(1, Math.sqrt(h)));
    }

    /**
     * Returns the distance between two points in decimal degrees.
     *
     * @param lat1 Latitude of the first point.
     * @param lon1 Longitude of the first point.
     * @param lat2 Latitude of the second point.
     * @param lon2 Longitude of the second point.
     * @return distance in kilometers.
     */
    public static double haversinFast(double lat1, double lon1, double lat2, double lon2) {
        double x1 = lat1 * TO_RADIANS;
        double x2 = lat2 * TO_RADIANS;
        double h1 = (1 - cos(x1 - x2)) / 2;
        double h2 = (1 - cos((lon1 - lon2) * TO_RADIANS)) / 2;
        double h = h1 + cos(x1) * cos(x2) * h2;
        return TO_KILOMETERS * 2 * asin(Math.min(1, Math.sqrt(h)));
    }

    /**
     * Returns the trigonometric cosine of an angle.
     * <p/>
     * Error is around 1E-15.
     * <p/>
     * Special cases:
     * <ul>
     * <li>If the argument is {@code NaN} or an infinity, then the result is {@code NaN}.
     * </ul>
     *
     * @param a an angle, in radians.
     * @return the cosine of the argument.
     * @see Math#cos(double)
     */
    public static double cos(double a) {
        if (a < 0.0) {
            a = -a;
        }
        if (a > SIN_COS_MAX_VALUE_FOR_INT_MODULO) {
            return Math.cos(a);
        }
        // index: possibly outside tables range.
        int index = (int) (a * SIN_COS_INDEXER + 0.5);
        double delta = (a - index * SIN_COS_DELTA_HI) - index * SIN_COS_DELTA_LO;
        // Making sure index is within tables range.
        // Last value of each table is the same than first, so we ignore it (tabs size minus one) for modulo.
        index &= (SIN_COS_TABS_SIZE - 2); // index % (SIN_COS_TABS_SIZE-1)
        double indexCos = cosTab[index];
        double indexSin = sinTab[index];
        return indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
    }

    /**
     * Returns the arc sine of a value.
     * <p/>
     * The returned angle is in the range <i>-pi</i>/2 through <i>pi</i>/2.
     * Error is around 1E-7.
     * <p/>
     * Special cases:
     * <ul>
     * <li>If the argument is {@code NaN} or its absolute value is greater than 1, then the result is {@code NaN}.
     * </ul>
     *
     * @param a the value whose arc sine is to be returned.
     * @return arc sine of the argument
     * @see Math#asin(double)
     */
// because asin(-x) = -asin(x), asin(x) only needs to be computed on [0,1].
// ---> we only have to compute asin(x) on [0,1].
// For values not close to +-1, we use look-up tables;
// for values near +-1, we use code derived from fdlibm.
    public static double asin(double a) {
        boolean negateResult;
        if (a < 0.0) {
            a = -a;
            negateResult = true;
        } else {
            negateResult = false;
        }
        if (a <= ASIN_MAX_VALUE_FOR_TABS) {
            int index = (int) (a * ASIN_INDEXER + 0.5);
            double delta = a - index * ASIN_DELTA;
            double result = asinTab[index] + delta * (asinDer1DivF1Tab[index] + delta * (asinDer2DivF2Tab[index] + delta * (asinDer3DivF3Tab[index] + delta * asinDer4DivF4Tab[index])));
            return negateResult ? -result : result;
        } else { // value > ASIN_MAX_VALUE_FOR_TABS, or value is NaN
            // This part is derived from fdlibm.
            if (a < 1.0) {
                double t = (1.0 - a) * 0.5;
                double p = t * (ASIN_PS0 + t * (ASIN_PS1 + t * (ASIN_PS2 + t * (ASIN_PS3 + t * (ASIN_PS4 + t * ASIN_PS5)))));
                double q = 1.0 + t * (ASIN_QS1 + t * (ASIN_QS2 + t * (ASIN_QS3 + t * ASIN_QS4)));
                double s = Math.sqrt(t);
                double z = s + s * (p / q);
                double result = ASIN_PIO2_HI - ((z + z) - ASIN_PIO2_LO);
                return negateResult ? -result : result;
            } else { // value >= 1.0, or value is NaN
                if (a == 1.0) {
                    return negateResult ? -Math.PI / 2 : Math.PI / 2;
                } else {
                    return Double.NaN;
                }
            }
        }
    }

    // haversin
    private static final double TO_RADIANS = Math.PI / 180D;
    private static final double TO_KILOMETERS = 6371.0087714D;

    // cos/asin
    private static final double ONE_DIV_F2 = 1 / 2.0;
    private static final double ONE_DIV_F3 = 1 / 6.0;
    private static final double ONE_DIV_F4 = 1 / 24.0;

    private static final double PIO2_HI = Double.longBitsToDouble(0x3FF921FB54400000L); // 1.57079632673412561417e+00 first 33 bits of pi/2
    private static final double PIO2_LO = Double.longBitsToDouble(0x3DD0B4611A626331L); // 6.07710050650619224932e-11 pi/2 - PIO2_HI
    private static final double TWOPI_HI = 4 * PIO2_HI;
    private static final double TWOPI_LO = 4 * PIO2_LO;
    private static final int SIN_COS_TABS_SIZE = (1 << 11) + 1;
    private static final double SIN_COS_DELTA_HI = TWOPI_HI / (SIN_COS_TABS_SIZE - 1);
    private static final double SIN_COS_DELTA_LO = TWOPI_LO / (SIN_COS_TABS_SIZE - 1);
    private static final double SIN_COS_INDEXER = 1 / (SIN_COS_DELTA_HI + SIN_COS_DELTA_LO);
    private static final double[] sinTab = new double[SIN_COS_TABS_SIZE];
    private static final double[] cosTab = new double[SIN_COS_TABS_SIZE];

    // Max abs value for fast modulo, above which we use regular angle normalization.
// This value must be < (Integer.MAX_VALUE / SIN_COS_INDEXER), to stay in range of int type.
// The higher it is, the higher the error, but also the faster it is for lower values.
// If you set it to ((Integer.MAX_VALUE / SIN_COS_INDEXER) * 0.99), worse accuracy on double range is about 1e-10.
    static final double SIN_COS_MAX_VALUE_FOR_INT_MODULO = ((Integer.MAX_VALUE >> 9) / SIN_COS_INDEXER) * 0.99;

    // Supposed to be >= sin(77.2deg), as fdlibm code is supposed to work with values > 0.975,
// but seems to work well enough as long as value >= sin(25deg).
    private static final double ASIN_MAX_VALUE_FOR_TABS = StrictMath.sin(Math.toRadians(73.0));

    private static final int ASIN_TABS_SIZE = (1 << 13) + 1;
    private static final double ASIN_DELTA = ASIN_MAX_VALUE_FOR_TABS / (ASIN_TABS_SIZE - 1);
    private static final double ASIN_INDEXER = 1 / ASIN_DELTA;
    private static final double[] asinTab = new double[ASIN_TABS_SIZE];
    private static final double[] asinDer1DivF1Tab = new double[ASIN_TABS_SIZE];
    private static final double[] asinDer2DivF2Tab = new double[ASIN_TABS_SIZE];
    private static final double[] asinDer3DivF3Tab = new double[ASIN_TABS_SIZE];
    private static final double[] asinDer4DivF4Tab = new double[ASIN_TABS_SIZE];

    private static final double ASIN_PIO2_HI = Double.longBitsToDouble(0x3FF921FB54442D18L); // 1.57079632679489655800e+00
    private static final double ASIN_PIO2_LO = Double.longBitsToDouble(0x3C91A62633145C07L); // 6.12323399573676603587e-17
    private static final double ASIN_PS0 = Double.longBitsToDouble(0x3fc5555555555555L); //  1.66666666666666657415e-01
    private static final double ASIN_PS1 = Double.longBitsToDouble(0xbfd4d61203eb6f7dL); // -3.25565818622400915405e-01
    private static final double ASIN_PS2 = Double.longBitsToDouble(0x3fc9c1550e884455L); //  2.01212532134862925881e-01
    private static final double ASIN_PS3 = Double.longBitsToDouble(0xbfa48228b5688f3bL); // -4.00555345006794114027e-02
    private static final double ASIN_PS4 = Double.longBitsToDouble(0x3f49efe07501b288L); //  7.91534994289814532176e-04
    private static final double ASIN_PS5 = Double.longBitsToDouble(0x3f023de10dfdf709L); //  3.47933107596021167570e-05
    private static final double ASIN_QS1 = Double.longBitsToDouble(0xc0033a271c8a2d4bL); // -2.40339491173441421878e+00
    private static final double ASIN_QS2 = Double.longBitsToDouble(0x40002ae59c598ac8L); //  2.02094576023350569471e+00
    private static final double ASIN_QS3 = Double.longBitsToDouble(0xbfe6066c1b8d0159L); // -6.88283971605453293030e-01
    private static final double ASIN_QS4 = Double.longBitsToDouble(0x3fb3b8c5b12e9282L); //  7.70381505559019352791e-02

    /** Initializes look-up tables. */
    static {
        // sin and cos
        final int SIN_COS_PI_INDEX = (SIN_COS_TABS_SIZE - 1) / 2;
        final int SIN_COS_PI_MUL_2_INDEX = 2 * SIN_COS_PI_INDEX;
        final int SIN_COS_PI_MUL_0_5_INDEX = SIN_COS_PI_INDEX / 2;
        final int SIN_COS_PI_MUL_1_5_INDEX = 3 * SIN_COS_PI_INDEX / 2;
        for (int i = 0; i < SIN_COS_TABS_SIZE; i++) {
            // angle: in [0,2*PI].
            double angle = i * SIN_COS_DELTA_HI + i * SIN_COS_DELTA_LO;
            double sinAngle = StrictMath.sin(angle);
            double cosAngle = StrictMath.cos(angle);
            // For indexes corresponding to null cosine or sine, we make sure the value is zero
            // and not an epsilon. This allows for a much better accuracy for results close to zero.
            if (i == SIN_COS_PI_INDEX) {
                sinAngle = 0.0;
            } else if (i == SIN_COS_PI_MUL_2_INDEX) {
                sinAngle = 0.0;
            } else if (i == SIN_COS_PI_MUL_0_5_INDEX) {
                cosAngle = 0.0;
            } else if (i == SIN_COS_PI_MUL_1_5_INDEX) {
                cosAngle = 0.0;
            }
            sinTab[i] = sinAngle;
            cosTab[i] = cosAngle;
        }

        // asin
        for (int i = 0; i < ASIN_TABS_SIZE; i++) {
            // x: in [0,ASIN_MAX_VALUE_FOR_TABS].
            double x = i * ASIN_DELTA;
            asinTab[i] = StrictMath.asin(x);
            double oneMinusXSqInv = 1.0 / (1 - x * x);
            double oneMinusXSqInv0_5 = StrictMath.sqrt(oneMinusXSqInv);
            double oneMinusXSqInv1_5 = oneMinusXSqInv0_5 * oneMinusXSqInv;
            double oneMinusXSqInv2_5 = oneMinusXSqInv1_5 * oneMinusXSqInv;
            double oneMinusXSqInv3_5 = oneMinusXSqInv2_5 * oneMinusXSqInv;
            asinDer1DivF1Tab[i] = oneMinusXSqInv0_5;
            asinDer2DivF2Tab[i] = (x * oneMinusXSqInv1_5) * ONE_DIV_F2;
            asinDer3DivF3Tab[i] = ((1 + 2 * x * x) * oneMinusXSqInv2_5) * ONE_DIV_F3;
            asinDer4DivF4Tab[i] = ((5 + 2 * x * (2 + x * (5 - 2 * x))) * oneMinusXSqInv3_5) * ONE_DIV_F4;
        }
    }

}

Awesome! I was already thinking of doing a Caliper benchmark after seeing the claims made on the JIRA.