BoostBeta.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.  */

  18. //  Copyright John Maddock 2006.
  19. //  Use, modification and distribution are subject to the
  20. //  Boost Software License, Version 1.0. (See accompanying file
  21. //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

  23. package org.apache.commons.numbers.gamma;

  25. import java.util.function.DoubleSupplier;
  26. import java.util.function.Supplier;
  27. import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction;
  28. import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;

  30. /**
  31.  * Implementation of the
  32.  * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">regularized beta functions</a> and
  33.  * <a href="https://mathworld.wolfram.com/IncompleteBetaFunction.html">incomplete beta functions</a>.
  34.  *
  35.  * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
  36.  * {@code c++} implementation {@code <boost/math/special_functions/beta.hpp>}.
  37.  * All work is copyright to the original authors and subject to the Boost Software License.
  38.  *
  39.  * @see
  40.  * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta.html">
  41.  * Boost C++ beta functions</a>
  42.  */
  43. final class BoostBeta {
  44.     //
  45.     // Code ported from Boost 1.77.0
  46.     //
  47.     // boost/math/special_functions/beta.hpp
  48.     //
  49.     // Original code comments are preserved.
  50.     //
  51.     // Changes to the Boost implementation:
  52.     // - Update method names to replace underscores with camel case
  53.     // - Remove checks for under/overflow. In this implementation no error is raised
  54.     //   for overflow (infinity is returned) or underflow (sub-normal or zero is returned).
  55.     //   This follows the conventions in java.lang.Math for the same conditions.
  56.     // - Removed the pointer p_derivative in the betaIncompleteImp. This is used
  57.     //   in the Boost code for the beta_inv functions for a derivative
  58.     //   based inverse function. This is currently not supported.
  59.     // - Altered series generators to use integer counters added to the double term
  60.     //   replacing directly incrementing a double term. When the term is large it cannot
  61.     //   be incremented: 1e16 + 1 == 1e16.
  62.     // - Added the use of the classic continued fraction representation for cases
  63.     //   where the Boost implementation detects sub-normal terms and does not evaluate.
  64.     // - Updated the method used to compute the binomialCoefficient. This can use a
  65.     //   series evaluation when n > max factorial given that n - k < 40.
  66.     // - Changed convergence criteria for betaSmallBLargeASeries to stop when r has
  67.     //   no effect on the sum. The Boost code uses machine epsilon (ignoring the policy eps).

  69.     /** Default epsilon value for relative error.
  70.      * This is equal to the Boost constant {@code boost::math::EPSILON}. */
  71.     private static final double EPSILON = 0x1.0p-52;
  72.     /** Approximate value for ln(Double.MAX_VALUE).
  73.      * This is equal to the Boost constant {@code boost::math::tools::log_max_value<double>()}.
  74.      * No term {@code x} should be used in {@code exp(x)} if {@code x > LOG_MAX_VALUE} to avoid
  75.      * overflow. */
  76.     private static final int LOG_MAX_VALUE = 709;
  77.     /** Approximate value for ln(Double.MIN_VALUE).
  78.      * This is equal to the Boost constant {@code boost::math::tools::log_min_value<double>()}.
  79.      * No term {@code x} should be used in {@code exp(x)} if {@code x < LOG_MIN_VALUE} to avoid
  80.      * underflow to sub-normal or zero. */
  81.     private static final int LOG_MIN_VALUE = -708;
  82.     /** pi/2. */
  83.     private static final double HALF_PI = Math.PI / 2;
  84.     /** The largest factorial that can be represented as a double.
  85.      * This is equal to the Boost constant {@code boost::math::max_factorial<double>::value}. */
  86.     private static final int MAX_FACTORIAL = 170;
  87.     /** Size of the table of Pn's.
  88.      * Equal to {@code Pn_size<double>} suitable for 16-20 digit accuracy. */
  89.     private static final int PN_SIZE = 30;
  90.     /** 2^53. Used to scale sub-normal values. */
  91.     private static final double TWO_POW_53 = 0x1.0p53;
  92.     /** 2^-53. Used to rescale values. */
  93.     private static final double TWO_POW_M53 = 0x1.0p-53;

  95.     /** Private constructor. */
  96.     private BoostBeta() {
  97.         // intentionally empty.
  98.     }

  100.     /**
  101.      * Beta function.
  102.      * <p>\[ B(p, q) = \frac{\Gamma(p) \Gamma(q)}{\Gamma(p+q)} \]
  103.      *
  104.      * @param p Argument p
  105.      * @param q Argument q
  106.      * @return beta value
  107.      */
  108.     static double beta(double p, double q) {
  109.         if (!(p > 0 && q > 0)) {
  110.             // Domain error
  111.             return Double.NaN;
  112.         }

  114.         final double c = p + q;

  116.         // Special cases:
  117.         if (c == p && q < EPSILON) {
  118.             return 1 / q;
  119.         } else if (c == q && p < EPSILON) {
  120.             return 1 / p;
  121.         }
  122.         if (q == 1) {
  123.             return 1 / p;
  124.         } else if (p == 1) {
  125.             return 1 / q;
  126.         } else if (c < EPSILON) {
  127.             return (c / p) / q;
  128.         }

  130.         // Create input a > b
  131.         final double a = p < q ? q : p;
  132.         final double b = p < q ? p : q;

  134.         // Lanczos calculation:
  135.         final double agh = a + BoostGamma.Lanczos.GMH;
  136.         final double bgh = b + BoostGamma.Lanczos.GMH;
  137.         final double cgh = c + BoostGamma.Lanczos.GMH;
  138.         double result = BoostGamma.Lanczos.lanczosSumExpGScaled(a) *
  139.                 (BoostGamma.Lanczos.lanczosSumExpGScaled(b) / BoostGamma.Lanczos.lanczosSumExpGScaled(c));
  140.         final double ambh = a - 0.5f - b;
  141.         if (Math.abs(b * ambh) < cgh * 100 && a > 100) {
  142.             // Special case where the base of the power term is close to 1
  143.             // compute (1+x)^y instead:
  144.             result *= Math.exp(ambh * Math.log1p(-b / cgh));
  145.         } else {
  146.             result *= Math.pow(agh / cgh, ambh);
  147.         }

  149.         if (cgh > 1e10f) {
  150.             // this avoids possible overflow, but appears to be marginally less accurate:
  151.             result *= Math.pow((agh / cgh) * (bgh / cgh), b);
  152.         } else {
  153.             result *= Math.pow((agh * bgh) / (cgh * cgh), b);
  154.         }
  155.         result *= Math.sqrt(Math.E / bgh);

  157.         return result;
  158.     }

  160.     /**
  161.      * Derivative of the regularised incomplete beta.
  162.      * <p>\[ \frac{\delta}{\delta x} I_x(a, b) = \frac{(1-x)^{b-1} x^{a-1}}{\B(a, b)} \]
  163.      *
  164.      * <p>Adapted from {@code boost::math::ibeta_derivative}.
  165.      *
  166.      * @param a Argument a
  167.      * @param b Argument b
  168.      * @param x Argument x
  169.      * @return ibeta derivative
  170.      */
  171.     static double ibetaDerivative(double a, double b, double x) {
  172.         //
  173.         // start with the usual error checks:
  174.         //
  175.         if (!(a > 0 && b > 0) || !(x >= 0 && x <= 1)) {
  176.             // Domain error
  177.             return Double.NaN;
  178.         }
  179.         //
  180.         // Now the corner cases:
  181.         //
  182.         if (x == 0) {
  183.             if (a > 1) {
  184.                 return 0;
  185.             }
  186.             // a == 1 : return 1 / beta(a, b) == b
  187.             return a == 1 ? b : Double.POSITIVE_INFINITY;
  188.         } else if (x == 1) {
  189.             if (b > 1) {
  190.                 return 0;
  191.             }
  192.             // b == 1 : return 1 / beta(a, b) == a
  193.             return b == 1 ? a : Double.POSITIVE_INFINITY;
  194.         }

  196.         // Update with extra edge cases
  197.         if (b == 1) {
  198.             // ibeta = x^a
  199.             return a * Math.pow(x, a - 1);
  200.         }
  201.         if (a == 1) {
  202.             // ibeta = 1 - (1-x)^b
  203.             if (x >= 0.5) {
  204.                 return b * Math.pow(1 - x, b - 1);
  205.             }
  206.             return b * Math.exp(Math.log1p(-x) * (b - 1));
  207.         }

  209.         //
  210.         // Now the regular cases:
  211.         //
  212.         final double y = (1 - x) * x;
  213.         return ibetaPowerTerms(a, b, x, 1 - x, true, 1 / y);
  214.     }

  216.     /**
  217.      * Compute the leading power terms in the incomplete Beta.
  218.      *
  219.      * <p>Utility function to call
  220.      * {@link #ibetaPowerTerms(double, double, double, double, boolean, double)}
  221.      * using a multiplication prefix of {@code 1.0}.
  222.      *
  223.      * @param a Argument a
  224.      * @param b Argument b
  225.      * @param x Argument x
  226.      * @param y Argument 1-x
  227.      * @param normalised true to divide by beta(a, b)
  228.      * @return incomplete beta power terms
  229.      */
  230.     private static double ibetaPowerTerms(double a, double b, double x,
  231.             double y, boolean normalised) {
  232.         return ibetaPowerTerms(a, b, x, y, normalised, 1);
  233.     }

  235.     /**
  236.      * Compute the leading power terms in the incomplete Beta.
  237.      *
  238.      * <pre>
  239.      * (x^a)(y^b)/Beta(a,b) when normalised, and
  240.      * (x^a)(y^b) otherwise.
  241.      * </pre>
  242.      *
  243.      * <p>Almost all of the error in the incomplete beta comes from this function:
  244.      * particularly when a and b are large. Computing large powers are *hard*
  245.      * though, and using logarithms just leads to horrendous cancellation errors.
  246.      *
  247.      * @param a Argument a
  248.      * @param b Argument b
  249.      * @param x Argument x
  250.      * @param y Argument 1-x
  251.      * @param normalised true to divide by beta(a, b)
  252.      * @param prefix Prefix to multiply by the result
  253.      * @return incomplete beta power terms
  254.      */
  255.     private static double ibetaPowerTerms(double a, double b, double x,
  256.             double y, boolean normalised, double prefix) {
  257.         if (!normalised) {
  258.             // can we do better here?
  259.             return Math.pow(x, a) * Math.pow(y, b);
  260.         }

  262.         double result;

  264.         final double c = a + b;

  266.         // combine power terms with Lanczos approximation:
  267.         final double agh = a + BoostGamma.Lanczos.GMH;
  268.         final double bgh = b + BoostGamma.Lanczos.GMH;
  269.         final double cgh = c + BoostGamma.Lanczos.GMH;
  270.         result = BoostGamma.Lanczos.lanczosSumExpGScaled(c) /
  271.                 (BoostGamma.Lanczos.lanczosSumExpGScaled(a) * BoostGamma.Lanczos.lanczosSumExpGScaled(b));
  272.         result *= prefix;
  273.         // combine with the leftover terms from the Lanczos approximation:
  274.         result *= Math.sqrt(bgh / Math.E);
  275.         result *= Math.sqrt(agh / cgh);

  277.         // l1 and l2 are the base of the exponents minus one:
  278.         double l1 = (x * b - y * agh) / agh;
  279.         double l2 = (y * a - x * bgh) / bgh;
  280.         if (Math.min(Math.abs(l1), Math.abs(l2)) < 0.2) {
  281.             // when the base of the exponent is very near 1 we get really
  282.             // gross errors unless extra care is taken:
  283.             if (l1 * l2 > 0 || Math.min(a, b) < 1) {
  284.                 //
  285.                 // This first branch handles the simple cases where either:
  286.                 //
  287.                 // * The two power terms both go in the same direction
  288.                 // (towards zero or towards infinity). In this case if either
  289.                 // term overflows or underflows, then the product of the two must
  290.                 // do so also.
  291.                 // * Alternatively if one exponent is less than one, then we
  292.                 // can't productively use it to eliminate overflow or underflow
  293.                 // from the other term. Problems with spurious overflow/underflow
  294.                 // can't be ruled out in this case, but it is *very* unlikely
  295.                 // since one of the power terms will evaluate to a number close to 1.
  296.                 //
  297.                 if (Math.abs(l1) < 0.1) {
  298.                     result *= Math.exp(a * Math.log1p(l1));
  299.                 } else {
  300.                     result *= Math.pow((x * cgh) / agh, a);
  301.                 }
  302.                 if (Math.abs(l2) < 0.1) {
  303.                     result *= Math.exp(b * Math.log1p(l2));
  304.                 } else {
  305.                     result *= Math.pow((y * cgh) / bgh, b);
  306.                 }
  307.             } else if (Math.max(Math.abs(l1), Math.abs(l2)) < 0.5) {
  308.                 //
  309.                 // Both exponents are near one and both the exponents are
  310.                 // greater than one and further these two
  311.                 // power terms tend in opposite directions (one towards zero,
  312.                 // the other towards infinity), so we have to combine the terms
  313.                 // to avoid any risk of overflow or underflow.
  314.                 //
  315.                 // We do this by moving one power term inside the other, we have:
  316.                 //
  317.                 //   (1 + l1)^a * (1 + l2)^b
  318.                 // = ((1 + l1)*(1 + l2)^(b/a))^a
  319.                 // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1
  320.                 //                                = exp((b/a) * log(1 + l2)) - 1
  321.                 //
  322.                 // The tricky bit is deciding which term to move inside :-)
  323.                 // By preference we move the larger term inside, so that the
  324.                 // size of the largest exponent is reduced. However, that can
  325.                 // only be done as long as l3 (see above) is also small.
  326.                 //
  327.                 final boolean smallA = a < b;
  328.                 final double ratio = b / a;
  329.                 if ((smallA && ratio * l2 < 0.1) || (!smallA && l1 / ratio > 0.1)) {
  330.                     double l3 = Math.expm1(ratio * Math.log1p(l2));
  331.                     l3 = l1 + l3 + l3 * l1;
  332.                     l3 = a * Math.log1p(l3);
  333.                     result *= Math.exp(l3);
  334.                 } else {
  335.                     double l3 = Math.expm1(Math.log1p(l1) / ratio);
  336.                     l3 = l2 + l3 + l3 * l2;
  337.                     l3 = b * Math.log1p(l3);
  338.                     result *= Math.exp(l3);
  339.                 }
  340.             } else if (Math.abs(l1) < Math.abs(l2)) {
  341.                 // First base near 1 only:
  342.                 double l = a * Math.log1p(l1) + b * Math.log((y * cgh) / bgh);
  343.                 if (l <= LOG_MIN_VALUE || l >= LOG_MAX_VALUE) {
  344.                     l += Math.log(result);
  345.                     // Update: Allow overflow to infinity
  346.                     result = Math.exp(l);
  347.                 } else {
  348.                     result *= Math.exp(l);
  349.                 }
  350.             } else {
  351.                 // Second base near 1 only:
  352.                 double l = b * Math.log1p(l2) + a * Math.log((x * cgh) / agh);
  353.                 if (l <= LOG_MIN_VALUE || l >= LOG_MAX_VALUE) {
  354.                     l += Math.log(result);
  355.                     // Update: Allow overflow to infinity
  356.                     result = Math.exp(l);
  357.                 } else {
  358.                     result *= Math.exp(l);
  359.                 }
  360.             }
  361.         } else {
  362.             // general case:
  363.             final double b1 = (x * cgh) / agh;
  364.             final double b2 = (y * cgh) / bgh;
  365.             l1 = a * Math.log(b1);
  366.             l2 = b * Math.log(b2);
  367.             if (l1 >= LOG_MAX_VALUE || l1 <= LOG_MIN_VALUE || l2 >= LOG_MAX_VALUE || l2 <= LOG_MIN_VALUE) {
  368.                 // Oops, under/overflow, sidestep if we can:
  369.                 if (a < b) {
  370.                     final double p1 = Math.pow(b2, b / a);
  371.                     final double l3 = a * (Math.log(b1) + Math.log(p1));
  372.                     if (l3 < LOG_MAX_VALUE && l3 > LOG_MIN_VALUE) {
  373.                         result *= Math.pow(p1 * b1, a);
  374.                     } else {
  375.                         l2 += l1 + Math.log(result);
  376.                         // Update: Allow overflow to infinity
  377.                         result = Math.exp(l2);
  378.                     }
  379.                 } else {
  380.                     final double p1 = Math.pow(b1, a / b);
  381.                     final double l3 = (Math.log(p1) + Math.log(b2)) * b;
  382.                     if (l3 < LOG_MAX_VALUE && l3 > LOG_MIN_VALUE) {
  383.                         result *= Math.pow(p1 * b2, b);
  384.                     } else {
  385.                         l2 += l1 + Math.log(result);
  386.                         // Update: Allow overflow to infinity
  387.                         result = Math.exp(l2);
  388.                     }
  389.                 }
  390.             } else {
  391.                 // finally the normal case:
  392.                 result *= Math.pow(b1, a) * Math.pow(b2, b);
  393.             }
  394.         }

  396.         return result;
  397.     }

  399.     /**
  400.      * Full incomplete beta.
  401.      * <p>\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
  402.      *
  403.      * @param a Argument a
  404.      * @param b Argument b
  405.      * @param x Argument x
  406.      * @return lower beta value
  407.      */
  408.     static double beta(double a, double b, double x) {
  409.         return betaIncompleteImp(a, b, x, Policy.getDefault(), false, false);
  410.     }

  412.     /**
  413.      * Full incomplete beta.
  414.      * <p>\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
  415.      *
  416.      * @param a Argument a
  417.      * @param b Argument b
  418.      * @param x Argument x
  419.      * @param policy Function evaluation policy
  420.      * @return lower beta value
  421.      */
  422.     static double beta(double a, double b, double x, Policy policy) {
  423.         return betaIncompleteImp(a, b, x, policy, false, false);
  424.     }

  426.     /**
  427.      * Complement of the incomplete beta.
  428.      * <p>\[ betac(a, b, x) = beta(b, a, 1-x) = B_{1-x}(b,a) \]
  429.      *
  430.      * @param a Argument a
  431.      * @param b Argument b
  432.      * @param x Argument x
  433.      * @return upper beta value
  434.      */
  435.     static double betac(double a, double b, double x) {
  436.         return betaIncompleteImp(a, b, x, Policy.getDefault(), false, true);
  437.     }

  439.     /**
  440.      * Complement of the incomplete beta.
  441.      * <p>\[ betac(a, b, x) = beta(b, a, 1-x) = B_{1-x}(b,a) \]
  442.      *
  443.      * @param a Argument a
  444.      * @param b Argument b
  445.      * @param x Argument x
  446.      * @param policy Function evaluation policy
  447.      * @return upper beta value
  448.      */
  449.     static double betac(double a, double b, double x, Policy policy) {
  450.         return betaIncompleteImp(a, b, x, policy, false, true);
  451.     }

  453.     /**
  454.      * Regularised incomplete beta.
  455.      * <p>\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
  456.      *
  457.      * @param a Argument a
  458.      * @param b Argument b
  459.      * @param x Argument x
  460.      * @return p
  461.      */
  462.     static double ibeta(double a, double b, double x) {
  463.         return betaIncompleteImp(a, b, x, Policy.getDefault(), true, false);
  464.     }

  466.     /**
  467.      * Regularised incomplete beta.
  468.      * <p>\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
  469.      *
  470.      * @param a Argument a
  471.      * @param b Argument b
  472.      * @param x Argument x
  473.      * @param policy Function evaluation policy
  474.      * @return p
  475.      */
  476.     static double ibeta(double a, double b, double x, Policy policy) {
  477.         return betaIncompleteImp(a, b, x, policy, true, false);
  478.     }

  480.     /**
  481.      * Complement of the regularised incomplete beta.
  482.      * <p>\[ ibetac(a, b, x) = 1 - I_x(a,b) = I_{1-x}(b, a) \]
  483.      *
  484.      * @param a Argument a
  485.      * @param b Argument b
  486.      * @param x Argument x
  487.      * @return q
  488.      */
  489.     static double ibetac(double a, double b, double x) {
  490.         return betaIncompleteImp(a, b, x, Policy.getDefault(), true, true);
  491.     }

  493.     /**
  494.      * Complement of the regularised incomplete beta.
  495.      * <p>\[ ibetac(a, b, x) = 1 - I_x(a,b) = I_{1-x}(b, a) \]
  496.      *
  497.      * @param a Argument a
  498.      * @param b Argument b
  499.      * @param x Argument x
  500.      * @param policy Function evaluation policy
  501.      * @return q
  502.      */
  503.     static double ibetac(double a, double b, double x, Policy policy) {
  504.         return betaIncompleteImp(a, b, x, policy, true, true);
  505.     }

  507.     /**
  508.      * Main incomplete beta entry point, handles all four incomplete betas.
  509.      * Adapted from {@code boost::math::detail::ibeta_imp}.
  510.      *
  511.      * <p>The Boost code has a pointer {@code p_derivative} that can be set to the
  512.      * value of the derivative. This is used for the inverse incomplete
  513.      * beta functions. It is not required for the forward evaluation functions.
  514.      *
  515.      * @param a Argument a
  516.      * @param b Argument b
  517.      * @param x Argument x
  518.      * @param pol Function evaluation policy
  519.      * @param normalised true to compute the regularised value
  520.      * @param inv true to compute the complement value
  521.      * @return incomplete beta value
  522.      */
  523.     private static double betaIncompleteImp(double a, double b, double x,
  524.             Policy pol, boolean normalised, boolean inv) {
  525.         //
  526.         // The incomplete beta function implementation:
  527.         // This is just a big bunch of spaghetti code to divide up the
  528.         // input range and select the right implementation method for
  529.         // each domain:
  530.         //

  532.         if (!(x >= 0 && x <= 1)) {
  533.             // Domain error
  534.             return Double.NaN;
  535.         }

  537.         if (normalised) {
  538.             if (!(a >= 0 && b >= 0)) {
  539.                 // Domain error
  540.                 return Double.NaN;
  541.             }
  542.             // extend to a few very special cases:
  543.             if (a == 0) {
  544.                 if (b == 0) {
  545.                     // a and b cannot both be zero
  546.                     return Double.NaN;
  547.                 }
  548.                 // Assume b > 0
  549.                 return inv ? 0 : 1;
  550.             } else if (b == 0) {
  551.                 // assume a > 0
  552.                 return inv ? 1 : 0;
  553.             }
  554.         } else {
  555.             if (!(a > 0 && b > 0)) {
  556.                 // Domain error
  557.                 return Double.NaN;
  558.             }
  559.         }

  561.         if (x == 0) {
  562.             if (inv) {
  563.                 return normalised ? 1 : beta(a, b);
  564.             }
  565.             return 0;
  566.         }
  567.         if (x == 1) {
  568.             if (!inv) {
  569.                 return normalised ? 1 : beta(a, b);
  570.             }
  571.             return 0;
  572.         }

  574.         if (a == 0.5f && b == 0.5f) {
  575.             // We have an arcsine distribution:
  576.             final double z = inv ? 1 - x : x;
  577.             final double asin = Math.asin(Math.sqrt(z));
  578.             return normalised ? asin / HALF_PI : 2 * asin;
  579.         }

  581.         boolean invert = inv;
  582.         double y = 1 - x;
  583.         if (a == 1) {
  584.             // swap(a, b)
  585.             double tmp = a;
  586.             a = b;
  587.             b = tmp;
  588.             // swap(x, y)
  589.             tmp = x;
  590.             x = y;
  591.             y = tmp;
  592.             invert = !invert;
  593.         }
  594.         if (b == 1) {
  595.             //
  596.             // Special case see:
  597.             // http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
  598.             //
  599.             if (a == 1) {
  600.                 return invert ? y : x;
  601.             }

  603.             double p;
  604.             if (y < 0.5) {
  605.                 p = invert ? -Math.expm1(a * Math.log1p(-y)) : Math.exp(a * Math.log1p(-y));
  606.             } else {
  607.                 p = invert ? -BoostMath.powm1(x, a) : Math.pow(x, a);
  608.             }
  609.             if (!normalised) {
  610.                 p /= a;
  611.             }
  612.             return p;
  613.         }

  615.         double fract;
  616.         if (Math.min(a, b) <= 1) {
  617.             if (x > 0.5) {
  618.                 // swap(a, b)
  619.                 double tmp = a;
  620.                 a = b;
  621.                 b = tmp;
  622.                 // swap(x, y)
  623.                 tmp = x;
  624.                 x = y;
  625.                 y = tmp;
  626.                 invert = !invert;
  627.             }
  628.             if (Math.max(a, b) <= 1) {
  629.                 // Both a,b < 1:
  630.                 if (a >= Math.min(0.2, b) || Math.pow(x, a) <= 0.9) {
  631.                     if (invert) {
  632.                         fract = -(normalised ? 1 : beta(a, b));
  633.                         invert = false;
  634.                         fract = -ibetaSeries(a, b, x, fract, normalised, pol);
  635.                     } else {
  636.                         fract = ibetaSeries(a, b, x, 0, normalised, pol);
  637.                     }
  638.                 } else {
  639.                     // swap(a, b)
  640.                     double tmp = a;
  641.                     a = b;
  642.                     b = tmp;
  643.                     // swap(x, y)
  644.                     tmp = x;
  645.                     x = y;
  646.                     y = tmp;
  647.                     invert = !invert;
  648.                     if (y >= 0.3) {
  649.                         if (invert) {
  650.                             fract = -(normalised ? 1 : beta(a, b));
  651.                             invert = false;
  652.                             fract = -ibetaSeries(a, b, x, fract, normalised, pol);
  653.                         } else {
  654.                             fract = ibetaSeries(a, b, x, 0, normalised, pol);
  655.                         }
  656.                     } else {
  657.                         // Sidestep on a, and then use the series representation:
  658.                         final double prefix;
  659.                         if (normalised) {
  660.                             prefix = 1;
  661.                         } else {
  662.                             prefix = risingFactorialRatio(a + b, a, 20);
  663.                         }
  664.                         fract = ibetaAStep(a, b, x, y, 20, normalised);
  665.                         if (invert) {
  666.                             fract -= normalised ? 1 : beta(a, b);
  667.                             invert = false;
  668.                             fract = -betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
  669.                         } else {
  670.                             fract = betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
  671.                         }
  672.                     }
  673.                 }
  674.             } else {
  675.                 // One of a, b < 1 only:
  676.                 if (b <= 1 || (x < 0.1 && Math.pow(b * x, a) <= 0.7)) {
  677.                     if (invert) {
  678.                         fract = -(normalised ? 1 : beta(a, b));
  679.                         invert = false;
  680.                         fract = -ibetaSeries(a, b, x, fract, normalised, pol);
  681.                     } else {
  682.                         fract = ibetaSeries(a, b, x, 0, normalised, pol);
  683.                     }
  684.                 } else {
  685.                     // swap(a, b)
  686.                     double tmp = a;
  687.                     a = b;
  688.                     b = tmp;
  689.                     // swap(x, y)
  690.                     tmp = x;
  691.                     x = y;
  692.                     y = tmp;
  693.                     invert = !invert;

  695.                     if (y >= 0.3) {
  696.                         if (invert) {
  697.                             fract = -(normalised ? 1 : beta(a, b));
  698.                             invert = false;
  699.                             fract = -ibetaSeries(a, b, x, fract, normalised, pol);
  700.                         } else {
  701.                             fract = ibetaSeries(a, b, x, 0, normalised, pol);
  702.                         }
  703.                     } else if (a >= 15) {
  704.                         if (invert) {
  705.                             fract = -(normalised ? 1 : beta(a, b));
  706.                             invert = false;
  707.                             fract = -betaSmallBLargeASeries(a, b, x, y, fract, 1, pol, normalised);
  708.                         } else {
  709.                             fract = betaSmallBLargeASeries(a, b, x, y, 0, 1, pol, normalised);
  710.                         }
  711.                     } else {
  712.                         // Sidestep to improve errors:
  713.                         final double prefix;
  714.                         if (normalised) {
  715.                             prefix = 1;
  716.                         } else {
  717.                             prefix = risingFactorialRatio(a + b, a, 20);
  718.                         }
  719.                         fract = ibetaAStep(a, b, x, y, 20, normalised);
  720.                         if (invert) {
  721.                             fract -= normalised ? 1 : beta(a, b);
  722.                             invert = false;
  723.                             fract = -betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
  724.                         } else {
  725.                             fract = betaSmallBLargeASeries(a + 20, b, x, y, fract, prefix, pol, normalised);
  726.                         }
  727.                     }
  728.                 }
  729.             }
  730.         } else {
  731.             // Both a,b >= 1:
  732.             // Note:
  733.             // median ~ (a - 1/3) / (a + b - 2/3) ~ a / (a + b)
  734.             // if x > a / (a + b) => a - (a + b) * x < 0
  735.             final double lambda;
  736.             if (a < b) {
  737.                 lambda = a - (a + b) * x;
  738.             } else {
  739.                 lambda = (a + b) * y - b;
  740.             }
  741.             if (lambda < 0) {
  742.                 // swap(a, b)
  743.                 double tmp = a;
  744.                 a = b;
  745.                 b = tmp;
  746.                 // swap(x, y)
  747.                 tmp = x;
  748.                 x = y;
  749.                 y = tmp;
  750.                 invert = !invert;
  751.             }

  753.             if (b < 40) {
  754.                 // Note: y != 1 check is required for non-zero x < epsilon
  755.                 if (Math.rint(a) == a && Math.rint(b) == b && a < (Integer.MAX_VALUE - 100) && y != 1) {
  756.                     // Here: a in [2, 2^31 - 102] && b in [2, 39]

  758.                     // relate to the binomial distribution and use a finite sum:
  759.                     final int k = (int) (a - 1);
  760.                     final int n = (int) (b + k);
  761.                     fract = binomialCCdf(n, k, x, y);
  762.                     if (!normalised) {
  763.                         fract *= beta(a, b);
  764.                     }
  765.                 } else if (b * x <= 0.7) {
  766.                     if (invert) {
  767.                         fract = -(normalised ? 1 : beta(a, b));
  768.                         invert = false;
  769.                         fract = -ibetaSeries(a, b, x, fract, normalised, pol);
  770.                     } else {
  771.                         fract = ibetaSeries(a, b, x, 0, normalised, pol);
  772.                     }
  773.                 } else if (a > 15) {
  774.                     // sidestep so we can use the series representation:
  775.                     int n = (int) b;
  776.                     if (n == b) {
  777.                         --n;
  778.                     }
  779.                     final double bbar = b - n;
  780.                     final double prefix;
  781.                     if (normalised) {
  782.                         prefix = 1;
  783.                     } else {
  784.                         prefix = risingFactorialRatio(a + bbar, bbar, n);
  785.                     }
  786.                     fract = ibetaAStep(bbar, a, y, x, n, normalised);
  787.                     fract = betaSmallBLargeASeries(a, bbar, x, y, fract, 1, pol, normalised);
  788.                     fract /= prefix;
  789.                 } else if (normalised) {
  790.                     // The formula here for the non-normalised case is tricky to figure
  791.                     // out (for me!!), and requires two pochhammer calculations rather
  792.                     // than one, so leave it for now and only use this in the normalized case....
  793.                     int n = (int) Math.floor(b);
  794.                     double bbar = b - n;
  795.                     if (bbar <= 0) {
  796.                         --n;
  797.                         bbar += 1;
  798.                     }
  799.                     fract = ibetaAStep(bbar, a, y, x, n, normalised);
  800.                     fract += ibetaAStep(a, bbar, x, y, 20, normalised);
  801.                     if (invert) {
  802.                         // Note this line would need changing if we ever enable this branch in
  803.                         // non-normalized case
  804.                         fract -= 1;
  805.                     }
  806.                     fract = betaSmallBLargeASeries(a + 20, bbar, x, y, fract, 1, pol, normalised);
  807.                     if (invert) {
  808.                         fract = -fract;
  809.                         invert = false;
  810.                     }
  811.                 } else {
  812.                     fract = ibetaFraction2(a, b, x, y, pol, normalised);
  813.                 }
  814.             } else {
  815.                 fract = ibetaFraction2(a, b, x, y, pol, normalised);
  816.             }
  817.         }
  818.         if (invert) {
  819.             return (normalised ? 1 : beta(a, b)) - fract;
  820.         }
  821.         return fract;
  822.     }

  824.     /**
  825.      * Series approximation to the incomplete beta.
  826.      *
  827.      * @param a Argument a
  828.      * @param b Argument b
  829.      * @param x Argument x
  830.      * @param s0 Initial sum for the series
  831.      * @param normalised true to compute the regularised value
  832.      * @param pol Function evaluation policy
  833.      * @return incomplete beta series
  834.      */
  835.     private static double ibetaSeries(double a, double b, double x, double s0, boolean normalised, Policy pol) {
  836.         double result;

  838.         if (normalised) {
  839.             final double c = a + b;

  841.             // incomplete beta power term, combined with the Lanczos approximation:
  842.             final double agh = a + BoostGamma.Lanczos.GMH;
  843.             final double bgh = b + BoostGamma.Lanczos.GMH;
  844.             final double cgh = c + BoostGamma.Lanczos.GMH;
  845.             result = BoostGamma.Lanczos.lanczosSumExpGScaled(c) /
  846.                     (BoostGamma.Lanczos.lanczosSumExpGScaled(a) * BoostGamma.Lanczos.lanczosSumExpGScaled(b));

  848.             final double l1 = Math.log(cgh / bgh) * (b - 0.5f);
  849.             final double l2 = Math.log(x * cgh / agh) * a;
  850.             //
  851.             // Check for over/underflow in the power terms:
  852.             //
  853.             if (l1 > LOG_MIN_VALUE && l1 < LOG_MAX_VALUE && l2 > LOG_MIN_VALUE && l2 < LOG_MAX_VALUE) {
  854.                 if (a * b < bgh * 10) {
  855.                     result *= Math.exp((b - 0.5f) * Math.log1p(a / bgh));
  856.                 } else {
  857.                     result *= Math.pow(cgh / bgh, b - 0.5f);
  858.                 }
  859.                 result *= Math.pow(x * cgh / agh, a);
  860.                 result *= Math.sqrt(agh / Math.E);
  861.             } else {
  862.                 //
  863.                 // Oh dear, we need logs, and this *will* cancel:
  864.                 //
  865.                 result = Math.log(result) + l1 + l2 + (Math.log(agh) - 1) / 2;
  866.                 result = Math.exp(result);
  867.             }
  868.         } else {
  869.             // Non-normalised, just compute the power:
  870.             result = Math.pow(x, a);
  871.         }
  872.         double rescale = 1.0;
  873.         if (result < Double.MIN_NORMAL) {
  874.             // Safeguard: series can't cope with denorms.

  876.             // Update:
  877.             // The entire series is only based on the magnitude of 'result'.
  878.             // If the first term can be added to s0 (e.g. if s0 == 0) then
  879.             // scale s0 and result, compute the series and then rescale.

  881.             // Intentional floating-point equality check.
  882.             if (s0 + result / a == s0) {
  883.                 return s0;
  884.             }
  885.             s0 *= TWO_POW_53;
  886.             result *= TWO_POW_53;
  887.             rescale = TWO_POW_M53;
  888.         }

  890.         final double eps = pol.getEps();
  891.         final int maxIterations = pol.getMaxIterations();

  893.         // Create effectively final 'result' for initialisation
  894.         final double result1 = result;
  895.         final DoubleSupplier gen = new DoubleSupplier() {
  896.             /** Next result term. */
  897.             private double result = result1;
  898.             /** Pochhammer term. */
  899.             private final double poch = -b;
  900.             /** Iteration. */
  901.             private int n;

  903.             @Override
  904.             public double getAsDouble() {
  905.                 final double r = result / (a + n);
  906.                 n++;
  907.                 result *= (n + poch) * x / n;
  908.                 return r;
  909.             }
  910.         };

  912.         return BoostTools.sumSeries(gen, eps, maxIterations, s0) * rescale;
  913.     }

  915.     /**
  916.      * Rising factorial ratio.
  917.      * This function is only needed for the non-regular incomplete beta,
  918.      * it computes the delta in:
  919.      * <pre>
  920.      * beta(a,b,x) = prefix + delta * beta(a+k,b,x)
  921.      * </pre>
  922.      * <p>It is currently only called for small k.
  923.      *
  924.      * @param a Argument a
  925.      * @param b Argument b
  926.      * @param k Argument k
  927.      * @return the rising factorial ratio
  928.      */
  929.     private static double risingFactorialRatio(double a, double b, int k) {
  930.         // calculate:
  931.         // (a)(a+1)(a+2)...(a+k-1)
  932.         // _______________________
  933.         // (b)(b+1)(b+2)...(b+k-1)

  935.         // This is only called with small k, for large k
  936.         // it is grossly inefficient, do not use outside it's
  937.         // intended purpose!!!
  938.         double result = 1;
  939.         for (int i = 0; i < k; ++i) {
  940.             result *= (a + i) / (b + i);
  941.         }
  942.         return result;
  943.     }

  945.     /**
  946.      * Binomial complementary cdf.
  947.      * For integer arguments we can relate the incomplete beta to the
  948.      * complement of the binomial distribution cdf and use this finite sum.
  949.      *
  950.      * @param n Argument n (called with {@code [2, 39] + k})
  951.      * @param k Argument k (called with {@code k in [1, Integer.MAX_VALUE - 102]})
  952.      * @param x Argument x
  953.      * @param y Argument 1-x
  954.      * @return Binomial complementary cdf
  955.      */
  956.     private static double binomialCCdf(int n, int k, double x, double y) {
  957.         double result = Math.pow(x, n);

  959.         if (result > Double.MIN_NORMAL) {
  960.             double term = result;
  961.             for (int i = n - 1; i > k; --i) {
  962.                 term *= ((i + 1) * y) / ((n - i) * x);
  963.                 result += term;
  964.             }
  965.         } else {
  966.             // First term underflows so we need to start at the mode of the
  967.             // distribution and work outwards:
  968.             int start = (int) (n * x);
  969.             if (start <= k + 1) {
  970.                 start = k + 2;
  971.             }
  972.             // Update:
  973.             // Carefully compute this term to guard against extreme parameterisation
  974.             result = binomialTerm(n, start, x, y);
  975.             if (result == 0) {
  976.                 // OK, starting slightly above the mode didn't work,
  977.                 // we'll have to sum the terms the old fashioned way:
  978.                 for (int i = start - 1; i > k; --i) {
  979.                     result += binomialTerm(n, i, x, y);
  980.                 }
  981.             } else {
  982.                 double term = result;
  983.                 final double startTerm = result;
  984.                 for (int i = start - 1; i > k; --i) {
  985.                     term *= ((i + 1) * y) / ((n - i) * x);
  986.                     result += term;
  987.                 }
  988.                 term = startTerm;
  989.                 for (int i = start + 1; i <= n; ++i) {
  990.                     term *= (n - i + 1) * x / (i * y);
  991.                     result += term;
  992.                 }
  993.             }
  994.         }

  996.         return result;
  997.     }

  999.     /**
  1000.      * Compute the binomial term.
  1001.      * <pre>
  1002.      * x^k * (1-x)^(n-k) * binomial(n, k)
  1003.      * </pre>
  1004.      * <p>This is a helper function used to guard against extreme values generated
  1005.      * in the term which can produce NaN from zero multiplied by infinity.
  1006.      *
  1007.      * @param n Argument n
  1008.      * @param k Argument k
  1009.      * @param x Argument x
  1010.      * @param y Argument 1-x
  1011.      * @return Binomial term
  1012.      */
  1013.     private static double binomialTerm(int n, int k, double x, double y) {
  1014.         // This function can be called with the following extreme that will overflow:
  1015.         // binomial(2147483545 + 39, 38) ~ 7.84899e309
  1016.         // Guard this as the power functions are very likely to be zero with large n and k.

  1018.         final double binom = binomialCoefficient(n, k);
  1019.         if (!Double.isFinite(binom)) {
  1020.             // Product of the power functions will be zero with large n and k
  1021.             return 0;
  1022.         }

  1024.         // The power terms are below 1.
  1025.         // Multiply by the largest so that any sub-normal term is used last
  1026.         // This method is called where x^n is sub-normal so assume the other term is larger.
  1027.         return binom * Math.pow(y, n - k) * Math.pow(x, k);
  1028.     }

  1030.     /**
  1031.      * Computes the binomial coefficient.
  1032.      *
  1033.      * <p>Adapted from
  1034.      * {@code org.apache.commons.numbers.combinatorics.BinomialCoefficientDouble}.
  1035.      *
  1036.      * <p>Note: This does not use {@code BinomialCoefficientDouble}
  1037.      * to avoid a circular dependency as the combinatorics depends on the gamma package.
  1038.      * No checks are made on the arguments.
  1039.      *
  1040.      * @param n Size of the set (must be positive).
  1041.      * @param k Size of the subsets to be counted (must be in [0, n]).
  1042.      * @return {@code n choose k}.
  1043.      */
  1044.     static double binomialCoefficient(int n, int k) {
  1045.         // Assume: n >= 0; 0 <= k <= n
  1046.         // Use symmetry
  1047.         final int m = Math.min(k, n - k);

  1049.         // This function is typically called with m <= 3 so handle these special cases
  1050.         if (m == 0) {
  1051.             return 1;
  1052.         }
  1053.         if (m == 1) {
  1054.             return n;
  1055.         }
  1056.         if (m == 2) {
  1057.             return 0.5 * n * (n - 1);
  1058.         }
  1059.         if (m == 3) {
  1060.             // Divide by 3 at the end to avoid using an 1/6 inexact initial multiplier.
  1061.             return 0.5 * n * (n - 1) * (n - 2) / 3;
  1062.         }

  1064.         double result;
  1065.         if (n <= MAX_FACTORIAL) {
  1066.             // Use fast table lookup:
  1067.             result = BoostGamma.uncheckedFactorial(n);
  1068.             // Smaller m will have a more accurate factorial
  1069.             result /= BoostGamma.uncheckedFactorial(m);
  1070.             result /= BoostGamma.uncheckedFactorial(n - m);
  1071.         } else {
  1072.             // Updated:
  1073.             // Do not use the beta function as per <boost/math/special_functions/binomial.hpp>,
  1074.             // e.g. result = 1 / (m * beta(m, n - m + 1)) == gamma(n+1) / (m * gamma(m) * gamma(n-m+1))

  1076.             // Use a summation as per BinomialCoefficientDouble which is more accurate
  1077.             // than using the beta function. This is only used up to m = 39 so
  1078.             // may overflow *only* on the final term (i.e. m << n when overflow can occur).

  1080.             result = 1;
  1081.             for (int i = 1; i < m; i++) {
  1082.                 result *= n - m + i;
  1083.                 result /= i;
  1084.             }
  1085.             // Final term may cause overflow.
  1086.             if (result * n > Double.MAX_VALUE) {
  1087.                 result /= m;
  1088.                 result *= n;
  1089.             } else {
  1090.                 result *= n;
  1091.                 result /= m;
  1092.             }
  1093.         }
  1094.         // convert to nearest integer:
  1095.         return Math.ceil(result - 0.5f);
  1096.     }

  1098.     /**
  1099.      * Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x).
  1100.      *
  1101.      * @param a Argument a
  1102.      * @param b Argument b
  1103.      * @param x Argument x
  1104.      * @param y Argument 1-x
  1105.      * @param k Argument k
  1106.      * @param normalised true to compute the regularised value
  1107.      * @return ibeta difference
  1108.      */
  1109.     private static double ibetaAStep(double a, double b, double x, double y, int k, boolean normalised) {
  1110.         double prefix = ibetaPowerTerms(a, b, x, y, normalised);
  1111.         prefix /= a;
  1112.         if (prefix == 0) {
  1113.             return prefix;
  1114.         }
  1115.         double sum = 1;
  1116.         double term = 1;
  1117.         // series summation from 0 to k-1:
  1118.         for (int i = 0; i < k - 1; ++i) {
  1119.             term *= (a + b + i) * x / (a + i + 1);
  1120.             sum += term;
  1121.         }
  1122.         prefix *= sum;

  1124.         return prefix;
  1125.     }

  1127.     /**
  1128.      * Computes beta(a, b, x) using small b and large a.
  1129.      *
  1130.      * @param a Argument a
  1131.      * @param b Argument b
  1132.      * @param x Argument x
  1133.      * @param y Argument 1-x
  1134.      * @param s0 Initial sum for the series
  1135.      * @param mult Multiplication prefix factor
  1136.      * @param pol Function evaluation policy
  1137.      * @param normalised true to compute the regularised value
  1138.      * @return beta series
  1139.      */
  1140.     private static double betaSmallBLargeASeries(double a, double b, double x, double y, double s0, double mult,
  1141.             Policy pol, boolean normalised) {
  1142.         //
  1143.         // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
  1144.         //
  1145.         // Some values we'll need later, these are Eq 9.1:
  1146.         //
  1147.         final double bm1 = b - 1;
  1148.         final double t = a + bm1 / 2;
  1149.         final double lx;
  1150.         if (y < 0.35) {
  1151.             lx = Math.log1p(-y);
  1152.         } else {
  1153.             lx = Math.log(x);
  1154.         }
  1155.         final double u = -t * lx;
  1156.         // and from from 9.2:
  1157.         final double h = BoostGamma.regularisedGammaPrefix(b, u);
  1158.         if (h <= Double.MIN_NORMAL) {
  1159.             // Update:
  1160.             // Boost returns s0.
  1161.             // If this is zero then compute an expected sub-normal value
  1162.             // using the classic continued fraction representation.
  1163.             if (s0 == 0) {
  1164.                 return ibetaFraction(a, b, x, y, pol, normalised);
  1165.             }
  1166.             return s0;
  1167.         }
  1168.         double prefix;
  1169.         if (normalised) {
  1170.             prefix = h / GammaRatio.delta(a, b);
  1171.             prefix /= Math.pow(t, b);
  1172.         } else {
  1173.             prefix = BoostGamma.fullIgammaPrefix(b, u) / Math.pow(t, b);
  1174.         }
  1175.         prefix *= mult;
  1176.         //
  1177.         // now we need the quantity Pn, unfortunately this is computed
  1178.         // recursively, and requires a full history of all the previous values
  1179.         // so no choice but to declare a big table and hope it's big enough...
  1180.         //
  1181.         final double[] p = new double[PN_SIZE];
  1182.         // see 9.3.
  1183.         p[0] = 1;
  1184.         //
  1185.         // Now an initial value for J, see 9.6:
  1186.         //
  1187.         double j = BoostGamma.gammaQ(b, u, pol) / h;
  1188.         //
  1189.         // Now we can start to pull things together and evaluate the sum in Eq 9:
  1190.         //
  1191.         // Value at N = 0
  1192.         double sum = s0 + prefix * j;
  1193.         // some variables we'll need:
  1194.         // 2*N+1
  1195.         int tnp1 = 1;
  1196.         double lx2 = lx / 2;
  1197.         lx2 *= lx2;
  1198.         double lxp = 1;
  1199.         final double t4 = 4 * t * t;
  1200.         double b2n = b;

  1202.         for (int n = 1; n < PN_SIZE; ++n) {
  1203.             //
  1204.             // begin by evaluating the next Pn from Eq 9.4:
  1205.             //
  1206.             tnp1 += 2;
  1207.             p[n] = 0;
  1208.             int tmp1 = 3;
  1209.             for (int m = 1; m < n; ++m) {
  1210.                 final double mbn = m * b - n;
  1211.                 p[n] += mbn * p[n - m] / BoostGamma.uncheckedFactorial(tmp1);
  1212.                 tmp1 += 2;
  1213.             }
  1214.             p[n] /= n;
  1215.             p[n] += bm1 / BoostGamma.uncheckedFactorial(tnp1);
  1216.             //
  1217.             // Now we want Jn from Jn-1 using Eq 9.6:
  1218.             //
  1219.             j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
  1220.             lxp *= lx2;
  1221.             b2n += 2;
  1222.             //
  1223.             // pull it together with Eq 9:
  1224.             //
  1225.             final double r = prefix * p[n] * j;
  1226.             final double previous = sum;
  1227.             sum += r;
  1228.             // Update:
  1229.             // This continues until convergence at machine epsilon
  1230.             // |r| < eps * |sum|; r < 1
  1231.             // |r| / eps < |sum|; r > 1
  1232.             if (sum == previous) {
  1233.                 break;
  1234.             }
  1235.         }
  1236.         return sum;
  1237.     }

  1239.     /**
  1240.      * Evaluate the incomplete beta via a continued fraction representation.
  1241.      *
  1242.      * <p>Note: This is not a generic continued fraction. The formula is from <a
  1243.      * href="https://dl.acm.org/doi/10.1145/131766.131776">Didonato and Morris</a> and is
  1244.      * only used when a and b are above 1. See <a
  1245.      * href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta/ibeta_function.html">Incomplete
  1246.      * Beta Function: Implementation</a>.
  1247.      *
  1248.      * @param a Argument a
  1249.      * @param b Argument b
  1250.      * @param x Argument x
  1251.      * @param y Argument 1-x
  1252.      * @param pol Function evaluation policy
  1253.      * @param normalised true to compute the regularised value
  1254.      * @return incomplete beta
  1255.      */
  1256.     static double ibetaFraction2(double a, double b, double x, double y, Policy pol, boolean normalised) {
  1257.         final double result = ibetaPowerTerms(a, b, x, y, normalised);
  1258.         if (result == 0) {
  1259.             return result;
  1260.         }

  1262.         final double eps = pol.getEps();
  1263.         final int maxIterations = pol.getMaxIterations();

  1265.         final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
  1266.             /** Iteration. */
  1267.             private int m;

  1269.             @Override
  1270.             public Coefficient get() {
  1271.                 double aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x;
  1272.                 final double denom = a + 2 * m - 1;
  1273.                 aN /= denom * denom;

  1275.                 double bN = m;
  1276.                 bN += (m * (b - m) * x) / (a + 2 * m - 1);
  1277.                 bN += ((a + m) * (a * y - b * x + 1 + m * (2 - x))) / (a + 2 * m + 1);

  1279.                 ++m;
  1280.                 return Coefficient.of(aN, bN);
  1281.             }
  1282.         };

  1284.         // Note: The first generated term a0 is discarded
  1285.         final double fract = GeneralizedContinuedFraction.value(gen, eps, maxIterations);
  1286.         return result / fract;
  1287.     }

  1289.     /**
  1290.      * Evaluate the incomplete beta via the classic continued fraction representation.
  1291.      *
  1292.      * <p>Note: This is not part of the Boost C++ implementation.
  1293.      * This is a generic continued fraction applicable to all arguments.
  1294.      * It is used when alternative methods are unsuitable due to the presence of sub-normal
  1295.      * terms and the result is expected to be sub-normal. In this case the original Boost
  1296.      * implementation would return zero.
  1297.      *
  1298.      * <p>This continued fraction was the evaluation method used in Commons Numbers 1.0.
  1299.      * This method has been updated to use the {@code ibetaPowerTerms} function to compute
  1300.      * the power terms. Reversal of the arguments to call {@code 1 - ibetaFraction(b, a, 1 - x)}
  1301.      * is not performed as the result is expected to be very small and this optimisation
  1302.      * for accuracy is not required.
  1303.      *
  1304.      * @param a Argument a
  1305.      * @param b Argument b
  1306.      * @param x Argument x
  1307.      * @param y Argument 1-x
  1308.      * @param pol Function evaluation policy
  1309.      * @param normalised true to compute the regularised value
  1310.      * @return incomplete beta
  1311.      */
  1312.     static double ibetaFraction(double a, double b, double x, double y, Policy pol, boolean normalised) {
  1313.         final double result = ibetaPowerTerms(a, b, x, y, normalised);
  1314.         if (result == 0) {
  1315.             return result;
  1316.         }

  1318.         final double eps = pol.getEps();
  1319.         final int maxIterations = pol.getMaxIterations();

  1321.         final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
  1322.             /** Iteration. */
  1323.             private int n;

  1325.             @Override
  1326.             public Coefficient get() {
  1327.                 // https://functions.wolfram.com/GammaBetaErf/Beta3/10/0001/

  1329.                 final int m = n;
  1330.                 final int k = m / 2;
  1331.                 final double aN;
  1332.                 if ((m & 0x1) == 0) {
  1333.                     // even
  1334.                     // m = 2k
  1335.                     aN = (k * (b - k) * x) /
  1336.                         ((a + m - 1) * (a + m));
  1337.                 } else {
  1338.                     // odd
  1339.                     // k = (m - 1) / 2 due to integer truncation
  1340.                     // m - 1 = 2k
  1341.                     aN = -((a + k) * (a + b + k) * x) /
  1342.                         ((a + m - 1) * (a + m));
  1343.                 }

  1345.                 n = m + 1;
  1346.                 return Coefficient.of(aN, 1);
  1347.             }
  1348.         };

  1350.         // Note: The first generated term a0 is discarded
  1351.         final double fract = GeneralizedContinuedFraction.value(gen, eps, maxIterations);
  1352.         return (result / a) / fract;
  1353.     }
  1354. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Numbers 1.2