| /* |
| * 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. */ |
| 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. |
| * <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 {@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. |
| * |
| * @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 |
| * {@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. |
| * <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; |
| } |
| } |
| } |
| } |
| |
| /** |
| * Convert to degrees. |
| * @param radians radians to convert to degrees |
| * @return degrees |
| * @deprecated Use {@link Math#toDegrees(double)} on Java 9+. |
| */ |
| @Deprecated |
| public static double toDegrees(final double radians) { |
| return radians * TO_DEGREES; |
| } |
| |
| /** |
| * Convert to radians. |
| * @param degrees degrees to convert to radians |
| * @return radians |
| * @deprecated Use {@link Math#toRadians(double)} on Java 9+. |
| */ |
| @Deprecated |
| 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. */ |
| 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; |
| } |
| } |
| } |