GeoMath.java

Go to the documentation of this file.

package net.sf.geographiclib;

public class GeoMath {
  public static final int digits = 53;
  public static final double epsilon = Math.pow(0.5, digits - 1);
  public static final double min = Math.pow(0.5, 1022);

  public static double sq(double x) { return x * x; }

  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://doi.org/10.1147/rd.276.0577
    // and A. A. Dubrulle (1983) https://doi.org/10.1147/rd.276.0582
  }

  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;
  }

  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;
  }

  public static double copysign(double x, double y) {
    return Math.abs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
  }

  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);
  }

  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);
  }

  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 converts -0 to
    // +0; however tiny negative numbers get converted to -0.
    final double z = 1/16.0;
    if (x == 0) return 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 ? -y : y;
  }

  public static double AngNormalize(double x) {
    x = x % 360.0;
    return x <= -180 ? x + 360 : (x <= 180 ? x : x - 360);
  }

  public static double LatFix(double x) {
    return Math.abs(x) > 90 ? Double.NaN : x;
  }

  public static Pair AngDiff(double x, double y) {
    double d, t;
    {
      Pair r = sum(AngNormalize(-x), AngNormalize(y));
      d = AngNormalize(r.first); t = r.second;
    }
    return sum(d == 180 && t > 0 ? -180 : d, t);
  }

  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 = -s; break;
    case  2: sinx = -s; cosx = -c; break;
    default: sinx = -c; cosx =  s; break; // case 3
    }
    if (x != 0) { sinx += 0.0; cosx += 0.0; }
    return new Pair(sinx, cosx);
  }

  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;
  }

  public static boolean isfinite(double x) {
    return Math.abs(x) <= Double.MAX_VALUE;
  }

  private GeoMath() {}
}