commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From aherb...@apache.org
Subject [commons-numbers] branch master updated: Updated ContinuedFraction and the test.
Date Mon, 06 Apr 2020 10:53:23 GMT
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 a2ce6c7  Updated ContinuedFraction and the test.
a2ce6c7 is described below

commit a2ce6c77d5fa90385b34b9b95bf46a7e7c98efc9
Author: aherbert <aherbert@apache.org>
AuthorDate: Mon Apr 6 11:53:18 2020 +0100

    Updated ContinuedFraction and the test.
    
    Re-factor the close-to-zero logic into a method.
    
    Increase test coverage to 100% with the exception of an approach to zero
    test.
    
    Document the representation of A and B coefficients.
    This is currently the inverse of the original paper and the wolfram
    reference.
---
 .../numbers/fraction/ContinuedFraction.java        | 60 ++++++++++-----
 .../numbers/fraction/ContinuedFractionTest.java    | 88 +++++++++++++++++++++-
 2 files changed, 128 insertions(+), 20 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 cd98311..4801318 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
@@ -23,9 +23,28 @@ import org.apache.commons.numbers.core.Precision;
  * <a href="https://mathworld.wolfram.com/ContinuedFraction.html">continued fractions</a>.
  * Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
  * coefficients to evaluate the continued fraction.
+ *
+ * <p>The fraction uses the following form for the {@code a} and {@code b} coefficients:
+ * <pre>
+ *              b1
+ * a0 + ------------------
+ *      a1 +      b2
+ *           -------------
+ *           a2 +    b3
+ *                --------
+ *                a3 + ...
+ * </pre>
  */
 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;
+
+    /**
      * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
      * {@code n}-th "a" coefficient</a> of the continued fraction.
      *
@@ -65,14 +84,16 @@ public abstract class ContinuedFraction {
      * Evaluates the continued fraction.
      * <p>
      * The implementation of this method is based on the modified Lentz algorithm as described
-     * on page 18 ff. in:
+     * on page 508 in:
      * </p>
      *
      * <ul>
      *   <li>
-     *   I. J. Thompson,  A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments
and Order."
-     *   <a target="_blank" href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
-     *   http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
+     *   I. J. Thompson,  A. R. Barnett (1986).
+     *   "Coulomb and Bessel Functions of Complex Arguments and Order."
+     *   Journal of Computational Physics 64, 490-509.
+     *   <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
+     *   https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
      *   </li>
      * </ul>
      *
@@ -85,13 +106,7 @@ public abstract class ContinuedFraction {
      * before the expected convergence is achieved.
      */
     public double evaluate(double x, double epsilon, int maxIterations) {
-        final double small = 1e-50;
-        double hPrev = getA(0, x);
-
-        // use the value of small as epsilon criteria for zero checks
-        if (Precision.equals(hPrev, 0.0, small)) {
-            hPrev = small;
-        }
+        double hPrev = updateIfCloseToZero(getA(0, x));
 
         int n = 1;
         double dPrev = 0.0;
@@ -102,14 +117,8 @@ public abstract class ContinuedFraction {
             final double a = getA(n, x);
             final double b = getB(n, x);
 
-            double dN = a + b * dPrev;
-            if (Precision.equals(dN, 0.0, small)) {
-                dN = small;
-            }
-            double cN = a + b / cPrev;
-            if (Precision.equals(cN, 0.0, small)) {
-                cN = small;
-            }
+            double dN = updateIfCloseToZero(a + b * dPrev);
+            final double cN = updateIfCloseToZero(a + b / cPrev);
 
             dN = 1 / dN;
             final double deltaN = cN * dN;
@@ -136,4 +145,17 @@ public abstract class ContinuedFraction {
 
         throw new FractionException("maximal count ({0}) exceeded", maxIterations);
     }
+
+    /**
+     * 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;
+    }
 }
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 a6fe628..ebb03ae 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
@@ -16,6 +16,7 @@
  */
 package org.apache.commons.numbers.fraction;
 
+import java.util.Locale;
 import org.junit.jupiter.api.Assertions;
 import org.junit.jupiter.api.Test;
 
@@ -45,6 +46,41 @@ public class ContinuedFractionTest {
     }
 
     @Test
+    public void test415Over93() throws Exception {
+        // https://en.wikipedia.org/wiki/Continued_fraction
+        // 415             1
+        // ---  = 4 + ---------
+        //  93        2 +   1
+        //                -----
+        //                6 + 1
+        //                    -
+        //                    7
+        //      = [4; 2, 6, 7]
+
+        ContinuedFraction cf = new ContinuedFraction() {
+            @Override
+            public double getA(int n, double x) {
+                switch (n) {
+                    case 0: return 4;
+                    case 1: return 2;
+                    case 2: return 6;
+                    case 3: return 7;
+                    default: return 1;
+                }
+            }
+
+            @Override
+            public double getB(int n, double x) {
+                return n <= 3 ? 1 : 0;
+            }
+        };
+
+        final double eps = 1e-8;
+        double gr = cf.evaluate(0, eps, 5);
+        Assertions.assertEquals(415.0 / 93.0, gr, eps);
+    }
+
+    @Test
     public void testMaxIterationsThrows() throws Exception {
         ContinuedFraction cf = new ContinuedFraction() {
             @Override
@@ -60,7 +96,57 @@ public class ContinuedFractionTest {
 
         final double eps = 1e-8;
         final int maxIterations = 3;
-        Assertions.assertThrows(FractionException.class, () -> cf.evaluate(0, eps, maxIterations));
+        final Throwable t = Assertions.assertThrows(FractionException.class,
+            () -> cf.evaluate(0, eps, maxIterations));
+        assertExceptionMessageContains(t, "max");
+    }
+
+    @Test
+    public void testNaNThrows() throws Exception {
+        // Create a NaN during the iteration
+        ContinuedFraction cf = new ContinuedFraction() {
+            @Override
+            public double getA(int n, double x) {
+                return n == 0 ? 1 : Double.NaN;
+            }
+
+            @Override
+            public double getB(int n, double x) {
+                return 1;
+            }
+        };
+
+        final double eps = 1e-8;
+        final Throwable t = Assertions.assertThrows(FractionException.class,
+            () -> cf.evaluate(0, eps, 5));
+        assertExceptionMessageContains(t, "nan");
+    }
+
+    @Test
+    public void testInfThrows() throws Exception {
+        // Create an infinity during the iteration:
+        // b / cPrev  => b_1 / a_0 => Double.MAX_VALUE / 0.5
+        ContinuedFraction cf = new ContinuedFraction() {
+            @Override
+            public double getA(int n, double x) {
+                return 0.5;
+            }
+
+            @Override
+            public double getB(int n, double x) {
+                return n == 0 ? 1 : Double.MAX_VALUE;
+            }
+        };
+
+        final double eps = 1e-8;
+        final Throwable t = Assertions.assertThrows(FractionException.class,
+            () -> cf.evaluate(0, eps, 5));
+        assertExceptionMessageContains(t, "infinity");
+    }
+
+    private static void assertExceptionMessageContains(Throwable t, String text) {
+        Assertions.assertTrue(t.getMessage().toLowerCase(Locale.ROOT).contains(text),
+            () -> "Missing '" + text + "' from exception message: " + t.getMessage());
     }
 
     // NUMBERS-46


Mime
View raw message