net.sf.geographiclib.GeoMath Maven / Gradle / Ivy
Show all versions of GeographicLib-Java Show documentation
/**
* Implementation of the net.sf.geographiclib.GeoMath class
*
* Copyright (c) Charles Karney (2013-2015) and licensed
* under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/
**********************************************************************/
package net.sf.geographiclib;
/**
* Mathematical functions needed by GeographicLib.
*
* Define mathematical functions and constants so that any version of Java
* can be used.
**********************************************************************/
public class GeoMath {
/**
* The number of binary digits in the fraction of a double precision
* number (equivalent to C++'s {@code numeric_limits::digits}).
**********************************************************************/
public static final int digits = 53;
/**
* Equivalent to C++'s {@code numeric_limits::epsilon()}. In Java
* version 1.5 and later, Math.ulp(1.0) can be used.
**********************************************************************/
public static final double epsilon = Math.pow(0.5, digits - 1);
/**
* Equivalent to C++'s {@code numeric_limits::min()}. In Java
* version 1.6 and later, Double.MIN_NORMAL can be used.
**********************************************************************/
public static final double min = Math.pow(0.5, 1022);
/**
* Square a number.
*
* @param x the argument.
* @return x2.
**********************************************************************/
public static double sq(double x) { return x * x; }
/**
* The hypotenuse function avoiding underflow and overflow. In Java version
* 1.5 and later, Math.hypot can be used.
*
* @param x the first argument.
* @param y the second argument.
* @return sqrt(x2 + y2).
**********************************************************************/
public static double hypot(double x, double y) {
x = Math.abs(x); y = Math.abs(y);
double a = Math.max(x, y), b = Math.min(x, y) / (a != 0 ? a : 1);
return a * Math.sqrt(1 + b * b);
// For an alternative method see
// C. Moler and D. Morrision (1983) https://dx.doi.org/10.1147/rd.276.0577
// and A. A. Dubrulle (1983) https://dx.doi.org/10.1147/rd.276.0582
}
/**
* log(1 + x) accurate near x = 0. In Java version 1.5 and
* later, Math.log1p can be used.
*
* This is taken from D. Goldberg,
* What every computer
* scientist should know about floating-point arithmetic (1991),
* Theorem 4. See also, N. J. Higham, Accuracy and Stability of Numerical
* Algorithms, 2nd Edition (SIAM, 2002), Answer to Problem 1.5, p 528.
*
* @param x the argument.
* @return log(1 + x).
**********************************************************************/
public static double log1p(double x) {
double
y = 1 + x,
z = y - 1;
// Here's the explanation for this magic: y = 1 + z, exactly, and z
// approx x, thus log(y)/z (which is nearly constant near z = 0) returns
// a good approximation to the true log(1 + x)/x. The multiplication x *
// (log(y)/z) introduces little additional error.
return z == 0 ? x : x * Math.log(y) / z;
}
/**
* The inverse hyperbolic tangent function. This is defined in terms of
* GeoMath.log1p(x) in order to maintain accuracy near x = 0.
* In addition, the odd parity of the function is enforced.
*
* @param x the argument.
* @return atanh(x).
**********************************************************************/
public static double atanh(double x) {
double y = Math.abs(x); // Enforce odd parity
y = Math.log1p(2 * y/(1 - y))/2;
return x < 0 ? -y : y;
}
/**
* The cube root function. In Java version 1.5 and later, Math.cbrt can be
* used.
*
* @param x the argument.
* @return the real cube root of x.
**********************************************************************/
public static double cbrt(double x) {
double y = Math.pow(Math.abs(x), 1/3.0); // Return the real cube root
return x < 0 ? -y : y;
}
public static Pair norm(double sinx, double cosx) {
double r = hypot(sinx, cosx);
return new Pair(sinx/r, cosx/r);
}
/**
* The error-free sum of two numbers.
*
* @param u the first number in the sum.
* @param v the second number in the sum.
* @return Pair(s, t) with s = round(u +
* v) and t = u + v - s.
*
* See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.
**********************************************************************/
public static Pair sum(double u, double v) {
double s = u + v;
double up = s - v;
double vpp = s - up;
up -= u;
vpp -= v;
double t = -(up + vpp);
// u + v = s + t
// = round(u + v) + t
return new Pair(s, t);
}
/**
* Evaluate a polynomial.
*
* @param N the order of the polynomial.
* @param p the coefficient array (of size N + s + 1 or more).
* @param s starting index for the array.
* @param x the variable.
* @return the value of the polynomial.
*
* Evaluate y = ∑n=0..N
* ps+n
* xN−n. Return 0 if N < 0.
* Return ps, if N = 0 (even if x is
* infinite or a nan). The evaluation uses Horner's method.
**********************************************************************/
public static double polyval(int N, double p[], int s, double x) {
double y = N < 0 ? 0 : p[s++];
while (--N >= 0) y = y * x + p[s++];
return y;
}
public static double AngRound(double x) {
// The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
// for reals = 0.7 pm on the earth if x is an angle in degrees. (This
// is about 1000 times more resolution than we get with angles around 90
// degrees.) We use this to avoid having to deal with near singular
// cases when x is non-zero but tiny (e.g., 1.0e-200). This also converts
// -0 to +0.
final double z = 1/16.0;
double y = Math.abs(x);
// The compiler mustn't "simplify" z - (z - y) to y
y = y < z ? z - (z - y) : y;
return x < 0 ? 0 - y : y;
}
/**
* Normalize an angle (restricted input range).
*
* @param x the angle in degrees.
* @return the angle reduced to the range [−180°, 180°).
*
* The range of x is unrestricted.
**********************************************************************/
public static double AngNormalize(double x) {
x = x % 360.0;
return x < -180 ? x + 360 : (x < 180 ? x : x - 360);
}
/**
* Normalize a latitude.
*
* @param x the angle in degrees.
* @return x if it is in the range [−90°, 90°], otherwise
* return NaN.
**********************************************************************/
public static double LatFix(double x) {
return Math.abs(x) > 90 ? Double.NaN : x;
}
/**
* Difference of two angles reduced to [−180°, 180°]
*
* @param x the first angle in degrees.
* @param y the second angle in degrees.
* @return y − x, reduced to the range [−180°,
* 180°].
*
* x and y must both lie in [−180°, 180°]. The
* result is equivalent to computing the difference exactly, reducing it to
* (−180°, 180°] and rounding the result. Note that this
* prescription allows −180° to be returned (e.g., if x is
* tiny and negative and y = 180°).
**********************************************************************/
public static double AngDiff(double x, double y) {
double d, t;
{
Pair r = sum(AngNormalize(x), AngNormalize(-y));
d = - AngNormalize(r.first); t = r.second;
}
return (d == 180 && t < 0 ? -180 : d) - t;
}
/**
* Evaluate the sine and cosine function with the argument in degrees
*
* @param x in degrees.
* @return Pair(s, t) with s = sin(x) and
* c = cos(x).
*
* The results obey exactly the elementary properties of the trigonometric
* functions, e.g., sin 9° = cos 81° = − sin 123456789°.
**********************************************************************/
public static Pair sincosd(double x) {
// In order to minimize round-off errors, this function exactly reduces
// the argument to the range [-45, 45] before converting it to radians.
double r; int q;
r = x % 360.0;
q = (int)Math.floor(r / 90 + 0.5);
r -= 90 * q;
// now abs(r) <= 45
r = Math.toRadians(r);
// Possibly could call the gnu extension sincos
double s = Math.sin(r), c = Math.cos(r);
double sinx, cosx;
switch (q & 3) {
case 0: sinx = s; cosx = c; break;
case 1: sinx = c; cosx = 0 - s; break;
case 2: sinx = 0 - s; cosx = 0 - c; break;
default: sinx = 0 - c; cosx = s; break; // case 3
}
return new Pair(sinx, cosx);
}
/**
* Evaluate the atan2 function with the result in degrees
*
* @param y
* @param x
* @return atan2(y, x) in degrees.
*
* The result is in the range [−180° 180°). N.B.,
* atan2d(±0, −1) = −180°; atan2d(+ε,
* −1) = +180°, for ε positive and tiny;
* atan2d(±0, 1) = ±0°.
**********************************************************************/
public static double atan2d(double y, double x) {
// In order to minimize round-off errors, this function rearranges the
// arguments so that result of atan2 is in the range [-pi/4, pi/4] before
// converting it to degrees and mapping the result to the correct
// quadrant.
int q = 0;
if (Math.abs(y) > Math.abs(x)) { double t; t = x; x = y; y = t; q = 2; }
if (x < 0) { x = -x; ++q; }
// here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
double ang = Math.toDegrees(Math.atan2(y, x));
switch (q) {
// Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
// atan2d will not be called with y = -0. If need be, include
//
// case 0: ang = 0 + ang; break;
//
// and handle mpfr as in AngRound.
case 1: ang = (y > 0 ? 180 : -180) - ang; break;
case 2: ang = 90 - ang; break;
case 3: ang = -90 + ang; break;
}
return ang;
}
/**
* Test for finiteness.
*
* @param x the argument.
* @return true if number is finite, false if NaN or infinite.
**********************************************************************/
public static boolean isfinite(double x) {
return Math.abs(x) <= Double.MAX_VALUE;
}
private GeoMath() {}
}