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 99c9ae2 NUMBERS-171: InverseErf to support the full range of [0, 2]
99c9ae2 is described below
commit 99c9ae29ae49634c7c777863f92078a8fad9b41e
Author: Alex Herbert <aherb...@apache.org>
AuthorDate: Thu Oct 14 08:21:15 2021 +0100
NUMBERS-171: InverseErf to support the full range of [0, 2]
Ported the Boost C++ function boost::math::erfc_inv. This supports
evaluation from double min value, superceding the support from 2^-53 in
the previous implementation.
Added a JMH benchmark for the Boost implementation verses the previous
version.
Documented special case return values for the inverse error functions
public API.
---
LICENSE | 7 +
NOTICE | 6 +
commons-numbers-examples/examples-jmh/pom.xml | 5 +
.../numbers/examples/jmh/gamma/ErfPerformance.java | 720 +++++++++++++++++++++
.../numbers/examples/jmh/gamma/package-info.java | 22 +
commons-numbers-gamma/LICENSE | 39 ++
commons-numbers-gamma/NOTICE | 4 +
.../org/apache/commons/numbers/gamma/BoostErf.java | 372 +++++++++++
.../apache/commons/numbers/gamma/InverseErf.java | 109 +---
.../apache/commons/numbers/gamma/InverseErfc.java | 20 +-
.../apache/commons/numbers/gamma/BoostErfTest.java | 433 +++++++++++++
.../apache/commons/numbers/gamma/TestUtils.java | 256 ++++++++
.../commons/numbers/gamma/TestUtilsTest.java | 168 +++++
.../apache/commons/numbers/gamma/erf_inv_data.csv | 122 ++++
.../commons/numbers/gamma/erf_inv_limit_data.csv | 129 ++++
.../commons/numbers/gamma/erfc_inv_big_data.csv | 71 ++
.../apache/commons/numbers/gamma/erfc_inv_data.csv | 122 ++++
.../commons/numbers/gamma/erfc_inv_limit_data.csv | 126 ++++
pom.xml | 5 +
.../checkstyle/checkstyle-suppressions.xml | 4 +
src/main/resources/pmd/pmd-ruleset.xml | 13 +-
21 files changed, 2644 insertions(+), 109 deletions(-)
diff --git a/LICENSE b/LICENSE
index daa5921..f859d57 100644
--- a/LICENSE
+++ b/LICENSE
@@ -215,6 +215,13 @@ component:
* Part of Complex is licensed under the Freely Distributable Math
Library License and the Boost Software License Version 1.0
+For the Commons Numbers Gamma (commons-numbers-gamma-x.x.x.jar)
+component:
+
+ * Part of the gamma package is licensed under the
+ Boost Software License Version 1.0
+
+
Freely Distributable Math Library License
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
diff --git a/NOTICE b/NOTICE
index baedbfa..8a308cd 100644
--- a/NOTICE
+++ b/NOTICE
@@ -14,3 +14,9 @@ The org.apache.commons.numbers.complex package contains
derivative
work originating from the "Boost C++ Libraries" <boost/math/complex>.
https://www.boost.org/
Copyright 2005 John Maddock 2005.
+
+For portions of the Commons Numbers Gamma component
+The org.apache.commons.numbers.gamma package contains derivative
+work originating from the "Boost C++ Libraries" <boost/math/special_functions>.
+https://www.boost.org/
+Copyright 2006-7 John Maddock 2006-7.
diff --git a/commons-numbers-examples/examples-jmh/pom.xml
b/commons-numbers-examples/examples-jmh/pom.xml
index 2e1facc..cc6fe35 100644
--- a/commons-numbers-examples/examples-jmh/pom.xml
+++ b/commons-numbers-examples/examples-jmh/pom.xml
@@ -53,6 +53,11 @@
<dependency>
<groupId>org.apache.commons</groupId>
+ <artifactId>commons-numbers-gamma</artifactId>
+ </dependency>
+
+ <dependency>
+ <groupId>org.apache.commons</groupId>
<artifactId>commons-math3</artifactId>
</dependency>
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
new file mode 100644
index 0000000..07726c7
--- /dev/null
+++
b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/ErfPerformance.java
@@ -0,0 +1,720 @@
+/*
+ * 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.examples.jmh.gamma;
+
+import java.util.SplittableRandom;
+import java.util.concurrent.ThreadLocalRandom;
+import java.util.concurrent.TimeUnit;
+import java.util.function.DoubleSupplier;
+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.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.openjdk.jmh.annotations.Benchmark;
+import org.openjdk.jmh.annotations.BenchmarkMode;
+import org.openjdk.jmh.annotations.Fork;
+import org.openjdk.jmh.annotations.Measurement;
+import org.openjdk.jmh.annotations.Mode;
+import org.openjdk.jmh.annotations.OutputTimeUnit;
+import org.openjdk.jmh.annotations.Param;
+import org.openjdk.jmh.annotations.Scope;
+import org.openjdk.jmh.annotations.Setup;
+import org.openjdk.jmh.annotations.State;
+import org.openjdk.jmh.annotations.Warmup;
+import org.openjdk.jmh.infra.Blackhole;
+
+/**
+ * Executes a benchmark to estimate the speed of error function operations.
+ */
+@BenchmarkMode(Mode.AverageTime)
+@OutputTimeUnit(TimeUnit.NANOSECONDS)
+@Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
+@Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
+@State(Scope.Benchmark)
+@Fork(value = 1, jvmArgs = {"-server", "-Xms512M", "-Xmx512M"})
+public class ErfPerformance {
+ /** The threshold value for returning the extreme value. */
+ private static final double EXTREME_VALUE_BOUND = 40;
+ /** Commons Numbers 1.0 implementation. */
+ private static final String IMP_NUMBERS_1_0 = "Numbers 1.0";
+ /** Commons Numbers 1.1 implementation. */
+ private static final String IMP_NUMBERS_1_1 = "Boost";
+ /** Uniform numbers in the appropriate domain of the function. */
+ private static final String NUM_UNIFORM = "uniform";
+ /** Uniform numbers in the domain of the error function result, [1, 1] or
[0, 2]. */
+ private static final String NUM_INVERSE_UNIFORM = "inverse uniform";
+ /** Message prefix for an unknown parameter. */
+ private static final String UNKNOWN = "unknown parameter: ";
+ /** Message prefix for a erf domain error. */
+ private static final String ERF_DOMAIN_ERROR = "erf domain error: ";
+ /** Message prefix for a erf domain error. */
+ private static final String ERFC_DOMAIN_ERROR = "erfc domain error: ";
+
+ /**
+ * The seed for random number generation. Ensures the same numbers are
generated
+ * for each implementation of the function.
+ *
+ * <p>Note: Numbers will not be the same for the error function and
complementary error
+ * function. The domain is shifted from [-1, 1] to [0, 2]. However number
generation
+ * may wish to target specific regions of the domain where there is
limited precision
+ * as {@code |x| -> 1} or high precision as {@code |x| -> 0}.
+ */
+ private static final long SEED = ThreadLocalRandom.current().nextLong();
+
+ /**
+ * Contains an array of numbers.
+ */
+ public abstract static class NumberData {
+ /** The size of the data. */
+ @Param({"1000"})
+ private int size;
+
+ /** The numbers. */
+ private double[] numbers;
+
+ /**
+ * Gets the size of the array of numbers.
+ *
+ * @return the size
+ */
+ public int getSize() {
+ return size;
+ }
+
+ /**
+ * Gets the numbers.
+ *
+ * @return the numbers
+ */
+ public double[] getNumbers() {
+ return numbers;
+ }
+
+ /**
+ * Create the numbers.
+ */
+ @Setup
+ public void setup() {
+ numbers = createNumbers(new SplittableRandom(SEED));
+ }
+
+ /**
+ * Creates the numbers.
+ *
+ * @param rng Random number generator.
+ * @return the random numbers
+ * @see #getSize()
+ */
+ protected abstract double[] createNumbers(SplittableRandom rng);
+ }
+
+ /**
+ * Contains an array of numbers. This is used to test the JMH overhead for
calling
+ * a function on each number in the array.
+ */
+ @State(Scope.Benchmark)
+ public static class BaseData extends NumberData {
+ /** {@inheritDoc} */
+ @Override
+ protected double[] createNumbers(SplittableRandom rng) {
+ return rng.doubles().limit(getSize()).toArray();
+ }
+ }
+
+ /**
+ * Contains an array of numbers and the method to compute the error
function.
+ */
+ public abstract static class FunctionData extends NumberData {
+
+ /** The function. */
+ private DoubleUnaryOperator function;
+
+ /**
+ * The implementation of the function.
+ */
+ @Param({IMP_NUMBERS_1_0, IMP_NUMBERS_1_1})
+ private String implementation;
+
+ /**
+ * Gets the implementation.
+ *
+ * @return the implementation
+ */
+ public String getImplementation() {
+ return implementation;
+ }
+
+ /**
+ * Gets the function.
+ *
+ * @return the function
+ */
+ public DoubleUnaryOperator getFunction() {
+ return function;
+ }
+
+ /**
+ * Create the numbers and the function.
+ */
+ @Override
+ @Setup
+ public void setup() {
+ super.setup();
+ function = createFunction();
+ verify();
+ }
+
+ /**
+ * Creates the function from the implementation name.
+ *
+ * @return the inverse error function
+ * @see #getImplementation()
+ */
+ protected abstract DoubleUnaryOperator createFunction();
+
+ /**
+ * Verify the numbers for the function. This is called after the
numbers and function
+ * have been created.
+ *
+ * @see #getNumbers()
+ * @see #getFunction()
+ */
+ protected abstract void verify();
+ }
+
+ /**
+ * Contains an array of numbers in the range for the error function.
+ */
+ @State(Scope.Benchmark)
+ public static class ErfData extends FunctionData {
+ /** The type of the data. */
+ @Param({NUM_UNIFORM, NUM_INVERSE_UNIFORM})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected double[] createNumbers(SplittableRandom rng) {
+ DoubleSupplier generator;
+ if (NUM_INVERSE_UNIFORM.equals(type)) {
+ // p range: [-1, 1)
+ // The final value is generated using the inverse erf function.
+ generator = () -> InverseErf.value(makeSignedDouble(rng));
+ } else if (NUM_UNIFORM.equals(type)) {
+ // range [-6, 6)
+ // Note: Values are not distinguishable from +/-1 when |x| > 6
+ generator = () -> makeSignedDouble(rng) * 6;
+ } else {
+ throw new IllegalStateException(UNKNOWN + type);
+ }
+ return DoubleStream.generate(generator).limit(getSize()).toArray();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected DoubleUnaryOperator createFunction() {
+ final String impl = getImplementation();
+ if (IMP_NUMBERS_1_0.equals(impl)) {
+ return ErfPerformance::erf;
+ } else if (IMP_NUMBERS_1_1.equals(impl)) {
+ return Erf::value;
+ } else {
+ throw new IllegalStateException(UNKNOWN + impl);
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected void verify() {
+ final DoubleUnaryOperator function = getFunction();
+ final double relativeEps = 1e-6;
+ for (final double x : getNumbers()) {
+ final double p = function.applyAsDouble(x);
+ assert -1 <= p & p <= 1 : ERF_DOMAIN_ERROR + p;
+
+ // Implementations may not compute a round-trip
+ // to a suitable accuracy as:
+ // |p| -> 0 : x -> 0
+ // |p| -> 1 : x -> +/-big
+ if (p < 1e-10 || Math.abs(p - 1) < 1e-10) {
+ continue;
+ }
+ assertEquals(x, InverseErf.value(p), Math.abs(x) * relativeEps,
+ () -> getImplementation() + " inverse erf " + p);
+ }
+ }
+ }
+
+ /**
+ * Contains an array of numbers in the range for the complementary error
function.
+ */
+ @State(Scope.Benchmark)
+ public static class ErfcData extends FunctionData {
+ /** The type of the data. */
+ @Param({NUM_UNIFORM, NUM_INVERSE_UNIFORM})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected double[] createNumbers(SplittableRandom rng) {
+ DoubleSupplier generator;
+ if (NUM_INVERSE_UNIFORM.equals(type)) {
+ // q range: [0, 2)
+ // The final value is generated using the inverse erfc
function.
+ generator = () -> InverseErfc.value(rng.nextDouble() * 2);
+ } else if (NUM_UNIFORM.equals(type)) {
+ // range [-6, 28)
+ // Note: Values are not distinguishable from 2 when x < -6
+ // Shift the range [-17, 17) to [-6, 28)
+ generator = () -> makeSignedDouble(rng) * 17 + 11;
+ } else {
+ throw new IllegalStateException(UNKNOWN + type);
+ }
+ return DoubleStream.generate(generator).limit(getSize()).toArray();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected DoubleUnaryOperator createFunction() {
+ final String impl = getImplementation();
+ if (IMP_NUMBERS_1_0.equals(impl)) {
+ return ErfPerformance::erfc;
+ } else if (IMP_NUMBERS_1_1.equals(impl)) {
+ return Erfc::value;
+ } else {
+ throw new IllegalStateException(UNKNOWN + impl);
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected void verify() {
+ final DoubleUnaryOperator function = getFunction();
+ final double relativeEps = 1e-6;
+ for (final double x : getNumbers()) {
+ final double q = function.applyAsDouble(x);
+ assert 0 <= q && q <= 2 : ERFC_DOMAIN_ERROR + q;
+
+ // Implementations may not compute a round-trip
+ // to a suitable accuracy as:
+ // q -> 0 : x -> big
+ // |q| -> 1 : x -> 0
+ // q -> 2 : x -> -big
+ if (q < 1e-10 || Math.abs(q - 1) < 1e-10 || q > 2 - 1e-10) {
+ continue;
+ }
+ assertEquals(x, InverseErfc.value(q), Math.abs(x) *
relativeEps,
+ () -> getImplementation() + " inverse erfc " + q);
+ }
+ }
+ }
+
+ /**
+ * Contains an array of numbers in the range [-1, 1] for the inverse error
function.
+ */
+ @State(Scope.Benchmark)
+ public static class InverseErfData extends FunctionData {
+ /**
+ * The type of the data.
+ */
+ @Param({NUM_UNIFORM})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected double[] createNumbers(SplittableRandom rng) {
+ DoubleSupplier generator;
+ if (NUM_UNIFORM.equals(type)) {
+ // range [-1, 1)
+ generator = () -> makeSignedDouble(rng);
+ } else {
+ throw new IllegalStateException(UNKNOWN + type);
+ }
+ return DoubleStream.generate(generator).limit(getSize()).toArray();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected DoubleUnaryOperator createFunction() {
+ final String impl = getImplementation();
+ if (IMP_NUMBERS_1_0.equals(impl)) {
+ return ErfPerformance::inverseErf;
+ } else if (IMP_NUMBERS_1_1.equals(impl)) {
+ return InverseErf::value;
+ } else {
+ throw new IllegalStateException(UNKNOWN + impl);
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected void verify() {
+ final DoubleUnaryOperator function = getFunction();
+ final double relativeEps = 1e-12;
+ for (final double x : getNumbers()) {
+ assert -1 <= x && x <= 1 : ERF_DOMAIN_ERROR + x;
+
+ // Implementations may not compute a round-trip
+ // to a suitable accuracy as:
+ // |x| -> 0 : t -> 0
+ // |x| -> 1 : t -> +/-big
+ if (x < 1e-10 || Math.abs(x - 1) < 1e-10) {
+ continue;
+ }
+ final double t = function.applyAsDouble(x);
+ assertEquals(x, Erf.value(t), Math.abs(x) * relativeEps,
+ () -> getImplementation() + " erf " + t);
+ }
+ }
+ }
+
+ /**
+ * Contains an array of numbers in the range [0, 2] for the inverse
complementary error function.
+ */
+ @State(Scope.Benchmark)
+ public static class InverseErfcData extends FunctionData {
+ /**
+ * The type of the data.
+ */
+ @Param({NUM_UNIFORM})
+ private String type;
+
+ /** {@inheritDoc} */
+ @Override
+ protected double[] createNumbers(SplittableRandom rng) {
+ DoubleSupplier generator;
+ if (NUM_UNIFORM.equals(type)) {
+ // range [0, 2)
+ generator = () -> rng.nextDouble() * 2;
+ } else {
+ throw new IllegalStateException(UNKNOWN + type);
+ }
+ return DoubleStream.generate(generator).limit(getSize()).toArray();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected DoubleUnaryOperator createFunction() {
+ final String impl = getImplementation();
+ if (IMP_NUMBERS_1_0.equals(impl)) {
+ return ErfPerformance::inverseErfc;
+ } else if (IMP_NUMBERS_1_1.equals(impl)) {
+ return InverseErfc::value;
+ } else {
+ throw new IllegalStateException(UNKNOWN + impl);
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected void verify() {
+ final DoubleUnaryOperator function = getFunction();
+ final double relativeEps = 1e-12;
+ for (final double x : getNumbers()) {
+ assert 0 <= x && x <= 2 : ERFC_DOMAIN_ERROR + x;
+
+ // Implementations may not compute a round-trip
+ // to a suitable accuracy as:
+ // x -> 0 : t -> big
+ // |x| -> 1 : t -> 0
+ // x -> 2 : t -> -big
+ if (x < 1e-10 || Math.abs(x - 1) < 1e-10 || x > 2 - 1e-10) {
+ continue;
+ }
+ final double t = function.applyAsDouble(x);
+ assertEquals(x, Erfc.value(t), Math.abs(x) * relativeEps,
+ () -> getImplementation() + " erfc " + t);
+ }
+ }
+ }
+
+ /**
+ * Make a signed double in the range [-1, 1).
+ *
+ * @param rng Random generator
+ * @return u in [-1, 1)
+ */
+ private static double makeSignedDouble(SplittableRandom rng) {
+ // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but
using a signed
+ // shift of 10 in place of an unsigned shift of 11.
+ // Use the upper 54 bits on the assumption they are more random.
+ // The sign bit is maintained by the signed shift.
+ // The next 53 bits generates a magnitude in the range [0, 2^53) or
[-2^53, 0).
+ return (rng.nextLong() >> 10) * 0x1.0p-53;
+ }
+
+ /**
+ * Returns the inverse complementary error function.
+ *
+ * <p>This is the implementation in Commons Numbers 1.0.
+ *
+ * @param x Value in [0, 2].
+ * @return t such that {@code x = erfc(t)}
+ */
+ private static double inverseErfc(final double x) {
+ return inverseErf(1 - x);
+ }
+
+ /**
+ * Returns the inverse error function.
+ *
+ * <p>This implementation is described in the paper:
+ * <a
href="http://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf";>Approximating
+ * the erfinv function</a> by Mike Giles, Oxford-Man Institute of
Quantitative Finance,
+ * which was published in GPU Computing Gems, volume 2, 2010.
+ * The source code is available <a
href="http://gpucomputing.net/?q=node/1828";>here</a>.
+ * </p>
+ *
+ * <p>This is the implementation in Commons Numbers 1.0.
+ *
+ * @param x Value in [-1, 1].
+ * @return t such that {@code x = erf(t)}
+ */
+ private static double inverseErf(final double x) {
+ // Beware that the logarithm argument must be
+ // computed as (1 - x) * (1 + x),
+ // it must NOT be simplified as 1 - x * x as this
+ // would induce rounding errors near the boundaries +/-1
+ double w = -Math.log((1 - x) * (1 + x));
+ double p;
+
+ if (w < 6.25) {
+ w -= 3.125;
+ p = -3.6444120640178196996e-21;
+ p = -1.685059138182016589e-19 + p * w;
+ p = 1.2858480715256400167e-18 + p * w;
+ p = 1.115787767802518096e-17 + p * w;
+ p = -1.333171662854620906e-16 + p * w;
+ p = 2.0972767875968561637e-17 + p * w;
+ p = 6.6376381343583238325e-15 + p * w;
+ p = -4.0545662729752068639e-14 + p * w;
+ p = -8.1519341976054721522e-14 + p * w;
+ p = 2.6335093153082322977e-12 + p * w;
+ p = -1.2975133253453532498e-11 + p * w;
+ p = -5.4154120542946279317e-11 + p * w;
+ p = 1.051212273321532285e-09 + p * w;
+ p = -4.1126339803469836976e-09 + p * w;
+ p = -2.9070369957882005086e-08 + p * w;
+ p = 4.2347877827932403518e-07 + p * w;
+ p = -1.3654692000834678645e-06 + p * w;
+ p = -1.3882523362786468719e-05 + p * w;
+ p = 0.0001867342080340571352 + p * w;
+ p = -0.00074070253416626697512 + p * w;
+ p = -0.0060336708714301490533 + p * w;
+ p = 0.24015818242558961693 + p * w;
+ p = 1.6536545626831027356 + p * w;
+ } else if (w < 16.0) {
+ w = Math.sqrt(w) - 3.25;
+ p = 2.2137376921775787049e-09;
+ p = 9.0756561938885390979e-08 + p * w;
+ p = -2.7517406297064545428e-07 + p * w;
+ p = 1.8239629214389227755e-08 + p * w;
+ p = 1.5027403968909827627e-06 + p * w;
+ p = -4.013867526981545969e-06 + p * w;
+ p = 2.9234449089955446044e-06 + p * w;
+ p = 1.2475304481671778723e-05 + p * w;
+ p = -4.7318229009055733981e-05 + p * w;
+ p = 6.8284851459573175448e-05 + p * w;
+ p = 2.4031110387097893999e-05 + p * w;
+ p = -0.0003550375203628474796 + p * w;
+ p = 0.00095328937973738049703 + p * w;
+ p = -0.0016882755560235047313 + p * w;
+ p = 0.0024914420961078508066 + p * w;
+ p = -0.0037512085075692412107 + p * w;
+ p = 0.005370914553590063617 + p * w;
+ p = 1.0052589676941592334 + p * w;
+ p = 3.0838856104922207635 + p * w;
+ } else if (w < Double.POSITIVE_INFINITY) {
+ w = Math.sqrt(w) - 5;
+ p = -2.7109920616438573243e-11;
+ p = -2.5556418169965252055e-10 + p * w;
+ p = 1.5076572693500548083e-09 + p * w;
+ p = -3.7894654401267369937e-09 + p * w;
+ p = 7.6157012080783393804e-09 + p * w;
+ p = -1.4960026627149240478e-08 + p * w;
+ p = 2.9147953450901080826e-08 + p * w;
+ p = -6.7711997758452339498e-08 + p * w;
+ p = 2.2900482228026654717e-07 + p * w;
+ p = -9.9298272942317002539e-07 + p * w;
+ p = 4.5260625972231537039e-06 + p * w;
+ p = -1.9681778105531670567e-05 + p * w;
+ p = 7.5995277030017761139e-05 + p * w;
+ p = -0.00021503011930044477347 + p * w;
+ p = -0.00013871931833623122026 + p * w;
+ p = 1.0103004648645343977 + p * w;
+ p = 4.8499064014085844221 + p * w;
+ } else if (w == Double.POSITIVE_INFINITY) {
+ // this branch does not appears in the original code, it
+ // was added because the previous branch does not handle
+ // x = +/-1 correctly. In this case, w is positive infinity
+ // and as the first coefficient (-2.71e-11) is negative.
+ // Once the first multiplication is done, p becomes negative
+ // infinity and remains so throughout the polynomial evaluation.
+ // So the branch above incorrectly returns negative infinity
+ // instead of the correct positive infinity.
+ p = Double.POSITIVE_INFINITY;
+ } else {
+ // this branch does not appears in the original code, it
+ // occurs when the input is NaN or not in the range [-1, 1].
+ return Double.NaN;
+ }
+
+ return p * x;
+ }
+
+ /**
+ * Returns the complementary error function.
+ *
+ * <p>This implementation computes erfc(x) using the
+ * {@link RegularizedGamma.Q#value(double, double, double, int)
regularized gamma function},
+ * following <a href="http://mathworld.wolfram.com/Erf.html";>Erf</a>,
equation (3).
+ *
+ * <p>This is the implementation in Commons Numbers 1.0.
+ *
+ * @param x Value in [0, 2].
+ * @return t such that {@code x = erfc(t)}
+ */
+ private static double erfc(final double x) {
+ if (Math.abs(x) > EXTREME_VALUE_BOUND) {
+ return x > 0 ? 0 : 2;
+ }
+ final double ret = RegularizedGamma.Q.value(0.5, x * x, 1e-15, 10000);
+ return x < 0 ? 2 - ret : ret;
+ }
+
+ /**
+ * Returns the error function.
+ *
+ * <p>This implementation computes erf(x) using the
+ * {@link RegularizedGamma.P#value(double, double, double, int)
regularized gamma function},
+ * following <a href="http://mathworld.wolfram.com/Erf.html";> Erf</a>,
equation (3)
+ *
+ * <p>This is the implementation in Commons Numbers 1.0.
+ *
+ * @param x Value in [-1, 1].
+ * @return t such that {@code x = erf(t)}
+ */
+ private static double erf(final double x) {
+ if (Math.abs(x) > EXTREME_VALUE_BOUND) {
+ return x > 0 ? 1 : -1;
+ }
+ final double ret = RegularizedGamma.P.value(0.5, x * x, 1e-15, 10000);
+ return x < 0 ? -ret : ret;
+ }
+
+ /**
+ * Assert the values are equal to the given epsilon, else throw an
AssertionError.
+ *
+ * @param x the x
+ * @param y the y
+ * @param eps the max epsilon for equality
+ * @param msg the message upon failure
+ */
+ static void assertEquals(double x, double y, double eps, Supplier<String>
msg) {
+ if (!Precision.equalsIncludingNaN(x, y, eps)) {
+ throw new AssertionError(msg.get() + ": " + x + " != " + y);
+ }
+ }
+
+ /**
+ * Apply the function to all the numbers.
+ *
+ * @param numbers Numbers.
+ * @param fun Function.
+ * @param bh Data sink.
+ */
+ private static void apply(double[] numbers, DoubleUnaryOperator fun,
Blackhole bh) {
+ for (int i = 0; i < numbers.length; i++) {
+ bh.consume(fun.applyAsDouble(numbers[i]));
+ }
+ }
+
+ /**
+ * Identity function. This can be used to measure the JMH overhead of
calling a function
+ * on an array of numbers.
+ *
+ * @param z Number.
+ * @return the number
+ */
+ private static double identity(double z) {
+ return z;
+ }
+
+ // Benchmark methods.
+ // Benchmarks use function references to perform different operations on
the numbers.
+
+ /**
+ * Baseline the JMH overhead for all the benchmarks that evaluate a
function of
+ * an array of numbers. All other methods are expected to be slower than
this.
+ *
+ * @param numbers Numbers.
+ * @param bh Data sink.
+ */
+ @Benchmark
+ public void baseline(BaseData numbers, Blackhole bh) {
+ apply(numbers.getNumbers(), ErfPerformance::identity, bh);
+ }
+
+ /**
+ * Benchmark the error function.
+ *
+ * @param data Test data.
+ * @param bh Data sink.
+ */
+ @Benchmark
+ public void erf(ErfData data, Blackhole bh) {
+ apply(data.getNumbers(), data.getFunction(), bh);
+ }
+
+ /**
+ * Benchmark the complementary error function.
+ *
+ * @param data Test data.
+ * @param bh Data sink.
+ */
+ @Benchmark
+ public void erfc(ErfcData data, Blackhole bh) {
+ apply(data.getNumbers(), data.getFunction(), bh);
+ }
+
+ /**
+ * Benchmark the inverse error function.
+ *
+ * @param data Test data.
+ * @param bh Data sink.
+ */
+ @Benchmark
+ public void inverseErf(InverseErfData data, Blackhole bh) {
+ apply(data.getNumbers(), data.getFunction(), bh);
+ }
+
+ /**
+ * Benchmark the inverse complementary error function.
+ *
+ * @param data Test data.
+ * @param bh Data sink.
+ */
+ @Benchmark
+ public void inverseErfc(InverseErfcData data, Blackhole bh) {
+ apply(data.getNumbers(), data.getFunction(), bh);
+ }
+}
diff --git
a/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/package-info.java
b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/package-info.java
new file mode 100644
index 0000000..57911bd
--- /dev/null
+++
b/commons-numbers-examples/examples-jmh/src/main/java/org/apache/commons/numbers/examples/jmh/gamma/package-info.java
@@ -0,0 +1,22 @@
+/*
+ * 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.
+ */
+
+/**
+ * Benchmarks for the {@code org.apache.commons.numbers.gamma} components.
+ */
+
+package org.apache.commons.numbers.examples.jmh.gamma;
diff --git a/commons-numbers-gamma/LICENSE b/commons-numbers-gamma/LICENSE
index 261eeb9..d8ef8fb 100644
--- a/commons-numbers-gamma/LICENSE
+++ b/commons-numbers-gamma/LICENSE
@@ -199,3 +199,42 @@
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.
+
+
+
+APACHE COMMONS NUMBERS SUBCOMPONENTS:
+
+Apache Commons Numbers includes a number of subcomponents with separate
+copyright notices and license terms. Your use of these subcomponents is
+subject to the terms and conditions of the following licenses.
+
+
+For the Commons Numbers Gamma (commons-numbers-gamma-x.x.x.jar)
+component:
+
+ * Part of the gamma package is licensed under the
+ Boost Software License Version 1.0
+
+Boost Software License - Version 1.0 - August 17th, 2003
+
+Permission is hereby granted, free of charge, to any person or organization
+obtaining a copy of the software and accompanying documentation covered by
+this license (the "Software") to use, reproduce, display, distribute,
+execute, and transmit the Software, and to prepare derivative works of the
+Software, and to permit third-parties to whom the Software is furnished to
+do so, all subject to the following:
+
+The copyright notices in the Software and this entire statement, including
+the above license grant, this restriction and the following disclaimer,
+must be included in all copies of the Software, in whole or in part, and
+all derivative works of the Software, unless such copies or derivative
+works are solely in the form of machine-executable object code generated by
+a source language processor.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
+SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
+FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+DEALINGS IN THE SOFTWARE.
diff --git a/commons-numbers-gamma/NOTICE b/commons-numbers-gamma/NOTICE
index c1a3017..5bbc4dd 100644
--- a/commons-numbers-gamma/NOTICE
+++ b/commons-numbers-gamma/NOTICE
@@ -4,3 +4,7 @@ Copyright 2001-2021 The Apache Software Foundation
This product includes software developed at
The Apache Software Foundation (http://www.apache.org/).
+The org.apache.commons.numbers.gamma package contains derivative
+work originating from the "Boost C++ Libraries" <boost/math/special_functions>.
+https://www.boost.org/
+Copyright 2006-7 John Maddock 2006-7.
diff --git
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
new file mode 100644
index 0000000..08845fc
--- /dev/null
+++
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
@@ -0,0 +1,372 @@
+/*
+ * 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.
+ */
+
+// (C) Copyright John Maddock 2006.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+package org.apache.commons.numbers.gamma;
+
+/**
+ * Implementation of the inverse of the <a
href="http://mathworld.wolfram.com/Erf.html";>error function</a>.
+ *
+ * <p>This code has been adapted from the <a
href="https://www.boost.org/";>Boost</a>
+ * {@code c++} implementation {@code <boost/math/cspecial_functions/erf.hpp>}.
+ * All work is copyright John Maddock 2006 and subject to the Boost Software
License.
+ *
+ * @see
+ * <a
href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_erf/error_function.html";>
+ * Boost C++ Error functions</a>
+ */
+final class BoostErf {
+ /** Private constructor. */
+ private BoostErf() {
+ // intentionally empty.
+ }
+
+ // Code ported from Boost 1.77.0
+ //
+ // boost/math/special_functions/erf.hpp
+ // boost/math/special_functions/detail/erf_inv.hpp
+ //
+ // Original code comments, including measured deviations, are preserved.
+ //
+ // Changes to the Boost implementation:
+ // - Update method names to replace underscores with camel case
+ // - Explicitly inline the polynomial function evaluation
+ // using Horner's method
(https://en.wikipedia.org/wiki/Horner%27s_method)
+ // - Support odd function for f(0.0) = -f(-0.0)
+ // Inverse erf:
+ // - Change inverse erf edge case detection to include NaN
+ //
+ // Note:
+ // Constants using the 'f' suffix are machine
+ // representable as a float, e.g.
+ // assert 0.0891314744949340820313f == 0.0891314744949340820313;
+ // The values are unchanged from the Boost reference.
+
+ /**
+ * Returns the inverse complementary error function.
+ *
+ * @param z Value (in {@code [0, 2]}).
+ * @return t such that {@code z = erfc(t)}
+ */
+ static double erfcInv(double z) {
+ //
+ // Begin by testing for domain errors, and other special cases:
+ //
+ if (z < 0 || z > 2 || Double.isNaN(z)) {
+ // Argument outside range [0,2] in inverse erfc function
+ return Double.NaN;
+ }
+ // Domain bounds must be detected as the implementation computes NaN.
+ // (log(q=0) creates infinity and the rational number is
+ // infinity / infinity)
+ if (z == (int) z) {
+ // z return
+ // 2 -inf
+ // 1 0
+ // 0 inf
+ return z == 1 ? 0 : (1 - z) * Double.POSITIVE_INFINITY;
+ }
+
+ //
+ // Normalise the input, so it's in the range [0,1], we will
+ // negate the result if z is outside that range. This is a simple
+ // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
+ //
+ double p;
+ double q;
+ double s;
+ if (z > 1) {
+ q = 2 - z;
+ p = 1 - q;
+ s = -1;
+ } else {
+ p = 1 - z;
+ q = z;
+ s = 1;
+ }
+
+ //
+ // And get the result, negating where required:
+ //
+ return s * erfInvImp(p, q);
+ }
+
+ /**
+ * Returns the inverse error function.
+ *
+ * @param z Value (in {@code [-1, 1]}).
+ * @return t such that {@code z = erf(t)}
+ */
+ static double erfInv(double z) {
+ //
+ // Begin by testing for domain errors, and other special cases:
+ //
+ if (z < -1 || z > 1 || Double.isNaN(z)) {
+ // Argument outside range [-1, 1] in inverse erf function
+ return Double.NaN;
+ }
+ // Domain bounds must be detected as the implementation computes NaN.
+ // (log(q=0) creates infinity and the rational number is
+ // infinity / infinity)
+ if (z == (int) z) {
+ // z return
+ // -1 -inf
+ // -0 -0
+ // 0 0
+ // 1 inf
+ return z == 0 ? z : z * Double.POSITIVE_INFINITY;
+ }
+
+ //
+ // Normalise the input, so it's in the range [0,1], we will
+ // negate the result if z is outside that range. This is a simple
+ // application of the erf reflection formula: erf(-z) = -erf(z)
+ //
+ double p;
+ double q;
+ double s;
+ if (z < 0) {
+ p = -z;
+ q = 1 - p;
+ s = -1;
+ } else {
+ p = z;
+ q = 1 - z;
+ s = 1;
+ }
+ //
+ // And get the result, negating where required:
+ //
+ return s * erfInvImp(p, q);
+ }
+
+ /**
+ * Common implementation for inverse erf and erfc functions.
+ *
+ * @param p P-value
+ * @param q Q-value (1-p)
+ * @return the inverse
+ */
+ private static double erfInvImp(double p, double q) {
+ double result = 0;
+
+ if (p <= 0.5) {
+ //
+ // Evaluate inverse erf using the rational approximation:
+ //
+ // x = p(p+10)(Y+R(p))
+ //
+ // Where Y is a constant, and R(p) is optimised for a low
+ // absolute error compared to |Y|.
+ //
+ // double: Max error found: 2.001849e-18
+ // long double: Max error found: 1.017064e-20
+ // Maximum Deviation Found (actual error term at infinite
precision) 8.030e-21
+ //
+ final float Y = 0.0891314744949340820313f;
+ double P;
+ P = -0.00538772965071242932965;
+ P = 0.00822687874676915743155 + P * p;
+ P = 0.0219878681111168899165 + P * p;
+ P = -0.0365637971411762664006 + P * p;
+ P = -0.0126926147662974029034 + P * p;
+ P = 0.0334806625409744615033 + P * p;
+ P = -0.00836874819741736770379 + P * p;
+ P = -0.000508781949658280665617 + P * p;
+ double Q;
+ Q = 0.000886216390456424707504;
+ Q = -0.00233393759374190016776 + Q * p;
+ Q = 0.0795283687341571680018 + Q * p;
+ Q = -0.0527396382340099713954 + Q * p;
+ Q = -0.71228902341542847553 + Q * p;
+ Q = 0.662328840472002992063 + Q * p;
+ Q = 1.56221558398423026363 + Q * p;
+ Q = -1.56574558234175846809 + Q * p;
+ Q = -0.970005043303290640362 + Q * p;
+ Q = 1.0 + Q * p;
+ final double g = p * (p + 10);
+ final double r = P / Q;
+ result = g * Y + g * r;
+ } else if (q >= 0.25) {
+ //
+ // Rational approximation for 0.5 > q >= 0.25
+ //
+ // x = sqrt(-2*log(q)) / (Y + R(q))
+ //
+ // Where Y is a constant, and R(q) is optimised for a low
+ // absolute error compared to Y.
+ //
+ // double : Max error found: 7.403372e-17
+ // long double : Max error found: 6.084616e-20
+ // Maximum Deviation Found (error term) 4.811e-20
+ //
+ final float Y = 2.249481201171875f;
+ final double xs = q - 0.25f;
+ double P;
+ P = -3.67192254707729348546;
+ P = 21.1294655448340526258 + P * xs;
+ P = 17.445385985570866523 + P * xs;
+ P = -44.6382324441786960818 + P * xs;
+ P = -18.8510648058714251895 + P * xs;
+ P = 17.6447298408374015486 + P * xs;
+ P = 8.37050328343119927838 + P * xs;
+ P = 0.105264680699391713268 + P * xs;
+ P = -0.202433508355938759655 + P * xs;
+ double Q;
+ Q = 1.72114765761200282724;
+ Q = -22.6436933413139721736 + Q * xs;
+ Q = 10.8268667355460159008 + Q * xs;
+ Q = 48.5609213108739935468 + Q * xs;
+ Q = -20.1432634680485188801 + Q * xs;
+ Q = -28.6608180499800029974 + Q * xs;
+ Q = 3.9713437953343869095 + Q * xs;
+ Q = 6.24264124854247537712 + Q * xs;
+ Q = 1.0 + Q * xs;
+ final double g = Math.sqrt(-2 * Math.log(q));
+ final double r = P / Q;
+ result = g / (Y + r);
+ } else {
+ //
+ // For q < 0.25 we have a series of rational approximations all
+ // of the general form:
+ //
+ // let: x = sqrt(-log(q))
+ //
+ // Then the result is given by:
+ //
+ // x(Y+R(x-B))
+ //
+ // where Y is a constant, B is the lowest value of x for which
+ // the approximation is valid, and R(x-B) is optimised for a low
+ // absolute error compared to Y.
+ //
+ // Note that almost all code will really go through the first
+ // or maybe second approximation. After than we're dealing with
very
+ // small input values indeed.
+ //
+ // Limit for a double: Math.sqrt(-Math.log(Double.MIN_VALUE)) =
27.28...
+ // Branches for x >= 44 (supporting 80 and 128 bit long double)
have been removed.
+ final double x = Math.sqrt(-Math.log(q));
+ if (x < 3) {
+ // Max error found: 1.089051e-20
+ final float Y = 0.807220458984375f;
+ final double xs = x - 1.125f;
+ double P;
+ P = -0.681149956853776992068e-9;
+ P = 0.285225331782217055858e-7 + P * xs;
+ P = -0.679465575181126350155e-6 + P * xs;
+ P = 0.00214558995388805277169 + P * xs;
+ P = 0.0290157910005329060432 + P * xs;
+ P = 0.142869534408157156766 + P * xs;
+ P = 0.337785538912035898924 + P * xs;
+ P = 0.387079738972604337464 + P * xs;
+ P = 0.117030156341995252019 + P * xs;
+ P = -0.163794047193317060787 + P * xs;
+ P = -0.131102781679951906451 + P * xs;
+ double Q;
+ Q = 0.01105924229346489121;
+ Q = 0.152264338295331783612 + Q * xs;
+ Q = 0.848854343457902036425 + Q * xs;
+ Q = 2.59301921623620271374 + Q * xs;
+ Q = 4.77846592945843778382 + Q * xs;
+ Q = 5.38168345707006855425 + Q * xs;
+ Q = 3.46625407242567245975 + Q * xs;
+ Q = 1.0 + Q * xs;
+ final double R = P / Q;
+ result = Y * x + R * x;
+ } else if (x < 6) {
+ // Max error found: 8.389174e-21
+ final float Y = 0.93995571136474609375f;
+ final double xs = x - 3;
+ double P;
+ P = 0.266339227425782031962e-11;
+ P = -0.230404776911882601748e-9 + P * xs;
+ P = 0.460469890584317994083e-5 + P * xs;
+ P = 0.000157544617424960554631 + P * xs;
+ P = 0.00187123492819559223345 + P * xs;
+ P = 0.00950804701325919603619 + P * xs;
+ P = 0.0185573306514231072324 + P * xs;
+ P = -0.00222426529213447927281 + P * xs;
+ P = -0.0350353787183177984712 + P * xs;
+ double Q;
+ Q = 0.764675292302794483503e-4;
+ Q = 0.00263861676657015992959 + Q * xs;
+ Q = 0.0341589143670947727934 + Q * xs;
+ Q = 0.220091105764131249824 + Q * xs;
+ Q = 0.762059164553623404043 + Q * xs;
+ Q = 1.3653349817554063097 + Q * xs;
+ Q = 1.0 + Q * xs;
+ final double R = P / Q;
+ result = Y * x + R * x;
+ } else if (x < 18) {
+ // Max error found: 1.481312e-19
+ final float Y = 0.98362827301025390625f;
+ final double xs = x - 6;
+ double P;
+ P = 0.99055709973310326855e-16;
+ P = -0.281128735628831791805e-13 + P * xs;
+ P = 0.462596163522878599135e-8 + P * xs;
+ P = 0.449696789927706453732e-6 + P * xs;
+ P = 0.149624783758342370182e-4 + P * xs;
+ P = 0.000209386317487588078668 + P * xs;
+ P = 0.00105628862152492910091 + P * xs;
+ P = -0.00112951438745580278863 + P * xs;
+ P = -0.0167431005076633737133 + P * xs;
+ double Q;
+ Q = 0.282243172016108031869e-6;
+ Q = 0.275335474764726041141e-4 + Q * xs;
+ Q = 0.000964011807005165528527 + Q * xs;
+ Q = 0.0160746087093676504695 + Q * xs;
+ Q = 0.138151865749083321638 + Q * xs;
+ Q = 0.591429344886417493481 + Q * xs;
+ Q = 1.0 + Q * xs;
+ final double R = P / Q;
+ result = Y * x + R * x;
+ } else {
+ // x < 44
+ // Max error found: 5.697761e-20
+ final float Y = 0.99714565277099609375f;
+ final double xs = x - 18;
+ double P;
+ P = -0.116765012397184275695e-17;
+ P = 0.145596286718675035587e-11 + P * xs;
+ P = 0.411632831190944208473e-9 + P * xs;
+ P = 0.396341011304801168516e-7 + P * xs;
+ P = 0.162397777342510920873e-5 + P * xs;
+ P = 0.254723037413027451751e-4 + P * xs;
+ P = -0.779190719229053954292e-5 + P * xs;
+ P = -0.0024978212791898131227 + P * xs;
+ double Q;
+ Q = 0.509761276599778486139e-9;
+ Q = 0.144437756628144157666e-6 + Q * xs;
+ Q = 0.145007359818232637924e-4 + Q * xs;
+ Q = 0.000690538265622684595676 + Q * xs;
+ Q = 0.0169410838120975906478 + Q * xs;
+ Q = 0.207123112214422517181 + Q * xs;
+ Q = 1.0 + Q * xs;
+ final double R = P / Q;
+ result = Y * x + R * x;
+ }
+ }
+ return result;
+ }
+}
+
diff --git
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErf.java
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErf.java
index 802262a..3ffd2f4 100644
---
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErf.java
+++
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErf.java
@@ -18,13 +18,6 @@ package org.apache.commons.numbers.gamma;
/**
* Inverse of the <a href="http://mathworld.wolfram.com/Erf.html";>error
function</a>.
- * <p>
- * This implementation is described in the paper:
- * <a
href="http://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf";>Approximating
- * the erfinv function</a> by Mike Giles, Oxford-Man Institute of Quantitative
Finance,
- * which was published in GPU Computing Gems, volume 2, 2010.
- * The source code is available <a
href="http://gpucomputing.net/?q=node/1828";>here</a>.
- * </p>
*/
public final class InverseErf {
/** Private constructor. */
@@ -35,99 +28,19 @@ public final class InverseErf {
/**
* Returns the inverse error function.
*
- * @param x Value.
+ * <p>Special cases:
+ * <ul>
+ * <li>If the argument is 0, then the result is 0.
+ * <li>If the argument is 1, then the result is positive infinity.
+ * <li>If the argument is -1, then the result is negative infinity.
+ * <li>If the argument is outside the interval {@code [-1, 1]}, then the
result is nan.
+ * <li>If the argument is nan, then the result is nan.
+ * </ul>
+ *
+ * @param x Value (in {@code [-1, 1]})
* @return t such that {@code x =} {@link Erf#value(double) Erf.value(t)}.
*/
public static double value(final double x) {
- // Beware that the logarithm argument must be
- // computed as (1 - x) * (1 + x),
- // it must NOT be simplified as 1 - x * x as this
- // would induce rounding errors near the boundaries +/-1
- double w = -Math.log((1 - x) * (1 + x));
- double p;
-
- if (w < 6.25) {
- w -= 3.125;
- p = -3.6444120640178196996e-21;
- p = -1.685059138182016589e-19 + p * w;
- p = 1.2858480715256400167e-18 + p * w;
- p = 1.115787767802518096e-17 + p * w;
- p = -1.333171662854620906e-16 + p * w;
- p = 2.0972767875968561637e-17 + p * w;
- p = 6.6376381343583238325e-15 + p * w;
- p = -4.0545662729752068639e-14 + p * w;
- p = -8.1519341976054721522e-14 + p * w;
- p = 2.6335093153082322977e-12 + p * w;
- p = -1.2975133253453532498e-11 + p * w;
- p = -5.4154120542946279317e-11 + p * w;
- p = 1.051212273321532285e-09 + p * w;
- p = -4.1126339803469836976e-09 + p * w;
- p = -2.9070369957882005086e-08 + p * w;
- p = 4.2347877827932403518e-07 + p * w;
- p = -1.3654692000834678645e-06 + p * w;
- p = -1.3882523362786468719e-05 + p * w;
- p = 0.0001867342080340571352 + p * w;
- p = -0.00074070253416626697512 + p * w;
- p = -0.0060336708714301490533 + p * w;
- p = 0.24015818242558961693 + p * w;
- p = 1.6536545626831027356 + p * w;
- } else if (w < 16.0) {
- w = Math.sqrt(w) - 3.25;
- p = 2.2137376921775787049e-09;
- p = 9.0756561938885390979e-08 + p * w;
- p = -2.7517406297064545428e-07 + p * w;
- p = 1.8239629214389227755e-08 + p * w;
- p = 1.5027403968909827627e-06 + p * w;
- p = -4.013867526981545969e-06 + p * w;
- p = 2.9234449089955446044e-06 + p * w;
- p = 1.2475304481671778723e-05 + p * w;
- p = -4.7318229009055733981e-05 + p * w;
- p = 6.8284851459573175448e-05 + p * w;
- p = 2.4031110387097893999e-05 + p * w;
- p = -0.0003550375203628474796 + p * w;
- p = 0.00095328937973738049703 + p * w;
- p = -0.0016882755560235047313 + p * w;
- p = 0.0024914420961078508066 + p * w;
- p = -0.0037512085075692412107 + p * w;
- p = 0.005370914553590063617 + p * w;
- p = 1.0052589676941592334 + p * w;
- p = 3.0838856104922207635 + p * w;
- } else if (w < Double.POSITIVE_INFINITY) {
- w = Math.sqrt(w) - 5;
- p = -2.7109920616438573243e-11;
- p = -2.5556418169965252055e-10 + p * w;
- p = 1.5076572693500548083e-09 + p * w;
- p = -3.7894654401267369937e-09 + p * w;
- p = 7.6157012080783393804e-09 + p * w;
- p = -1.4960026627149240478e-08 + p * w;
- p = 2.9147953450901080826e-08 + p * w;
- p = -6.7711997758452339498e-08 + p * w;
- p = 2.2900482228026654717e-07 + p * w;
- p = -9.9298272942317002539e-07 + p * w;
- p = 4.5260625972231537039e-06 + p * w;
- p = -1.9681778105531670567e-05 + p * w;
- p = 7.5995277030017761139e-05 + p * w;
- p = -0.00021503011930044477347 + p * w;
- p = -0.00013871931833623122026 + p * w;
- p = 1.0103004648645343977 + p * w;
- p = 4.8499064014085844221 + p * w;
- } else if (w == Double.POSITIVE_INFINITY) {
- // this branch does not appears in the original code, it
- // was added because the previous branch does not handle
- // x = +/-1 correctly. In this case, w is positive infinity
- // and as the first coefficient (-2.71e-11) is negative.
- // Once the first multiplication is done, p becomes negative
- // infinity and remains so throughout the polynomial evaluation.
- // So the branch above incorrectly returns negative infinity
- // instead of the correct positive infinity.
- p = Double.POSITIVE_INFINITY;
- } else {
- // this branch does not appears in the original code, it
- // occurs when the input is NaN or not in the range [-1, 1].
- return Double.NaN;
- }
-
- return p * x;
+ return BoostErf.erfInv(x);
}
}
-
diff --git
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErfc.java
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErfc.java
index 60a8e7e..26ea355 100644
---
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErfc.java
+++
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/InverseErfc.java
@@ -18,13 +18,6 @@ package org.apache.commons.numbers.gamma;
/**
* Inverse of the <a
href="http://mathworld.wolfram.com/Erfc.html";>complementary error function</a>.
- * <p>
- * This implementation is described in the paper:
- * <a
href="http://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf";>Approximating
- * the erfinv function</a> by Mike Giles, Oxford-Man Institute of Quantitative
Finance,
- * which was published in GPU Computing Gems, volume 2, 2010.
- * The source code is available <a
href="http://gpucomputing.net/?q=node/1828";>here</a>.
- * </p>
*/
public final class InverseErfc {
/** Private constructor. */
@@ -35,10 +28,19 @@ public final class InverseErfc {
/**
* Returns the inverse complementary error function.
*
- * @param x Value.
+ * <p>Special cases:
+ * <ul>
+ * <li>If the argument is 1, then the result is 0.
+ * <li>If the argument is 0, then the result is positive infinity.
+ * <li>If the argument is 2, then the result is negative infinity.
+ * <li>If the argument is outside the interval {@code [0, 2]}, then the
result is nan.
+ * <li>If the argument is nan, then the result is nan.
+ * </ul>
+ *
+ * @param x Value (in {@code [0, 2]})
* @return t such that {@code x =} {@link Erfc#value(double)
Erfc.value(t)}.
*/
public static double value(double x) {
- return InverseErf.value(1 - x);
+ return BoostErf.erfcInv(x);
}
}
diff --git
a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/BoostErfTest.java
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/BoostErfTest.java
new file mode 100644
index 0000000..8f2a6f0
--- /dev/null
+++
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/BoostErfTest.java
@@ -0,0 +1,433 @@
+/*
+ * 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 java.math.BigDecimal;
+import java.util.function.DoubleUnaryOperator;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.MethodOrderer;
+import org.junit.jupiter.api.Order;
+import org.junit.jupiter.api.Test;
+import org.junit.jupiter.api.TestMethodOrder;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvFileSource;
+import org.junit.jupiter.params.provider.CsvSource;
+
+/**
+ * Tests for {@link BoostErf}.
+ *
+ * <p>Note: Some resource data files used in these tests have been extracted
+ * from the Boost test files for the Erf functions.
+ */
+@TestMethodOrder(MethodOrderer.OrderAnnotation.class)
+class BoostErfTest {
+ // Test order:
+ // Any methods not annotated with @Order will be executed as if using int
max value.
+ // This class uses the @Order annotation to build the RMS values before
asserting
+ // the max RMS. For convenience when using reporting the ulp errors this
is done
+ // in the order of the class methods.
+
+ /**
+ * Define the function to test.
+ * This is a utility to identify if the function is an odd function: f(x)
= -f(-x).
+ */
+ private enum TestFunction {
+ ERF_INV(BoostErf::erfInv, true),
+ ERFC_INV(BoostErf::erfcInv, false);
+
+ /** The function. */
+ private final DoubleUnaryOperator fun;
+
+ /** Odd function flag: f(x) = -f(-x). */
+ private final boolean odd;
+
+ /**
+ * @param fun function to test
+ * @param odd true if an odd function
+ */
+ TestFunction(DoubleUnaryOperator fun, boolean odd) {
+ this.fun = fun;
+ this.odd = odd;
+ }
+
+ /**
+ * @return the function
+ */
+ DoubleUnaryOperator getFunction() {
+ return fun;
+ }
+
+ /**
+ * @return true if an odd function: f(x) = -f(-x)
+ */
+ boolean isOdd() {
+ return odd;
+ }
+ }
+
+ /**
+ * Define the test cases for each resource file.
+ * This encapsulates the function to test, the expected maximum and RMS
error, and
+ * methods to accumulate error and compute RMS error.
+ *
+ * <h2>Note on accuracy</h2>
+ *
+ * <p>The Boost functions use the default policy of internal promotion
+ * of double to long double if it offers more precision. Code comments
+ * in the implementations for the maximum error are using the defaults with
+ * promotion enabled where the error is 'effectively zero'. Java does not
+ * support long double computation. Tolerances have been set to allow
tests to
+ * pass. Spot checks on larger errors have been verified against the
reference
+ * implementation compiled with promotion of double
<strong>disabled</strong>.
+ *
+ * @see <a
href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/relative_error.html";>Relative
error</a>
+ * @see <a
href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/pol_tutorial/policy_tut_defaults.html";>Policy
defaults</a>
+ */
+ private enum TestCase {
+ /** Inverse Erf Boost data. */
+ ERF_INV(TestFunction.ERF_INV, 1.8, 0.65),
+ /** Inverse Erfc Boost data. */
+ ERFC_INV(TestFunction.ERFC_INV, 1.8, 0.65),
+ /** Inverse Erfc Boost big data. */
+ ERFC_INV_BIG(TestFunction.ERFC_INV, 1.8, 0.5),
+ /** Inverse Erf limit data. */
+ ERF_INV_LIMIT(TestFunction.ERF_INV, 1.4, 0.5),
+ /** Inverse Erfc limit data. */
+ ERFC_INV_LIMIT(TestFunction.ERFC_INV, 1.2, 0.5);
+
+ /** Sum of the squared ULP error and count n. */
+ private final RMS rms = new RMS();
+
+ /** The function. */
+ private final TestFunction fun;
+
+ /** The maximum allowed ulp. */
+ private final double maxUlp;
+
+ /** The maximum allowed RMS ulp. */
+ private final double rmsUlp;
+
+ /**
+ * @param fun function to test
+ * @param maxUlp maximum allowed ulp
+ * @param rmsUlp maximum allowed RMS ulp
+ */
+ TestCase(TestFunction fun, double maxUlp, double rmsUlp) {
+ this.fun = fun;
+ this.maxUlp = maxUlp;
+ this.rmsUlp = rmsUlp;
+ }
+
+ /**
+ * @return function to test
+ */
+ DoubleUnaryOperator getFunction() {
+ return fun.getFunction();
+ }
+
+ /**
+ * @return true if an odd function
+ */
+ boolean isOdd() {
+ return fun.isOdd();
+ }
+
+ /**
+ * @param ulp error in ulp
+ */
+ void addError(double ulp) {
+ rms.add(ulp);
+ }
+
+ /**
+ * @return Root Mean Squared measured error
+ */
+ double getRMSError() {
+ return rms.getRMS();
+ }
+
+ /**
+ * @return maximum measured error
+ */
+ double getMaxError() {
+ return rms.getMax();
+ }
+
+ /**
+ * @return maximum allowed error
+ */
+ double getTolerance() {
+ return maxUlp;
+ }
+
+ /**
+ * @return maximum allowed RMS error
+ */
+ double getRmsTolerance() {
+ return rmsUlp;
+ }
+ }
+
+ /**
+ * Class to compute the root mean squared error (RMS).
+ * @see <a
href="https://en.wikipedia.org/wiki/Root_mean_square";>Wikipedia: RMS</a>
+ */
+ private static class RMS {
+ private double ss;
+ private double max;
+ private int n;
+
+ /**
+ * @param x Value (assumed to be positive)
+ */
+ void add(double x) {
+ // Overflow is not supported.
+ // Assume the expected and actual are quite close when measuring
the RMS.
+ ss += x * x;
+ n++;
+ // Assume absolute error when detecting the maximum
+ max = max < x ? x : max;
+ }
+
+ /**
+ * Gets the maximum error.
+ *
+ * <p>This is not used for assertions. It can be used to set maximum
ULP thresholds
+ * for test data if the TestUtils.assertEquals method is used with a
large maxUlps
+ * to measure the ulp (and effectively ignore failures) and the
maximum reported
+ * as the end of testing.
+ *
+ * @return maximum error
+ */
+ double getMax() {
+ return max;
+ }
+
+ /**
+ * Gets the root mean squared error (RMS).
+ *
+ * <p> Note: If no data has been added this will return 0/0 = nan.
+ * This prevents using in assertions without adding data.
+ *
+ * @return root mean squared error (RMS)
+ */
+ double getRMS() {
+ return Math.sqrt(ss / n);
+ }
+ }
+
+ @ParameterizedTest
+ @CsvSource({
+ "0, 0",
+ // Odd function
+ "-0.0, -0.0",
+ "1, Infinity",
+ "-1, -Infinity",
+ "NaN, NaN",
+ // Domain errors do not throw
+ "-1.1, NaN",
+ "1.1, NaN",
+ })
+ void testInverseErfEdgeCases(double p, double z) {
+ Assertions.assertEquals(z, BoostErf.erfInv(p));
+ }
+
+ @ParameterizedTest
+ @CsvSource({
+ "1, 0",
+ "0, Infinity",
+ "2, -Infinity",
+ "NaN, NaN",
+ // Domain errors do not throw
+ "-0.1, NaN",
+ "2.1, NaN",
+ })
+ void testInverseErfcEdgeCases(double p, double z) {
+ Assertions.assertEquals(z, BoostErf.erfcInv(p));
+ }
+
+ @ParameterizedTest
+ @Order(1)
+ @CsvFileSource(resources = "erf_inv_data.csv")
+ void testInverseErf(double p, BigDecimal z) {
+ assertErf(TestCase.ERF_INV, p, z);
+ }
+
+ @Test
+ @Order(2010)
+ void testInverseErfRMS() {
+ assertRms(TestCase.ERF_INV);
+ }
+
+ @ParameterizedTest
+ @Order(1)
+ @CsvFileSource(resources = "erfc_inv_data.csv")
+ void testInverseErfc(double q, BigDecimal z) {
+ assertErf(TestCase.ERFC_INV, q, z);
+ }
+
+ @Test
+ @Order(2020)
+ void testInverseErfcRMS() {
+ assertRms(TestCase.ERFC_INV);
+ }
+
+ @ParameterizedTest
+ @Order(1)
+ @CsvFileSource(resources = "erfc_inv_big_data.csv")
+ void testInverseErfcBig(double q, BigDecimal z) {
+ assertErf(TestCase.ERFC_INV_BIG, q, z);
+ }
+
+ @Test
+ @Order(2030)
+ void testInverseErfcBigRMS() {
+ assertRms(TestCase.ERFC_INV_BIG);
+ }
+
+ @ParameterizedTest
+ @Order(1)
+ @CsvFileSource(resources = "erf_inv_limit_data.csv")
+ void testInverseErfLimit(double p, BigDecimal z) {
+ assertErf(TestCase.ERF_INV_LIMIT, p, z);
+ }
+
+ @Test
+ @Order(2040)
+ void testInverseErfLimitRMS() {
+ assertRms(TestCase.ERF_INV_LIMIT);
+ }
+
+ @ParameterizedTest
+ @Order(1)
+ @CsvFileSource(resources = "erfc_inv_limit_data.csv")
+ void testInverseErfcLimit(double q, BigDecimal z) {
+ assertErf(TestCase.ERFC_INV_LIMIT, q, z);
+ }
+
+ @Test
+ @Order(2050)
+ void testInverseErfcLimitRMS() {
+ assertRms(TestCase.ERFC_INV_LIMIT);
+ }
+
+ /**
+ * Round-trip test of the inverse erf and then the erf to return to the
original value.
+ */
+ @Test
+ void testErfRoundTrip() {
+ assertRoundTrip("erf(erfInv(x))",
+ x -> Erf.value(BoostErf.erfInv(x)),
+ // Inverse Erf domain: [-1, 1]
+ -0.95, 1, 0.125,
+ 5L, 2.99);
+ }
+
+ /**
+ * Round-trip test of the inverse erfc and then the erfc to return to the
original value.
+ */
+ @Test
+ void testErfcRoundTrip() {
+ assertRoundTrip("erfc(erfcInv(x))",
+ x -> Erfc.value(BoostErf.erfcInv(x)),
+ // Inverse Erfc domain: [0, 2]
+ 0.125, 2, 0.125,
+ 15L, 3.99);
+ }
+
+ /**
+ * Test a round-trip function (foward and reverse) returns to the original
value.
+ *
+ * @param name Test name
+ * @param fun Round-trip function
+ * @param low Low bound to test
+ * @param high Upper bound to test
+ * @param increment Increment between bounds
+ * @param tolerance Maximum ULP tolerance
+ * @param rmsUlp Maximum RMS ULP
+ */
+ private static void assertRoundTrip(String name,
+ DoubleUnaryOperator fun,
+ double low, double high, double increment,
+ long tolerance, double rmsUlp) {
+ final RMS data = new RMS();
+ for (double p = low; p <= high; p += increment) {
+ final double pp = fun.applyAsDouble(p);
+ TestUtils.assertEquals(p, pp, tolerance, ulp -> data.add(ulp), ()
-> name);
+ }
+ assertRms(name, data, rmsUlp);
+ }
+
+ /**
+ * Assert the function using extended precision.
+ *
+ * @param tc Test case
+ * @param x Input value
+ * @param expected Expected value
+ */
+ private static void assertErf(TestCase tc, double x, BigDecimal expected) {
+ final double actual = tc.getFunction().applyAsDouble(x);
+ TestUtils.assertEquals(expected, actual, tc.getTolerance(),
tc::addError,
+ () -> tc + " x=" + x);
+ if (tc.isOdd()) {
+ // Use a delta of 0.0 to allow -0.0 == 0.0
+ Assertions.assertEquals(actual,
-tc.getFunction().applyAsDouble(-x), 0.0, "odd function: f(x) = -f(-x)");
+ }
+ }
+
+ /**
+ * Assert the Root Mean Square (RMS) error of the function is below the
allowed maximum.
+ *
+ * @param tc Test case
+ */
+ private static void assertRms(TestCase tc) {
+ final double rms = tc.getRMSError();
+ debugRms(tc.toString(), tc.getMaxError(), rms);
+ Assertions.assertTrue(rms <= tc.getRmsTolerance(),
+ () -> String.format("%s RMS %s < %s", tc, rms,
tc.getRmsTolerance()));
+ }
+
+ /**
+ * Assert the Root Mean Square (RMS) error of the function is below the
allowed maximum.
+ *
+ * @param name Test name
+ * @param data Test data
+ * @param rmsTolerance RMS tolerance
+ */
+ private static void assertRms(String name, RMS data, double rmsTolerance) {
+ final double rms = data.getRMS();
+ debugRms(name, data.getMax(), rms);
+ Assertions.assertTrue(rms <= rmsTolerance,
+ () -> String.format("%s RMS %s < %s", name, rms, rmsTolerance));
+ }
+
+ /**
+ * Output the maximum and RMS ulp for the named test.
+ * Used for reporting the errors and setting appropriate test tolerances.
+ * This is relevant across different JDK implementations where the
java.util.Math
+ * functions used in BoostErf may compute to different accuracy.
+ *
+ * @param name Test name
+ * @param maxUlp Maximum ulp
+ * @param rmsUlp RMS ulp
+ */
+ private static void debugRms(String name, double maxUlp, double rmsUlp) {
+ // CHECKSTYLE: stop regexp
+ // Debug output of max and RMS error.
+ //System.out.printf("%-25s max %10.6g RMS %10.6g%n", name, maxUlp,
rmsUlp);
+ }
+}
diff --git
a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TestUtils.java
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TestUtils.java
new file mode 100644
index 0000000..5c4a269
--- /dev/null
+++
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TestUtils.java
@@ -0,0 +1,256 @@
+/*
+ * 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 java.math.BigDecimal;
+import java.math.MathContext;
+import java.util.function.DoubleConsumer;
+import java.util.function.LongConsumer;
+import java.util.function.Supplier;
+import org.junit.jupiter.api.Assertions;
+
+/**
+ * Test utilities.
+ */
+final class TestUtils {
+ /** Positive zero bits. */
+ private static final long POSITIVE_ZERO_DOUBLE_BITS =
Double.doubleToRawLongBits(+0.0);
+ /** Negative zero bits. */
+ private static final long NEGATIVE_ZERO_DOUBLE_BITS =
Double.doubleToRawLongBits(-0.0);
+ /** Set this to true to report all deviations to System out when the
maximum ULPs is negative. */
+ private static boolean reportAllDeviations = false;
+
+ /** Private constructor. */
+ private TestUtils() {
+ // intentionally empty.
+ }
+
+ /**
+ * Assert the two numbers are equal within the provided units of least
precision.
+ * The maximum count of numbers allowed between the two values is {@code
maxUlps - 1}.
+ *
+ * <p>The values -0.0 and 0.0 are considered equal.
+ *
+ * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and
ignore
+ * failures.
+ *
+ * @param expected expected value
+ * @param actual actual value
+ * @param maxUlps maximum units of least precision between the two values
+ * @return ulp difference between the values (always positive; may be
truncated to Long.MAX_VALUE)
+ */
+ static long assertEquals(double expected, double actual, long maxUlps) {
+ return assertEquals(expected, actual, maxUlps, null,
(Supplier<String>) null);
+ }
+
+ /**
+ * Assert the two numbers are equal within the provided units of least
+ * precision. The maximum count of numbers allowed between the two values
is
+ * {@code maxUlps - 1}.
+ *
+ * <p>The values -0.0 and 0.0 are considered equal.
+ *
+ * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and
+ * ignore failures.
+ *
+ * @param expected expected value
+ * @param actual actual value
+ * @param maxUlps maximum units of least precision between the two values
+ * @param error Consumer for the ulp difference between the values (always
positive)
+ * @param msg failure message
+ * @return ulp difference between the values (always positive; may be
truncated to Long.MAX_VALUE)
+ */
+ static long assertEquals(double expected, double actual, long maxUlps,
+ LongConsumer error, Supplier<String> msg) {
+ final long e = Double.doubleToLongBits(expected);
+ final long a = Double.doubleToLongBits(actual);
+
+ // Code adapted from Precision#equals(double, double, int) so we
maintain the delta
+ // for the message and return it for reporting.
+
+ long delta;
+ boolean equal;
+ if (e == a) {
+ // Binary equal
+ equal = true;
+ delta = 0;
+ } else if ((a ^ e) < 0L) {
+ // The difference is the count of numbers between each and zero.
+ // This makes -0.0 and 0.0 equal.
+ long d1;
+ long d2;
+ if (a < e) {
+ d1 = e - POSITIVE_ZERO_DOUBLE_BITS;
+ d2 = a - NEGATIVE_ZERO_DOUBLE_BITS;
+ } else {
+ d1 = a - POSITIVE_ZERO_DOUBLE_BITS;
+ d2 = e - NEGATIVE_ZERO_DOUBLE_BITS;
+ }
+ // This may overflow but we report it using an unsigned formatter.
+ delta = d1 + d2;
+ if (delta < 0) {
+ // Overflow
+ equal = false;
+ } else {
+ // Allow input of a negative maximum ULPs
+ equal = delta <= ((maxUlps < 0) ? (-maxUlps - 1) : maxUlps);
+ }
+ } else {
+ delta = Math.abs(e - a);
+ // Allow input of a negative maximum ULPs
+ equal = delta <= ((maxUlps < 0) ? (-maxUlps - 1) : maxUlps);
+ }
+
+ // DEBUG:
+ if (maxUlps < 0) {
+ // CHECKSTYLE: stop Regex
+ if (!equal || reportAllDeviations) {
+ System.out.printf("%sexpected <%s> != actual <%s> (ulps=%s)%n",
+ prefix(msg), expected, actual,
Long.toUnsignedString(delta));
+ }
+ // CHECKSTYLE: resume Regex
+ } else if (!equal) {
+ Assertions.fail(String.format("%sexpected <%s> != actual <%s>
(ulps=%s)",
+ prefix(msg), expected, actual, Long.toUnsignedString(delta)));
+ }
+
+ // This may have overflowed.
+ delta = delta < 0 ? Long.MAX_VALUE : delta;
+ if (error != null) {
+ error.accept(delta);
+ }
+ return delta;
+ }
+
+ /**
+ * Assert the two numbers are equal within the provided units of least
precision.
+ *
+ * <p>This method is for values that can be computed to arbitrary
precision.
+ * It raises an exception when the actual value is not finite and the
expected value
+ * has a non-infinite representation; or the actual value is finite and
the expected
+ * value has a infinite representation. In this case the computed ulp
difference
+ * is infinite.
+ *
+ * <p>This method expresses the error relative the units in the last place
of the
+ * expected value when converted to a {@code double} type
+ * (see {@link #assertEquals(BigDecimal, double, double, DoubleConsumer,
Supplier)} for details).
+ *
+ * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and
ignore
+ * failures.
+ *
+ * @param expected expected value
+ * @param actual actual value
+ * @param maxUlps maximum units of least precision between the two values
+ * @return ulp difference between the values (always positive)
+ */
+ static double assertEquals(BigDecimal expected, double actual, double
maxUlps) {
+ return assertEquals(expected, actual, maxUlps, null,
(Supplier<String>) null);
+ }
+
+ /**
+ * Assert the two numbers are equal within the provided units of least
+ * precision.
+ *
+ * <p>This method is for values that can be computed to arbitrary
precision. It
+ * raises an exception when the actual value is not finite and the expected
+ * value has a non-infinite representation; or the actual value is finite
and
+ * the expected value has a infinite representation. In this case the
computed
+ * ulp difference is infinite.
+ *
+ * <p>This method expresses the error relative the units in the last place
(ulp)
+ * of the expected value when converted to a {@code double} type. If the
actual
+ * value equals the expected value the error is 0. Otherwise the error is
+ * computed relative to the ulp of the expected value. The minimum non-zero
+ * value for the error is 0.5. A ulp of < 1.0 indicates the value is the
closest
+ * value to the result that is not exact.
+ *
+ * <pre>
+ * ulp -1 0 1 2
+ * --------------|---------------|---------------|---------------|--------
+ * ^
+ * |
+ * expected
+ *
+ * a b c
+ *
+ * a = 0.75 ulp
+ * b = 0 ulp
+ * c = 1.25 ulp
+ * </pre>
+ *
+ * <p>Set {@code maxUlps} to negative to report the ulps to the stdout and
+ * ignore failures.
+ *
+ * @param expected expected value
+ * @param actual actual value
+ * @param maxUlps maximum units of least precision between the two values
+ * @param error Consumer for the ulp difference between the values (always
positive)
+ * @param msg failure message
+ * @return ulp difference between the values (always positive)
+ */
+ static double assertEquals(BigDecimal expected, double actual, double
maxUlps,
+ DoubleConsumer error, Supplier<String> msg) {
+ final double e = expected.doubleValue();
+
+ double delta;
+ boolean equal;
+ if (e == actual) {
+ // Binary equal. This will match infinity if expected is very
large.
+ equal = true;
+ delta = 0;
+ } else if (!Double.isFinite(e) || !Double.isFinite(actual)) {
+ // No representable delta between infinite and non-infinite values
+ equal = false;
+ delta = Double.POSITIVE_INFINITY;
+ } else {
+ // Two finite numbers
+ delta = expected.subtract(new BigDecimal(actual)).abs()
+ .divide(new BigDecimal(Math.ulp(e)),
MathContext.DECIMAL64).doubleValue();
+ // Allow input of a negative maximum ULPs
+ equal = delta <= Math.abs(maxUlps);
+ }
+
+ if (error != null) {
+ error.accept(delta);
+ }
+
+ // DEBUG:
+ if (maxUlps < 0) {
+ // CHECKSTYLE: stop Regex
+ if (!equal || reportAllDeviations) {
+ System.out.printf("%sexpected <%s> != actual <%s> (ulps=%s)%n",
+ prefix(msg), expected, actual, delta);
+ }
+ // CHECKSTYLE: resume Regex
+ } else if (!equal) {
+ Assertions.fail(String.format("%sexpected <%s> != actual <%s>
(ulps=%s)",
+ prefix(msg), expected, actual, delta));
+ }
+
+ return delta;
+ }
+
+ /**
+ * Get the prefix for the message.
+ *
+ * @param msg Message supplier
+ * @return the prefix
+ */
+ private static String prefix(Supplier<String> msg) {
+ return msg == null ? "" : msg.get() + ": ";
+ }
+}
diff --git
a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TestUtilsTest.java
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TestUtilsTest.java
new file mode 100644
index 0000000..6447735
--- /dev/null
+++
b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/TestUtilsTest.java
@@ -0,0 +1,168 @@
+/*
+ * 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 java.math.BigDecimal;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+/**
+ * Tests for {@link TestUtils}.
+ * This verifies the custom ULP assertions function as expected.
+ */
+class TestUtilsTest {
+ @Test
+ void testAssertEquals() {
+ final double x = 1.23;
+ final double ulp = Math.ulp(x);
+ Assertions.assertEquals(0, TestUtils.assertEquals(x, x, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(x, x, 1));
+ Assertions.assertEquals(0, TestUtils.assertEquals(x, x, 30));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(x, x + ulp, 0));
+ Assertions.assertEquals(1, TestUtils.assertEquals(x, x + ulp, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(x, x - ulp, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(x, x + ulp, 30));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(x, x + 30 * ulp, 29));
+ Assertions.assertEquals(30, TestUtils.assertEquals(x, x + 30 * ulp,
30));
+ Assertions.assertEquals(30, TestUtils.assertEquals(x, x - 30 * ulp,
30));
+
+ // Check order and sign
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(-x - ulp, x, 0));
+ Assertions.assertEquals(1, TestUtils.assertEquals(-x - ulp, -x, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(x - ulp, x, 1));
+ Assertions.assertEquals(2, TestUtils.assertEquals(x - ulp, x + ulp,
2));
+
+ // Opposite signs
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(-x, x, 500000));
+ }
+
+ @Test
+ void testAssertEqualsZero() {
+ // These are not binary equal but for numeric purposes they are
treated as equal
+ final double[] zero = {-0.0, 0.0};
+ for (final double a : zero) {
+ for (final double b : zero) {
+ Assertions.assertEquals(0, TestUtils.assertEquals(a, b, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(a, b, 1));
+ Assertions.assertEquals(0, TestUtils.assertEquals(a, b, 30));
+ }
+ }
+
+ // Difference from zero
+ final double x = Double.MIN_VALUE;
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(0.0, x, 0));
+ Assertions.assertEquals(1, TestUtils.assertEquals(0.0, x, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(-0.0, x, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(0.0, -x, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(-0.0, -x, 1));
+ Assertions.assertEquals(2, TestUtils.assertEquals(0.0, 2 * x, 2));
+ Assertions.assertEquals(2, TestUtils.assertEquals(-0.0, 2 * x, 2));
+ }
+
+ @Test
+ void testAssertEqualsEdgeCases() {
+ final double nan = Double.NaN;
+ final double inf = Double.POSITIVE_INFINITY;
+ final double max = Double.MAX_VALUE;
+ Assertions.assertEquals(0, TestUtils.assertEquals(nan, nan, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(inf, inf, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(max, max, 0));
+
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(nan, inf, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(nan, max, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(inf, max, 0));
+
+ Assertions.assertEquals(1, TestUtils.assertEquals(inf, max, 1));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(inf, -inf, 0));
+ }
+
+ @Test
+ void testBigDecimalAssertEquals() {
+ final double x = 1.23;
+ final double ulp = Math.ulp(x);
+
+ final BigDecimal bdx = new BigDecimal(x);
+
+ Assertions.assertEquals(0, TestUtils.assertEquals(bdx, x, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(bdx, x, 1));
+ Assertions.assertEquals(0, TestUtils.assertEquals(bdx, x, 30));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdx, x + ulp, 0));
+ Assertions.assertEquals(1, TestUtils.assertEquals(bdx, x + ulp, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(bdx, x - ulp, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(bdx, x + ulp, 30));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdx, x + 30 * ulp, 29));
+ Assertions.assertEquals(30, TestUtils.assertEquals(bdx, x + 30 * ulp,
30));
+ Assertions.assertEquals(30, TestUtils.assertEquals(bdx, x - 30 * ulp,
30));
+
+ // Opposite signs
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdx.negate(), x, 500000));
+
+ // High precision: Add 1/4 ulp to the value:
+ // ---|----------|----------|-----
+ // -1 x | +1
+ // y
+ final BigDecimal bdy = bdx.add(new
BigDecimal(ulp).divide(BigDecimal.valueOf(4)));
+ Assertions.assertEquals(0, TestUtils.assertEquals(bdy, x, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdy, x + ulp, 0));
+ Assertions.assertEquals(0.75, TestUtils.assertEquals(bdy, x + ulp, 1));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdy, x - ulp, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdy, x - ulp, 1));
+ Assertions.assertEquals(1.25, TestUtils.assertEquals(bdy, x - ulp,
1.5), 1e-16);
+ }
+
+ @Test
+ void testBigDecimalAssertEqualsZero() {
+ // These are not binary equal but for numeric purposes they are
treated as equal
+ final double[] zero = {-0.0, 0.0};
+ final BigDecimal a = BigDecimal.ZERO;
+ for (final double b : zero) {
+ Assertions.assertEquals(0, TestUtils.assertEquals(a, b, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(a, b, 1));
+ Assertions.assertEquals(0, TestUtils.assertEquals(a, b, 30));
+ }
+
+ // Difference from zero
+ final double x = Double.MIN_VALUE;
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(a, x, 0));
+ Assertions.assertEquals(1, TestUtils.assertEquals(a, x, 1));
+ Assertions.assertEquals(1, TestUtils.assertEquals(a, -x, 1));
+ Assertions.assertEquals(2, TestUtils.assertEquals(a, 2 * x, 2));
+ Assertions.assertEquals(2, TestUtils.assertEquals(a, -2 * x, 2));
+ }
+
+ @Test
+ void testBigDecimalAssertEqualsEdgeCases() {
+ final double nan = Double.NaN;
+ final double inf = Double.POSITIVE_INFINITY;
+ final double max = Double.MAX_VALUE;
+
+ final BigDecimal bdinf = new BigDecimal("1e310");
+ final BigDecimal bdmax = new BigDecimal(max);
+
+ Assertions.assertEquals(0, TestUtils.assertEquals(bdinf, inf, 0));
+ Assertions.assertEquals(0, TestUtils.assertEquals(bdmax, max, 0));
+
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdinf, max, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdinf, nan, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdmax, inf, 0));
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdmax, nan, 0));
+
+ final double x = Math.nextDown(max);
+ Assertions.assertThrows(AssertionError.class, () ->
TestUtils.assertEquals(bdmax, x, 0));
+ Assertions.assertEquals(1, TestUtils.assertEquals(bdmax, x, 1));
+ }
+}
diff --git
a/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erf_inv_data.csv
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erf_inv_data.csv
new file mode 100644
index 0000000..483b8b1
--- /dev/null
+++
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erf_inv_data.csv
@@ -0,0 +1,122 @@
+# 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.
+
+# (C) Copyright John Maddock 2006-7.
+# Use, modification and distribution are subject to the
+# Boost Software License, Version 1.0. (See accompanying file
+# LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# Extracted from: boost/libs/math/test/erf_inv_data.ipp
+
+-0.990433037281036376953125,-1.832184533179510927322805923563700329767
+-0.936334311962127685546875,-1.311339282092737086640055105484822812599
+-0.931107819080352783203125,-1.286316461685184857889337272829270644576
+-0.928576648235321044921875,-1.274755355308702535979544636706295190547
+-0.92711746692657470703125,-1.268242390461126936128446743938528753699
+-0.907657206058502197265625,-1.190178701802009872651528128114629153367
+-0.89756715297698974609375,-1.154826929034364581020764782706068245251
+-0.80573642253875732421875,-0.917873491512746130598510739795171676962
+-0.804919183254241943359375,-0.9161942446913841771580833184012189161559
+-0.780276477336883544921875,-0.867806431357551134410815525153998291827
+-0.775070965290069580078125,-0.8580919924152284867224297625328485999768
+-0.7496345043182373046875,-0.8127924290810407931222432768905514928867
+-0.74820673465728759765625,-0.8103475955423936417307157186107256388648
+-0.74602639675140380859375,-0.806632688757205819001462030577472890805
+-0.72904598712921142578125,-0.7784313874823551106598826200232666844639
+-0.7162272930145263671875,-0.7579355487144890063429377586715787838945
+-0.701772034168243408203125,-0.7355614373173595299453301005770428516518
+-0.68477380275726318359375,-0.7101588399949052872532446197423489724803
+-0.657626628875732421875,-0.6713881533266128640126408255047881638389
+-0.652269661426544189453125,-0.6639738263692763456669198307149427734317
+-0.6262547969818115234375,-0.6289572925573171740836428308746584439674
+-0.62323606014251708984375,-0.6249936843093662471116097431474787933967
+-0.57958185672760009765625,-0.5697131213589784467617578394703976041604
+-0.576151371002197265625,-0.5655172244109430153330195150336141500139
+-0.5579319000244140625,-0.5435569422159360687847790186563654276103
+-0.446154057979583740234375,-0.4186121208546731033057205459902879301199
+-0.44300353527069091796875,-0.4152898953738801984047941692529271391195
+-0.40594112873077392578125,-0.3768620801611051992528860948080812212023
+-0.396173775196075439453125,-0.3669220210390825311547962776125822899061
+-0.38366591930389404296875,-0.3542977152760563782151668726041057557165
+-0.36689913272857666015625,-0.3375493847053488720470432821496358279516
+-0.365801036357879638671875,-0.3364591774366710656954166264654945559873
+-0.277411997318267822265625,-0.2510244067029671790889794981353227476998
+-0.236883103847503662109375,-0.213115119330839975829499967157244997714
+-0.215545952320098876953125,-0.1934073617841803428235669261097060281642
+-0.202522933483123779296875,-0.1814532246720147926398927046057793150106
+-0.18253767490386962890625,-0.1632073953550647568421821286058243218715
+-0.156477451324462890625,-0.1395756320903277910768376053314442757507
+-0.1558246612548828125,-0.1389857795955756484965030151195660030168
+-0.12251126766204833984375,-0.109002961098867662134935094105847496074
+-0.1088275909423828125,-0.09674694516640724629590870677194632943569
+-0.08402168750762939453125,-0.07460044047654119877070700345343119035515
+-0.05048263072967529296875,-0.04476895818328636312384562686715995170129
+-0.029248714447021484375,-0.0259268064334840921104659134138093242797
+-0.02486217021942138671875,-0.02203709146986755832638577823832075055744
+-0.02047121524810791015625,-0.01814413302702029459097557481591610553903
+-0.018821895122528076171875,-0.01668201759439857888105181293763417899072
+0.0073254108428955078125,0.006492067534753215749601670217642082465642
+0.09376299381256103515625,0.08328747254794857150987333986733043734817
+0.0944411754608154296875,0.08389270963798942778622198997355058545872
+0.264718532562255859375,0.2390787735821979528028028789569770109829
+0.27952671051025390625,0.2530214201700340392837551955289041822603
+0.29262602329254150390625,0.2654374523135675523971788948011709467352
+0.3109557628631591796875,0.282950508503826367238408926581528085458
+0.31148135662078857421875,0.2834552014554130441860525970673030809536
+0.32721102237701416015625,0.2986277427848421570858990348074985028421
+0.3574702739715576171875,0.3282140305634627092431945088114761850208
+0.362719058990478515625,0.3334035993712283467959295804559099468454
+0.3896572589874267578125,0.3603304982893212173104266596348905175268
+0.4120922088623046875,0.3831602323665075533579267768785894144888
+0.41872966289520263671875,0.3899906753567599452444107492361433402154
+0.45167791843414306640625,0.4244594733907945411184647153213164209335
+0.48129451274871826171875,0.4563258063707025027210352963461819167707
+0.4862649440765380859375,0.4617640058971775089811390737537561779898
+0.50937330722808837890625,0.4874174763856674076219106695373814892182
+0.5154802799224853515625,0.4943041993872143888987628020569772222018
+0.52750003337860107421875,0.5079978091910991117615000459548117088362
+0.53103363513946533203125,0.5120597685873370942783226077302881881069
+0.58441460132598876953125,0.5756584292527058478710392476034273328569
+0.5879499912261962890625,0.5800336103175463592377907341030447077804
+0.59039986133575439453125,0.5830784871670823806198622501806646778319
+0.59455978870391845703125,0.588273673825686998734497652983815773459
+0.59585726261138916015625,0.5899005483108011364541949539839185473259
+0.5962116718292236328125,0.5903454775096607218832535637355431851718
+0.6005609035491943359375,0.5958247243549040349587326482492767206448
+0.6150619983673095703125,0.6143583249050861028039832921829036722514
+0.62944734096527099609375,0.6331707263097125575937994856370309207836
+0.64380657672882080078125,0.6524069265890823819975498133014027247554
+0.6469156742095947265625,0.656635855345815020063728463464436343698
+0.67001712322235107421875,0.6888269167957872563013714303376548038671
+0.6982586383819580078125,0.7302336318927408409119676651737758401138
+0.74485766887664794921875,0.8046505193013635090578266413458426260098
+0.75686132907867431640625,0.8253191678260588578995203396384711816647
+0.81158387660980224609375,0.9300427626888758122211127950646282789481
+0.826751708984375,0.9629665092443368464606966822833571908852
+0.83147108554840087890625,0.9736479209913771931387473923084901789046
+0.84174954891204833984375,0.997713670556719074960678197806799852186
+0.8679864406585693359375,1.065050516333636716777334376076374184102
+0.90044414997100830078125,1.164612422633086435501625591693259387477
+0.91433393955230712890625,1.215315881176612875682446995412738776976
+0.91501367092132568359375,1.217962731073139868794942852653058932976
+0.918984889984130859375,1.233778505900771488542027767896521427575
+0.92977702617645263671875,1.28019542575660930623179572273596558907
+0.93538987636566162109375,1.306695301483797253764522033930388453334
+0.93773555755615234375,1.318335478463913327121670503572736587296
+0.94118559360504150390625,1.33613349872692113073358883961598631154
+0.96221935749053955078125,1.468821071545234761861756248744372345584
+0.98576259613037109375,1.733272259459038694476413373595347034928
+0.9881370067596435546875,1.77921769652839903464038407684397479173
+0.99292266368865966796875,1.904368122482929779094714951471938518496
diff --git
a/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erf_inv_limit_data.csv
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erf_inv_limit_data.csv
new file mode 100644
index 0000000..1cf6883
--- /dev/null
+++
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erf_inv_limit_data.csv
@@ -0,0 +1,129 @@
+# 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.
+
+# Erf limits: [-1, 1]
+# Odd function (f(x) = -f(-x)) so test the range [0, 1] and
+# its mirror in the unit test.
+# Test the critical regions of the range using x -> 0, x -> 0.5 and x -> 1.
+#
+# Data generated using the following matlab script:
+#
+# f = fopen('erf_inv_data2.csv', 'wt');
+#
+# % Uses variable precision arithmetic (vpa) to exactly approach 0, 0.5 and 1.
+# digits(30);
+# x = [1:10];
+# % anonymous function writes the argument and the inverse erf function value
to file
+# pr = @(p) fprintf(f, "%.17g,%s\n", p, vpa(erfinv(p)));
+# for i = [1074:-1:1022]
+# pr(sym(2^-i));
+# end
+# for i = flip(x)
+# pr(sym(1/2) - sym(i*2^-54));
+# end
+# for i = x
+# pr(sym(1/2) + sym(i*2^-53));
+# end
+# for i = flip(x)
+# pr(sym(1) - sym(i*2^-53));
+# end
+#
+# clear pr;
+# fclose(f);
+
+4.9406564584124654e-324,4.37854278285719143435811695173e-324
+9.8813129168249309e-324,8.75708556571438286871623390345e-324
+1.9762625833649862e-323,1.75141711314287657374324678069e-323
+3.9525251667299724e-323,3.50283422628575314748649356138e-323
+7.9050503334599447e-323,7.00566845257150629497298712276e-323
+1.5810100666919889e-322,1.40113369051430125899459742455e-322
+3.1620201333839779e-322,2.8022673810286025179891948491e-322
+6.3240402667679558e-322,5.60453476205720503597838969821e-322
+1.2648080533535912e-321,1.12090695241144100719567793964e-321
+2.5296161067071823e-321,2.24181390482288201439135587928e-321
+5.0592322134143646e-321,4.48362780964576402878271175857e-321
+1.0118464426828729e-320,8.96725561929152805756542351713e-321
+2.0236928853657458e-320,1.79345112385830561151308470343e-320
+4.0473857707314917e-320,3.58690224771661122302616940685e-320
+8.0947715414629834e-320,7.17380449543322244605233881371e-320
+1.6189543082925967e-319,1.43476089908664448921046776274e-319
+3.2379086165851934e-319,2.86952179817328897842093552548e-319
+6.4758172331703867e-319,5.73904359634657795684187105097e-319
+1.2951634466340773e-318,1.14780871926931559136837421019e-318
+2.5903268932681547e-318,2.29561743853863118273674842039e-318
+5.1806537865363094e-318,4.59123487707726236547349684077e-318
+1.0361307573072619e-317,9.18246975415452473094699368155e-318
+2.0722615146145237e-317,1.83649395083090494618939873631e-317
+4.1445230292290475e-317,3.67298790166180989237879747262e-317
+8.289046058458095e-317,7.34597580332361978475759494524e-317
+1.657809211691619e-316,1.46919516066472395695151898905e-316
+3.315618423383238e-316,2.93839032132944791390303797809e-316
+6.631236846766476e-316,5.87678064265889582780607595619e-316
+1.3262473693532952e-315,1.17535612853177916556121519124e-315
+2.6524947387065904e-315,2.35071225706355833112243038248e-315
+5.3049894774131808e-315,4.70142451412711666224486076495e-315
+1.0609978954826362e-314,9.4028490282542333244897215299e-315
+2.1219957909652723e-314,1.88056980565084666489794430598e-314
+4.2439915819305446e-314,3.76113961130169332979588861196e-314
+8.4879831638610893e-314,7.52227922260338665959177722392e-314
+1.6975966327722179e-313,1.50445584452067733191835544478e-313
+3.3951932655444357e-313,3.00891168904135466383671088957e-313
+6.7903865310888714e-313,6.01782337808270932767342177914e-313
+1.3580773062177743e-312,1.20356467561654186553468435583e-312
+2.7161546124355486e-312,2.40712935123308373106936871166e-312
+5.4323092248710971e-312,4.81425870246616746213873742331e-312
+1.0864618449742194e-311,9.62851740493233492427747484662e-312
+2.1729236899484389e-311,1.92570348098646698485549496932e-311
+4.3458473798968777e-311,3.85140696197293396971098993865e-311
+8.6916947597937554e-311,7.7028139239458679394219798773e-311
+1.7383389519587511e-310,1.54056278478917358788439597546e-310
+3.4766779039175022e-310,3.08112556957834717576879195092e-310
+6.9533558078350043e-310,6.16225113915669435153758390184e-310
+1.3906711615670009e-309,1.23245022783133887030751678037e-309
+2.7813423231340017e-309,2.46490045566267774061503356073e-309
+5.5626846462680035e-309,4.92980091132535548123006712147e-309
+1.1125369292536007e-308,9.85960182265071096246013424294e-309
+2.2250738585072014e-308,1.97192036453014219249202684859e-308
+0.49999999999999944,0.476936276204469255772776918108
+0.4999999999999995,0.476936276204469317533641061662
+0.49999999999999956,0.476936276204469379294505205215
+0.49999999999999961,0.476936276204469441055369348769
+0.49999999999999967,0.476936276204469502816233492322
+0.49999999999999972,0.476936276204469564577097635876
+0.49999999999999978,0.476936276204469626337961779429
+0.49999999999999983,0.476936276204469688098825922983
+0.49999999999999989,0.476936276204469749859690066536
+0.49999999999999994,0.47693627620446981162055421009
+0.50000000000000011,0.47693627620446999690314664075
+0.50000000000000022,0.476936276204470120424874927857
+0.50000000000000033,0.476936276204470243946603214964
+0.50000000000000044,0.476936276204470367468331502071
+0.50000000000000056,0.476936276204470490990059789178
+0.50000000000000067,0.476936276204470614511788076286
+0.50000000000000078,0.476936276204470738033516363393
+0.50000000000000089,0.4769362762044708615552446505
+0.500000000000001,0.476936276204470985076972937607
+0.50000000000000111,0.476936276204471108598701224714
+0.99999999999999889,5.66676502191633181061441214145
+0.999999999999999,5.67591573974471317882020497616
+0.99999999999999911,5.68612844131039098005160412756
+0.99999999999999922,5.69768514275317414947187188334
+0.99999999999999933,5.71099812575505447299174850294
+0.99999999999999944,5.72670523026542455395142564004
+0.99999999999999956,5.74587239219118047033348138881
+0.99999999999999967,5.77049185355333287054782202161
+0.99999999999999978,5.80501868319345330018125827039
+0.99999999999999989,5.86358474875516792720766258583
+
diff --git
a/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_big_data.csv
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_big_data.csv
new file mode 100644
index 0000000..a1fb022
--- /dev/null
+++
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_big_data.csv
@@ -0,0 +1,71 @@
+# 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.
+
+# (C) Copyright John Maddock 2006-7.
+# Use, modification and distribution are subject to the
+# Boost Software License, Version 1.0. (See accompanying file
+# LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# Extracted from: boost/libs/math/test/erf_inv_big_data.ipp
+
+0.1040364789689237257780791904879173249974e-287,25.67656975866313564221220720210665888424
+0.2450719147172718682892603854216257656211e-282,25.43473899116627259523415573506406682191
+0.1540635014300999959290867576988859893076e-280,25.35326739382267978364899192426764093558
+0.1611739209855438499046222515484436760851e-273,25.03273169525729329624232590800413859205
+0.2781078901633922872827030068281452871724e-257,24.27512000889987237167959986465966890736
+0.8837820681846044382839818307987545000847e-230,22.93495456229819340317355998004446065241
+0.2845511846316072211634522373099423142668e-217,22.29887672117465208331735160799934601323
+0.662011186385949636511655220013184768187e-212,22.02033498263166080103335807537366542004
+0.1113036584088835972703336205837326051999e-204,21.63966022176108818726656867477709317659
+0.891278753475701910699424577762349195273e-195,21.10677883438545215904960900943812596651
+0.1243375606125932867518669565282002251187e-191,20.93474611896499049089384377937974962862
+0.963481566993839758726177317091390488753e-188,20.72000418866581389344212197908962674052
+0.1401265332225531798315780528328964001988e-167,19.56909048260682534539216053987194710384
+0.7560947579095629715718115752684784187336e-154,18.74494560773490096872968649468667611647
+0.2179858592534951368829782528927445425064e-151,18.59346814776643084908600446605555382182
+0.5208405693419352209961246987938629746382e-131,17.28776816936911830963159975455855394785
+0.5376896865772724579278898126622880802843e-127,17.01882418428649428645068852272615582504
+0.1449951754876351411423239477096839160613e-114,16.15763261758311680195806763260744183267
+0.3189212836315126918184908007667112228414e-108,15.70012576642621709878560027430710336585
+0.1297369886602648502231884634897739722895e-107,15.65546665388935777226018284593032382368
+0.2602821327589283546034644516083144727197e-103,15.33647743152077165227354683273671003413
+0.1590035594081776210615338926776093483999e-96,14.81946110818452708436758045284808001241
+0.2720908074515874075065365442027770927174e-96,14.80136596875698966601900922593357151598
+0.3706750818860597777102228512408751197629e-87,14.07473190298865115531360325271322276549
+0.3489539664837006993042769801460156282803e-83,13.74669364199157990777768669458685265472
+0.1947620269151566018960312643171165602404e-79,13.43009987431500967622734624289673694188
+0.5129832297533069185312436951335270626736e-73,12.86957639029561313119037589037764922829
+0.9272206234111032518003166213125699095688e-70,12.57574048509944911276422466059543556873
+0.1460449453933976163637836336160553873706e-68,12.46599637375040876166599138465603014844
+0.8768204308820258161647995326423963318019e-68,12.39412855962299302779918858089973570208
+0.3176564512338672251415481009303329209995e-61,11.77127316747622053532997731227084711173
+0.1172937774279091349820752991146331290881e-54,11.11297358832626209702547835299261897395
+0.3822700066028246575792241970441049035743e-52,10.85058828114725350216779162566256390403
+0.4793602326185941396320885530096873382082e-41,9.607343723542133339487829511083057069273
+0.106201511009181306109371950926034372112e-40,9.566077632118694585964030456246300276422
+0.1026229395621872573137267432392849554594e-24,7.413111396559360630306454710337113823294
+0.1433834013674021626308791056277965507134e-19,6.574536501026477531602186421374903506051
+0.2287484604245551701383736767135709586086e-12,5.183672646329402197657171322814515766894
+0.4054382118091501052221361988694253843776e-11,4.903978376537783549083010816442468785255
+0.5517250582467272934968166422809728429772e-9,4.386634719241854622420456791837035899963
+0.1708920553446964553659243704866375034024e-5,3.383585107532389209677350898502297644278
+0.2892930668406028955896934416060055101499e-5,3.308042643013470022433136596357284983766
+0.7029998339557698614154475319387477445787e-5,3.17687395056466498458108882620439294841
+0.3973417837559780922949941034264323413971e-4,2.90551585834893682957467771364140403613
+0.04246334897368917726269942070072005435577,1.434684541881674882915735690265484394337
+0.159920364968560744996532371753339418774,0.9937250495458864822841753804433979428041
+0.2425390104583346613758947085681683120129,0.826370276571483139432927467391341431105
+0.7954739362140872788414606986417965117653,0.1832884885283639734589195991287031281603
+0.9236374866673857278562449757419727802699,0.06777816075457032123722158048171767182127
diff --git
a/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_data.csv
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_data.csv
new file mode 100644
index 0000000..1c7b2a9
--- /dev/null
+++
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_data.csv
@@ -0,0 +1,122 @@
+# 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.
+
+# (C) Copyright John Maddock 2006-7.
+# Use, modification and distribution are subject to the
+# Boost Software License, Version 1.0. (See accompanying file
+# LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# Extracted from: boost/libs/math/test/erfc_inv_data.ipp
+
+0.00956696830689907073974609375,1.832184391051582711731256541599359331735
+0.063665688037872314453125,1.311339282092737086640055105484822812599
+0.068892158567905426025390625,1.286316565305373898738127895195338854718
+0.071423359215259552001953125,1.274755321776344058215704428086211324587
+0.0728825032711029052734375,1.268242522387371561738393687518023984868
+0.09234277904033660888671875,1.190178756246535875259766567441510867604
+0.102432854473590850830078125,1.154826903977823078146497880706118113588
+0.1942635476589202880859375,0.9178735528443878579995511780412810469667
+0.19508080184459686279296875,0.9161942752629032646404767631869277618212
+0.21972350776195526123046875,0.8678064594007161661713535461829067693456
+0.224929034709930419921875,0.8580919924152284867224297625328485999768
+0.2503655254840850830078125,0.8127923779477598070926819995714374417663
+0.25179326534271240234375,0.8103475955423936417307157186107256388648
+0.2539736330509185791015625,0.8066326381314558738773191719462921998277
+0.27095401287078857421875,0.7784313874823551106598826200232666844639
+0.2837726771831512451171875,0.7579355956263440770864088550908621329508
+0.298227965831756591796875,0.7355614373173595299453301005770428516518
+0.3152261674404144287109375,0.7101588837290742217667270852502075077668
+0.342373371124267578125,0.6713881533266128640126408255047881638389
+0.347730338573455810546875,0.6639738263692763456669198307149427734317
+0.3737452030181884765625,0.6289572925573171740836428308746584439674
+0.37676393985748291015625,0.6249936843093662471116097431474787933967
+0.42041814327239990234375,0.5697131213589784467617578394703976041604
+0.4238486588001251220703125,0.5655171880456876504494070613171955224472
+0.4420680999755859375,0.5435569422159360687847790186563654276103
+0.553845942020416259765625,0.4186121208546731033057205459902879301199
+0.55699646472930908203125,0.4152898953738801984047941692529271391195
+0.59405887126922607421875,0.3768620801611051992528860948080812212023
+0.603826224803924560546875,0.3669220210390825311547962776125822899061
+0.61633408069610595703125,0.3542977152760563782151668726041057557165
+0.63310086727142333984375,0.3375493847053488720470432821496358279516
+0.634198963642120361328125,0.3364591774366710656954166264654945559873
+0.722588002681732177734375,0.2510244067029671790889794981353227476998
+0.763116896152496337890625,0.213115119330839975829499967157244997714
+0.784454047679901123046875,0.1934073617841803428235669261097060281642
+0.797477066516876220703125,0.1814532246720147926398927046057793150106
+0.81746232509613037109375,0.1632073953550647568421821286058243218715
+0.843522548675537109375,0.1395756320903277910768376053314442757507
+0.8441753387451171875,0.1389857795955756484965030151195660030168
+0.87748873233795166015625,0.109002961098867662134935094105847496074
+0.8911724090576171875,0.09674694516640724629590870677194632943569
+0.91597831249237060546875,0.07460044047654119877070700345343119035515
+0.94951736927032470703125,0.04476895818328636312384562686715995170129
+0.970751285552978515625,0.0259268064334840921104659134138093242797
+0.97513782978057861328125,0.02203709146986755832638577823832075055744
+0.97952878475189208984375,0.01814413302702029459097557481591610553903
+0.981178104877471923828125,0.01668201759439857888105181293763417899072
+1.0073254108428955078125,-0.006492067534753215749601670217642082465642
+1.09376299381256103515625,-0.08328747254794857150987333986733043734817
+1.0944411754608154296875,-0.08389270963798942778622198997355058545872
+1.264718532562255859375,-0.2390787735821979528028028789569770109829
+1.27952671051025390625,-0.2530214201700340392837551955289041822603
+1.29262602329254150390625,-0.2654374523135675523971788948011709467352
+1.3109557628631591796875,-0.282950508503826367238408926581528085458
+1.31148135662078857421875,-0.2834552014554130441860525970673030809536
+1.32721102237701416015625,-0.2986277427848421570858990348074985028421
+1.3574702739715576171875,-0.3282140305634627092431945088114761850208
+1.362719058990478515625,-0.3334035993712283467959295804559099468454
+1.3896572589874267578125,-0.3603304982893212173104266596348905175268
+1.4120922088623046875,-0.3831602323665075533579267768785894144888
+1.41872966289520263671875,-0.3899906753567599452444107492361433402154
+1.45167791843414306640625,-0.4244594733907945411184647153213164209335
+1.48129451274871826171875,-0.4563258063707025027210352963461819167707
+1.4862649440765380859375,-0.4617640058971775089811390737537561779898
+1.50937330722808837890625,-0.4874174763856674076219106695373814892182
+1.5154802799224853515625,-0.4943041993872143888987628020569772222018
+1.52750003337860107421875,-0.5079978091910991117615000459548117088362
+1.53103363513946533203125,-0.5120597685873370942783226077302881881069
+1.58441460132598876953125,-0.5756584292527058478710392476034273328569
+1.5879499912261962890625,-0.5800336103175463592377907341030447077804
+1.59039986133575439453125,-0.5830784871670823806198622501806646778319
+1.59455978870391845703125,-0.588273673825686998734497652983815773459
+1.59585726261138916015625,-0.5899005483108011364541949539839185473259
+1.5962116718292236328125,-0.5903454775096607218832535637355431851718
+1.6005609035491943359375,-0.5958247243549040349587326482492767206448
+1.6150619983673095703125,-0.6143583249050861028039832921829036722514
+1.62944734096527099609375,-0.6331707263097125575937994856370309207836
+1.64380657672882080078125,-0.6524069265890823819975498133014027247554
+1.6469156742095947265625,-0.656635855345815020063728463464436343698
+1.67001712322235107421875,-0.6888269167957872563013714303376548038671
+1.6982586383819580078125,-0.7302336318927408409119676651737758401138
+1.74485766887664794921875,-0.8046505193013635090578266413458426260098
+1.75686132907867431640625,-0.8253191678260588578995203396384711816647
+1.81158387660980224609375,-0.9300427626888758122211127950646282789481
+1.826751708984375,-0.9629665092443368464606966822833571908852
+1.83147108554840087890625,-0.9736479209913771931387473923084901789046
+1.84174954891204833984375,-0.997713670556719074960678197806799852186
+1.8679864406585693359375,-1.065050516333636716777334376076374184102
+1.90044414997100830078125,-1.164612422633086435501625591693259387477
+1.91433393955230712890625,-1.215315881176612875682446995412738776976
+1.91501367092132568359375,-1.217962731073139868794942852653058932976
+1.918984889984130859375,-1.233778505900771488542027767896521427575
+1.92977702617645263671875,-1.28019542575660930623179572273596558907
+1.93538987636566162109375,-1.306695301483797253764522033930388453334
+1.93773555755615234375,-1.318335478463913327121670503572736587296
+1.94118559360504150390625,-1.33613349872692113073358883961598631154
+1.96221935749053955078125,-1.468821071545234761861756248744372345584
+1.98576259613037109375,-1.733272259459038694476413373595347034928
+1.9881370067596435546875,-1.77921769652839903464038407684397479173
+1.99292266368865966796875,-1.904368122482929779094714951471938518496
diff --git
a/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_limit_data.csv
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_limit_data.csv
new file mode 100644
index 0000000..e11d63b
--- /dev/null
+++
b/commons-numbers-gamma/src/test/resources/org/apache/commons/numbers/gamma/erfc_inv_limit_data.csv
@@ -0,0 +1,126 @@
+# 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.
+
+# Erfc limits: [2, 0]
+# Test the critical regions of the range using x -> 2, x -> 1 and x -> 0.
+#
+# Data generated using the following matlab script:
+#
+# f = fopen('erfc_inv_data2.csv', 'wt');
+#
+# % Uses variable precision arithmetic (vpa) to exactly approach 2, 1, and 0.
+# digits(30);
+# x = [1:10];
+# % anonymous function writes the argument and the inverse erfc function value
to file
+# pr = @(p) fprintf(f, "%.17g,%s\n", p, vpa(erfcinv(p)));
+# for i = x
+# pr(sym(2) - sym(i*2^-52));
+# end
+# for i = flip(x)
+# pr(sym(1) + sym(i*2^-52));
+# end
+# for i = x
+# pr(sym(1) - sym(i*2^-53));
+# end
+# for i = [1022:1074]
+# pr(sym(2^-i));
+# end
+#
+# clear pr;
+# fclose(f);
+
+1.9999999999999998,-5.80501868319345330018125827039
+1.9999999999999996,-5.74587239219118047033348138881
+1.9999999999999993,-5.71099812575505447299174850294
+1.9999999999999991,-5.68612844131039098005160412756
+1.9999999999999989,-5.66676502191633181061441214145
+1.9999999999999987,-5.65089613274800609228436271922
+1.9999999999999984,-5.63744534527960023748265465926
+1.9999999999999982,-5.62576850990669864999311406674
+1.999999999999998,-5.61544924946836459402473243698
+1.9999999999999978,-5.60620273614591880997628748599
+1.0000000000000022,-0.00000000000000196781907536082825960007390103
+1.000000000000002,-0.00000000000000177103716782474543364006651093
+1.0000000000000018,-0.00000000000000157425526028866260768005912082
+1.0000000000000016,-0.00000000000000137747335275257978172005173072
+1.0000000000000013,-0.00000000000000118069144521649695576004434062
+1.0000000000000011,-0.000000000000000983909537680414129800036950514
+1.0000000000000009,-0.000000000000000787127630144331303840029560411
+1.0000000000000007,-0.000000000000000590345722608248477880022170308
+1.0000000000000004,-0.000000000000000393563815072165651920014780205
+1.0000000000000002,-0.000000000000000196781907536082825960007390103
+0.99999999999999989,0.0000000000000000983909537680414129800036950513
+0.99999999999999978,0.000000000000000196781907536082825960007390103
+0.99999999999999967,0.000000000000000295172861304124238940011085154
+0.99999999999999956,0.000000000000000393563815072165651920014780205
+0.99999999999999944,0.000000000000000491954768840207064900018475257
+0.99999999999999933,0.000000000000000590345722608248477880022170308
+0.99999999999999922,0.00000000000000068873667637628989086002586536
+0.99999999999999911,0.000000000000000787127630144331303840029560411
+0.999999999999999,0.000000000000000885518583912372716820033255462
+0.99999999999999889,0.000000000000000983909537680414129800036950514
+2.2250738585072014e-308,26.5432584542509813824702197142
+1.1125369292536007e-308,26.5563029415321441679153405313
+5.5626846462680035e-309,26.5693410335354631524963692932
+2.7813423231340017e-309,26.5823727396534961059275476332
+1.3906711615670009e-309,26.5953980692558387594638242111
+6.9533558078350043e-310,26.6084170316892033030771306907
+3.4766779039175022e-310,26.6214296362774965380042408218
+1.7383389519587511e-310,26.6344358923218976865134392447
+8.6916947597937554e-311,26.6474358091009358607255378148
+4.3458473798968777e-311,26.6604293958705671923131736698
+2.1729236899484389e-311,26.6734166618642516248908042397
+1.0864618449742194e-311,26.6863976162930293708963792246
+5.4323092248710971e-312,26.699372268345597034754317556
+2.7161546124355486e-312,26.7123406271883834040981478159
+1.3580773062177743e-312,26.7253027019656249108199828507
+6.7903865310888714e-313,26.7382585017994407637028926996
+3.3951932655444357e-313,26.7512080357899077543812138002
+1.6975966327722179e-313,26.7641513130151347383628860778
+8.4879831638610893e-314,26.7770883425313367928370423088
+4.2439915819305446e-314,26.7900191333729090529792854372
+2.1219957909652723e-314,26.8029436945525002284563786595
+1.0609978954826362e-314,26.815862035061085801821439453
+5.3049894774131808e-315,26.8287741638680409104801716724
+2.6524947387065904e-315,26.8416800899212129138981887478
+1.3262473693532952e-315,26.8545798221469936477090752797
+6.631236846766476e-316,26.8674733694503913663725033119
+3.315618423383238e-316,26.8803607407151023760214626828
+1.657809211691619e-316,26.8932419448035823591274814881
+8.289046058458095e-317,26.9061169905571173926026022553
+4.1445230292290475e-317,26.9189858867958946609468413249
+2.0722615146145237e-317,26.9318486423190728660398925827
+1.0361307573072619e-317,26.9447052659048523351659414994
+5.1806537865363094e-318,26.9575557663105448288506308428
+2.5903268932681547e-318,26.9704001522726430500794648573
+1.2951634466340773e-318,26.9832384325068898564572535936
+6.4758172331703867e-319,26.9960706157083471768585828536
+3.2379086165851934e-319,27.0088967105514646341097473482
+1.6189543082925967e-319,27.0217167256901478752331045806
+8.0947715414629834e-320,27.0345306697578266107753941403
+4.0473857707314917e-320,27.0473385513675223647322209666
+2.0236928853657458e-320,27.0601403791119159365716211893
+1.0118464426828729e-320,27.0729361615634145768504148458
+5.0592322134143646e-321,27.0857259072742188779079005788
+2.5296161067071823e-321,27.0985096247763893811123628229
+1.2648080533535912e-321,27.1112873225819129021268414685
+6.3240402667679558e-322,27.1240590091827685756516570398
+3.1620201333839779e-322,27.1368246930509936210922905338
+1.5810100666919889e-322,27.1495843826387488305923857306
+7.9050503334599447e-323,27.1623380863783837808628725138
+3.9525251667299724e-323,27.175085812682501770229502028
+1.9762625833649862e-323,27.1878275699440244823124378705
+9.8813129168249309e-324,27.2005633665362563777429614682
+4.9406564584124654e-324,27.2132932108129488153138238647
diff --git a/pom.xml b/pom.xml
index 2cf2e8a..cbc990d 100644
--- a/pom.xml
+++ b/pom.xml
@@ -156,6 +156,11 @@
</dependency>
<dependency>
<groupId>org.apache.commons</groupId>
+ <artifactId>commons-numbers-gamma</artifactId>
+ <version>${project.version}</version>
+ </dependency>
+ <dependency>
+ <groupId>org.apache.commons</groupId>
<artifactId>commons-numbers-arrays</artifactId>
<version>${project.version}</version>
</dependency>
diff --git a/src/main/resources/checkstyle/checkstyle-suppressions.xml
b/src/main/resources/checkstyle/checkstyle-suppressions.xml
index 4bb2be9..4937eac 100644
--- a/src/main/resources/checkstyle/checkstyle-suppressions.xml
+++ b/src/main/resources/checkstyle/checkstyle-suppressions.xml
@@ -24,6 +24,10 @@
<suppress checks="FileLengthCheck" files=".*[/\\]Complex(Test)?\.java" />
<suppress checks="FileLengthCheck"
files=".*jmh[/\\]core[/\\]LinearCombinations.java" />
<suppress checks="MethodLength" files=".*[/\\]Complex\.java" />
+ <!-- Naming and method logic preserved from the original Boost C++ code -->
+ <suppress checks="MethodLength" files=".*[/\\]BoostErf\.java" />
+ <suppress checks="LocalVariableName" files=".*[/\\]BoostErf\.java" />
+ <suppress checks="LocalFinalVariableName" files=".*[/\\]BoostErf\.java" />
<!-- Method "hashCode()" is defined in the "Angle" parent class; subclasses
have no additional fields to include in the hash. -->
<suppress checks="EqualsHashCode" files=".*[/\\]Angle\.java" />
diff --git a/src/main/resources/pmd/pmd-ruleset.xml
b/src/main/resources/pmd/pmd-ruleset.xml
index d86342e..4ba91b2 100644
--- a/src/main/resources/pmd/pmd-ruleset.xml
+++ b/src/main/resources/pmd/pmd-ruleset.xml
@@ -83,6 +83,13 @@
<property name="utilityClassPattern" value="[A-Z][a-zA-Z0-9]*" />
</properties>
</rule>
+ <rule ref="category/java/codestyle.xml/LocalVariableNamingConventions">
+ <properties>
+ <!-- Naming conventions preserved from the original Boost C++ code -->
+ <property name="violationSuppressXPath"
+ value="//ClassOrInterfaceDeclaration[@SimpleName='BoostErf']"/>
+ </properties>
+ </rule>
<rule ref="category/java/design.xml/CyclomaticComplexity">
<properties>
@@ -122,7 +129,8 @@
or @SimpleName='LogGammaSum'
or @SimpleName='InverseErf'
or @SimpleName='InvGamma1pm1'
- or @SimpleName='RegularizedBeta']"/>
+ or @SimpleName='RegularizedBeta'
+ or @SimpleName='BoostErf']"/>
</properties>
</rule>
@@ -141,7 +149,8 @@
<property name="violationSuppressXPath"
value="//ClassOrInterfaceDeclaration[@SimpleName='Complex'
or @SimpleName='BigFraction'
- or @SimpleName='Fraction']"/>
+ or @SimpleName='Fraction'
+ or @SimpleName='BoostErf']"/>
</properties>
</rule>
