This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
The following commit(s) were added to refs/heads/master by this push:
new 29388f2 Add Commons Numbers 1.0 implementation to ErfPerformance
benchmark
29388f2 is described below
commit 29388f2a26621539891ef6e2068919ca77fb5cb6
Author: aherbert <aherb...@apache.org>
AuthorDate: Tue Nov 30 11:30:01 2021 +0000
Add Commons Numbers 1.0 implementation to ErfPerformance benchmark
Previously the code called RegularizedGamma.P or Q. Since these have
been updated the 1.0 version has been preserved for a reference
benchmark.
---
.../numbers/examples/jmh/gamma/ErfPerformance.java | 167 ++++++++++++++++++++-
1 file changed, 166 insertions(+), 1 deletion(-)
diff --git
a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/ErfPerformance.java
b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/ErfPerformance.java
index 07726c7..3df24b7 100644
---
a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/ErfPerformance.java
+++
b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/ErfPerformance.java
@@ -25,11 +25,12 @@ import java.util.function.DoubleUnaryOperator;
import java.util.function.Supplier;
import java.util.stream.DoubleStream;
import org.apache.commons.numbers.core.Precision;
+import org.apache.commons.numbers.fraction.ContinuedFraction;
import org.apache.commons.numbers.gamma.Erf;
import org.apache.commons.numbers.gamma.Erfc;
import org.apache.commons.numbers.gamma.InverseErf;
import org.apache.commons.numbers.gamma.InverseErfc;
-import org.apache.commons.numbers.gamma.RegularizedGamma;
+import org.apache.commons.numbers.gamma.LogGamma;
import org.openjdk.jmh.annotations.Benchmark;
import org.openjdk.jmh.annotations.BenchmarkMode;
import org.openjdk.jmh.annotations.Fork;
@@ -622,6 +623,170 @@ public class ErfPerformance {
}
/**
+ * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";>
+ * Regularized Gamma functions</a>.
+ *
+ * <p>This is the Commons Numbers 1.0 implementation. Later versions of
+ * RegularizedGamma changed to compute using more than the continued
fraction
+ * representation (Q) or lower gamma series representation (P).
+ *
+ * <p>The ContinuedFraction and LogGamma class use the current version
+ * and are not preserved from Commons Numbers 1.0.
+ */
+ private static final class RegularizedGamma {
+ /** Private constructor. */
+ private RegularizedGamma() {
+ // intentionally empty.
+ }
+
+ /**
+ * \( P(a, x) \) <a
href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";>
+ * regularized Gamma function</a>.
+ *
+ * Class is immutable.
+ */
+ static final class P {
+ /** Prevent instantiation. */
+ private P() {}
+
+ /**
+ * Computes the regularized gamma function \( P(a, x) \).
+ *
+ * The implementation of this method is based on:
+ * <ul>
+ * <li>
+ * <a
href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";>
+ * Regularized Gamma Function</a>, equation (1)
+ * </li>
+ * <li>
+ * <a
href="http://mathworld.wolfram.com/IncompleteGammaFunction.html";>
+ * Incomplete Gamma Function</a>, equation (4).
+ * </li>
+ * <li>
+ * <a
href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html";>
+ * Confluent Hypergeometric Function of the First Kind</a>,
equation (1).
+ * </li>
+ * </ul>
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @param epsilon Tolerance in continued fraction evaluation.
+ * @param maxIterations Maximum number of iterations in continued
fraction evaluation.
+ * @return \( P(a, x) \).
+ * @throws ArithmeticException if the continued fraction fails to
converge.
+ */
+ static double value(double a,
+ double x,
+ double epsilon,
+ int maxIterations) {
+ if (Double.isNaN(a) ||
+ Double.isNaN(x) ||
+ a <= 0 ||
+ x < 0) {
+ return Double.NaN;
+ } else if (x == 0) {
+ return 0;
+ } else if (x >= a + 1) {
+ // Q should converge faster in this case.
+ return 1 - RegularizedGamma.Q.value(a, x, epsilon,
maxIterations);
+ } else {
+ // Series.
+ double n = 0; // current element index
+ double an = 1 / a; // n-th element in the series
+ double sum = an; // partial sum
+ while (Math.abs(an / sum) > epsilon &&
+ n < maxIterations &&
+ sum < Double.POSITIVE_INFINITY) {
+ // compute next element in the series
+ n += 1;
+ an *= x / (a + n);
+
+ // update partial sum
+ sum += an;
+ }
+ if (n >= maxIterations) {
+ throw new ArithmeticException("Max iterations
exceeded: " + maxIterations);
+ } else if (Double.isInfinite(sum)) {
+ return 1;
+ } else {
+ // Ensure result is in the range [0, 1]
+ final double result = Math.exp(-x + (a * Math.log(x))
- LogGamma.value(a)) * sum;
+ return result > 1.0 ? 1.0 : result;
+ }
+ }
+ }
+ }
+
+ /**
+ * Creates the \( Q(a, x) \equiv 1 - P(a, x) \) <a
href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";>
+ * regularized Gamma function</a>.
+ *
+ * Class is immutable.
+ */
+ static final class Q {
+ /** Prevent instantiation. */
+ private Q() {}
+
+ /**
+ * Computes the regularized gamma function \( Q(a, x) = 1 - P(a,
x) \).
+ *
+ * The implementation of this method is based on:
+ * <ul>
+ * <li>
+ * <a
href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";>
+ * Regularized Gamma Function</a>, equation (1).
+ * </li>
+ * <li>
+ * <a
href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/";>
+ * Regularized incomplete gamma function: Continued fraction
representations
+ * (formula 06.08.10.0003)</a>
+ * </li>
+ * </ul>
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @param epsilon Tolerance in continued fraction evaluation.
+ * @param maxIterations Maximum number of iterations in continued
fraction evaluation.
+ * @throws ArithmeticException if the continued fraction fails to
converge.
+ * @return \( Q(a, x) \).
+ */
+ static double value(final double a,
+ double x,
+ double epsilon,
+ int maxIterations) {
+ if (Double.isNaN(a) ||
+ Double.isNaN(x) ||
+ a <= 0 ||
+ x < 0) {
+ return Double.NaN;
+ } else if (x == 0) {
+ return 1;
+ } else if (x < a + 1) {
+ // P should converge faster in this case.
+ return 1 - RegularizedGamma.P.value(a, x, epsilon,
maxIterations);
+ } else {
+ final ContinuedFraction cf = new ContinuedFraction() {
+ /** {@inheritDoc} */
+ @Override
+ protected double getA(int n, double x) {
+ return n * (a - n);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected double getB(int n, double x) {
+ return ((2 * n) + 1) - a + x;
+ }
+ };
+
+ return Math.exp(-x + (a * Math.log(x)) -
LogGamma.value(a)) /
+ cf.evaluate(x, epsilon, maxIterations);
+ }
+ }
+ }
+ }
+
+ /**
* Assert the values are equal to the given epsilon, else throw an
AssertionError.
*
* @param x the x
