BoostTools.java
/*
* 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;
import java.util.function.DoubleSupplier;
/**
* Utility tools used by the Boost functions.
*
* <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
* {@code c++} implementations in {@code <boost/math/tools/>}.
* All work is copyright John Maddock 2006 and subject to the Boost Software License.
*/
final class BoostTools {
/**
* The minimum epsilon value for relative error in the summation.
* Equal to Math.ulp(1.0) or 2^-52.
*
* <h2>Note</h2>
*
* <p>The summation will terminate when any additional terms are too small to
* change the sum. Assuming additional terms are reducing in magnitude this
* occurs when the term is 1 ULP from the sum:
* <pre>
* ulp(sum) >= term
* </pre>
*
* <p>The epsilon is used to set a configurable threshold using:
* <pre>
* sum * eps >= term
* </pre>
* <p>The minimum epsilon is the smallest value that will create a value
* {@code >= 1} ULP and {@code < 2} ULP of the sum. For any normal number the ULP
* of all values with the same exponent b is scalb(1.0, b - 52). This can
* be achieved by multiplication by 2^-52.
*/
private static final double EPSILON = 0x1.0p-52;
/**
* The minimum epsilon value for relative error in the Kahan summation.
* This can be lower than {@link #EPSILON}. Set to 2^-62.
*
* <h2>Note</h2>
*
* <p>The Kahan summation uses a carry term to extend the precision of the sum.
* This extends the 53-bit mantissa by adding more bits to hold round-off error.
* Thus the term may effect the sum when it has a magnitude smaller than 1 ULP
* of the current sum. The epsilon can be lowered from 2^-52 to include the
* extra bits in the convergence criteria. The lower limit for the epsilon is
* 2^-104. Boost uses an epsilon specified as a number of bits of accuracy. Each
* lowering of the epsilon by a factor of 2 adds a guard digit to the sum.
* Slower converging series will benefit from a lower epsilon. This uses 2^-62
* to add 10 guard digits and allow testing of different thresholds when the
* Kahan summation is used, for example in the gamma function. Lowering the
* epsilon further is only of benefit if the terms can be computed exactly.
* Otherwise the rounding errors of the early terms affect the final result as
* much as the inclusion of extra guard digits.
*/
private static final double KAHAN_EPSILON = 0x1.0p-62;
/** Message for failure to converge. */
private static final String MSG_FAILED_TO_CONVERGE = "Failed to converge within %d iterations";
/** Private constructor. */
private BoostTools() {
// intentionally empty.
}
/**
* Sum the series.
*
* <p>Adapted from {@code boost/math/tools/series.hpp}.
*
* @param func Series generator
* @param epsilon Maximum relative error allowed
* @param maxTerms Maximum number of terms
* @return result
*/
static double sumSeries(DoubleSupplier func, double epsilon, int maxTerms) {
return sumSeries(func, epsilon, maxTerms, 0);
}
/**
* Sum the series.
*
* <p>Adapted from {@code boost/math/tools/series.hpp}.
*
* @param func Series generator
* @param epsilon Maximum relative error allowed
* @param maxTerms Maximum number of terms
* @param initValue Initial value
* @return result
*/
static double sumSeries(DoubleSupplier func, double epsilon, int maxTerms, double initValue) {
// Note:
// The Boost code requires eps to be non-zero. It is created in the
// <boost/math/policies/policy.hpp> as a non-zero relative error term.
// An alternative termination condition with a divide is:
// (eps < Math.abs(nextTerm / result))
//
// Here the argument is checked against the minimum epsilon for a double
// to provide functional equivalence with the Boost policy.
// In the min eps case the loop terminates if the most recently added term is
// 0 or 1 ulp of the result. This condition is acceptable if the next
// computed term will be at most half of the most recent term (thus
// cannot be added to the current result).
final double eps = getEpsilon(epsilon, EPSILON);
int counter = maxTerms;
double result = initValue;
double nextTerm;
do {
nextTerm = func.getAsDouble();
result += nextTerm;
} while (Math.abs(eps * result) < Math.abs(nextTerm) && --counter > 0);
if (counter <= 0) {
throw new ArithmeticException(
String.format(MSG_FAILED_TO_CONVERGE, maxTerms));
}
return result;
}
/**
* Sum the series using Kahan summation.
*
* <p>Adapted from {@code boost/math/tools/series.hpp}.
*
* @param func Series generator
* @param epsilon Maximum relative error allowed
* @param maxTerms Maximum number of terms
* @return result
*/
static double kahanSumSeries(DoubleSupplier func, double epsilon, int maxTerms) {
return kahanSumSeries(func, epsilon, maxTerms, 0);
}
/**
* Sum the series using Kahan summation.
*
* <p>Adapted from {@code boost/math/tools/series.hpp}.
*
* @param func Series generator
* @param epsilon Maximum relative error allowed
* @param maxTerms Maximum number of terms
* @param initValue Initial value
* @return result
*/
static double kahanSumSeries(DoubleSupplier func, double epsilon, int maxTerms, double initValue) {
final double eps = getEpsilon(epsilon, KAHAN_EPSILON);
int counter = maxTerms;
// Kahan summation:
// https://en.wikipedia.org/wiki/Kahan_summation_algorithm
// This summation is accurate if the term is smaller in magnitude
// than the current sum. This is a condition required for the
// series termination thus the extended precision sum need not
// check magnitudes of terms to compute the carry.
double result = initValue;
double carry = 0;
double nextTerm;
do {
nextTerm = func.getAsDouble();
final double y = nextTerm - carry;
final double t = result + y;
carry = t - result;
carry -= y;
result = t;
} while (Math.abs(eps * result) < Math.abs(nextTerm) && --counter > 0);
if (counter <= 0) {
throw new ArithmeticException(
String.format(MSG_FAILED_TO_CONVERGE, maxTerms));
}
return result;
}
/**
* Gets the epsilon ensuring it satisfies the minimum allowed value.
*
* <p>This is returning the maximum of the two arguments.
* Do not use Math.max as it returns NaN if either value is NaN.
* In this case the desired result in the default minEpsilon, not NaN.
* Math.max will also check ordering when terms are equal to support
* -0.0 and 0.0. This does not apply here and a single conditional
* returns the desired result.
*
* @param epsilon Configured epsilon
* @param minEpsilon Minimum allowed epsilon
* @return the epsilon
*/
private static double getEpsilon(double epsilon, double minEpsilon) {
return epsilon > minEpsilon ? epsilon : minEpsilon;
}
/**
* Evaluate the polynomial using Horner's method.
* The coefficients are used in descending order, for example a polynomial of order
* 3 requires 4 coefficients:
* <pre>
* f(x) = c[3] * x^3 + c[2] * x^2 + c[1] * x + c[0]
* </pre>
*
* @param c Polynomial coefficients (must have length > 0)
* @param x Argument x
* @return polynomial value
*/
static double evaluatePolynomial(double[] c, double x) {
final int count = c.length;
double sum = c[count - 1];
for (int i = count - 2; i >= 0; --i) {
sum *= x;
sum += c[i];
}
return sum;
}
}