/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.lucene.util;


/* some code derived from jodk: http://code.google.com/p/jodk/ (apache 2.0)
 * asin() derived from fdlibm: http://www.netlib.org/fdlibm/e_asin.c (public domain):
 * =============================================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * =============================================================================
 */

Math functions that trade off accuracy for speed. /** Math functions that trade off accuracy for speed. */
public class SloppyMath {
  
  
Returns the Haversine distance in meters between two points
specified in decimal degrees (latitude/longitude).  This works correctly
even if the dateline is between the two points.

Error is at most 4E-1 (40cm) from the actual haversine distance, but is typically
much smaller for reasonable distances: around 1E-5 (0.01mm) for distances less than
1000km.
Params:
lat1 – Latitude of the first point.
lon1 – Longitude of the first point.
lat2 – Latitude of the second point.
lon2 – Longitude of the second point.
Returns:
distance in meters./**
   * Returns the Haversine distance in meters between two points
   * specified in decimal degrees (latitude/longitude).  This works correctly
   * even if the dateline is between the two points.
   * <p>
   * Error is at most 4E-1 (40cm) from the actual haversine distance, but is typically
   * much smaller for reasonable distances: around 1E-5 (0.01mm) for distances less than
   * 1000km.
   *
   * @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 meters.
   */
  public static double haversinMeters(double lat1, double lon1, double lat2, double lon2) {
    return haversinMeters(haversinSortKey(lat1, lon1, lat2, lon2));
  }
  
  
Returns the Haversine distance in meters between two points given the previous result from haversinSortKey(double, double, double, double) 
Returns:
distance in meters./**
   * Returns the Haversine distance in meters between two points
   * given the previous result from {@link #haversinSortKey(double, double, double, double)}
   * @return distance in meters.
   */
  public static double haversinMeters(double sortKey) {
    return TO_METERS * 2 * asin(Math.min(1, Math.sqrt(sortKey * 0.5)));
  }
  
  
Returns the Haversine distance in kilometers between two points
specified in decimal degrees (latitude/longitude).  This works correctly
even if the dateline is between the two points.
Params:
lat1 – Latitude of the first point.
lon1 – Longitude of the first point.
lat2 – Latitude of the second point.
lon2 – Longitude of the second point.
Returns:
distance in kilometers.
Deprecated:
Use instead/**
   * Returns the Haversine distance in kilometers between two points
   * specified in decimal degrees (latitude/longitude).  This works correctly
   * even if the dateline is between the two points.
   *
   * @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.
   * @deprecated Use {@link #haversinMeters(double, double, double, double) instead}
   */
  @Deprecated
  public static double haversinKilometers(double lat1, double lon1, double lat2, double lon2) {
    double h = haversinSortKey(lat1, lon1, lat2, lon2);
    return TO_KILOMETERS * 2 * asin(Math.min(1, Math.sqrt(h * 0.5)));
  }
  
  
Returns a sort key for distance. This is less expensive to compute than haversinMeters(double, double, double, double), but it always compares the same. This can be converted into an actual distance with haversinMeters(double), which effectively does the second half of the computation. /**
   * Returns a sort key for distance. This is less expensive to compute than 
   * {@link #haversinMeters(double, double, double, double)}, but it always compares the same.
   * This can be converted into an actual distance with {@link #haversinMeters(double)}, which
   * effectively does the second half of the computation.
   */
  public static double haversinSortKey(double lat1, double lon1, double lat2, double lon2) {
    double x1 = lat1 * TO_RADIANS;
    double x2 = lat2 * TO_RADIANS;
    double h1 = 1 - cos(x1 - x2);
    double h2 = 1 - cos((lon1 - lon2) * TO_RADIANS);
    double h = h1 + cos(x1) * cos(x2) * h2;
    // clobber crazy precision so subsequent rounding does not create ties.
    return Double.longBitsToDouble(Double.doubleToRawLongBits(h) & 0xFFFFFFFFFFFFFFF8L);
  }

  
Returns the trigonometric cosine of an angle.

Error is around 1E-15.

Special cases:

 
If the argument is NaN or an infinity, then the result is NaN. 
Params:
a – an angle, in radians.
See Also:
Math.cos(double)
Returns:
the cosine of the argument./**
   * 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.

The returned angle is in the range -pi/2 through pi/2. 
Error is around 1E-7.

Special cases:

 
If the argument is NaN or its absolute value is greater than 1, then the result is NaN. 
Params:
a – the value whose arc sine is to be returned.
See Also:
Math.asin(double)
Returns:
arc sine of the argument/**
   * 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;
        }
      }
    }
  }

  
Convert to degrees.
Params:
radians – radians to convert to degrees
Returns:
degrees/**
   * Convert to degrees.
   * @param radians radians to convert to degrees
   * @return degrees
   */
  public static double toDegrees(final double radians) {
    return radians * TO_DEGREES;
  }
  
  
Convert to radians.
Params:
degrees – degrees to convert to radians
Returns:
radians/**
   * Convert to radians.
   * @param degrees degrees to convert to radians
   * @return radians
   */
  public static double toRadians(final double degrees) {
    return degrees * TO_RADIANS;
  }
  
  // haversin
  // TODO: remove these for java 9, they fixed Math.toDegrees()/toRadians() to work just like this.
  public static final double TO_RADIANS = Math.PI / 180D;
  public static final double TO_DEGREES = 180D / Math.PI;

  // Earth's mean radius, in meters and kilometers; see http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
  private static final double TO_METERS = 6_371_008.7714D; // equatorial radius
  private static final double TO_KILOMETERS = 6_371.0087714D; // equatorial radius
  
  // 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(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. /** 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;
    }
  }
}