LogBeta.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.numbers.gamma;

  19. /**
  20.  * Computes \( log_e B(p, q) \).
  21.  * <p>
  22.  * This class is immutable.
  23.  * </p>
  24.  */
  25. public final class LogBeta {
  26.     /** The threshold value of 10 where the series expansion of the \( \Delta \) function applies. */
  27.     private static final double TEN = 10;
  28.     /** The threshold value of 2 for algorithm switch. */
  29.     private static final double TWO = 2;
  30.     /** The threshold value of 1000 for algorithm switch. */
  31.     private static final double THOUSAND = 1000;

  33.     /** The constant value of ½log 2π. */
  34.     private static final double HALF_LOG_TWO_PI = 0.9189385332046727;

  36.     /**
  37.      * The coefficients of the series expansion of the \( \Delta \) function.
  38.      * This function is defined as follows:
  39.      * \[
  40.      *  \Delta(x) = \log \Gamma(x) - (x - \frac{1}{2}) \log a + a - \frac{1}{2} \log 2\pi,
  41.      * \]
  42.      * <p>
  43.      * See equation (23) in Didonato and Morris (1992). The series expansion,
  44.      * which applies for \( x \geq 10 \), reads
  45.      * </p>
  46.      * \[
  47.      *  \Delta(x) = \frac{1}{x} \sum_{n = 0}^{14} d_n (\frac{10}{x})^{2 n}
  48.      * \]
  49.      */
  50.     private static final double[] DELTA = {
  51.         .833333333333333333333333333333E-01,
  52.         -.277777777777777777777777752282E-04,
  53.         .793650793650793650791732130419E-07,
  54.         -.595238095238095232389839236182E-09,
  55.         .841750841750832853294451671990E-11,
  56.         -.191752691751854612334149171243E-12,
  57.         .641025640510325475730918472625E-14,
  58.         -.295506514125338232839867823991E-15,
  59.         .179643716359402238723287696452E-16,
  60.         -.139228964661627791231203060395E-17,
  61.         .133802855014020915603275339093E-18,
  62.         -.154246009867966094273710216533E-19,
  63.         .197701992980957427278370133333E-20,
  64.         -.234065664793997056856992426667E-21,
  65.         .171348014966398575409015466667E-22
  66.     };

  68.     /** Private constructor. */
  69.     private LogBeta() {
  70.         // intentionally empty.
  71.     }

  73.     /**
  74.      * Returns the value of \( \Delta(b) - \Delta(a + b) \),
  75.      * with \( 0 \leq a \leq b \) and \( b \geq 10 \).
  76.      * Based on equations (26), (27) and (28) in Didonato and Morris (1992).
  77.      *
  78.      * @param a First argument.
  79.      * @param b Second argument.
  80.      * @return the value of \( \Delta(b) - \Delta(a + b) \)
  81.      * @throws IllegalArgumentException if {@code a < 0} or {@code a > b}
  82.      * @throws IllegalArgumentException if {@code b < 10}
  83.      */
  84.     private static double deltaMinusDeltaSum(final double a,
  85.                                              final double b) {
  86.         // Assumptions:
  87.         // Never called with a < 0; a > b; or b < 10

  89.         final double h = a / b;
  90.         final double p = h / (1 + h);
  91.         final double q = 1 / (1 + h);
  92.         final double q2 = q * q;
  93.         /*
  94.          * s[i] = 1 + q + ... - q**(2 * i)
  95.          */
  96.         final double[] s = new double[DELTA.length];
  97.         s[0] = 1;
  98.         for (int i = 1; i < s.length; i++) {
  99.             s[i] = 1 + (q + q2 * s[i - 1]);
  100.         }
  101.         /*
  102.          * w = Delta(b) - Delta(a + b)
  103.          */
  104.         final double sqrtT = 10 / b;
  105.         final double t = sqrtT * sqrtT;
  106.         double w = DELTA[DELTA.length - 1] * s[s.length - 1];
  107.         for (int i = DELTA.length - 2; i >= 0; i--) {
  108.             w = t * w + DELTA[i] * s[i];
  109.         }
  110.         return w * p / b;
  111.     }

  113.     /**
  114.      * Returns the value of \( \Delta(p) + \Delta(q) - \Delta(p + q) \),
  115.      * with \( p, q \geq 10 \).
  116.      * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
  117.      * {@code DBCORR}.
  118.      *
  119.      * @param p First argument.
  120.      * @param q Second argument.
  121.      * @return the value of \( \Delta(p) + \Delta(q) - \Delta(p + q) \).
  122.      * @throws IllegalArgumentException if {@code p < 10} or {@code q < 10}.
  123.      */
  124.     private static double sumDeltaMinusDeltaSum(final double p,
  125.                                                 final double q) {

  127.         if (p < TEN) {
  128.             throw new GammaException(GammaException.OUT_OF_RANGE, p, TEN, Double.POSITIVE_INFINITY);
  129.         }
  130.         if (q < TEN) {
  131.             throw new GammaException(GammaException.OUT_OF_RANGE, q, TEN, Double.POSITIVE_INFINITY);
  132.         }

  134.         final double a = Math.min(p, q);
  135.         final double b = Math.max(p, q);
  136.         final double sqrtT = 10 / a;
  137.         final double t = sqrtT * sqrtT;
  138.         double z = DELTA[DELTA.length - 1];
  139.         for (int i = DELTA.length - 2; i >= 0; i--) {
  140.             z = t * z + DELTA[i];
  141.         }
  142.         return z / a + deltaMinusDeltaSum(a, b);
  143.     }

  145.     /**
  146.      * Returns the value of \( \log B(p, q) \) for \( 0 \leq x \leq 1 \) and \( p, q &gt; 0 \).
  147.      * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
  148.      * {@code DBETLN}.
  149.      *
  150.      * @param p First argument.
  151.      * @param q Second argument.
  152.      * @return the value of \( \log B(p, q) \), or {@code NaN} if
  153.      * {@code p <= 0} or {@code q <= 0}.
  154.      */
  155.     public static double value(double p,
  156.                                double q) {
  157.         if (Double.isNaN(p) ||
  158.             Double.isNaN(q) ||
  159.             p <= 0 ||
  160.             q <= 0) {
  161.             return Double.NaN;
  162.         }

  164.         final double a = Math.min(p, q);
  165.         final double b = Math.max(p, q);
  166.         if (a >= TEN) {
  167.             final double w = sumDeltaMinusDeltaSum(a, b);
  168.             final double h = a / b;
  169.             final double c = h / (1 + h);
  170.             final double u = -(a - 0.5) * Math.log(c);
  171.             final double v = b * Math.log1p(h);
  172.             if (u <= v) {
  173.                 return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
  174.             }
  175.             return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
  176.         } else if (a > TWO) {
  177.             if (b > THOUSAND) {
  178.                 final int n = (int) Math.floor(a - 1);
  179.                 double prod = 1;
  180.                 double ared = a;
  181.                 for (int i = 0; i < n; i++) {
  182.                     ared -= 1;
  183.                     prod *= ared / (1 + ared / b);
  184.                 }
  185.                 return (Math.log(prod) - n * Math.log(b)) +
  186.                         (LogGamma.value(ared) +
  187.                          logGammaMinusLogGammaSum(ared, b));
  188.             }
  189.             double prod1 = 1;
  190.             double ared = a;
  191.             while (ared > TWO) {
  192.                 ared -= 1;
  193.                 final double h = ared / b;
  194.                 prod1 *= h / (1 + h);
  195.             }
  196.             if (b < TEN) {
  197.                 double prod2 = 1;
  198.                 double bred = b;
  199.                 while (bred > TWO) {
  200.                     bred -= 1;
  201.                     prod2 *= bred / (ared + bred);
  202.                 }
  203.                 return Math.log(prod1) +
  204.                        Math.log(prod2) +
  205.                        (LogGamma.value(ared) +
  206.                        (LogGamma.value(bred) -
  207.                         LogGammaSum.value(ared, bred)));
  208.             }
  209.             return Math.log(prod1) +
  210.                    LogGamma.value(ared) +
  211.                    logGammaMinusLogGammaSum(ared, b);
  212.         } else if (a >= 1) {
  213.             if (b > TWO) {
  214.                 if (b < TEN) {
  215.                     double prod = 1;
  216.                     double bred = b;
  217.                     while (bred > TWO) {
  218.                         bred -= 1;
  219.                         prod *= bred / (a + bred);
  220.                     }
  221.                     return Math.log(prod) +
  222.                            (LogGamma.value(a) +
  223.                             (LogGamma.value(bred) -
  224.                              LogGammaSum.value(a, bred)));
  225.                 }
  226.                 return LogGamma.value(a) +
  227.                     logGammaMinusLogGammaSum(a, b);
  228.             }
  229.             return LogGamma.value(a) +
  230.                    LogGamma.value(b) -
  231.                    LogGammaSum.value(a, b);
  232.         } else {
  233.             if (b >= TEN) {
  234.                 return LogGamma.value(a) +
  235.                        logGammaMinusLogGammaSum(a, b);
  236.             }
  237.             // The original NSWC implementation was
  238.             //   LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
  239.             // but the following command turned out to be more accurate.
  240.             // Note: Check for overflow that occurs if a and/or b are tiny.
  241.             final double beta = Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b);
  242.             if (Double.isFinite(beta)) {
  243.                 return Math.log(beta);
  244.             }
  245.             return LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
  246.         }
  247.     }

  249.     /**
  250.      * Returns the value of \( \log ( \Gamma(b) / \Gamma(a + b) ) \)
  251.      * for \( a \geq 0 \) and \( b \geq 10 \).
  252.      * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
  253.      * {@code DLGDIV}.
  254.      *
  255.      * <p>This method assumes \( a \leq b \).
  256.      *
  257.      * @param a First argument.
  258.      * @param b Second argument.
  259.      * @return the value of \( \log(\Gamma(b) / \Gamma(a + b) \).
  260.      * @throws IllegalArgumentException if {@code a < 0} or {@code b < 10}.
  261.      */
  262.     private static double logGammaMinusLogGammaSum(double a,
  263.                                                    double b) {
  264.         if (a < 0) {
  265.             throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, Double.POSITIVE_INFINITY);
  266.         }
  267.         if (b < TEN) {
  268.             throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY);
  269.         }

  271.         /*
  272.          * d = a + b - 0.5
  273.          */
  274.         // Assumption: always called with a <= b.
  275.         final double d = b + (a - 0.5);
  276.         final double w = deltaMinusDeltaSum(a, b);

  278.         final double u = d * Math.log1p(a / b);
  279.         final double v = a * (Math.log(b) - 1);

  281.         return u <= v ?
  282.             (w - u) - v :
  283.             (w - v) - u;
  284.     }
  285. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Numbers 1.2