Compute the hyperbolic tangent of a number.
@param x number on which evaluation is done
@return hyperbolic tangent of x
2665
public static double tanh(double x) {
      boolean negate = false;

      if (x != x) {
          return x;
      }

      // tanh[z] = sinh[z] / cosh[z]
      // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
      // = (exp(2x) - 1) / (exp(2x) + 1)

      // for magnitude > 20, sinh[z] == cosh[z] in double precision

      if (x > 20.0) {
          return 1.0;
      }

      if (x < -20) {
          return -1.0;
      }

      if (x == 0) {
          return x;
      }

      if (x < 0.0) {
          x = -x;
          negate = true;
      }

      double result;
      if (x >= 0.5) {
          double hiPrec[] = new double[2];
          // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
          exp(x*2.0, 0.0, hiPrec);

          double ya = hiPrec[0] + hiPrec[1];
          double yb = -(ya - hiPrec[0] - hiPrec[1]);

          /* Numerator */
          double na = -1.0 + ya;
          double nb = -(na + 1.0 - ya);
          double temp = na + yb;
          nb += -(temp - na - yb);
          na = temp;

          /* Denominator */
          double da = 1.0 + ya;
          double db = -(da - 1.0 - ya);
          temp = da + yb;
          db += -(temp - da - yb);
          da = temp;

          temp = da * HEX_40000000;
          double daa = da + temp - temp;
          double dab = da - daa;

          // ratio = na/da
          double ratio = na/da;
          temp = ratio * HEX_40000000;
          double ratioa = ratio + temp - temp;
          double ratiob = ratio - ratioa;

          // Correct for rounding in division
          ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;

          // Account for nb
          ratiob += nb / da;
          // Account for db
          ratiob += -db * na / da / da;

          result = ratioa + ratiob;
      }
      else {
          double hiPrec[] = new double[2];
          // tanh(x) = expm1(2x) / (expm1(2x) + 2)
          expm1(x*2.0, hiPrec);

          double ya = hiPrec[0] + hiPrec[1];
          double yb = -(ya - hiPrec[0] - hiPrec[1]);

          /* Numerator */
          double na = ya;
          double nb = yb;

          /* Denominator */
          double da = 2.0 + ya;
          double db = -(da - 2.0 - ya);
          double temp = da + yb;
          db += -(temp - da - yb);
          da = temp;

          temp = da * HEX_40000000;
          double daa = da + temp - temp;
          double dab = da - daa;

          // ratio = na/da
          double ratio = na/da;
          temp = ratio * HEX_40000000;
          double ratioa = ratio + temp - temp;
          double ratiob = ratio - ratioa;

          // Correct for rounding in division
          ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;

          // Account for nb
          ratiob += nb / da;
          // Account for db
          ratiob += -db * na / da / da;

          result = ratioa + ratiob;
      }

      if (negate) {
          result = -result;
      }

      return result;
    }
684475
public static double cbrt(double x) {
      /* Convert input double to bits */
      long inbits = Double.doubleToLongBits(x);
      int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
      boolean subnormal = false;

      if (exponent == -1023) {
          if (x == 0) {
              return x;
          }

          /* Subnormal, so normalize */
          subnormal = true;
          x *= 1.8014398509481984E16;  // 2^54
          inbits = Double.doubleToLongBits(x);
          exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
      }

      if (exponent == 1024) {
          // Nan or infinity.  Don't care which.
          return x;
      }

      /* Divide the exponent by 3 */
      int exp3 = exponent / 3;

      /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
      double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
                                          (long)(((exp3 + 1023) & 0x7ff)) << 52);

      /* This will be a number between 1 and 2 */
      final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);

      /* Estimate the cube root of mant by polynomial */
      double est = -0.010714690733195933;
      est = est * mant + 0.0875862700108075;
      est = est * mant + -0.3058015757857271;
      est = est * mant + 0.7249995199969751;
      est = est * mant + 0.5039018405998233;

      est *= CBRTTWO[exponent % 3 + 2];

      // est should now be good to about 15 bits of precision.   Do 2 rounds of
      // Newton's method to get closer,  this should get us full double precision
      // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
      final double xs = x / (p2*p2*p2);
      est += (xs - est*est*est) / (3*est*est);
      est += (xs - est*est*est) / (3*est*est);

      // Do one round of Newton's method in extended precision to get the last bit right.
      double temp = est * HEX_40000000;
      double ya = est + temp - temp;
      double yb = est - ya;

      double za = ya * ya;
      double zb = ya * yb * 2.0 + yb * yb;
      temp = za * HEX_40000000;
      double temp2 = za + temp - temp;
      zb += za - temp2;
      za = temp2;

      zb = za * yb + ya * zb + zb * yb;
      za = za * ya;

      double na = xs - za;
      double nb = -(na - xs + za);
      nb -= zb;

      est += (na+nb)/(3*est*est);

      /* Scale by a power of two, so this is exact. */
      est *= p2;

      if (subnormal) {
          est *= 3.814697265625E-6;  // 2^-18
      }

      return est;
    }
925304
public static double exp(double x) {
    double val = Math.min(MAX, Math.exp(x));
//    if (val == MAX) Log.warn("Exp overflow: exp(" + x + ") truncated to " + MAX);
    return val;
  }
1339497
public static double Y0(double x) {
        if (x < 8.0) {
            double y = x * x;

            double ans1 = -2957821389.0 + y * (7062834065.0 + y * (-512359803.6
                    + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))));
            double ans2 = 40076544269.0 + y * (745249964.8 + y * (7189466.438
                    + y * (47447.26470 + y * (226.1030244 + y * 1.0))));

            return (ans1 / ans2) + 0.636619772 * J0(x) * Math.log(x);
        } else {
            double z = 8.0 / x;
            double y = z * z;
            double xx = x - 0.785398164;

            double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
            double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                    + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                    + y * (-0.934945152e-7))));
            return Math.sqrt(0.636619772 / x) *
                    (Math.sin(xx) * ans1 + z * Math.cos(xx) * ans2);
        }
    }
778092
public static double atan2(double y, double x) {
        if (x != x || y != y) {
            return Double.NaN;
        }

        if (y == 0) {
            final double result = x * y;
            final double invx = 1d / x;
            final double invy = 1d / y;

            if (invx == 0) { // X is infinite
                if (x > 0) {
                    return y; // return +/- 0.0
                } else {
                    return copySign(Math.PI, y);
                }
            }

            if (x < 0 || invx < 0) {
                if (y < 0 || invy < 0) {
                    return -Math.PI;
                } else {
                    return Math.PI;
                }
            } else {
                return result;
            }
        }

        // y cannot now be zero

        if (y == Double.POSITIVE_INFINITY) {
            if (x == Double.POSITIVE_INFINITY) {
                return Math.PI * F_1_4;
            }

            if (x == Double.NEGATIVE_INFINITY) {
                return Math.PI * F_3_4;
            }

            return Math.PI * F_1_2;
        }

        if (y == Double.NEGATIVE_INFINITY) {
            if (x == Double.POSITIVE_INFINITY) {
                return -Math.PI * F_1_4;
            }

            if (x == Double.NEGATIVE_INFINITY) {
                return -Math.PI * F_3_4;
            }

            return -Math.PI * F_1_2;
        }

        if (x == Double.POSITIVE_INFINITY) {
            if (y > 0 || 1 / y > 0) {
                return 0d;
            }

            if (y < 0 || 1 / y < 0) {
                return -0d;
            }
        }

        if (x == Double.NEGATIVE_INFINITY)
        {
            if (y > 0.0 || 1 / y > 0.0) {
                return Math.PI;
            }

            if (y < 0 || 1 / y < 0) {
                return -Math.PI;
            }
        }

        // Neither y nor x can be infinite or NAN here

        if (x == 0) {
            if (y > 0 || 1 / y > 0) {
                return Math.PI * F_1_2;
            }

            if (y < 0 || 1 / y < 0) {
                return -Math.PI * F_1_2;
            }
        }

        // Compute ratio r = y/x
        final double r = y / x;
        if (Double.isInfinite(r)) { // bypass calculations that can create NaN
            return atan(r, 0, x < 0);
        }

        double ra = doubleHighPart(r);
        double rb = r - ra;

        // Split x
        final double xa = doubleHighPart(x);
        final double xb = x - xa;

        rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;

        final double temp = ra + rb;
        rb = -(temp - ra - rb);
        ra = temp;

        if (ra == 0) { // Fix up the sign so atan works correctly
            ra = copySign(0d, y);
        }

        // Call atan
        final double result = atan(ra, rb, x < 0);

        return result;
    }
272044
@Override 
    public double fromReal(double x) {
        long sign = POSITIVE;
        if (x < 0) {
            sign = NEGATIVE;
            x = -x;  
        }
        return compact(sign, FastMath.log(x));
    }
13174
public double get(int x, int y) {
    return (x == y) ? 0 : (x < y) //
        ? matrix[MatrixParadigm.triangleSize(y) + x] //
        : matrix[MatrixParadigm.triangleSize(x) + y];
  }
598560
public static double digamma(double x) {
        final double C7[][] = {
            {
                1.3524999667726346383e4, 4.5285601699547289655e4,
                4.5135168469736662555e4, 1.8529011818582610168e4,
                3.3291525149406935532e3, 2.4068032474357201831e2,
                5.1577892000139084710, 6.2283506918984745826e-3
            },
            {
                6.9389111753763444376e-7, 1.9768574263046736421e4,
                4.1255160835353832333e4, 2.9390287119932681918e4,
                9.0819666074855170271e3, 1.2447477785670856039e3,
                6.7429129516378593773e1, 1.0
            }
        };

        final double C4[][] = {
            {
                -2.728175751315296783e-15, -6.481571237661965099e-1,
                -4.486165439180193579, -7.016772277667586642,
                -2.129404451310105168
            },
            {
                7.777885485229616042, 5.461177381032150702e1,
                8.929207004818613702e1, 3.227034937911433614e1,
                1.0
            }
        };

        double prodPj = 0.0;
        double prodQj = 0.0;
        double digX = 0.0;

        if (x >= 3.0) {
            double x2 = 1.0 / (x * x);
            for (int j = 4; j >= 0; j--) {
                prodPj = prodPj * x2 + C4[0][j];
                prodQj = prodQj * x2 + C4[1][j];
            }
            digX = Math.log(x) - (0.5 / x) + (prodPj / prodQj);

        } else if (x >= 0.5) {
            final double X0 = 1.46163214496836234126;
            for (int j = 7; j >= 0; j--) {
                prodPj = x * prodPj + C7[0][j];
                prodQj = x * prodQj + C7[1][j];
            }
            digX = (x - X0) * (prodPj / prodQj);

        } else {
            double f = (1.0 - x) - Math.floor(1.0 - x);
            digX = digamma(1.0 - x) + Math.PI / Math.tan(Math.PI * f);
        }

        return digX;
    }
1017609
public static double sinh(double x) {
      boolean negate = false;
      if (x != x) {
          return x;
      }

      // sinh[z] = (exp(z) - exp(-z) / 2

      // for values of z larger than about 20,
      // exp(-z) can be ignored in comparison with exp(z)

      if (x > 20) {
          if (x >= LOG_MAX_VALUE) {
              // Avoid overflow (MATH-905).
              final double t = exp(0.5 * x);
              return (0.5 * t) * t;
          } else {
              return 0.5 * exp(x);
          }
      } else if (x < -20) {
          if (x <= -LOG_MAX_VALUE) {
              // Avoid overflow (MATH-905).
              final double t = exp(-0.5 * x);
              return (-0.5 * t) * t;
          } else {
              return -0.5 * exp(-x);
          }
      }

      if (x == 0) {
          return x;
      }

      if (x < 0.0) {
          x = -x;
          negate = true;
      }

      double result;

      if (x > 0.25) {
          double hiPrec[] = new double[2];
          exp(x, 0.0, hiPrec);

          double ya = hiPrec[0] + hiPrec[1];
          double yb = -(ya - hiPrec[0] - hiPrec[1]);

          double temp = ya * HEX_40000000;
          double yaa = ya + temp - temp;
          double yab = ya - yaa;

          // recip = 1/y
          double recip = 1.0/ya;
          temp = recip * HEX_40000000;
          double recipa = recip + temp - temp;
          double recipb = recip - recipa;

          // Correct for rounding in division
          recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
          // Account for yb
          recipb += -yb * recip * recip;

          recipa = -recipa;
          recipb = -recipb;

          // y = y + 1/y
          temp = ya + recipa;
          yb += -(temp - ya - recipa);
          ya = temp;
          temp = ya + recipb;
          yb += -(temp - ya - recipb);
          ya = temp;

          result = ya + yb;
          result *= 0.5;
      }
      else {
          double hiPrec[] = new double[2];
          expm1(x, hiPrec);

          double ya = hiPrec[0] + hiPrec[1];
          double yb = -(ya - hiPrec[0] - hiPrec[1]);

          /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
          double denom = 1.0 + ya;
          double denomr = 1.0 / denom;
          double denomb = -(denom - 1.0 - ya) + yb;
          double ratio = ya * denomr;
          double temp = ratio * HEX_40000000;
          double ra = ratio + temp - temp;
          double rb = ratio - ra;

          temp = denom * HEX_40000000;
          double za = denom + temp - temp;
          double zb = denom - za;

          rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;

          // Adjust for yb
          rb += yb*denomr;                        // numerator
          rb += -ya * denomb * denomr * denomr;   // denominator

          // y = y - 1/y
          temp = ya + ra;
          yb += -(temp - ya - ra);
          ya = temp;
          temp = ya + rb;
          yb += -(temp - ya - rb);
          ya = temp;

          result = ya + yb;
          result *= 0.5;
      }

      if (negate) {
          result = -result;
      }

      return result;
    }
629647
public static double max(double[] x) {
        double m = Double.NEGATIVE_INFINITY;

        for (double n : x) {
            if (n > m) {
                m = n;
            }
        }

        return m;
    }
1217742
public static double log1p(final double x) {
        if (x == -1) {
            return Double.NEGATIVE_INFINITY;
        }

        if (x == Double.POSITIVE_INFINITY) {
            return Double.POSITIVE_INFINITY;
        }

        if (x > 1e-6 ||
            x < -1e-6) {
            final double xpa = 1 + x;
            final double xpb = -(xpa - 1 - x);

            final double[] hiPrec = new double[2];
            final double lores = log(xpa, hiPrec);
            if (Double.isInfinite(lores)) { // Don't allow this to be converted to NaN
                return lores;
            }

            // Do a taylor series expansion around xpa:
            //   f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
            final double fx1 = xpb / xpa;
            final double epsilon = 0.5 * fx1 + 1;
            return epsilon * fx1 + hiPrec[1] + hiPrec[0];
        } else {
            // Value is small |x| < 1e6, do a Taylor series centered on 1.
            final double y = (x * F_1_3 - F_1_2) * x + 1;
            return y * x;
        }
    }
664730
public static double Function(double x) {
        double[] P =
                {
                        1.60119522476751861407E-4,
                        1.19135147006586384913E-3,
                        1.04213797561761569935E-2,
                        4.76367800457137231464E-2,
                        2.07448227648435975150E-1,
                        4.94214826801497100753E-1,
                        9.99999999999999996796E-1
                };
        double[] Q =
                {
                        -2.31581873324120129819E-5,
                        5.39605580493303397842E-4,
                        -4.45641913851797240494E-3,
                        1.18139785222060435552E-2,
                        3.58236398605498653373E-2,
                        -2.34591795718243348568E-1,
                        7.14304917030273074085E-2,
                        1.00000000000000000320E0
                };

        double p, z;

        double q = Math.abs(x);

        if (q > 33.0) {
            if (x < 0.0) {
                p = Math.floor(q);

                if (p == q) {
                    try {
                        throw new ArithmeticException("Overflow");
                    } catch (Exception e) {
                        e.printStackTrace();
                    }
                }

                z = q - p;
                if (z > 0.5) {
                    p += 1.0;
                    z = q - p;
                }
                z = q * Math.sin(Math.PI * z);

                if (z == 0.0) {
                    try {
                        throw new ArithmeticException("Overflow");
                    } catch (Exception e) {
                        e.printStackTrace();
                    }
                }

                z = Math.abs(z);
                z = Math.PI / (z * Stirling(q));

                return -z;
            } else {
                return Stirling(x);
            }
        }

        z = 1.0;
        while (x >= 3.0) {
            x -= 1.0;
            z *= x;
        }

        while (x < 0.0) {
            if (x == 0.0) {
                throw new ArithmeticException();
            } else if (x > -1.0E-9) {
                return (z / ((1.0 + 0.5772156649015329 * x) * x));
            }
            z /= x;
            x += 1.0;
        }

        while (x < 2.0) {
            if (x == 0.0) {
                throw new ArithmeticException();
            } else if (x < 1.0E-9) {
                return (z / ((1.0 + 0.5772156649015329 * x) * x));
            }

            z /= x;
            x += 1.0;
        }

        if ((x == 2.0) || (x == 3.0)) return z;

        x -= 2.0;
        p = Special.Polevl(x, P, 6);
        q = Special.Polevl(x, Q, 7);
        return z * p / q;

    }
250080
@Pure
	public double toGlobalX(double x) {
		final double adjustedX = x - this.translationX.get();
		return this.invertX.get() ? -adjustedX : adjustedX;
	}
800240
@SuppressWarnings("WeakerAccess")
    public static double pow2(double x) {
        final double exponent = Math.floor(x) - 1;
        x = x - exponent;
        double m = 3 - sqrt(7 - 3 * x);
        return Math.pow(2, exponent + 1) * m;
    }
1010103
public static double factorial(int x) {
    double result = 1.0;
    for (int i=x; i>1; i--) {
      result *= i;
    }
    return result;
  }
354785
public static BigFloat atanh(BigFloat x) {
		return x.context.valueOf(BigDecimalMath.atanh(x.value, x.context.mathContext));
	}
281662
public double interpolate(double... x) {
        if (x.length != this.x[0].length) {
            throw new IllegalArgumentException(String.format("Invalid input vector size: %d, expected: %d", x.length, this.x[0].length));
        }

        double sum = 0.0, sumw = 0.0;
        for (int i = 0; i < this.x.length; i++) {
            double f = rbf.f(Math.distance(x, this.x[i]));
            sumw += w[i] * f;
            sum += f;
        }

        return normalized ? sumw / sum : sumw;
    }