You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/11/29 17:00:00 UTC
[commons-numbers] 04/04: NUMBERS-176: ContinuedFraction delegates to GeneralizedContinuedFraction
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 23ba4f8ea6109c9bac7b4ec43e17dcaaed358bd9
Author: aherbert <ah...@apache.org>
AuthorDate: Wed Nov 24 14:15:11 2021 +0000
NUMBERS-176: ContinuedFraction delegates to GeneralizedContinuedFraction
---
.../numbers/fraction/ContinuedFraction.java | 111 ++++++++++-----------
.../fraction/GeneralizedContinuedFraction.java | 8 +-
.../numbers/fraction/ContinuedFractionTest.java | 29 +++++-
src/changes/changes.xml | 11 +-
4 files changed, 90 insertions(+), 69 deletions(-)
diff --git a/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/ContinuedFraction.java b/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/ContinuedFraction.java
index 38d1708..3001ac6 100644
--- a/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/ContinuedFraction.java
+++ b/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/ContinuedFraction.java
@@ -16,7 +16,8 @@
*/
package org.apache.commons.numbers.fraction;
-import org.apache.commons.numbers.core.Precision;
+import java.util.function.Supplier;
+import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;
/**
* Provides a generic means to evaluate
@@ -36,19 +37,15 @@ import org.apache.commons.numbers.core.Precision;
*
* <p>Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
* coefficients to evaluate the continued fraction.
+ *
+ * <p>This class allows evaluation of the fraction for a specified evaluation point {@code x};
+ * the point can be used to express the values of the coefficients.
+ * Evaluation of a continued fraction from a generator of the coefficients can be performed using
+ * {@link GeneralizedContinuedFraction}. This may be preferred if the coefficients can be computed
+ * with updates to the previous coefficients.
*/
public abstract class ContinuedFraction {
/**
- * The value for any number close to zero.
- *
- * <p>"The parameter small should be some non-zero number less than typical values of
- * eps * |b_n|, e.g., 1e-50".
- */
- private static final double SMALL = 1e-50;
- /** Minimum allowed relative error epsilon. Equal to Math.ulp(1.0). */
- private static final double MIN_EPSILON = 0x1.0p-52;
-
- /**
* Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
* {@code n}-th "a" coefficient</a> of the continued fraction.
*
@@ -81,7 +78,7 @@ public abstract class ContinuedFraction {
* @see #evaluate(double,double,int)
*/
public double evaluate(double x, double epsilon) {
- return evaluate(x, epsilon, Integer.MAX_VALUE);
+ return evaluate(x, epsilon, GeneralizedContinuedFraction.DEFAULT_ITERATIONS);
}
/**
@@ -110,56 +107,52 @@ public abstract class ContinuedFraction {
* before the expected convergence is achieved.
*/
public double evaluate(double x, double epsilon, int maxIterations) {
- // Relative error epsilon must not be zero.
- // Do not use Math.max as NaN would be returned for a NaN epsilon.
- final double eps = epsilon > MIN_EPSILON ? epsilon : MIN_EPSILON;
-
- double hPrev = updateIfCloseToZero(getB(0, x));
-
- int n = 1;
- double dPrev = 0.0;
- double cPrev = hPrev;
- double hN;
-
- while (n <= maxIterations) {
- final double a = getA(n, x);
- final double b = getB(n, x);
-
- double dN = updateIfCloseToZero(b + a * dPrev);
- final double cN = updateIfCloseToZero(b + a / cPrev);
-
- dN = 1 / dN;
- final double deltaN = cN * dN;
- hN = hPrev * deltaN;
-
- if (!Double.isFinite(hN)) {
- throw new FractionException(
- "Continued fraction diverged to %s for value %s", hN, x);
- }
-
- if (Math.abs(deltaN - 1) < eps) {
- return hN;
+ // Delegate to GeneralizedContinuedFraction
+
+ // Get the first coefficient
+ final double b0 = getB(0, x);
+
+ // Generate coefficients from (a1,b1)
+ final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
+ private int n;
+ @Override
+ public Coefficient get() {
+ n++;
+ final double a = getA(n, x);
+ final double b = getB(n, x);
+ return Coefficient.of(a, b);
}
-
- dPrev = dN;
- cPrev = cN;
- hPrev = hN;
- ++n;
+ };
+
+ // Invoke appropriate method based on magnitude of first term.
+
+ // If b0 is too small or zero it is set to a non-zero small number to allow
+ // magnitude updates. Avoid this by adding b0 at the end if b0 is small.
+ //
+ // This handles the use case of a negligible initial term. If b1 is also small
+ // then the evaluation starting at b0 or b1 may converge poorly.
+ // One solution is to manually compute the convergent until it is not small
+ // and then evaluate the fraction from the next term:
+ // h1 = b0 + a1 / b1
+ // h2 = b0 + a1 / (b1 + a2 / b2)
+ // ...
+ // hn not 'small', start generator at (n+1):
+ // value = GeneralizedContinuedFraction.value(hn, gen)
+ // This solution is not implemented to avoid recursive complexity.
+
+ if (Math.abs(b0) < GeneralizedContinuedFraction.SMALL) {
+ // Updates from initial convergent b1 and computes:
+ // b0 + a1 / [ b1 + a2 / (b2 + ... ) ]
+ return GeneralizedContinuedFraction.value(b0, gen, epsilon, maxIterations);
}
- throw new FractionException("Maximum iterations (%d) exceeded", maxIterations);
- }
+ // Use the package-private evaluate method.
+ // Calling GeneralizedContinuedFraction.value(gen, epsilon, maxIterations)
+ // requires the generator to start from (a0,b0) and repeats computation of b0
+ // and wastes computation of a0.
- /**
- * Returns the value, or if close to zero returns a small epsilon.
- *
- * <p>This method is used in Thompson & Barnett to monitor both the numerator and denominator
- * ratios for approaches to zero.
- *
- * @param value the value
- * @return the value (or small epsilon)
- */
- private static double updateIfCloseToZero(double value) {
- return Precision.equals(value, 0.0, SMALL) ? SMALL : value;
+ // Updates from initial convergent b0:
+ // b0 + a1 / (b1 + ... )
+ return GeneralizedContinuedFraction.evaluate(b0, gen, epsilon, maxIterations);
}
}
diff --git a/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/GeneralizedContinuedFraction.java b/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/GeneralizedContinuedFraction.java
index 3619240..43dc48b 100644
--- a/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/GeneralizedContinuedFraction.java
+++ b/commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/GeneralizedContinuedFraction.java
@@ -59,7 +59,9 @@ public final class GeneralizedContinuedFraction {
* <p>"The parameter small should be some non-zero number less than typical values of
* eps * |b_n|, e.g., 1e-50".
*/
- private static final double SMALL = 1e-50;
+ static final double SMALL = 1e-50;
+ /** Default maximum number of iterations. */
+ static final int DEFAULT_ITERATIONS = Integer.MAX_VALUE;
/**
* Minimum relative error epsilon. Equal to 1 - Math.nextDown(1.0), or 2^-53.
*
@@ -86,8 +88,6 @@ public final class GeneralizedContinuedFraction {
/** Default absolute difference threshold for change in magnitude. Precomputed using MIN_EPSILON.
* Equal to {@code 1 / (1 - 2^-53) = 2^-52}. */
private static final double DEFAULT_EPS = 0x1.0p-52;
- /** Default maximum number of iterations. */
- private static final int DEFAULT_ITERATIONS = Integer.MAX_VALUE;
/**
* Defines the <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">
@@ -358,7 +358,7 @@ public final class GeneralizedContinuedFraction {
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number
* of iterations is reached before the expected convergence is achieved.
*/
- private static double evaluate(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) {
+ static double evaluate(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) {
// Relative error epsilon should not be zero to prevent drift in the event
// that the update ratio never achieves 1.0.
diff --git a/commons-numbers-fraction/src/test/java/org/apache/commons/numbers/fraction/ContinuedFractionTest.java b/commons-numbers-fraction/src/test/java/org/apache/commons/numbers/fraction/ContinuedFractionTest.java
index b9c21e0..3bfccc2 100644
--- a/commons-numbers-fraction/src/test/java/org/apache/commons/numbers/fraction/ContinuedFractionTest.java
+++ b/commons-numbers-fraction/src/test/java/org/apache/commons/numbers/fraction/ContinuedFractionTest.java
@@ -78,6 +78,31 @@ class ContinuedFractionTest {
}
/**
+ * Test a continued fraction where the leading term is zero.
+ * Evaluates the reciprocal of the golden ratio.
+ */
+ @Test
+ void testGoldenRatioReciprocal() throws Exception {
+ final double eps = 1e-8;
+ final ContinuedFraction cf = new ContinuedFraction() {
+ @Override
+ public double getA(int n, double x) {
+ // Check this is not called with n=0
+ Assertions.assertNotEquals(0, n, "a0 should never require evaluation");
+ return 1;
+ }
+
+ @Override
+ public double getB(int n, double x) {
+ // b0 = 0
+ return n == 0 ? 0 : 1;
+ }
+ };
+ final double gr = cf.evaluate(0, eps);
+ Assertions.assertEquals(1 / GOLDEN_RATIO, gr, eps / GOLDEN_RATIO);
+ }
+
+ /**
* Test an invalid epsilon (zero, negative or NaN).
* See NUMBERS-173.
*
@@ -191,7 +216,7 @@ class ContinuedFractionTest {
// NUMBERS-46
@Test
void testOneIteration() {
- final double eps = 10;
+ final double eps = 0.5;
final double gr = GoldenRatio.getInstance().evaluate(0, eps, 1);
// Expected: 1 + 1 / 1
Assertions.assertEquals(2.0, gr);
@@ -200,7 +225,7 @@ class ContinuedFractionTest {
// NUMBERS-46
@Test
void testTwoIterations() {
- final double eps = 0.5;
+ final double eps = 0.25;
final double gr = GoldenRatio.getInstance().evaluate(0, eps, 2);
// Expected: 1 + 1 / (1 + 1 / 1)
Assertions.assertEquals(1.5, gr);
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 04424f6..954c4f3 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -72,17 +72,20 @@ N.B. the Performance testing module requires Java 9+.
(The unit tests require Java 8+)
">
- <action dev="aherbert" type="update" issue="174">
+ <action dev="aherbert" type="update" issue="NUMBERS-176">
+ "ContinuedFraction": Update to use a shared implementation with GeneralizedContinuedFraction.
+ </action>
+ <action dev="aherbert" type="update" issue="NUMBERS-174">
"Gamma/LogGamma/RegularizedGamma": Update the gamma function implementations to
increase accuracy and support for extreme values.
Functionality is ported from the Boost C++ library.
</action>
- <action dev="aherbert" type="add" issue="175">
+ <action dev="aherbert" type="add" issue="NUMBERS-175">
"GeneralizedContinuedFraction": A continued fraction class to compute using a generator.
Allows evaluation of continued fractions from a regular series where coefficients can
be computed iteratively from the previous coefficients.
</action>
- <action dev="aherbert" type="fix" issue="173">
+ <action dev="aherbert" type="fix" issue="NUMBERS-173">
"ContinuedFraction": Set a minimum bound on the relative error epsilon. Prevents
an infinite loop when the epsilon is zero.
</action>
@@ -101,7 +104,7 @@ N.B. the Performance testing module requires Java 9+.
minimum value from 2^-53 to double min value. Execution speed is improved.
Functionality is ported from the Boost C++ library.
</action>
- <action dev="aherbert" type="fix" issue="170">
+ <action dev="aherbert" type="fix" issue="NUMBERS-170">
"RegularizedBeta": Detect edge cases for arguments that can be evaluated by
exploiting properties of the regularized beta function.
</action>