You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/11/29 16:59:58 UTC
[commons-numbers] 02/04: NUMBERS-174: Update gamma functions public API
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
commit 8987d677e472bbf1ff28c296d7fd1ce0deb695ab
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue Nov 16 16:42:16 2021 +0000
NUMBERS-174: Update gamma functions public API
Use the Boost gamma function implementations.
Add new IncompleteGamma class and tests.
Add LogGamma with a sign output.
Update the test data for Gamma and tighten the tolerances for the
improved function. The data is not duplicated from the Boost test data.
Remove obsolete tests for RegularizedGamma.
Add Gamma P and Q derivative.
---
.../org/apache/commons/numbers/gamma/Gamma.java | 89 +-
.../commons/numbers/gamma/IncompleteGamma.java | 121 +++
.../org/apache/commons/numbers/gamma/LogGamma.java | 72 +-
.../commons/numbers/gamma/RegularizedGamma.java | 181 ++--
.../apache/commons/numbers/gamma/GammaTest.java | 1000 ++++++++++----------
.../commons/numbers/gamma/IncompleteGammaTest.java | 95 ++
.../apache/commons/numbers/gamma/LogGammaTest.java | 366 ++++---
.../numbers/gamma/RegularizedGammaTest.java | 87 +-
8 files changed, 1071 insertions(+), 940 deletions(-)
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Gamma.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Gamma.java
index 5d4cf8b..37a72db 100644
--- a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Gamma.java
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Gamma.java
@@ -18,25 +18,27 @@ package org.apache.commons.numbers.gamma;
/**
- * <a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
- * function</a>.
- * <p>
- * The <a href="http://mathworld.wolfram.com/GammaFunction.html">gamma
+ * <a href="https://mathworld.wolfram.com/GammaFunction.html">Gamma
+ * function</a> \( \Gamma(x) \).
+ *
+ * <p>The <a href="https://mathworld.wolfram.com/GammaFunction.html">gamma
* function</a> can be seen to extend the factorial function to cover real and
* complex numbers, but with its argument shifted by {@code -1}. This
* implementation supports real numbers.
- * </p>
- * <p>
- * This class is immutable.
- * </p>
+ *
+ * <p>This code has been adapted from:
+ * <ul>
+ * <li>The <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
+ * <li>The <em>NSWC Library of Mathematics Subroutines</em> double
+ * precision implementation, {@code DGAMMA}.
+ * </ul>
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/tgamma.html">
+ * Boost C++ Gamma functions</a>
*/
public final class Gamma {
- /** The threshold value for choosing the Lanczos approximation. */
- private static final double LANCZOS_THRESHOLD = 20;
-
- /** √(2π). */
- private static final double SQRT_TWO_PI = 2.506628274631000502;
-
/** Private constructor. */
private Gamma() {
// intentionally empty.
@@ -44,68 +46,11 @@ public final class Gamma {
/**
* Computes the value of \( \Gamma(x) \).
- * <p>
- * Based on the <em>NSWC Library of Mathematics Subroutines</em> double
- * precision implementation, {@code DGAMMA}.
*
* @param x Argument.
* @return \( \Gamma(x) \)
*/
public static double value(final double x) {
-
- if (x == Math.rint(x) && x <= 0.0) {
- return Double.NaN;
- }
-
- final double absX = Math.abs(x);
- if (absX <= LANCZOS_THRESHOLD) {
- if (x >= 1) {
- /*
- * From the recurrence relation
- * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
- * then
- * Gamma(t) = 1 / [1 + InvGamma1pm1.value(t - 1)],
- * where t = x - n. This means that t must satisfy
- * -0.5 <= t - 1 <= 1.5.
- */
- double prod = 1;
- double t = x;
- while (t > 2.5) {
- t -= 1;
- prod *= t;
- }
- return prod / (1 + InvGamma1pm1.value(t - 1));
- }
- /*
- * From the recurrence relation
- * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
- * then
- * Gamma(x + n + 1) = 1 / [1 + InvGamma1pm1.value(x + n)],
- * which requires -0.5 <= x + n <= 1.5.
- */
- double prod = x;
- double t = x;
- while (t < -0.5) {
- t += 1;
- prod *= t;
- }
- return 1 / (prod * (1 + InvGamma1pm1.value(t)));
- }
- final double y = absX + LanczosApproximation.g() + 0.5;
- final double gammaAbs = SQRT_TWO_PI / absX *
- Math.pow(y, absX + 0.5) *
- Math.exp(-y) * LanczosApproximation.value(absX);
- if (x > 0) {
- return gammaAbs;
- }
- /*
- * From the reflection formula
- * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
- * and the recurrence relation
- * Gamma(1 - x) = -x * Gamma(-x),
- * it is found
- * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
- */
- return -Math.PI / (x * Math.sin(Math.PI * x) * gammaAbs);
+ return BoostGamma.tgamma(x);
}
}
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/IncompleteGamma.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/IncompleteGamma.java
new file mode 100644
index 0000000..e1f2fb7
--- /dev/null
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/IncompleteGamma.java
@@ -0,0 +1,121 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.gamma;
+
+import org.apache.commons.numbers.gamma.BoostGamma.Policy;
+
+/**
+ * <a href="https://mathworld.wolfram.com/IncompleteGammaFunction.html">
+ * Incomplete Gamma functions</a>.
+ *
+ * <p>By definition, the lower and upper incomplete gamma functions satisfy:
+ *
+ * <p>\[ \Gamma(a) = \gamma(a, x) + \Gamma(a, x) \]
+ *
+ * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html">
+ * Boost C++ Incomplete Gamma functions</a>
+ */
+public final class IncompleteGamma {
+ /** No instances. */
+ private IncompleteGamma() {}
+
+ /**
+ * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
+ * Lower incomplete Gamma function</a> \( \gamma(a, x) \).
+ *
+ * <p>\[ \gamma(a,x) = \int_0^x t^{a-1}\,e^{-t}\,dt \]
+ */
+ public static final class Lower {
+ /** No instances. */
+ private Lower() {}
+
+ /**
+ * Computes the lower incomplete gamma function \( \gamma(a, x) \).
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @return \( \gamma(a, x) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double value(double a,
+ double x) {
+ return BoostGamma.tgammaLower(a, x);
+ }
+
+ /**
+ * Computes the lower incomplete gamma function \( \gamma(a, x) \).
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
+ * @return \( \gamma(a, x) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double value(final double a,
+ double x,
+ double epsilon,
+ int maxIterations) {
+ return BoostGamma.tgammaLower(a, x, new Policy(epsilon, maxIterations));
+ }
+ }
+
+ /**
+ * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
+ * Upper incomplete Gamma function</a> \( \Gamma(a, x) \).
+ *
+ * <p>\[ \Gamma(a,x) = \int_x^{\infty} t^{a-1}\,e^{-t}\,dt \]
+ */
+ public static final class Upper {
+ /** No instances. */
+ private Upper() {}
+
+ /**
+ * Computes the upper incomplete gamma function \( \Gamma(a, x) \).
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @return \( \Gamma(a, x) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double value(double a,
+ double x) {
+ return BoostGamma.tgamma(a, x);
+ }
+
+ /**
+ * Computes the upper incomplete gamma function \( \Gamma(a, x) \).
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
+ * @return \( \Gamma(a, x) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double value(double a,
+ double x,
+ double epsilon,
+ int maxIterations) {
+ return BoostGamma.tgamma(a, x, new Policy(epsilon, maxIterations));
+ }
+ }
+}
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogGamma.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogGamma.java
index d7f4472..76c020f 100644
--- a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogGamma.java
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogGamma.java
@@ -17,59 +17,49 @@
package org.apache.commons.numbers.gamma;
/**
- * Function \( \ln \Gamma(x) \).
+ * Natural logarithm of the absolute value of \( \Gamma(x) \).
*
- * Class is immutable.
+ * <p>\[ \operatorname{lgamma}(z) = \ln \lvert \Gamma(x) \rvert \]
+ *
+ * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html">
+ * Boost C++ Log Gamma functions</a>
*/
public final class LogGamma {
- /** Lanczos constant. */
- private static final double LANCZOS_G = 607d / 128d;
- /** Performance. */
- private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
-
/** Private constructor. */
private LogGamma() {
// intentionally empty.
}
/**
- * Computes the function \( \ln \Gamma(x) \) for {@code x >= 0}.
- *
- * For {@code x <= 8}, the implementation is based on the double precision
- * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
- * {@code DGAMLN}. For {@code x >= 8}, the implementation is based on
- * <ul>
- * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
- * Function</a>, equation (28).</li>
- * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
- * Lanczos Approximation</a>, equations (1) through (5).</li>
- * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
- * the computation of the convergent Lanczos complex Gamma
- * approximation</a></li>
- * </ul>
+ * Computes the function \( \ln \lvert \Gamma(x) \rvert \), the natural
+ * logarithm of the absolute value of \( \Gamma(x) \).
*
* @param x Argument.
- * @return \( \ln \Gamma(x) \), or {@code NaN} if {@code x <= 0}.
+ * @return \( \ln \lvert \Gamma(x) \rvert \), or {@code NaN} if {@code x <= 0}
+ * and is an integer.
*/
public static double value(double x) {
- if (Double.isNaN(x) || x <= 0.0) {
- return Double.NaN;
- } else if (x < 0.5) {
- return LogGamma1p.value(x) - Math.log(x);
- } else if (x <= 2.5) {
- return LogGamma1p.value((x - 0.5) - 0.5);
- } else if (x <= 8.0) {
- final int n = (int) Math.floor(x - 1.5);
- double prod = 1.0;
- for (int i = 1; i <= n; i++) {
- prod *= x - i;
- }
- return LogGamma1p.value(x - (n + 1)) + Math.log(prod);
- } else {
- final double sum = LanczosApproximation.value(x);
- final double tmp = x + LANCZOS_G + .5;
- return ((x + .5) * Math.log(tmp)) - tmp +
- HALF_LOG_2_PI + Math.log(sum / x);
- }
+ return BoostGamma.lgamma(x);
+ }
+
+ /**
+ * Computes the function \( \ln \lvert \Gamma(x) \rvert \), the natural
+ * logarithm of the absolute value of \( \Gamma(x) \).
+ *
+ * <p>The sign output is set to 1 if the sign of gamma(x) is positive or zero;
+ * otherwise it is set to -1.
+ *
+ * @param x Argument.
+ * @param sign Sign output. If a non-zero length the first index {@code sign[0]} is
+ * set on output to the sign of gamma(z).
+ * @return \( \ln \lvert \Gamma(x) \rvert \), or {@code NaN} if {@code x <= 0}
+ * and is an integer.
+ */
+ public static double value(double x, int[] sign) {
+ return BoostGamma.lgamma(x, sign);
}
}
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedGamma.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedGamma.java
index f607fb9..956798f 100644
--- a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedGamma.java
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedGamma.java
@@ -16,37 +16,41 @@
*/
package org.apache.commons.numbers.gamma;
-import java.text.MessageFormat;
-
-import org.apache.commons.numbers.fraction.ContinuedFraction;
+import org.apache.commons.numbers.gamma.BoostGamma.Policy;
/**
- * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
+ * <a href="https://mathworld.wolfram.com/RegularizedGammaFunction.html">
* Regularized Gamma functions</a>.
*
- * Class is immutable.
+ * <p>By definition, the lower and upper regularized gamma functions satisfy:
+ *
+ * <p>\[ 1 = P(a, x) + Q(a, x) \]
+ *
+ * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html">
+ * Boost C++ Incomplete Gamma functions</a>
*/
public final class RegularizedGamma {
- /** Maximum allowed numerical error. */
- private static final double DEFAULT_EPSILON = 1e-15;
-
/** Private constructor. */
private RegularizedGamma() {
// intentionally empty.
}
/**
- * \( P(a, x) \) <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
- * regularized Gamma function</a>.
+ * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
+ * Lower regularized Gamma function</a> \( P(a, x) \).
*
- * Class is immutable.
+ * <p>\[ P(a,x) = 1 - Q(a,x) = \frac{\gamma(a,x)}{\Gamma(a)} = \frac{1}{\Gamma(a)} \int_0^x t^{a-1}\,e^{-t}\,dt \]
*/
public static final class P {
/** Prevent instantiation. */
private P() {}
/**
- * Computes the regularized gamma function \( P(a, x) \).
+ * Computes the lower regularized gamma function \( P(a, x) \).
*
* @param a Argument.
* @param x Argument.
@@ -55,156 +59,97 @@ public final class RegularizedGamma {
*/
public static double value(double a,
double x) {
- return value(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
+ return BoostGamma.gammaP(a, x);
}
/**
- * 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>
+ * Computes the lower regularized gamma function \( P(a, x) \).
*
* @param a Argument.
* @param x Argument.
- * @param epsilon Tolerance in continued fraction evaluation.
- * @param maxIterations Maximum number of iterations in continued fraction evaluation.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
* @return \( P(a, x) \).
- * @throws ArithmeticException if the continued fraction fails to converge.
+ * @throws ArithmeticException if the series evaluation fails to converge.
*/
public 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);
+ return BoostGamma.gammaP(a, x, new Policy(epsilon, maxIterations));
+ }
- // update partial sum
- sum += an;
- }
- if (n >= maxIterations) {
- throw new ArithmeticException(
- MessageFormat.format("Failed to converge within {0} iterations", 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;
- }
- }
+ /**
+ * Computes the derivative of the lower regularized gamma function \( P(a, x) \).
+ *
+ * <p>\[ \frac{\delta}{\delta x} P(a,x) = \frac{e^{-x} x^{a-1}}{\Gamma(a)} \]
+ *
+ * <p>This function has uses in some statistical distributions.
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @return derivative of \( P(a,x) \) with respect to x.
+ */
+ public static double derivative(double a,
+ double x) {
+ return BoostGamma.gammaPDerivative(a, x);
}
}
/**
- * Creates the \( Q(a, x) \equiv 1 - P(a, x) \) <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
- * regularized Gamma function</a>.
+ * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
+ * Upper regularized Gamma function</a> \( Q(a, x) \).
*
- * Class is immutable.
+ * <p>\[ Q(a,x) = 1 - P(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)} = \frac{1}{\Gamma(a)} \int_x^{\infty} t^{a-1}\,e^{-t}\,dt \]
*/
public static final class Q {
/** Prevent instantiation. */
private Q() {}
/**
- * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \).
+ * Computes the upper regularized gamma function \( Q(a, x) \).
*
* @param a Argument.
* @param x Argument.
* @return \( Q(a, x) \).
- * @throws ArithmeticException if the continued fraction fails to converge.
+ * @throws ArithmeticException if the series evaluation fails to converge.
*/
public static double value(double a,
double x) {
- return value(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
+ return BoostGamma.gammaQ(a, x);
}
/**
- * 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>
+ * Computes the upper regularized gamma function \( Q(a, x) \).
*
* @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.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
* @return \( Q(a, x) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
*/
public 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 BoostGamma.gammaQ(a, x, new Policy(epsilon, maxIterations));
+ }
- return Math.exp(-x + (a * Math.log(x)) - LogGamma.value(a)) /
- cf.evaluate(x, epsilon, maxIterations);
- }
+ /**
+ * Computes the derivative of the upper regularized gamma function \( Q(a, x) \).
+ *
+ * <p>\[ \frac{\delta}{\delta x} Q(a,x) = -\frac{e^{-x} x^{a-1}}{\Gamma(a)} \]
+ *
+ * <p>This function has uses in some statistical distributions.
+ *
+ * @param a Argument.
+ * @param x Argument.
+ * @return derivative of \( Q(a,x) \) with respect to x.
+ */
+ public static double derivative(double a,
+ double x) {
+ return -BoostGamma.gammaPDerivative(a, x);
}
}
}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/GammaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/GammaTest.java
index 3da8b6e..5e2a08c 100644
--- a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/GammaTest.java
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/GammaTest.java
@@ -18,15 +18,22 @@ package org.apache.commons.numbers.gamma;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvSource;
/**
* Tests for {@link Gamma}.
+ *
+ * <p>The class directly calls the methods in {@link BoostGamma}. This test ensures
+ * the arguments are passed through correctly. Accuracy of the function is tested
+ * in {@link BoostGammaTest}.
*/
class GammaTest {
/**
* Reference data for the {@link Gamma#value(double)} function. This
* data was generated with the following <a
* href="http://maxima.sourceforge.net/">Maxima</a> script.
+ * Note: This data is different from the data used for {@link BoostGammaTest}.
*
* <pre>
* kill(all);
@@ -35,506 +42,516 @@ class GammaTest {
*
* EPSILON : 10**(-fpprec + 1);
* isInteger(x) := abs(x - floor(x)) <= EPSILON * abs(x);
+ * str(x) := ssubst("e","b",string(x));
*
- * x : makelist(bfloat(i / 8), i, -160, 160);
- * x : append(x, makelist(bfloat(i / 2), i, 41, 200));
+ * x : makelist(bfloat(i / 8) + 0.0625, i, -160, 160);
+ * x : append(x, makelist(bfloat(i / 2) + 0.125, i, 41, 200));
*
* for i : 1 while i <= length(x) do if not(isInteger(x[i])) then
- * print("{", float(x[i]), ",", float(gamma(x[i])), "},");
- * </pre>
+ * printf(true, "\"~f,~a\",~%", x[i], str(gamma(x[i]))), fpprintprec:30;
+ *</pre>
*/
- private static final double[][] GAMMA_REF = {
- {-19.875, 4.920331854832504e-18},
- {-19.75, 3.879938752480031e-18},
- {-19.625, 4.323498423815027e-18},
- {-19.5, 5.811045977502237e-18},
- {-19.375, 9.14330910942125e-18},
- {-19.25, 1.735229114436739e-17},
- {-19.125, 4.653521565668223e-17},
- {-18.875, -9.779159561479603e-17},
- {-18.75, -7.662879036148061e-17},
- {-18.625, -8.484865656736991e-17},
- {-18.5, -1.133153965612936e-16},
- {-18.375, -1.771516139950367e-16},
- {-18.25, -3.340316045290721e-16},
- {-18.125, -8.899859994340475e-16},
- {-17.875, 1.845816367229275e-15},
- {-17.75, 1.436789819277761e-15},
- {-17.625, 1.580306228567265e-15},
- {-17.5, 2.096334836383932e-15},
- {-17.375, 3.255160907158799e-15},
- {-17.25, 6.096076782655566e-15},
- {-17.125, 1.613099623974211e-14},
- {-16.875, -3.29939675642233e-14},
- {-16.75, -2.550301929218027e-14},
- {-16.625, -2.785289727849803e-14},
- {-16.5, -3.66858596367188e-14},
- {-16.375, -5.655842076188414e-14},
- {-16.25, -1.051573245008085e-13},
- {-16.125, -2.762433106055837e-13},
- {-15.875, 5.567732026462681e-13},
- {-15.75, 4.271755731440195e-13},
- {-15.625, 4.630544172550298e-13},
- {-15.5, 6.053166840058604e-13},
- {-15.375, 9.261441399758529e-13},
- {-15.25, 1.708806523138138e-12},
- {-15.125, 4.454423383515037e-12},
- {-14.875, -8.838774592009505e-12},
- {-14.75, -6.728015277018307e-12},
- {-14.625, -7.235225269609841e-12},
- {-14.5, -9.382408602090835e-12},
- {-14.375, -1.423946615212874e-11},
- {-14.25, -2.605929947785661e-11},
- {-14.125, -6.737315367566492e-11},
- {-13.875, 1.314767720561414e-10},
- {-13.75, 9.923822533602004e-11},
- {-13.625, 1.058151695680439e-10},
- {-13.5, 1.360449247303171e-10},
- {-13.375, 2.046923259368506e-10},
- {-13.25, 3.713450175594567e-10},
- {-13.125, 9.516457956687671e-10},
- {-12.875, -1.8242402122789617e-9},
- {-12.75, -1.3645255983702756e-9},
- {-12.625, -1.4417316853645984e-9},
- {-12.5, -1.836606483859281e-9},
- {-12.375, -2.7377598594053765e-9},
- {-12.25, -4.9203214826628017e-9},
- {-12.125, -1.2490351068152569e-8},
- {-11.875, 2.3487092733091633e-8},
- {-11.75, 1.7397701379221012e-8},
- {-11.625, 1.8201862527728055e-8},
- {-11.5, 2.295758104824101e-8},
- {-11.375, 3.3879778260141535e-8},
- {-11.25, 6.027393816261931e-8},
- {-11.125, 1.5144550670134987e-7},
- {-10.875, -2.7890922620546316e-7},
- {-10.75, -2.044229912058469e-7},
- {-10.625, -2.1159665188483867e-7},
- {-10.5, -2.640121820547716e-7},
- {-10.375, -3.8538247770911e-7},
- {-10.25, -6.780818043294673e-7},
- {-10.125, -1.6848312620525174e-6},
- {-9.875, 3.0331378349844124e-6},
- {-9.75, 2.1975471554628537e-6},
- {-9.625, 2.2482144262764103e-6},
- {-9.5, 2.772127911575102e-6},
- {-9.375, 3.998343206232017e-6},
- {-9.25, 6.95033849437704e-6},
- {-9.125, 1.7058916528281737e-5},
- {-8.875, -2.9952236120471065e-5},
- {-8.75, -2.1426084765762826e-5},
- {-8.625, -2.163906385291045e-5},
- {-8.5, -2.633521515996347e-5},
- {-8.375, -3.748446755842515e-5},
- {-8.25, -6.429063107298763e-5},
- {-8.125, -1.5566261332057085e-4},
- {-7.875, 2.658260955691807e-4},
- {-7.75, 1.874782417004247e-4},
- {-7.625, 1.8663692573135265e-4},
- {-7.5, 2.238493288596895e-4},
- {-7.375, 3.1393241580181064e-4},
- {-7.25, 5.303977063521479e-4},
- {-7.125, .001264758733229638},
- {-6.875, -.002093380502607298},
- {-6.75, -.001452956373178292},
- {-6.625, -.001423106558701564},
- {-6.5, -.001678869966447671},
- {-6.375, -.002315251566538353},
- {-6.25, -.003845383371053072},
- {-6.125, -.009011405974261174},
- {-5.875, .01439199095542518},
- {-5.75, .009807455518953468},
- {-5.625, .009428080951397862},
- {-5.5, .01091265478190986},
- {-5.375, 0.014759728736682},
- {-5.25, 0.0240336460690817},
- {-5.125, .05519486159234969},
- {-4.875, -0.0845529468631229},
- {-4.75, -.05639286923398244},
- {-4.625, -.05303295535161297},
- {-4.5, -.06001960130050425},
- {-4.375, -.07933354195966577},
- {-4.25, -.1261766418626789},
- {-4.125, -.2828736656607921},
- {-3.875, .4121956159577241},
- {-3.75, .2678661288614166},
- {-3.625, 0.24527741850121},
- {-3.5, .2700882058522691},
- {-3.375, .3470842460735378},
- {-3.25, .5362507279163854},
- {-3.125, 1.166853870850768},
- {-2.875, -1.597258011836181},
- {-2.75, -1.004497983230312},
- {-2.625, -.8891306420668862},
- {-2.5, -.9453087204829419},
- {-2.375, -1.17140933049819},
- {-2.25, -1.742814865728253},
- {-2.125, -3.646418346408649},
- {-1.875, 4.59211678402902},
- {-1.75, 2.762369453883359},
- {-1.625, 2.333967935425576},
- {-1.5, 2.363271801207355},
- {-1.375, 2.782097159933201},
- {-1.25, 3.921333447888569},
- {-1.125, 7.748638986118379},
- {-0.875, -8.610218970054413},
- {-0.75, -4.834146544295877},
- {-0.625, -3.792697895066561},
- {-0.5, -3.544907701811032},
- {-0.375, -3.825383594908152},
- {-0.25, -4.901666809860711},
- {-0.125, -8.717218859383175},
- {0.125, 7.533941598797612},
- {0.25, 3.625609908221908},
- {0.375, 2.370436184416601},
- {0.5, 1.772453850905516},
- {0.625, 1.434518848090557},
- {0.75, 1.225416702465178},
- {0.875, 1.089652357422897},
- {1.0, 1.0},
- {1.125, .9417426998497015},
- {1.25, 0.906402477055477},
- {1.375, .8889135691562253},
- {1.5, 0.886226925452758},
- {1.625, 0.896574280056598},
- {1.75, .9190625268488832},
- {1.875, .9534458127450348},
- {2.0, 1.0},
- {2.125, 1.059460537330914},
- {2.25, 1.133003096319346},
- {2.375, 1.22225615758981},
- {2.5, 1.329340388179137},
- {2.625, 1.456933205091972},
- {2.75, 1.608359421985546},
- {2.875, 1.78771089889694},
- {3.0, 2.0},
- {3.125, 2.251353641828193},
- {3.25, 2.549256966718529},
- {3.375, 2.902858374275799},
- {3.5, 3.323350970447843},
- {3.625, 3.824449663366426},
- {3.75, 4.422988410460251},
- {3.875, 5.139668834328703},
- {4.0, 6.0},
- {4.125, 7.035480130713102},
- {4.25, 8.28508514183522},
- {4.375, 9.797147013180819},
- {4.5, 11.63172839656745},
- {4.625, 13.86363002970329},
- {4.75, 16.58620653922594},
- {4.875, 19.91621673302373},
- {5.0, 24.0},
- {5.125, 29.02135553919155},
- {5.25, 35.21161185279968},
- {5.375, 42.86251818266609},
- {5.5, 52.34277778455352},
- {5.625, 64.11928888737773},
- {5.75, 78.78448106132322},
- {5.875, 97.09155657349066},
- {6.0, 120.0},
- {6.125, 148.7344471383567},
- {6.25, 184.8609622271983},
- {6.375, 230.3860352318302},
- {6.5, 287.8852778150443},
- {6.625, 360.6709999914997},
- {6.75, 453.0107661026085},
- {6.875, 570.4128948692577},
- {7.0, 720.0},
- {7.125, 910.9984887224346},
- {7.25, 1155.38101391999},
- {7.375, 1468.710974602918},
- {7.5, 1871.254305797788},
- {7.625, 2389.445374943686},
- {7.75, 3057.822671192607},
- {7.875, 3921.588652226146},
- {8.0, 5040.0},
- {8.125, 6490.864232147346},
- {8.25, 8376.512350919926},
- {8.375, 10831.74343769652},
- {8.5, 14034.40729348341},
- {8.625, 18219.5209839456},
- {8.75, 23698.12570174271},
- {8.875, 30882.5106362809},
- {9.0, 40320.0},
- {9.125, 52738.27188619719},
- {9.25, 69106.22689508938},
- {9.375, 90715.85129070834},
- {9.5, 119292.461994609},
- {9.625, 157143.3684865308},
- {9.75, 207358.5998902487},
- {9.875, 274082.281896993},
- {10.0, 362880.0},
- {10.125, 481236.7309615494},
- {10.25, 639232.5987795768},
- {10.375, 850461.1058503906},
- {10.5, 1133278.388948786},
- {10.625, 1512504.921682859},
- {10.75, 2021746.348929925},
- {10.875, 2706562.533732806},
- {11.0, 3628800.0},
- {11.125, 4872521.900985687},
- {11.25, 6552134.137490662},
- {11.375, 8823533.973197803},
- {11.5, 1.1899423083962249e+7},
- {11.625, 1.6070364792880382e+7},
- {11.75, 2.1733773250996688e+7},
- {11.875, 2.943386755434427e+7},
- {12.0, 3.99168e+7},
- {12.125, 5.420680614846578e+7},
- {12.25, 7.371150904676994e+7},
- {12.375, 1.0036769894512501e+8},
- {12.5, 1.3684336546556586e+8},
- {12.625, 1.8681799071723443e+8},
- {12.75, 2.5537183569921107e+8},
- {12.875, 3.4952717720783816e+8},
- {13.0, 4.790016e+8},
- {13.125, 6.572575245501475e+8},
- {13.25, 9.029659858229319e+8},
- {13.375, 1.2420502744459219e+9},
- {13.5, 1.7105420683195732e+9},
- {13.625, 2.3585771328050845e+9},
- {13.75, 3.2559909051649416e+9},
- {13.875, 4.500162406550916e+9},
- {14.0, 6.2270208e+9},
- {14.125, 8.626505009720685e+9},
- {14.25, 1.196429931215385e+10},
- {14.375, 1.66124224207142e+10},
- {14.5, 2.309231792231424e+10},
- {14.625, 3.213561343446927e+10},
- {14.75, 4.476987494601794e+10},
- {14.875, 6.243975339089396e+10},
- {15.0, 8.71782912e+10},
- {15.125, 1.218493832623047e+11},
- {15.25, 1.704912651981923e+11},
- {15.375, 2.388035722977667e+11},
- {15.5, 3.348386098735565e+11},
- {15.625, 4.699833464791132e+11},
- {15.75, 6.603556554537646e+11},
- {15.875, 9.287913316895475e+11},
- {16.0, 1.307674368e+12},
- {16.125, 1.842971921842358e+12},
- {16.25, 2.599991794272433e+12},
- {16.375, 3.671604924078163e+12},
- {16.5, 5.189998453040126e+12},
- {16.625, 7.343489788736144e+12},
- {16.75, 1.040060157339679e+13},
- {16.875, 1.474456239057157e+13},
- {17.0, 2.0922789888e+13},
- {17.125, 2.971792223970803e+13},
- {17.25, 4.224986665692704e+13},
- {17.375, 6.012253063177992e+13},
- {17.5, 8.563497447516206e+13},
- {17.625, 1.220855177377384e+14},
- {17.75, 1.742100763543963e+14},
- {17.875, 2.488144903408952e+14},
- {18.0, 3.55687428096e+14},
- {18.125, 5.08919418355e+14},
- {18.25, 7.288101998319914e+14},
- {18.375, 1.044628969727176e+15},
- {18.5, 1.498612053315336e+15},
- {18.625, 2.151757250127639e+15},
- {18.75, 3.092228855290534e+15},
- {18.875, 4.447559014843502e+15},
- {19.0, 6.402373705728e+15},
- {19.125, 9.224164457684374e+15},
- {19.25, 1.330078614693384e+16},
- {19.375, 1.919505731873686e+16},
- {19.5, 2.772432298633372e+16},
- {19.625, 4.007647878362728e+16},
- {19.75, 5.797929103669752e+16},
- {19.875, 8.39476764051711e+16},
- {20.0, 1.21645100408832e+17},
- {20.5, 5.406242982335075e+17},
- {21.0, 2.43290200817664e+18},
- {21.5, 1.10827981137869e+19},
- {22.0, 5.109094217170944e+19},
- {22.5, 2.382801594464184e+20},
- {23.0, 1.124000727777608e+21},
- {23.5, 5.361303587544415e+21},
- {24.0, 2.585201673888498e+22},
- {24.5, 1.259906343072938e+23},
- {25.0, 6.204484017332395e+23},
- {25.5, 3.086770540528697e+24},
- {26.0, 1.551121004333099e+25},
- {26.5, 7.871264878348176e+25},
- {27.0, 4.032914611266056e+26},
- {27.5, 2.085885192762267e+27},
- {28.0, 1.088886945041835e+28},
- {28.5, 5.736184280096234e+28},
- {29.0, 3.048883446117139e+29},
- {29.5, 1.634812519827427e+30},
- {30.0, 8.841761993739702e+30},
- {30.5, 4.822696933490909e+31},
- {31.0, 2.65252859812191e+32},
- {31.5, 1.470922564714727e+33},
- {32.0, 8.222838654177922e+33},
- {32.5, 4.633406078851391e+34},
- {33.0, 2.631308369336935e+35},
- {33.5, 1.505856975626702e+36},
- {34.0, 8.683317618811885e+36},
- {34.5, 5.044620868349451e+37},
- {35.0, 2.952327990396041e+38},
- {35.5, 1.740394199580561e+39},
- {36.0, 1.033314796638614e+40},
- {36.5, 6.178399408510991e+40},
- {37.0, 3.719933267899013e+41},
- {37.5, 2.255115784106512e+42},
- {38.0, 1.376375309122634e+43},
- {38.5, 8.456684190399419e+43},
- {39.0, 5.230226174666011e+44},
- {39.5, 3.255823413303776e+45},
- {40.0, 2.039788208119745e+46},
- {40.5, 1.286050248254992e+47},
- {41.0, 8.159152832478976e+47},
- {41.5, 5.208503505432716e+48},
- {42.0, 3.345252661316381e+49},
- {42.5, 2.161528954754577e+50},
- {43.0, 1.40500611775288e+51},
- {43.5, 9.186498057706953e+51},
- {44.0, 6.041526306337383e+52},
- {44.5, 3.996126655102524e+53},
- {45.0, 2.658271574788449e+54},
- {45.5, 1.778276361520623e+55},
- {46.0, 1.196222208654802e+56},
- {46.5, 8.091157444918836e+56},
- {47.0, 5.502622159812088e+57},
- {47.5, 3.762388211887259e+58},
- {48.0, 2.586232415111682e+59},
- {48.5, 1.787134400646448e+60},
- {49.0, 1.241391559253607e+61},
- {49.5, 8.667601843135273e+61},
- {50.0, 6.082818640342675e+62},
- {50.5, 4.290462912351959e+63},
- {51.0, 3.041409320171338e+64},
- {51.5, 2.16668377073774e+65},
- {52.0, 1.551118753287382e+66},
- {52.5, 1.115842141929936e+67},
- {53.0, 8.065817517094388e+67},
- {53.5, 5.858171245132164e+68},
- {54.0, 4.274883284060025e+69},
- {54.5, 3.134121616145708e+70},
- {55.0, 2.308436973392413e+71},
- {55.5, 1.70809628079941e+72},
- {56.0, 1.269640335365828e+73},
- {56.5, 9.479934358436728e+73},
- {57.0, 7.109985878048635e+74},
- {57.5, 5.356162912516752e+75},
- {58.0, 4.052691950487721e+76},
- {58.5, 3.079793674697132e+77},
- {59.0, 2.350561331282878e+78},
- {59.5, 1.801679299697822e+79},
- {60.0, 1.386831185456898e+80},
- {60.5, 1.071999183320204e+81},
- {61.0, 8.320987112741391e+81},
- {61.5, 6.485595059087236e+82},
- {62.0, 5.075802138772249e+83},
- {62.5, 3.98864096133865e+84},
- {63.0, 3.146997326038794e+85},
- {63.5, 2.492900600836656e+86},
- {64.0, 1.98260831540444e+87},
- {64.5, 1.582991881531277e+88},
- {65.0, 1.268869321858841e+89},
- {65.5, 1.021029763587673e+90},
- {66.0, 8.247650592082472e+90},
- {66.5, 6.687744951499262e+91},
- {67.0, 5.443449390774431e+92},
- {67.5, 4.447350392747009e+93},
- {68.0, 3.647111091818868e+94},
- {68.5, 3.001961515104231e+95},
- {69.0, 2.48003554243683e+96},
- {69.5, 2.056343637846398e+97},
- {70.0, 1.711224524281413e+98},
- {70.5, 1.429158828303247e+99},
- {71.0, 1.19785716699699e+100},
- {71.5, 1.00755697395379e+101},
- {72.0, 8.50478588567862e+101},
- {72.5, 7.20403236376959e+102},
- {73.0, 6.12344583768861e+103},
- {73.5, 5.22292346373295e+104},
- {74.0, 4.47011546151268e+105},
- {74.5, 3.83884874584372e+106},
- {75.0, 3.30788544151939e+107},
- {75.5, 2.85994231565357e+108},
- {76.0, 2.48091408113954e+109},
- {76.5, 2.15925644831845e+110},
- {77.0, 1.88549470166605e+111},
- {77.5, 1.65183118296361e+112},
- {78.0, 1.45183092028286e+113},
- {78.5, 1.2801691667968e+114},
- {79.0, 1.13242811782063e+115},
- {79.5, 1.00493279593549e+116},
- {80.0, 8.94618213078298e+116},
- {80.5, 7.98921572768712e+117},
- {81.0, 7.15694570462638e+118},
- {81.5, 6.43131866078814e+119},
- {82.0, 5.79712602074737e+120},
- {82.5, 5.24152470854233e+121},
- {83.0, 4.75364333701284e+122},
- {83.5, 4.32425788454742e+123},
- {84.0, 3.94552396972066e+124},
- {84.5, 3.6107553335971e+125},
- {85.0, 3.31424013456535e+126},
- {85.5, 3.05108825688955e+127},
- {86.0, 2.81710411438055e+128},
- {86.5, 2.60868045964056e+129},
- {87.0, 2.42270953836727e+130},
- {87.5, 2.25650859758909e+131},
- {88.0, 2.10775729837953e+132},
- {88.5, 1.97444502289045e+133},
- {89.0, 1.85482642257398e+134},
- {89.5, 1.74738384525805e+135},
- {90.0, 1.65079551609085e+136},
- {90.5, 1.56390854150595e+137},
- {91.0, 1.48571596448176e+138},
- {91.5, 1.41533723006289e+139},
- {92.0, 1.3520015276784e+140},
- {92.5, 1.29503356550754e+141},
- {93.0, 1.24384140546413e+142},
- {93.5, 1.19790604809448e+143},
- {94.0, 1.15677250708164e+144},
- {94.5, 1.12004215496834e+145},
- {95.0, 1.08736615665674e+146},
- {95.5, 1.05843983644508e+147},
- {96.0, 1.03299784882391e+148},
- {96.5, 1.01081004380505e+149},
- {97.0, 9.9167793487095e+149},
- {97.5, 9.75431692271873e+150},
- {98.0, 9.61927596824821e+151},
- {98.5, 9.51045899965076e+152},
- {99.0, 9.42689044888325e+153},
- {99.5, 9.367802114656e+154},
- {100.0, 9.33262154439441e+155},
- };
-
- @Test
- void testGamma() {
- for (int i = 0; i < GAMMA_REF.length; i++) {
- final double[] ref = GAMMA_REF[i];
- final double x = ref[0];
- final double expected = ref[1];
- final double actual = Gamma.value(x);
- final double absX = Math.abs(x);
- final int ulps;
- if (absX <= 8.0) {
- ulps = 3;
- } else if (absX <= 20.0) {
- ulps = 5;
- } else if (absX <= 30.0) {
- ulps = 50;
- } else if (absX <= 50.0) {
- ulps = 180;
- } else {
- ulps = 500;
- }
- final double tol = ulps * Math.ulp(expected);
- Assertions.assertEquals(expected, actual, tol, Double.toString(x));
+ @ParameterizedTest
+ @CsvSource(value = {
+ "-19.9375,7.99347140318629900587250189872e-18",
+ "-19.8125,4.09142642752295504387006593351e-18",
+ "-19.6875,3.98176742212573192270008693047e-18",
+ "-19.5625,4.91270643046446517527599456012e-18",
+ "-19.4375,7.14422504112472659254903607729e-18",
+ "-19.3125,1.22455052925828737010057593176e-17",
+ "-19.1875,2.66094388924946694269102198385e-17",
+ "-19.0625,1.09937522735384800858081253897e-16",
+ "-18.9375,-1.59369836101026836429583006606e-16",
+ "-18.8125,-8.10613860952985468066756813076e-17",
+ "-18.6875,-7.83910461231003472281579614437e-17",
+ "-18.5625,-9.61048195459610999913366435823e-17",
+ "-18.4375,-1.38865874236861873142671888752e-16",
+ "-18.3125,-2.36491320963006748350673726821e-16",
+ "-18.1875,-5.10568608749741469628839843151e-16",
+ "-18.0625,-2.09568402714327276635717390241e-15",
+ "-17.9375,3.0180662711631957148852281876e-15",
+ "-17.8125,1.5249673259178039118005862546e-15",
+ "-17.6875,1.46493267442543773882620190448e-15",
+ "-17.5625,1.7839457128219029185891864465e-15",
+ "-17.4375,2.56033955624214078606801294887e-15",
+ "-17.3125,4.33074731513506107917171262241e-15",
+ "-17.1875,9.2859665716359229788745246473e-15",
+ "-17.0625,3.78532927402753643423264536122e-14",
+ "-16.9375,-5.4136563738989823135753780615e-14",
+ "-16.8125,-2.71634804929108821789479426601e-14",
+ "-16.6875,-2.59109966788999300054884461855e-14",
+ "-16.5625,-3.13305465814346700077225869666e-14",
+ "-16.4375,-4.4645921011972329957060975796e-14",
+ "-16.3125,-7.49760628932757449331602747755e-14",
+ "-16.1875,-1.59602550449992426199405892375e-13",
+ "-16.0625,-6.45871807380948404090945114759e-13",
+ "-15.9375,9.16938048329140129361829659167e-13",
+ "-15.8125,4.56686015787064206633562285972e-13",
+ "-15.6875,4.3238975707914258196658844572e-13",
+ "-15.5625,5.18912177755011722002905346634e-13",
+ "-15.4375,7.33867326634295173669189789646e-13",
+ "-15.3125,1.22304702594656058922217698228e-12",
+ "-15.1875,2.58356628540925239910288288283e-12",
+ "-15.0625,1.03743159060564837407108059058e-11",
+ "-14.9375,-1.4613700145245670811704160193e-11",
+ "-14.8125,-7.22134762463295276739320364694e-12",
+ "-14.6875,-6.78311431417904925460085624223e-12",
+ "-14.5625,-8.07557076631236992367021445699e-12",
+ "-14.4375,-1.13290768549169317435181173777e-11",
+ "-14.3125,-1.87279075848067090224645850411e-11",
+ "-14.1875,-3.9237912959653020811375033783e-11",
+ "-14.0625,-1.56263133334975786344456513956e-10",
+ "-13.9375,2.18292145919607207749830892882e-10",
+ "-13.8125,1.0696621168987561286701182902e-10",
+ "-13.6875,9.96269914895047859269500760578e-11",
+ "-13.5625,1.1760049928442388701344749803e-10",
+ "-13.4375,1.6356354709286320204704281964e-10",
+ "-13.3125,2.68043177307546022884024373401e-10",
+ "-13.1875,5.56687890115077232761383291796e-10",
+ "-13.0625,2.19745031252309699546891972751e-9",
+ "-12.9375,-3.04244678375452545801326806955e-9",
+ "-12.8125,-1.47747079896640690272560088834e-9",
+ "-12.6875,-1.36364444601259675737512916604e-9",
+ "-12.5625,-1.59495677154499896761988169203e-9",
+ "-12.4375,-2.19788516406034927750713788891e-9",
+ "-12.3125,-3.5683247979067064296435744709e-9",
+ "-12.1875,-7.34132155089258100704074216056e-9",
+ "-12.0625,-2.87041947073329545033127639406e-8",
+ "-11.9375,3.93616552648241731130466556498e-8",
+ "-11.8125,1.89300946117570884411717613819e-8",
+ "-11.6875,1.73012389087848213591969512941e-8",
+ "-11.5625,2.00366444425340495307247637561e-8",
+ "-11.4375,2.73361967280005941389950274934e-8",
+ "-11.3125,4.39349990742263229149865106729e-8",
+ "-11.1875,8.94723564015033310233090450818e-8",
+ "-11.0625,3.46244348657203763696210215033e-7",
+ "-10.9375,-4.69879759723838566536994451819e-7",
+ "-10.8125,-2.23611742601380607211341431323e-7",
+ "-10.6875,-2.0220822974642259963561436825e-7",
+ "-10.5625,-2.3167370136679994769900508093e-7",
+ "-10.4375,-3.12657750076506795464755626955e-7",
+ "-10.3125,-4.97014677027185277975784901988e-7",
+ "-10.1875,-1.00097198724181851582326994185e-6",
+ "-10.0625,-3.83032810702031663588932550381e-6",
+ "-9.9375,5.13930987197948432149837681677e-6",
+ "-9.8125,2.41780196687742781547262922619e-6",
+ "-9.6875,2.16110045541489153360562856067e-6",
+ "-9.5625,2.44705347068682444757074116733e-6",
+ "-9.4375,3.26336526642353967766338685634e-6",
+ "-9.3125,5.12546385684284817912528180175e-6",
+ "-9.1875,1.01974021200260261299495625326e-5",
+ "-9.0625,3.85426765768919361486363378821e-5",
+ "-8.9375,-5.10718918527961254448901196167e-5",
+ "-8.8125,-2.37246817999847604393251742819e-5",
+ "-8.6875,-2.09356606618317617318045266815e-5",
+ "-8.5625,-2.33999488134427587798952124126e-5",
+ "-8.4375,-3.07980097018721557079482134568e-5",
+ "-8.3125,-4.77308821668490236681041867788e-5",
+ "-8.1875,-9.36886319777391150689116057684e-5",
+ "-8.0625,-3.49293006478083171347016812056e-4",
+ "-7.9375,4.56455033434365371163705444074e-4",
+ "-7.8125,2.0907375836236570137155309836e-4",
+ "-7.6875,1.81878551999663430045051825546e-4",
+ "-7.5625,2.00362061715103622052852756283e-4",
+ "-7.4375,2.59858206859546313785813051041e-4",
+ "-7.3125,3.96762958011932509241116052598e-4",
+ "-7.1875,7.67075674317739004626713772229e-4",
+ "-7.0625,2.8161748647295455689853230472e-3",
+ "-6.9375,-3.62311182788527513361191196234e-3",
+ "-6.8125,-1.63338873720598204196525858093e-3",
+ "-6.6875,-1.39819136849741261847133590888e-3",
+ "-6.5625,-1.51523809172047114177469896939e-3",
+ "-6.4375,-1.93269541351787570878198456712e-3",
+ "-6.3125,-2.90132913046225647382566113463e-3",
+ "-6.1875,-5.5133564091587490957545052379e-3",
+ "-6.0625,-1.98892349821524155809588440209e-2",
+ "-5.9375,2.51353383059540962394326392387e-2",
+ "-5.8125,1.11274607722157526608883240826e-2",
+ "-5.6875,9.35040477682644688602705889066e-3",
+ "-5.5625,9.9437499769155918678964619866e-3",
+ "-5.4375,1.24417267245213248752840256508e-2",
+ "-5.3125,1.83146401360429939910244859123e-2",
+ "-5.1875,3.41138927816697600299810011595e-2",
+ "-5.0625,1.20578487079299019459562991877e-1",
+ "-4.9375,-1.4924107119160244642163129548e-1",
+ "-4.8125,-6.46783657385040623414133837302e-2",
+ "-4.6875,-5.31804271682004166642788974406e-2",
+ "-4.5625,-5.53121092465929797651740698005e-2",
+ "-4.4375,-6.76518890645847040093568894764e-2",
+ "-4.3125,-9.72965257227284055773175814092e-2",
+ "-4.1875,-1.76965818804911880155526443515e-1",
+ "-4.0625,-6.10428590838951286014037646375e-1",
+ "-3.9375,7.36877789008537079206804521432e-1",
+ "-3.8125,3.11264635116550800018051909202e-1",
+ "-3.6875,2.49283252350939453113807331753e-1",
+ "-3.5625,2.52361498437580470178606693465e-1",
+ "-3.4375,3.00205257724094624041521197052e-1",
+ "-3.3125,4.19591267179266249052182069827e-1",
+ "-3.1875,7.41044366245568498151266982218e-1",
+ "-3.0625,2.4798661502832395994320279384e0",
+ "-2.9375,-2.90145629422111474937679280314e0",
+ "-2.8125,-1.18669642138184992506882290383e0",
+ "-2.6875,-9.19231993044089233357164535839e-1",
+ "-2.5625,-8.99037838183880425011286345468e-1",
+ "-2.4375,-1.03195557342657527014272911486e0",
+ "-2.3125,-1.3898960725313194499853531063e0",
+ "-2.1875,-2.36207891740774958785716350582e0",
+ "-2.0625,-7.59459008524242127326058556134e0",
+ "-1.9375,8.52302786427452457629432885922e0",
+ "-1.8125,3.33758368513645291425606441703e0",
+ "-1.6875,2.47043598130598981464737969007e0",
+ "-1.5625,2.30378446034619358909142126026e0",
+ "-1.4375,2.51539171022727722097290221748e0",
+ "-1.3125,3.21413466772867622809112905833e0",
+ "-1.1875,5.16704763182945222343754516898e0",
+ "-1.0625,1.56638420508124938760999577203e1",
+ "-0.9375,-1.65133664870318913665702621647e1",
+ "-0.8125,-6.04937042930982090708911675586e0",
+ "-0.6875,-4.16886071845385781221745322699e0",
+ "-0.5625,-3.59966321929092748295534571916e0",
+ "-0.4375,-3.61587558345171100514854693763e0",
+ "-0.3125,-4.21855175139388754936960688905e0",
+ "-0.1875,-6.13586906279747451533208488817e0",
+ "-0.0625,-1.66428321789882747433562050778e1",
+ "0.0625,1.54812810815923981561596207794e1",
+ "0.1875,4.91511347381422948700990736413e0",
+ "0.3125,2.86609174393702724589949909355e0",
+ "0.4375,2.02481056085114670916238196703e0",
+ "0.5625,1.58194556776012356475248928521e0",
+ "0.6875,1.31829742231058985917800215283e0",
+ "0.8125,1.15047544927452647162476591653e0",
+ "0.9375,1.04017701118676717145976281736e0",
+ "1.0625,9.67580067599524884759976298715e-1",
+ "1.1875,9.21583776340168028814357630775e-1",
+ "1.3125,8.95653669980321014343593466736e-1",
+ "1.4375,8.85854620372376685258542110574e-1",
+ "1.5625,8.89844381865069505173275222933e-1",
+ "1.6875,9.0632947783853052818487648007e-1",
+ "1.8125,9.34761302535552758195122307182e-1",
+ "1.9375,9.75165947987594223243527641277e-1",
+ "2.0625,1.02805382182449519005747481739e0",
+ "2.1875,1.09438073440394953421704968655e0",
+ "2.3125,1.17554544184917133132596642509e0",
+ "2.4375,1.27341601678529148505915428395e0",
+ "2.5625,1.39038184666417110183324253583e0",
+ "2.6875,1.52943099385252026631197906012e0",
+ "2.8125,1.69425486084568937422865918177e0",
+ "2.9375,1.88938402422596380753433480497e0",
+ "3.0625,2.12036100751302132949354181086e0",
+ "3.1875,2.39395785650863960609979618932e0",
+ "3.3125,2.71844883427620870369129735802e0",
+ "3.4375,3.10395154091414799483168856713e0",
+ "3.5625,3.56285348207693844844768399807e0",
+ "3.6875,4.11034579597864821571344372407e0",
+ "3.8125,4.76509179612850136501810394872e0",
+ "3.9375,5.55006557116376868463210848961e0",
+ "4.0625,6.49360558550862782157397179575e0",
+ "4.1875,7.63074066762128874444310035345e0",
+ "4.3125,9.00486176353994133097742249845e0",
+ "4.4375,1.06698334218923837322339294495e1",
+ "4.5625,1.26926655298990932225948742431e1",
+ "4.6875,1.51569001226712652954433237325e1",
+ "4.8125,1.81669124727399114541315213045e1",
+ "4.9375,2.18533831864573391957389271778e1",
+ "5.0625,2.63802726911288005251442604202e1",
+ "5.1875,3.19537265456641466173554827301e1",
+ "5.3125,3.88334663552659969898401345246e1",
+ "5.4375,4.73473858096474528117880619322e1",
+ "5.5625,5.79102864801646128280891137343e1",
+ "5.6875,7.10479693250215560723905799961e1",
+ "5.8125,8.74282662750608238730079462779e1",
+ "5.9375,1.07901079483133112278960952941e2",
+ "6.0625,1.33550130498839552658542818377e2",
+ "6.1875,1.65759956455632760577531566662e2",
+ "6.3125,2.06302790012350609008525714662e2",
+ "6.4375,2.57451410339958024664097586756e2",
+ "6.5625,3.22125968545915658856245695147e2",
+ "6.6875,4.04085325536060100161721423728e2",
+ "6.8125,5.0817679772379103876185868774e2",
+ "6.9375,6.40662659431102854156330658085e2",
+ "7.0625,8.09647666149214787992415836413e2",
+ "7.1875,1.02563973056922770607347656872e3",
+ "7.3125,1.3022863619529632193663185738e3",
+ "7.4375,1.65734345406347978377512821474e3",
+ "7.5625,2.1139516685825715112441123744e3",
+ "7.6875,2.70232061452240191983151202118e3",
+ "7.8125,3.46195443449332645156516231023e3",
+ "7.9375,4.44459719980327605070954394046e3",
+ "8.0625,5.71813664217882944019643684467e3",
+ "8.1875,7.3717855634663241374031128377e3",
+ "8.3125,9.52296902178104354161620457093e3",
+ "8.4375,1.23264919395971308918275160972e4",
+ "8.5625,1.59867594936556970537835998314e4",
+ "8.6875,2.07740897241409647587047486628e4",
+ "8.8125,2.70465190194791129028528305487e4",
+ "8.9375,3.52789902734385036525070050274e4",
+ "9.0625,4.61024766775668123615837720601e4",
+ "9.1875,6.03564943008805288749879863586e4",
+ "9.3125,7.91596799935549244396847004958e4",
+ "9.4375,1.0400477574035079189979466707e5",
+ "9.5625,1.36886628164426906023022073556e5",
+ "9.6875,1.80474904478474631341247504008e5",
+ "9.8125,2.3834744885915968245639056921e5",
+ "9.9375,3.15305975568856626394281357433e5",
+ "10.0625,4.17803694890449237026852934295e5",
+ "10.1875,5.5452529138933985903895212467e5",
+ "10.3125,7.37174519939980233844563773368e5",
+ "10.4375,9.81545071049560598554312170471e5",
+ "10.5625,1.30897838182233228884514857838e6",
+ "10.6875,1.74835063713522299111833519508e6",
+ "10.8125,2.33878434193050438410333246037e6",
+ "10.9375,3.13335313221551272479317098949e6",
+ "11.0625,4.20414967983514544758270765134e6",
+ "11.1875,5.64922640602889981395932477008e6",
+ "11.3125,7.60211223688104616152206391285e6",
+ "11.4375,1.02448766790797887474106332793e7",
+ "11.5625,1.38260841579983848009268818592e7",
+ "11.6875,1.86854974343826957175772073974e7",
+ "11.8125,2.52881056971235786531172822278e7",
+ "11.9375,3.42710498836071704274253076975e7",
+ "12.0625,4.6508405833176296513883703393e7",
+ "12.1875,6.32007204174483166686699458652e7",
+ "12.3125,8.59988946797168347022183480141e7",
+ "12.4375,1.17175777016975083798509118132e8",
+ "12.5625,1.59864098076856324260717071497e8",
+ "12.6875,2.18386751264347756199183611457e8",
+ "12.8125,2.98715748547272272839947896316e8",
+ "12.9375,4.09110657985560596977389610639e8",
+ "13.0625,5.61007645362689076698722172178e8",
+ "13.1875,7.70258780087651359399414965232e8",
+ "13.3125,1.05886139074401352727106340992e9",
+ "13.4375,1.45737372664862760474395715677e9",
+ "13.5625,2.00829273209050757352525821068e9",
+ "13.6875,2.77078190666641215677714207036e9",
+ "13.8125,3.82729552826192599576183242155e9",
+ "13.9375,5.29286913768819022339497808764e9",
+ "14.0625,7.32816236755012606437705837407e9",
+ "14.1875,1.0157787662405902302079784854e10",
+ "14.3125,1.40960922642796800817960316446e10",
+ "14.4375,1.9583459451840933438746924294e10",
+ "14.5625,2.72374701789775089659363144823e10",
+ "14.6875,3.79250773474965163958871320881e10",
+ "14.8125,5.28645194841178528164603103226e10",
+ "14.9375,7.37693636065291512385675070965e10",
+ "15.0625,1.03052283293673647780302383385e11",
+ "15.1875,1.44113612460383738910756947616e11",
+ "15.3125,2.01750320532502921170705702914e11",
+ "15.4375,2.82736195835953476521908719495e11",
+ "15.5625,3.96645659481359974316447579649e11",
+ "15.6875,5.57024573541355084564592252544e11",
+ "15.8125,7.83055694858495694843818346654e11",
+ "15.9375,1.10192986887252919662610213725e12",
+ "16.0625,1.55222501711095931969080464974e12",
+ "16.1875,2.18872548924207803470712114192e12",
+ "16.3125,3.08930178315395098042643107586e12",
+ "16.4375,4.36474002321753179380696585721e12",
+ "16.5625,6.17279807567866460029971545828e12",
+ "16.6875,8.73832299743000788910704096179e12",
+ "16.8125,1.23820681749499631747178776065e13",
+ "16.9375,1.75620072851559340712285028125e13",
+ "17.0625,2.49326143373447840725335496865e13",
+ "17.1875,3.54299938571061381868215234848e13",
+ "17.3125,5.0394235337698825368206156925e13",
+ "17.4375,7.17454141316381788607020012778e13",
+ "17.5625,1.02236968128427882442464037278e14",
+ "17.6875,1.4582076501961325664947374605e14",
+ "17.8125,2.08173521191346255874944317259e14",
+ "17.9375,2.97456498392328633331432766387e14",
+ "18.0625,4.25412732130945378237603691526e14",
+ "18.1875,6.08953019419011750085994934896e14",
+ "18.3125,8.72450199283910914187069091765e14",
+ "18.4375,1.25106065892044074388349114728e15",
+ "18.5625,1.79553675275551468539577465469e15",
+ "18.6875,2.57920478128440947698756688326e15",
+ "18.8125,3.70809084622085518277244565117e15",
+ "18.9375,5.33562593991239486038257524706e15",
+ "19.0625,7.68401747411520089441671667818e15",
+ "19.1875,1.10753330406832762046890328784e16",
+ "19.3125,1.59767442743866186160507027429e16",
+ "19.4375,2.3066430898845626215351868028e16",
+ "19.5625,3.33296509730242413476590670277e16",
+ "19.6875,4.81988893502524021012051561308e16",
+ "19.8125,6.97584590445298381259066338127e16",
+ "19.9375,1.01043416237090977668495018741e17",
+ "20.0625,1.46476583100321017049818661678e17",
+ "20.625,7.86500896128685297200838843101e17",
+ "21.125,3.55029442322092508739165194691e18",
+ "21.625,1.6221580982654134254767301139e19",
+ "22.125,7.49999696905420424711486473785e19",
+ "22.625,3.5079168874989565325934288713e20",
+ "23.125,1.65937432940324268967416382325e21",
+ "23.625,7.93666195796638915499263282131e21",
+ "24.125,3.83730313674499871987150384127e22",
+ "24.625,1.87503638756955943786700950404e23",
+ "25.125,9.25749381739730941169000301705e23",
+ "25.625,4.61727710439004011574751090369e24",
+ "26.125,2.32594532162107398968711325803e25",
+ "26.625,1.18317725799994777966029966907e26",
+ "27.125,6.07653215273505579805758338661e26",
+ "27.625,3.1502094494248609633455478689e27",
+ "28.125,1.64825934642938388522311949362e28",
+ "28.625,8.70245360403617841124207598783e28",
+ "29.125,4.6357294118326421771900235758e29",
+ "29.625,2.49107734415535607021804425152e30",
+ "30.125,1.35015619119625703410659436645e31",
+ "30.625,7.37981663206024235802095609512e31",
+ "31.125,4.06734552597872431524611552894e32",
+ "31.625,2.26006884356844922214391780413e33",
+ "32.125,1.26596129496087794312035345838e34",
+ "32.625,7.14746771778522066503014005556e34",
+ "33.125,4.06690066006182039227413548505e35",
+ "33.625,2.33186134292742824196608319313e36",
+ "34.125,1.34716084364547800494080737942e37",
+ "34.625,7.84088376559347746361095473689e37",
+ "35.125,4.59718637894019369186050518228e38",
+ "35.625,2.71490600383674157177529307765e39",
+ "36.125,1.61476171560274303426600244528e40",
+ "36.625,9.67185263866839184944948158912e40",
+ "37.125,5.83332669761490921128593383356e41",
+ "37.625,3.54231602891229851486087263202e42",
+ "38.125,2.16562253648953504468990293571e43",
+ "38.625,1.3327964058782523162164033278e44",
+ "39.125,8.2564359203663523578802549424e44",
+ "39.625,5.14792611770474957138585785361e45",
+ "40.125,3.23033055384333536002064974621e46",
+ "40.625,2.03986572414050701766164617449e47",
+ "41.125,1.29617013472963831320828571067e48",
+ "41.625,8.28695450432080975925043758388e48",
+ "42.125,5.33049967907563756306907498512e49",
+ "42.625,3.44944481242353706228799464429e50",
+ "43.125,2.24547298981061232344284783748e51",
+ "43.625,1.47032585129553267280025771713e52",
+ "44.125,9.68360226855826564484728129914e52",
+ "44.625,6.41429652627676128509112429097e53",
+ "45.125,4.27288950100133471578886287325e54",
+ "45.625,2.86237982485100472347191421485e55",
+ "46.125,1.92814138732685229049972437155e56",
+ "46.625,1.30596079508827090508406086052e57",
+ "47.125,8.89355214904510618992997866379e57",
+ "47.625,6.08904220709906309495443376219e58",
+ "48.125,4.19108645023750629200450244531e59",
+ "48.625,2.89990635113092879897204907924e60",
+ "49.125,2.01696035417679990302716680181e61",
+ "49.625,1.41007946323741412850015886478e62",
+ "50.125,9.90831773989352952362095691387e62",
+ "50.625,6.99751933631566761268203836648e63",
+ "51.125,4.96654426712163167371500465308e64",
+ "51.625,3.54249416400980672892028192303e65",
+ "52.125,2.53914575656593419318679612889e66",
+ "52.625,1.82881261217006272380509554277e67",
+ "53.125,1.32352972560999319819861748218e68",
+ "53.625,9.6241263715449550840243152938e68",
+ "54.125,7.03125166730308886543015537409e69",
+ "54.625,5.1609377667409821638080390763e70",
+ "55.125,3.80566496492779684841407159623e71",
+ "55.625,2.81916225508226150698014134543e72",
+ "56.125,2.09787281191644801268825696742e73",
+ "56.625,1.5681590043895079632577036234e74",
+ "57.125,1.17743111568810644712128422296e75",
+ "57.625,8.87970036235558884194674676748e75",
+ "58.125,6.72607524836830807918033612368e76",
+ "58.625,5.11692733380740807017181282476e77",
+ "59.125,3.90953123811407907102357037189e78",
+ "59.625,2.99979864944459298113822526851e79",
+ "60.125,2.31151034453494925074268598238e80",
+ "60.625,1.78862994473133856500366681635e81",
+ "61.125,1.38979559465163823700903994691e82",
+ "61.625,1.08435690399337400503347300741e83",
+ "62.125,8.49512557230813872371775667546e83",
+ "62.625,6.68234942085916730601877740818e84",
+ "63.125,5.27759676179643118210965633463e85",
+ "63.625,4.18482132481305352539425935187e86",
+ "64.125,3.33148295588399718370672056124e87",
+ "64.625,2.66259256791230530553209751263e88",
+ "65.125,2.13631344546061319405193455989e89",
+ "65.625,1.72070044701332730370011801754e90",
+ "66.125,1.39127413135622434262632238213e91",
+ "66.625,1.12920966835249604305320244901e92",
+ "67.125,9.19980019359303346561655675184e92",
+ "67.625,7.52335941539850488684196131652e93",
+ "68.125,6.17536587994932371379511371967e94",
+ "68.625,5.0876718046632389297268763403e95",
+ "69.125,4.20696800571547678002292122153e96",
+ "69.625,3.49141477595014771552506888853e97",
+ "70.125,2.90806663395082332419084429438e98",
+ "70.625,2.43089753775529034693432921364e99",
+ "71.125,2.03928172705801485608882956143e100",
+ "71.625,1.71682138603967380752237000713e101",
+ "72.125,1.45043912837001306639318002557e102",
+ "72.625,1.22967331775091636463789751761e103",
+ "73.125,1.04612922133687192413608109344e104",
+ "73.625,8.93050247016603009818273072163e104",
+ "74.125,7.6498199310258759452450929958e105",
+ "74.625,6.5750824436597396597870354938e106",
+ "75.125,5.67042902387293054441292518313e107",
+ "75.625,4.90665527358108072111607523725e108",
+ "76.125,4.25990980418453907149021004383e109",
+ "76.625,3.71065805064569229534403189817e110",
+ "77.125,3.24285633843548036817192239586e111",
+ "77.625,2.84329173130726172130736444197e112",
+ "78.125,2.50105295101836423395259514781e113",
+ "78.625,2.20710520642726191116484164808e114",
+ "79.125,1.95394761798309705777546495923e115",
+ "79.625,1.7353364685534346776533567458e116",
+ "80.125,1.54606105272912554696483664899e117",
+ "80.625,1.38176166308567236208148530885e118",
+ "81.125,1.238781418499211844505575365e119",
+ "81.625,1.11404534086282334192819753026e120",
+ "82.125,1.00496142575748560885514801486e121",
+ "82.625,9.09339509479279552848891234073e121",
+ "83.125,8.25324570903335056272290307202e122",
+ "83.625,7.51341769707254730541396382152e123",
+ "84.125,6.86051049563397265526341317862e124",
+ "84.625,6.28309554917691768415242724575e125",
+ "85.125,5.77140445445207949624034633651e126",
+ "85.625,5.31706960849096659021399155672e127",
+ "86.125,4.91290804185233267117459481896e128",
+ "86.625,4.55274085227039014287073027044e129",
+ "87.125,4.23124205104532151304911978783e130",
+ "87.625,3.94381176327922546126177009677e131",
+ "88.125,3.68646963697323636824404561514e132",
+ "88.625,3.45576505757342131043062604729e133",
+ "89.125,3.24870136758266454951506519834e134",
+ "89.625,3.06267178227444463636914233441e135",
+ "90.125,2.89540509385804977975530185802e136",
+ "90.625,2.74491958486347100534584381722e137",
+ "91.125,2.60948384083956736400446579954e138",
+ "91.625,2.48758337378252059859467095935e139",
+ "92.125,2.37789214996505576044906945984e140",
+ "92.625,2.27924826622823449846236726651e141",
+ "93.125,2.19063314315530761931370523987e142",
+ "93.625,2.1111537065939022042007676806e143",
+ "94.125,2.04002711456338022048588800463e144",
+ "94.625,1.97656765779854093868296874096e145",
+ "95.125,1.92017552158278163253234208436e146",
+ "95.625,1.87032714619186936322875917114e147",
+ "96.125,1.82656696490562102794639040775e148",
+ "96.625,1.7885003335459750785875009574e149",
+ "97.125,1.75578749501552821311346777945e150",
+ "97.625,1.72813844728879841968517280009e151",
+ "98.125,1.70530860453383177698645558079e152",
+ "98.625,1.68709515916568945721764994609e153",
+ "99.125,1.67333406819882243116795953865e154",
+ "99.625,1.66389760072716122718090725933e155",
+ "100.125,1.65869239510208273489523989268e156",
+ })
+ void testGamma(double x, double expected) {
+ final double actual = Gamma.value(x);
+ int ulps;
+ if (Math.abs(x) <= 8.0) {
+ ulps = 2;
+ } else {
+ ulps = 4;
}
+ final double tol = ulps * Math.ulp(expected);
+ Assertions.assertEquals(expected, actual, tol, () -> Double.toString(x));
}
@Test
void testGammaNegativeInteger() {
for (int i = -100; i <= 0; i++) {
- Assertions.assertTrue(Double.isNaN(Gamma.value(i)), Integer.toString(i));
+ final int x = i;
+ Assertions.assertTrue(Double.isNaN(Gamma.value(x)), () -> Integer.toString(x));
}
}
@@ -545,9 +562,8 @@ class GammaTest {
double previousGamma = Gamma.value(-18.5);
for (double x = -19.5; x > -25; x -= 1.0) {
- double gamma = Gamma.value(x);
+ final double gamma = Gamma.value(x);
Assertions.assertEquals((int) Math.signum(previousGamma), -(int) Math.signum(gamma));
-
previousGamma = gamma;
}
}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/IncompleteGammaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/IncompleteGammaTest.java
new file mode 100644
index 0000000..e09f3df
--- /dev/null
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/IncompleteGammaTest.java
@@ -0,0 +1,95 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.gamma;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvFileSource;
+import org.junit.jupiter.params.provider.CsvSource;
+
+/**
+ * Tests for {@link IncompleteGamma}.
+ *
+ * <p>The class directly calls the methods in {@link BoostGamma}. This test ensures
+ * the arguments are passed through correctly. Accuracy of the function is tested
+ * in {@link BoostGammaTest}.
+ */
+class IncompleteGammaTest {
+ private static final double EPS = Math.ulp(1.0);
+
+ @ParameterizedTest
+ @CsvFileSource(resources = "igamma_med_data.csv")
+ void testIGammaUpper(double a, double x, double upper) {
+ TestUtils.assertEquals(upper, IncompleteGamma.Upper.value(a, x), 10);
+ }
+
+ @ParameterizedTest
+ @CsvFileSource(resources = "igamma_med_data.csv")
+ void testIGammaLower(double a, double x, double upper, double p, double lower) {
+ TestUtils.assertEquals(lower, IncompleteGamma.Lower.value(a, x), 7);
+ }
+
+ /**
+ * Test the incomplete gamma function uses the policy containing the epsilon and
+ * maximum iterations for series evaluations. The data targets each method computed
+ * using a series component to check the policy is not ignored.
+ *
+ * @see BoostGammaTest#testIGammaPolicy(double, double, double, double, double, double)
+ */
+ @ParameterizedTest
+ @CsvSource(value = {
+ // Method 2: x > 1.1, x - (1 / (3 * x)) < a
+ "5.0,2.5,21.38827245393963,0.8911780189141513,2.6117275460603704,0.10882198108584876",
+ // Method 4: a < 20, x > 1.1, x - (1 / (3 * x)) > a
+ "19.24400520324707,21.168405532836914,4.0308280447358675E15,0.3084240508178698,9.038282597080282E15,0.6915759491821302",
+ // Method 7: (x > 1000) && (a < x * 0.75f)
+ "664.0791015625,1328.158203125,Infinity,4.90100553385586E-91,Infinity,1.0",
+ // Method 2: 0.5 < x < 1.1, x * 0.75f < a
+ "0.9759566783905029,1.0735523700714111,0.33659577343416824,0.33179703084688433,0.6778671124302277,0.6682029691531157",
+ // Method 3: 0.5 < x < 1.1, x * 0.75f > a
+ "0.4912221431732178,0.9824442863464355,0.2840949896471149,0.1575143024618326,1.519518937513272,0.8424856975381674",
+ })
+ void testIGammaPolicy(double a, double x, double upper, double q, double lower, double p) {
+ // Low iterations should fail to converge
+ Assertions.assertThrows(ArithmeticException.class, () -> IncompleteGamma.Upper.value(a, x, EPS, 1), "upper");
+ Assertions.assertThrows(ArithmeticException.class, () -> IncompleteGamma.Lower.value(a, x, EPS, 1), "lower");
+
+ // Low epsilon should not be as accurate
+
+ // Innore infinite
+ if (Double.isFinite(upper)) {
+ final double u1 = IncompleteGamma.Upper.value(a, x);
+ final double u2 = IncompleteGamma.Upper.value(a, x, 1e-3, Integer.MAX_VALUE);
+ assertCloser("upper", upper, u1, u2);
+ }
+ if (Double.isFinite(lower)) {
+ final double l1 = IncompleteGamma.Lower.value(a, x);
+ final double l2 = IncompleteGamma.Lower.value(a, x, 1e-3, Integer.MAX_VALUE);
+ assertCloser("lower", lower, l1, l2);
+ }
+ }
+
+ /**
+ * Assert x is closer to the expected result than y.
+ */
+ private static void assertCloser(String msg, double expected, double x, double y) {
+ final double dx = Math.abs(expected - x);
+ final double dy = Math.abs(expected - y);
+ Assertions.assertTrue(dx < dy,
+ () -> String.format("%s %s : %s (%s) : %s (%s)", msg, expected, x, dx, y, dy));
+ }
+}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogGammaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogGammaTest.java
index 3f18410..4d4bb78 100644
--- a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogGammaTest.java
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogGammaTest.java
@@ -17,14 +17,27 @@
package org.apache.commons.numbers.gamma;
import org.junit.jupiter.api.Assertions;
-import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvSource;
+import org.junit.jupiter.params.provider.ValueSource;
/**
* Tests for {@link LogGamma}.
+ *
+ * <p>The class directly calls the methods in {@link BoostGamma}. This test ensures
+ * the arguments are passed through correctly. Accuracy of the function is tested
+ * in {@link BoostGammaTest}.
*/
class LogGammaTest {
+ @ParameterizedTest
+ @ValueSource(doubles = {0, -1, -2, Double.NaN})
+ void testLogGammaEdgeCases(double x) {
+ Assertions.assertEquals(Double.NaN, LogGamma.value(x));
+ Assertions.assertEquals(Double.NaN, LogGamma.value(x, null));
+ }
+
/**
- * Reference data for the {@link Gamma#logGamma(double)} function. This data
+ * Reference data for the {@link LogGamma#value(double)} function. This data
* was generated with the following <a
* href="http://maxima.sourceforge.net/">Maxima</a> script.
*
@@ -32,199 +45,176 @@ class LogGammaTest {
* kill(all);
*
* fpprec : 64;
- * gamln(x) := log(gamma(x));
* x : append(makelist(bfloat(i / 8), i, 1, 80),
* [0.8b0, 1b2, 1b3, 1b4, 1b5, 1b6, 1b7, 1b8, 1b9, 1b10]);
*
* for i : 1 while i <= length(x) do
- * print("{", float(x[i]), ",", float(gamln(x[i])), "},");
+ * print("\"", float(x[i]), ",", float(log_gamma(x[i])), "\",");
* </pre>
*/
- private static final double[][] LOG_GAMMA_REF = {
- {0.125, 2.019418357553796},
- {0.25, 1.288022524698077},
- {0.375, .8630739822706475},
- {0.5, .5723649429247001},
- {0.625, .3608294954889402},
- {0.75, .2032809514312954},
- {0.875, .08585870722533433},
- {0.890625, .07353860936979656},
- {0.90625, .06169536624059108},
- {0.921875, .05031670080005688},
- {0.9375, 0.0393909017345823},
- {0.953125, .02890678734595923},
- {0.96875, .01885367233441289},
- {0.984375, .009221337197578781},
- {1.0, 0.0},
- {1.015625, -0.00881970970573307},
- {1.03125, -.01724677500176807},
- {1.046875, -.02528981394675729},
- {1.0625, -.03295710029357782},
- {1.078125, -.04025658272400143},
- {1.09375, -.04719590272716985},
- {1.109375, -.05378241123619192},
- {1.125, -.06002318412603958},
- {1.25, -.09827183642181316},
- {1.375, -.1177552707410788},
- {1.5, -.1207822376352452},
- {1.625, -.1091741337567954},
- {1.75, -.08440112102048555},
- {1.875, -0.0476726853991883},
- {1.890625, -.04229320615532515},
- {1.90625, -.03674470657266143},
- {1.921875, -.03102893865389552},
- {1.9375, -.02514761940298887},
- {1.953125, -.01910243184040138},
- {1.96875, -.01289502598016741},
- {1.984375, -.006527019770560387},
- {2.0, 0.0},
- {2.015625, .006684476830232185},
- {2.03125, .01352488366498562},
- {2.046875, .02051972208453692},
- {2.0625, .02766752152285702},
- {2.078125, 0.0349668385135861},
- {2.09375, .04241625596251728},
- {2.109375, .05001438244545164},
- {2.125, .05775985153034387},
- {2.25, .1248717148923966},
- {2.375, .2006984603774558},
- {2.5, .2846828704729192},
- {2.625, .3763336820249054},
- {2.75, .4752146669149371},
- {2.875, .5809359740231859},
- {2.890625, .5946142560817441},
- {2.90625, .6083932548009232},
- {2.921875, .6222723333588501},
- {2.9375, .6362508628423761},
- {2.953125, .6503282221022278},
- {2.96875, .6645037976116387},
- {2.984375, 0.678776983328359},
- {3.0, .6931471805599453},
- {3.015625, .7076137978322324},
- {3.03125, .7221762507608962},
- {3.046875, .7368339619260166},
- {3.0625, 0.751586360749556},
- {3.078125, .7664328833756681},
- {3.09375, .7813729725537568},
- {3.109375, .7964060775242092},
- {3.125, 0.811531653906724},
- {3.25, .9358019311087253},
- {3.375, 1.06569589786406},
- {3.5, 1.200973602347074},
- {3.625, 1.341414578068493},
- {3.75, 1.486815578593417},
- {3.875, 1.6369886482725},
- {4.0, 1.791759469228055},
- {4.125, 1.950965937095089},
- {4.25, 2.114456927450371},
- {4.375, 2.282091222188554},
- {4.5, 2.453736570842442},
- {4.625, 2.62926886637513},
- {4.75, 2.808571418575736},
- {4.875, 2.99153431107781},
- {5.0, 3.178053830347946},
- {5.125, 3.368031956881733},
- {5.25, 3.561375910386697},
- {5.375, 3.757997741998131},
- {5.5, 3.957813967618717},
- {5.625, 4.160745237339519},
- {5.75, 4.366716036622286},
- {5.875, 4.57565441552762},
- {6.0, 4.787491742782046},
- {6.125, 5.002162481906205},
- {6.25, 5.219603986990229},
- {6.375, 5.439756316011858},
- {6.5, 5.662562059857142},
- {6.625, 5.887966185430003},
- {6.75, 6.115915891431546},
- {6.875, 6.346360475557843},
- {7.0, 6.579251212010101},
- {7.125, 6.814541238336996},
- {7.25, 7.05218545073854},
- {7.375, 7.292140407056348},
- {7.5, 7.534364236758733},
- {7.625, 7.778816557302289},
- {7.75, 8.025458396315983},
- {7.875, 8.274252119110479},
- {8.0, 8.525161361065415},
- {8.125, 8.77815096449171},
- {8.25, 9.033186919605123},
- {8.375, 9.290236309282232},
- {8.5, 9.549267257300997},
- {8.625, 9.810248879795765},
- {8.75, 10.07315123968124},
- {8.875, 10.33794530382217},
- {9.0, 10.60460290274525},
- {9.125, 10.87309669270751},
- {9.25, 11.14340011995171},
- {9.375, 11.41548738699336},
- {9.5, 11.68933342079727},
- {9.625, 11.96491384271319},
- {9.75, 12.24220494005076},
- {9.875, 12.52118363918365},
- {10.0, 12.80182748008147},
- {0.8, .1520596783998376},
- {100.0, 359.1342053695754},
- {1000.0, 5905.220423209181},
- {10000.0, 82099.71749644238},
- {100000.0, 1051287.708973657},
- {1000000.0, 1.2815504569147612e+7},
- {10000000.0, 1.511809493694739e+8},
- {1.e+8, 1.7420680661038346e+9},
- {1.e+9, 1.972326582750371e+10},
- {1.e+10, 2.202585092888106e+11},
- };
-
- @Test
- void testLogGammaNan() {
- testLogGamma(Double.NaN, Double.NaN);
- }
-
- @Test
- void testLogGammaNegative() {
- testLogGamma(Double.NaN, -1.0);
- }
-
- @Test
- void testLogGammaZero() {
- testLogGamma(Double.NaN, 0.0);
- }
-
- @Test
- void testLogGammaPositive() {
- testLogGamma(0.6931471805599457, 3.0);
- }
-
- @Test
- void testLogGamma() {
- final int ulps = 3;
- for (int i = 0; i < LOG_GAMMA_REF.length; i++) {
- final double[] data = LOG_GAMMA_REF[i];
- final double x = data[0];
- final double expected = data[1];
- final double actual = LogGamma.value(x);
- final double tol;
- if (expected == 0.0) {
- tol = 1E-15;
- } else {
- tol = ulps * Math.ulp(expected);
- }
- Assertions.assertEquals(expected, actual, tol, Double.toString(x));
- }
- }
-
- @Test
- void testLogGammaPrecondition1() {
- Assertions.assertTrue(Double.isNaN(LogGamma.value(0.0)));
- }
-
- @Test
- void testLogGammaPrecondition2() {
- Assertions.assertTrue(Double.isNaN(LogGamma.value(-1.0)));
+ @ParameterizedTest
+ @CsvSource(value = {
+ "0.125, 2.019418357553796",
+ "0.25, 1.288022524698077",
+ "0.375, .8630739822706475",
+ "0.5, .5723649429247001",
+ "0.625, .3608294954889402",
+ "0.75, .2032809514312954",
+ "0.875, .08585870722533433",
+ "0.890625, .07353860936979656",
+ "0.90625, .06169536624059108",
+ "0.921875, .05031670080005688",
+ "0.9375, 0.0393909017345823",
+ "0.953125, .02890678734595923",
+ "0.96875, .01885367233441289",
+ "0.984375, .009221337197578781",
+ "1.0, 0.0",
+ "1.015625, -0.00881970970573307",
+ "1.03125, -.01724677500176807",
+ "1.046875, -.02528981394675729",
+ "1.0625, -.03295710029357782",
+ "1.078125, -.04025658272400143",
+ "1.09375, -.04719590272716985",
+ "1.109375, -.05378241123619192",
+ "1.125, -.06002318412603958",
+ "1.25, -.09827183642181316",
+ "1.375, -.1177552707410788",
+ "1.5, -.1207822376352452",
+ "1.625, -.1091741337567954",
+ "1.75, -.08440112102048555",
+ "1.875, -0.0476726853991883",
+ "1.890625, -.04229320615532515",
+ "1.90625, -.03674470657266143",
+ "1.921875, -.03102893865389552",
+ "1.9375, -.02514761940298887",
+ "1.953125, -.01910243184040138",
+ "1.96875, -.01289502598016741",
+ "1.984375, -.006527019770560387",
+ "2.0, 0.0",
+ "2.015625, .006684476830232185",
+ "2.03125, .01352488366498562",
+ "2.046875, .02051972208453692",
+ "2.0625, .02766752152285702",
+ "2.078125, 0.0349668385135861",
+ "2.09375, .04241625596251728",
+ "2.109375, .05001438244545164",
+ "2.125, .05775985153034387",
+ "2.25, .1248717148923966",
+ "2.375, .2006984603774558",
+ "2.5, .2846828704729192",
+ "2.625, .3763336820249054",
+ "2.75, .4752146669149371",
+ "2.875, .5809359740231859",
+ "2.890625, .5946142560817441",
+ "2.90625, .6083932548009232",
+ "2.921875, .6222723333588501",
+ "2.9375, .6362508628423761",
+ "2.953125, .6503282221022278",
+ "2.96875, .6645037976116387",
+ "2.984375, 0.678776983328359",
+ "3.0, .6931471805599453",
+ "3.015625, .7076137978322324",
+ "3.03125, .7221762507608962",
+ "3.046875, .7368339619260166",
+ "3.0625, 0.751586360749556",
+ "3.078125, .7664328833756681",
+ "3.09375, .7813729725537568",
+ "3.109375, .7964060775242092",
+ "3.125, 0.811531653906724",
+ "3.25, .9358019311087253",
+ "3.375, 1.06569589786406",
+ "3.5, 1.200973602347074",
+ "3.625, 1.341414578068493",
+ "3.75, 1.486815578593417",
+ "3.875, 1.6369886482725",
+ "4.0, 1.791759469228055",
+ "4.125, 1.950965937095089",
+ "4.25, 2.114456927450371",
+ "4.375, 2.282091222188554",
+ "4.5, 2.453736570842442",
+ "4.625, 2.62926886637513",
+ "4.75, 2.808571418575736",
+ "4.875, 2.99153431107781",
+ "5.0, 3.178053830347946",
+ "5.125, 3.368031956881733",
+ "5.25, 3.561375910386697",
+ "5.375, 3.757997741998131",
+ "5.5, 3.957813967618717",
+ "5.625, 4.160745237339519",
+ "5.75, 4.366716036622286",
+ "5.875, 4.57565441552762",
+ "6.0, 4.787491742782046",
+ "6.125, 5.002162481906205",
+ "6.25, 5.219603986990229",
+ "6.375, 5.439756316011858",
+ "6.5, 5.662562059857142",
+ "6.625, 5.887966185430003",
+ "6.75, 6.115915891431546",
+ "6.875, 6.346360475557843",
+ "7.0, 6.579251212010101",
+ "7.125, 6.814541238336996",
+ "7.25, 7.05218545073854",
+ "7.375, 7.292140407056348",
+ "7.5, 7.534364236758733",
+ "7.625, 7.778816557302289",
+ "7.75, 8.025458396315983",
+ "7.875, 8.274252119110479",
+ "8.0, 8.525161361065415",
+ "8.125, 8.77815096449171",
+ "8.25, 9.033186919605123",
+ "8.375, 9.290236309282232",
+ "8.5, 9.549267257300997",
+ "8.625, 9.810248879795765",
+ "8.75, 10.07315123968124",
+ "8.875, 10.33794530382217",
+ "9.0, 10.60460290274525",
+ "9.125, 10.87309669270751",
+ "9.25, 11.14340011995171",
+ "9.375, 11.41548738699336",
+ "9.5, 11.68933342079727",
+ "9.625, 11.96491384271319",
+ "9.75, 12.24220494005076",
+ "9.875, 12.52118363918365",
+ "10.0, 12.80182748008147",
+ "0.8, .1520596783998376",
+ "100.0, 359.1342053695754",
+ "1000.0, 5905.220423209181",
+ "10000.0, 82099.71749644238",
+ "100000.0, 1051287.708973657",
+ "1000000.0, 1.2815504569147612e+7",
+ "10000000.0, 1.511809493694739e+8",
+ "1.e+8, 1.7420680661038346e+9",
+ "1.e+9, 1.972326582750371e+10",
+ "1.e+10, 2.202585092888106e+11",
+ })
+ void testLogGamma(double x, double expected) {
+ final double actual = LogGamma.value(x);
+ Assertions.assertEquals(expected, actual, 3 * Math.ulp(expected), () -> Double.toString(x));
+ final int[] sign = {0};
+ Assertions.assertEquals(actual, LogGamma.value(x, sign));
+ Assertions.assertEquals(1, sign[0]);
}
- private void testLogGamma(double expected, double x) {
- double actual = LogGamma.value(x);
- Assertions.assertEquals(expected, actual, 1e-15);
+ @ParameterizedTest
+ @CsvSource(value = {
+ "1, 0, 1",
+ "1, 0, 1",
+ "-0.125, 2.1653002489051702517540619481440174064962195287626, -1",
+ "-3.125, 0.1543111276840418242676072830970532952413339012367, 1",
+ "-52.0009765625, -149.43323093420259741100038126078721302600128285894, -1",
+ // Overflow close to negative poles creates signed infinity
+ "-171.99999999999997, Infinity, 1",
+ "-172.00000000000003, -Infinity, -1",
+ "-172.99999999999997, -Infinity, -1",
+ "-173.00000000000003, Infinity, 1",
+ })
+ void testLogGammaSign(double x, double expected, int expectedSign) {
+ final double actual = LogGamma.value(x);
+ Assertions.assertEquals(expected, actual, 3 * Math.ulp(expected), () -> Double.toString(x));
+ Assertions.assertEquals(actual, LogGamma.value(x, null));
+ final int[] signEmpty = {};
+ Assertions.assertEquals(actual, LogGamma.value(x, signEmpty));
+ final int[] sign = {0};
+ Assertions.assertEquals(actual, LogGamma.value(x, sign));
}
}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedGammaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedGammaTest.java
index 2261788..c8b3016 100644
--- a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedGammaTest.java
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedGammaTest.java
@@ -17,8 +17,8 @@
package org.apache.commons.numbers.gamma;
import org.junit.jupiter.api.Assertions;
-import org.junit.jupiter.api.Test;
import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvFileSource;
import org.junit.jupiter.params.provider.CsvSource;
import org.junit.jupiter.params.provider.ValueSource;
@@ -26,6 +26,8 @@ import org.junit.jupiter.params.provider.ValueSource;
* Tests for {@link RegularizedGamma}.
*/
class RegularizedGammaTest {
+ private static final double EPS = Math.ulp(1.0);
+
/**
* Test argument A cannot be NaN, negative or zero.
*
@@ -50,25 +52,6 @@ class RegularizedGammaTest {
assertRegularizedGamma(1.0, x, Double.NaN);
}
- @Test
- void testRegularizedGammaPWithACloseToZero() {
- // Creates a case where the regularized gamma P series is evaluated and the
- // result is outside the expected bounds of [0, 1]. This should be clipped to 1.0.
- final double a = 1e-18;
- // x must force use of the series in regularized gamma P using x < a + 1
- final double x = 0.5;
- assertRegularizedGamma(a, x, 1.0);
- }
-
- @Test
- void testRegularizedGammaPWithAVeryCloseToZero() {
- // Creates a case where the partial sum is infinite due to inclusion of 1 / a
- final double a = Double.MIN_VALUE;
- // x must force use of the series in regularized gamma P using x < a + 1
- final double x = 0.5;
- assertRegularizedGamma(a, x, 1.0);
- }
-
@ParameterizedTest
@CsvSource({
// Generated using matlab's gammainc function
@@ -137,16 +120,52 @@ class RegularizedGammaTest {
Assertions.assertEquals(q, actualQ, q * eps);
}
- @Test
- void testRegularizedGammaMaxIterationsExceededThrows() {
- final double a = 1.0;
- final double x = 1.0;
- // OK without
- Assertions.assertEquals(0.632120558828558, RegularizedGamma.P.value(a, x), 1e-15);
+ /**
+ * Test the incomplete gamma function uses the policy containing the epsilon and
+ * maximum iterations for series evaluations. The data targets each method computed
+ * using a series component to check the policy is not ignored.
+ *
+ * @see BoostGammaTest#testIGammaPolicy(double, double, double, double, double, double)
+ */
+ @ParameterizedTest
+ @CsvSource(value = {
+ // Method 2: x > 1.1, x - (1 / (3 * x)) < a
+ "5.0,2.5,21.38827245393963,0.8911780189141513,2.6117275460603704,0.10882198108584876",
+ // Method 4: a < 20, x > 1.1, x - (1 / (3 * x)) > a
+ "19.24400520324707,21.168405532836914,4.0308280447358675E15,0.3084240508178698,9.038282597080282E15,0.6915759491821302",
+ // Method 7: (x > 1000) && (a < x * 0.75f)
+ "664.0791015625,1328.158203125,Infinity,4.90100553385586E-91,Infinity,1.0",
+ // Method 2: 0.5 < x < 1.1, x * 0.75f < a
+ "0.9759566783905029,1.0735523700714111,0.33659577343416824,0.33179703084688433,0.6778671124302277,0.6682029691531157",
+ // Method 3: 0.5 < x < 1.1, x * 0.75f > a
+ "0.4912221431732178,0.9824442863464355,0.2840949896471149,0.1575143024618326,1.519518937513272,0.8424856975381674",
+ })
+ void testIGammaPolicy(double a, double x, double upper, double q, double lower, double p) {
+ // Low iterations should fail to converge
+ Assertions.assertThrows(ArithmeticException.class, () -> RegularizedGamma.P.value(a, x, EPS, 1), "p");
+ Assertions.assertThrows(ArithmeticException.class, () -> RegularizedGamma.Q.value(a, x, EPS, 1), "q");
+
+ // Low epsilon should not be as accurate
- final int maxIterations = 3;
- Assertions.assertThrows(ArithmeticException.class, () ->
- RegularizedGamma.P.value(a, x, 1e-15, maxIterations));
+ // Ignore 0 or 1
+ if ((int) p != p) {
+ final double p1 = RegularizedGamma.P.value(a, x);
+ final double p2 = RegularizedGamma.P.value(a, x, 1e-3, Integer.MAX_VALUE);
+ assertCloser("p", p, p1, p2);
+ }
+ if ((int) q != q) {
+ final double q1 = RegularizedGamma.Q.value(a, x);
+ final double q2 = RegularizedGamma.Q.value(a, x, 1e-3, Integer.MAX_VALUE);
+ assertCloser("q", q, q1, q2);
+ }
+ }
+
+ @ParameterizedTest
+ @CsvFileSource(resources = "igamma_med_data_p_derivative.csv")
+ void testGammaPDerivative(double a, double x, double expected) {
+ final double actual = RegularizedGamma.P.derivative(a, x);
+ TestUtils.assertEquals(expected, actual, 5);
+ Assertions.assertEquals(-actual, RegularizedGamma.Q.derivative(a, x));
}
/**
@@ -162,4 +181,14 @@ class RegularizedGammaTest {
Assertions.assertEquals(p, actualP);
Assertions.assertEquals(1 - p, actualQ);
}
+
+ /**
+ * Assert x is closer to the expected result than y.
+ */
+ private static void assertCloser(String msg, double expected, double x, double y) {
+ final double dx = Math.abs(expected - x);
+ final double dy = Math.abs(expected - y);
+ Assertions.assertTrue(dx < dy,
+ () -> String.format("%s %s : %s (%s) : %s (%s)", msg, expected, x, dx, y, dy));
+ }
}