View Javadoc
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  
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)
22  
23  package org.apache.commons.numbers.gamma;
24  
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;
29  
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).
68  
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;
94  
95      /** Private constructor. */
96      private BoostBeta() {
97          // intentionally empty.
98      }
99  
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         }
113 
114         final double c = p + q;
115 
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         }
129 
130         // Create input a > b
131         final double a = p < q ? q : p;
132         final double b = p < q ? p : q;
133 
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         }
148 
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);
156 
157         return result;
158     }
159 
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         }
195 
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         }
208 
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     }
215 
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     }
234 
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         }
261 
262         double result;
263 
264         final double c = a + b;
265 
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);
276 
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         }
395 
396         return result;
397     }
398 
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     }
411 
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     }
425 
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     }
438 
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     }
452 
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     }
465 
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     }
479 
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     }
492 
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     }
506 
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         //
531 
532         if (!(x >= 0 && x <= 1)) {
533             // Domain error
534             return Double.NaN;
535         }
536 
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         }
560 
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         }
573 
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         }
580 
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             }
602 
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         }
614 
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;
694 
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             }
752 
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]
757 
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     }
823 
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;
837 
838         if (normalised) {
839             final double c = a + b;
840 
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));
847 
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.
875 
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.
880 
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         }
889 
890         final double eps = pol.getEps();
891         final int maxIterations = pol.getMaxIterations();
892 
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;
902 
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         };
911 
912         return BoostTools.sumSeries(gen, eps, maxIterations, s0) * rescale;
913     }
914 
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)
934 
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     }
944 
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);
958 
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         }
995 
996         return result;
997     }
998 
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.
1017 
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         }
1023 
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     }
1029 
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);
1048 
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         }
1063 
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))
1075 
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).
1079 
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     }
1097 
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;
1123 
1124         return prefix;
1125     }
1126 
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;
1201 
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     }
1238 
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         }
1261 
1262         final double eps = pol.getEps();
1263         final int maxIterations = pol.getMaxIterations();
1264 
1265         final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
1266             /** Iteration. */
1267             private int m;
1268 
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;
1274 
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);
1278 
1279                 ++m;
1280                 return Coefficient.of(aN, bN);
1281             }
1282         };
1283 
1284         // Note: The first generated term a0 is discarded
1285         final double fract = GeneralizedContinuedFraction.value(gen, eps, maxIterations);
1286         return result / fract;
1287     }
1288 
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         }
1317 
1318         final double eps = pol.getEps();
1319         final int maxIterations = pol.getMaxIterations();
1320 
1321         final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
1322             /** Iteration. */
1323             private int n;
1324 
1325             @Override
1326             public Coefficient get() {
1327                 // https://functions.wolfram.com/GammaBetaErf/Beta3/10/0001/
1328 
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                 }
1344 
1345                 n = m + 1;
1346                 return Coefficient.of(aN, 1);
1347             }
1348         };
1349 
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 }