All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.carrot2.mahout.collections.Arithmetic Maven / Gradle / Ivy

/* Imported from Mahout. */package org.carrot2.mahout.collections;


public class Arithmetic extends Constants {
  // for method stirlingCorrection(...)
  private static final double[] stirlingCorrection = {
      0.0,
      8.106146679532726e-02, 4.134069595540929e-02,
      2.767792568499834e-02, 2.079067210376509e-02,
      1.664469118982119e-02, 1.387612882307075e-02,
      1.189670994589177e-02, 1.041126526197209e-02,
      9.255462182712733e-03, 8.330563433362871e-03,
      7.573675487951841e-03, 6.942840107209530e-03,
      6.408994188004207e-03, 5.951370112758848e-03,
      5.554733551962801e-03, 5.207655919609640e-03,
      4.901395948434738e-03, 4.629153749334029e-03,
      4.385560249232324e-03, 4.166319691996922e-03,
      3.967954218640860e-03, 3.787618068444430e-03,
      3.622960224683090e-03, 3.472021382978770e-03,
      3.333155636728090e-03, 3.204970228055040e-03,
      3.086278682608780e-03, 2.976063983550410e-03,
      2.873449362352470e-03, 2.777674929752690e-03,
  };

  // for method logFactorial(...)
  // log(k!) for k = 0, ..., 29
  private static final double[] logFactorials = {
      0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
      1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
      6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
      12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
      19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
      27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
      36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
      45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
      54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
      64.55753862700633106, 67.88974313718153498, 71.25703896716800901
  };

  // k! for k = 0, ..., 20
  private static final long[] longFactorials = {
      1L,
      1L,
      2L,
      6L,
      24L,
      120L,
      720L,
      5040L,
      40320L,
      362880L,
      3628800L,
      39916800L,
      479001600L,
      6227020800L,
      87178291200L,
      1307674368000L,
      20922789888000L,
      355687428096000L,
      6402373705728000L,
      121645100408832000L,
      2432902008176640000L
  };

  // k! for k = 21, ..., 170
  private static final double[] doubleFactorials = {
      5.109094217170944E19,
      1.1240007277776077E21,
      2.585201673888498E22,
      6.204484017332394E23,
      1.5511210043330984E25,
      4.032914611266057E26,
      1.0888869450418352E28,
      3.048883446117138E29,
      8.841761993739701E30,
      2.652528598121911E32,
      8.222838654177924E33,
      2.6313083693369355E35,
      8.68331761881189E36,
      2.952327990396041E38,
      1.0333147966386144E40,
      3.719933267899013E41,
      1.3763753091226346E43,
      5.23022617466601E44,
      2.0397882081197447E46,
      8.15915283247898E47,
      3.34525266131638E49,
      1.4050061177528801E51,
      6.041526306337384E52,
      2.6582715747884495E54,
      1.196222208654802E56,
      5.502622159812089E57,
      2.5862324151116827E59,
      1.2413915592536068E61,
      6.082818640342679E62,
      3.0414093201713376E64,
      1.5511187532873816E66,
      8.06581751709439E67,
      4.274883284060024E69,
      2.308436973392413E71,
      1.2696403353658264E73,
      7.109985878048632E74,
      4.052691950487723E76,
      2.350561331282879E78,
      1.386831185456898E80,
      8.32098711274139E81,
      5.075802138772246E83,
      3.146997326038794E85,
      1.9826083154044396E87,
      1.2688693218588414E89,
      8.247650592082472E90,
      5.443449390774432E92,
      3.6471110918188705E94,
      2.48003554243683E96,
      1.7112245242814127E98,
      1.1978571669969892E100,
      8.504785885678624E101,
      6.123445837688612E103,
      4.470115461512686E105,
      3.307885441519387E107,
      2.4809140811395404E109,
      1.8854947016660506E111,
      1.451830920282859E113,
      1.1324281178206295E115,
      8.94618213078298E116,
      7.15694570462638E118,
      5.797126020747369E120,
      4.7536433370128435E122,
      3.94552396972066E124,
      3.314240134565354E126,
      2.8171041143805494E128,
      2.4227095383672744E130,
      2.107757298379527E132,
      1.854826422573984E134,
      1.6507955160908465E136,
      1.4857159644817605E138,
      1.3520015276784033E140,
      1.2438414054641305E142,
      1.156772507081641E144,
      1.0873661566567426E146,
      1.0329978488239061E148,
      9.916779348709491E149,
      9.619275968248216E151,
      9.426890448883248E153,
      9.332621544394415E155,
      9.332621544394418E157,
      9.42594775983836E159,
      9.614466715035125E161,
      9.902900716486178E163,
      1.0299016745145631E166,
      1.0813967582402912E168,
      1.1462805637347086E170,
      1.2265202031961373E172,
      1.324641819451829E174,
      1.4438595832024942E176,
      1.5882455415227423E178,
      1.7629525510902457E180,
      1.974506857221075E182,
      2.2311927486598138E184,
      2.543559733472186E186,
      2.925093693493014E188,
      3.393108684451899E190,
      3.96993716080872E192,
      4.6845258497542896E194,
      5.574585761207606E196,
      6.689502913449135E198,
      8.094298525273444E200,
      9.875044200833601E202,
      1.2146304367025332E205,
      1.506141741511141E207,
      1.882677176888926E209,
      2.3721732428800483E211,
      3.0126600184576624E213,
      3.856204823625808E215,
      4.974504222477287E217,
      6.466855489220473E219,
      8.471580690878813E221,
      1.1182486511960037E224,
      1.4872707060906847E226,
      1.99294274616152E228,
      2.690472707318049E230,
      3.6590428819525483E232,
      5.0128887482749884E234,
      6.917786472619482E236,
      9.615723196941089E238,
      1.3462012475717523E241,
      1.8981437590761713E243,
      2.6953641378881633E245,
      3.8543707171800694E247,
      5.550293832739308E249,
      8.047926057471989E251,
      1.1749972043909107E254,
      1.72724589045464E256,
      2.5563239178728637E258,
      3.8089226376305687E260,
      5.7133839564458575E262,
      8.627209774233244E264,
      1.3113358856834527E267,
      2.0063439050956838E269,
      3.0897696138473515E271,
      4.789142901463393E273,
      7.471062926282892E275,
      1.1729568794264134E278,
      1.8532718694937346E280,
      2.946702272495036E282,
      4.714723635992061E284,
      7.590705053947223E286,
      1.2296942187394494E289,
      2.0044015765453032E291,
      3.287218585534299E293,
      5.423910666131583E295,
      9.003691705778434E297,
      1.5036165148649983E300,
      2.5260757449731988E302,
      4.2690680090047056E304,
      7.257415615308004E306
  };

  
  protected Arithmetic() {
  }

  
  public static double binomial(double n, long k) {
    if (k < 0) {
      return 0;
    }
    if (k == 0) {
      return 1;
    }
    if (k == 1) {
      return n;
    }

    // binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
    double a = n - k + 1;
    double b = 1;
    double binomial = 1;
    for (long i = k; i-- > 0;) {
      binomial *= (a++) / (b++);
    }
    return binomial;
  }

  
  public static double binomial(long n, long k) {
    if (k < 0) {
      return 0;
    }
    if (k == 0 || k == n) {
      return 1;
    }
    if (k == 1 || k == n - 1) {
      return n;
    }

    // try quick version and see whether we get numeric overflows.
    // factorial(..) is O(1); requires no loop; only a table lookup.
    if (n > k) {
      int max = longFactorials.length + doubleFactorials.length;
      if (n < max) { // if (n! < inf && k! < inf)
        double n_fac = factorial((int) n);
        double k_fac = factorial((int) k);
        double n_minus_k_fac = factorial((int) (n - k));
        double nk = n_minus_k_fac * k_fac;
        if (nk != Double.POSITIVE_INFINITY) { // no numeric overflow?
          // now this is completely safe and accurate
          return n_fac / nk;
        }
      }
      if (k > n / 2) {
        k = n - k;
      } // quicker
    }

    // binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
    long a = n - k + 1;
    long b = 1;
    double binomial = 1;
    for (long i = k; i-- > 0;) {
      binomial *= ((double) (a++)) / (b++);
    }
    return binomial;
  }

  
  public static long ceil(double value) {
    return Math.round(Math.ceil(value));
  }

  
  public static double chbevl(double x, double[] coef, int N) throws ArithmeticException {

    int p = 0;

    double b0 = coef[p++];
    double b1 = 0.0;
    int i = N - 1;

    double b2;
    do {
      b2 = b1;
      b1 = b0;
      b0 = x * b1 - b2 + coef[p++];
    } while (--i > 0);

    return (0.5 * (b0 - b2));
  }

  
  private static double factorial(int k) {
    if (k < 0) {
      throw new IllegalArgumentException();
    }

    int length1 = longFactorials.length;
    if (k < length1) {
      return longFactorials[k];
    }

    int length2 = doubleFactorials.length;
    if (k < length1 + length2) {
      return doubleFactorials[k - length1];
    } else {
      return Double.POSITIVE_INFINITY;
    }
  }

  
  public static long floor(double value) {
    return Math.round(Math.floor(value));
  }

  
  public static double log(double base, double value) {
    return Math.log(value) / Math.log(base);
  }

  
  public static double log10(double value) {
    // 1.0 / Math.log(10) == 0.43429448190325176
    return Math.log(value) * 0.43429448190325176;
  }

  
  public static double log2(double value) {
    // 1.0 / Math.log(2) == 1.4426950408889634
    return Math.log(value) * 1.4426950408889634;
  }

  
  public static double logFactorial(int k) {
    if (k >= 30) {

      double r = 1.0 / (double) k;
      double rr = r * r;
      double C7 = -5.95238095238095238e-04;
      double C5 = 7.93650793650793651e-04;
      double C3 = -2.77777777777777778e-03;
      double C1 = 8.33333333333333333e-02;
      double C0 = 9.18938533204672742e-01;
      return (k + 0.5) * Math.log(k) - k + C0 + r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
    } else {
      return logFactorials[k];
    }
  }

  
  public static long longFactorial(int k) throws IllegalArgumentException {
    if (k < 0) {
      throw new IllegalArgumentException("Negative k");
    }

    if (k < longFactorials.length) {
      return longFactorials[k];
    }
    throw new IllegalArgumentException("Overflow");
  }

  
  public static double stirlingCorrection(int k) {

    if (k > 30) {
      double r = 1.0 / (double) k;
      double rr = r * r;
      double C7 = -5.95238095238095238e-04;     //  -1/1680
      double C5 = 7.93650793650793651e-04;     //  +1/1260
      double C3 = -2.77777777777777778e-03;     //  -1/360
      double C1 = 8.33333333333333333e-02;     //  +1/12
      return r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
    } else {
      return stirlingCorrection[k];
    }
  }

}




© 2015 - 2025 Weber Informatics LLC | Privacy Policy