[
https://issues.apache.org/jira/browse/NUMBERS-175?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17445445#comment-17445445
]
Alex Herbert edited comment on NUMBERS-175 at 11/22/21, 2:36 PM:
-----------------------------------------------------------------
I have started on an implementation using the generator approach.
Work in progress here [PR
109|https://github.com/apache/commons-numbers/pull/109].
I found some issues with the current implementation.
h3. It cannot compute values approaching zero with a small first term
Consider the evaluation of tan(z):
{noformat}
z
tan(z) = ---------------
1 - z^2
-----------
3 - z^2
-------
5 - ...
{noformat}
Here the leading term is zero. The current implementation must set this to a
non-zero small number as an estimate of the first convergent h0. The remaining
convergents h1, h2, h3, etc are computed using an update to the first.
If the expected value is far from the initial non-zero number then this works.
Otherwise the result is incorrect:
{code:java}
@Test
void test() {
// Evaluate tan(x)
ContinuedFraction cf = new ContinuedFraction() {
@Override
protected double getA(int n, double x) {
return n == 1 ? x : -x * x;
}
@Override
protected double getB(int n, double x) {
return n == 0 ? 0 : 2 * n - 1;
}
};
for (double x : new double[] {0.1, 1e-2, 1e-4, 1e-6, 1e-50, 1e-100,
1e-200}) {
try {
System.out.printf("%25s %25s%n", Math.tan(x), cf.evaluate(x,
1e-10));
} catch (ArithmeticException ex) {
// failed
}
}
}
{code}
Outputs:
{noformat}
0.10033467208545055 0.10033467208545055
0.010000333346667207 0.010000333346667207
1.0000000033333334E-4 1.0000000033333334E-4
1.0000000000003333E-6 1.0000000000003333E-6
1.0E-50 2.0E-50
1.0E-100 1.0E-50
1.0E-200 1.0E-50
{noformat}
This is not a major issue as you can evaluate the fraction advanced by 1
iteration and then divide the original first numerator coefficient a1 by the
evaluated fraction:
{code:java}
@Test
void test() {
// Evaluate tan(x)
ContinuedFraction cf = new ContinuedFraction() {
@Override
protected double getA(int n, double x) {
return -x * x;
}
@Override
protected double getB(int n, double x) {
return 2 * n + 1;
}
};
for (double x : new double[] {0.1, 1e-2, 1e-4, 1e-6, 1e-50, 1e-100,
1e-200}) {
try {
System.out.printf("%25s %25s%n", Math.tan(x), x / cf.evaluate(x,
1e-10));
} catch (ArithmeticException ex) {
// failed
}
}
}
{code}
This is correct:
{noformat}
0.10033467208545055 0.10033467208545054
0.010000333346667207 0.010000333346667205
1.0000000033333334E-4 1.0000000033333334E-4
1.0000000000003333E-6 1.0000000000003333E-6
1.0E-50 1.0E-50
1.0E-100 1.0E-100
1.0E-200 1.0E-200
{noformat}
h3. It does not detect an update multiplier that is zero
Consider this very bad fraction using max value for the numerator after a few
terms:
{noformat}
1
0.5 + ---------------------
0.5 + max
---------------
0.5 + max
---------
0.5 + ...
{noformat}
This will infinite loop:
{code:java}
@Test
void test() {
// This will create an update deltaN of zero.
// The algorithm gets stuck in a continuous loop.
ContinuedFraction cf = new ContinuedFraction() {
@Override
protected double getA(int n, double x) {
return n < 2 ? 1 : Double.MAX_VALUE;
}
@Override
protected double getB(int n, double x) {
return 0.5;
}
};
final double x = 1; // Ignored
Assertions.assertTimeoutPreemptively(Duration.ofSeconds(5),
() -> cf.evaluate(x, 1e-10));
}
{code}
The infinite loop is not detected due to the use of:
{code:java}
while (n <= maxIterations) {
{code}
With maxIterations set to Integer.MAX_VALUE this will loop forever. This
condition is mentioned in NUMBERS-46 but it was not fixed so that the algorithm
can enter at least 1 iteration.
This is a bug. I've fixed it in the PR above for the new generalized continued
fraction. An appropriate fix should be implemented in the ContinuedFraction
class.
was (Author: alexherbert):
I have started on an implementation using the generator approach.
Work in progress here [PR
109|https://github.com/apache/commons-numbers/pull/109].
I found some issues with the current implementation.
h3. It cannot compute values approaching zero with a small first term
Consider the evaluation of tan(z):
{noformat}
z
tan(z) = ---------------
1 - z^2
-----------
3 - z^2
-------
5 - ...
{noformat}
Here the leading term is zero. The current implementation must set this to a
non-zero small number as an estimate of the first convergent h0. The remaining
convergents h1, h2, h3, etc are computed using an update to the first.
If the expected value is far from the initial non-zero number then this works.
Otherwise the result is incorrect:
{code:java}
@Test
void test() {
// Evaluate tan(x)
ContinuedFraction cf = new ContinuedFraction() {
@Override
protected double getA(int n, double x) {
return n == 1 ? x : -x * x;
}
@Override
protected double getB(int n, double x) {
return n == 0 ? 0 : 2 * n - 1;
}
};
for (double x : new double[] {0.1, 1e-2, 1e-4, 1e-6, 1e-50, 1e-100,
1e-200}) {
try {
System.out.printf("%25s %25s%n", Math.tan(x), cf.evaluate(x,
1e-10));
} catch (ArithmeticException ex) {
// failed
}
}
}
{code}
Outputs:
{noformat}
0.10033467208545055 0.10033467208545055
0.010000333346667207 0.010000333346667207
1.0000000033333334E-4 1.0000000033333334E-4
1.0000000000003333E-6 1.0000000000003333E-6
1.0E-50 2.0E-50
1.0E-100 1.0E-50
1.0E-200 1.0E-50
{noformat}
This is not a major issue as you can evaluate the fraction advanced by 1
iteration and then divide the original first numerator coefficient a1 by the
evaluated fraction:
{code:java}
@Test
void test() {
// Evaluate tan(x)
ContinuedFraction cf = new ContinuedFraction() {
@Override
protected double getA(int n, double x) {
return -x * x;
}
@Override
protected double getB(int n, double x) {
return 2 * n + 1;
}
};
for (double x : new double[] {0.1, 1e-2, 1e-4, 1e-6, 1e-50, 1e-100,
1e-200}) {
try {
System.out.printf("%25s %25s%n", Math.tan(x), x / cf.evaluate(x,
1e-10));
} catch (ArithmeticException ex) {
// failed
}
}
}
{code}
This is correct:
{noformat}
0.10033467208545055 0.10033467208545054
0.010000333346667207 0.010000333346667205
1.0000000033333334E-4 1.0000000033333334E-4
1.0000000000003333E-6 1.0000000000003333E-6
1.0E-50 1.0E-50
1.0E-100 1.0E-100
1.0E-200 1.0E-200
{noformat}
h3. It does not detect an update multiplier than is zero
Consider this very bad fraction using max value for the numerator after a few
terms:
{noformat}
1
0.5 + ---------------------
0.5 + max
---------------
0.5 + max
---------
0.5 + ...
{noformat}
This will infinite loop:
{code:java}
@Test
void test() {
// This will create an update deltaN of zero.
// The algorithm gets stuck in a continuous loop.
ContinuedFraction cf = new ContinuedFraction() {
@Override
protected double getA(int n, double x) {
return n < 2 ? 1 : Double.MAX_VALUE;
}
@Override
protected double getB(int n, double x) {
return 0.5;
}
};
final double x = 1; // Ignored
Assertions.assertTimeoutPreemptively(Duration.ofSeconds(5),
() -> cf.evaluate(x, 1e-10));
}
{code}
The infinite loop is not detected due to the use of:
{code:java}
while (n <= maxIterations) {
{code}
With maxIterations set to Integer.MAX_VALUE this will loop forever. This
condition is mentioned in NUMBERS-46 but it was not fixed so that the algorithm
can enter at least 1 iteration.
This is a bug. I've fixed it in the PR above for the new generalized continued
fraction. An appropriate fix should be implemented in the ContinuedFraction
class.
> Add continued fraction implementations using a generator of terms
> -----------------------------------------------------------------
>
> Key: NUMBERS-175
> URL: https://issues.apache.org/jira/browse/NUMBERS-175
> Project: Commons Numbers
> Issue Type: Improvement
> Components: fraction
> Affects Versions: 1.0
> Reporter: Alex Herbert
> Priority: Minor
> Time Spent: 10m
> Remaining Estimate: 0h
>
> The ContinuedFraction class allows computation of:
> {noformat}
> b0 + a1
> ------------------
> b1 + a2
> -------------
> b2 + a3
> --------
> b3 + ...{noformat}
> This is done using an abstract class that is extended to implement the
> methods to get the terms:
> {code:java}
> double getA(int n, double x);
> double getB(int n, double x);{code}
> This allows the fraction to be reused to generate results for different
> points to evaluate. However:
> * It does not lend itself to fractions where the terms can be computed using
> recursion:
> {noformat}
> b(n+1) = f( b(n) ) {noformat}
> * It requires two method calls to generate terms a_n and b_n for each
> iteration thus preventing optimisation of the computation using the input n
> for values shared between computation of a and b.
> An alternative method is to support a generator of the paired terms a_n and
> b_n:
> {code:java}
> static double continuedFraction(Supplier<double[]> gen);{code}
> To be used in a single evaluation as:
> {code:java}
> Supplier<double[]> goldenRatio = () -> return new double[] {1, 1};
> double gr = continuedFraction(goldenRatio);
> // gr = 1.61803398874...{code}
> An additional feature is to support a simple continued fraction where all
> partial numerators are 1:
> {noformat}
> b0 + 1
> ------------------
> b1 + 1
> -------------
> b2 + 1
> --------
> b3 + ...
> {noformat}
> E.g. using:
> {code:java}
> static double simpleContinuedFraction(DoubleSupplier gen); {code}
> h3. Addition
> In some situations it is an advantage to not evaluate the leading term b0.
> The term may not be part of a regular series, or may be zero:
> {noformat}
> a1
> ------------------
> b1 + a2
> -------------
> b2 + a3
> --------
> b3 + ...
> {noformat}
> h3. API
> Using the nomeclature from Wikipedia and Wolfram suggests the following:
> {code:java}
> public static final class GeneralizedContinuedFraction {
> /**
> * Evaluate the continued fraction from the generator for partial
> numerator
> * a and partial denominator b.
> * <pre>
> * b0 + a1
> * ------------------
> * b1 + a2
> * -------------
> * b2 + a3
> * --------
> * b3 + ...
> * </pre>
> *
> * <p>Note: The first generated partial numerator a0 is discarded.
> *
> * @param gen Generator
> * @return the continued fraction value
> */
> public static double value(Supplier<double[]> gen);
> /**
> * Evaluate the continued fraction from the generator for partial
> numerator
> * a and partial denominator b.
> * <pre>
> * a1
> * ------------------
> * b1 + a2
> * -------------
> * b2 + a3
> * --------
> * b3 + ...
> * </pre>
> *
> * <p>Note: Both of the first terms a and b are used.
> *
> * @param gen Generator
> * @return the continued fraction value
> */
> public static double valueA(Supplier<double[]> gen);
> }
> public static final class SimpleContinuedFraction {
> public static double value(DoubleSupplier gen);
> public static double valueA(DoubleSupplier gen);
> }
> {code}
> The API should support the optional parameters to control the convergence
> tolerance and the maximum number of iterations.
> The variant to evaluate without the leading b0 term is not essential. It can
> be evaluated by starting the generator at the next iteration using:
> {code:java}
> Supplier<double[]> gen = ...; // Start at terms (a1,b1)
> // Evaluation will discard a1;
> // this value (separately computed) is then used to compute the result
> double value = a1 / GeneralizedContinuedFraction.value(gen);{code}
> This variant is present in the Boost c++ library and used to evaluate terms
> for the gamma function.
> See:
> [https://en.wikipedia.org/wiki/Continued_fraction]
> [https://mathworld.wolfram.com/SimpleContinuedFraction.html]
> [https://mathworld.wolfram.com/GeneralizedContinuedFraction.html]
> [Boost Continued Fraction (v
> 1_77_0)|https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/internals/cf.html]
>
--
This message was sent by Atlassian Jira
(v8.20.1#820001)
