[
https://issues.apache.org/jira/browse/NUMBERS-175?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17448057#comment-17448057
]
Alex Herbert commented on NUMBERS-175:
--------------------------------------
I used the GeneralizedContinuedFraction implementation in the revised gamma
functions (see NUMBERS-174). The implementation led to different, and often
slightly worse, results. This was tracked to the convergence of the continued
fraction and drift of the result by 1 ULP.
The modified Lentz algorithm computes the value by evaluating changes to
successive convergents (h).
{noformat}
h_{n+1} = h_n * delta{noformat}
A convergent is a rational number (P/Q) that approximates the value of the
continued fraction. If successive updates (delta) to the convergent approach 1
then the convergent is approaching a stable limit. The convergent limit is the
value returned by the algorithm.
For a simple continued fraction Wikipedia states that "Even-numbered
convergents are smaller than the original number, while odd-numbered ones are
larger". So the update multiplier should alternate between above 1 and below 1
until the updates are exactly 1. This is observed for various simple continued
fractions in the test suite. Here is the golden ratio updates until convergence
error of 1e-6:
{noformat}
1.0
2.0
0.75
1.111111111111111
0.9600000000000002
1.015625
0.9940828402366865
1.0022675736961453
0.9991349480968859
1.0003305785123966
0.9998737533139754
1.000048225308642
0.9999815800622593
1.0000070358617874
0.9999973125503896
1.0000010265159331
0.9999996079060259
{noformat}
Generalized continued fractions may not alternate the update around 1. Here is
the result for a generalized fraction for tan(x=0.5) until convergence error of
2^-52:
{noformat}
1.0
0.9166666666666666
0.9984591679506933
0.9999887082204155
0.9999999546818229
0.9999999998848089
0.9999999999997977
0.9999999999999998{noformat}
Note that the iteration is converging with smaller and smaller updates but all
are below 1. The issue here is that the accuracy of numbers below 1 is twice as
high as those above 1. With the lowest possible tolerance the algorithm should
stop when the update adjusted the convergent by 1 ULP. The algorithm was coded
as:
{code:java}
if (Math.abs(deltaN - 1) < eps) {
// stop
}{code}
So the error tolerance should be 2^-53 when the update is below 1 and 2^-52
when the update is above 1. Previously the eps was set to a minimum of 2^-52.
This sets the stopping criteria as an increase of 0 ULP or decrease of 1 ULP:
{noformat}
0.9999999999999998 < deltaN < 1.0000000000000002
or
0.9999999999999999 <= deltaN <= 1.0{noformat}
This means that if the last update was computed as 1.0000000000000002 then
another iteration will be performed. Ideally the next update would be
0.9999999999999999 or 1. Unfortunately due to floating point error the next
update can be computed again as 1.0000000000000002 and the convergent drifts
upwards.
If the eps is set to 2^-53 the same can occur with downward drift. Here the
algorithm is evaluating the golden ratio and the update delta never reaches 1:
{noformat}
1.0
2.0
0.75
1.111111111111111
0.9600000000000002
1.015625
0.9940828402366865
1.0022675736961453
0.9991349480968859
1.0003305785123966
0.9998737533139754
...
1.0000000000000047
0.9999999999999981
1.0000000000000007
0.9999999999999999
0.9999999999999999
0.9999999999999999
0.9999999999999999
...{noformat}
So the algorithm should be reconfigured to stop at the epsilon, not below it:
{code:java}
// Change '<' to '<='
if (Math.abs(deltaN - 1) <= eps) {
// stop
} {code}
Supporting a minimum eps of 2^-53 allows convergence when the updates approach
1 from less than 1. But if the convergence is from above then 2^-52 is the
limit. A value of 2^-53 is too small and may never be achieved.
The eps can be considered the absolute deviation from 1 for convergence:
{noformat}
1 - eps <= deltaN <= 1.0 + eps{noformat}
Or it can be considered the relative error from 1. The low bound is (1 - eps) /
1 or simply 1 - eps. The upper bound can be set using the inverse.
{noformat}
1 - eps <= deltaN <= 1.0 / (1 - eps){noformat}
This sets the following bounds:
||eps||eps||1-eps||ulps from 1||1+eps||ulps from 1||1 / (1 - eps)||ulps from 1||
|2^-53|1.1102230246251565E-16|0.9999999999999999|-1|1.0|0|1.0000000000000002|1|
|2^-52|2.220446049250313E-16|0.9999999999999998|-2|1.0000000000000002|1|1.0000000000000002|1|
|2^-51|4.440892098500626E-16|0.9999999999999996|-4|1.0000000000000004|2|1.0000000000000004|2|
|2^-50|8.881784197001252E-16|0.9999999999999991|-8|1.0000000000000009|4|1.0000000000000009|4|
|2^-49|1.7763568394002505E-15|0.9999999999999982|-16|1.0000000000000018|8|1.0000000000000018|8|
|2^-48|3.552713678800501E-15|0.9999999999999964|-32|1.0000000000000036|16|1.0000000000000036|16|
|2^-47|7.105427357601002E-15|0.9999999999999929|-64|1.000000000000007|32|1.000000000000007|32|
|2^-46|1.4210854715202004E-14|0.9999999999999858|-128|1.0000000000000142|64|1.0000000000000142|64|
|2^-45|2.8421709430404007E-14|0.9999999999999716|-256|1.0000000000000284|128|1.0000000000000284|128|
|2^-44|5.6843418860808015E-14|0.9999999999999432|-512|1.0000000000000568|256|1.0000000000000568|256|
|2^-43|1.1368683772161603E-13|0.9999999999998863|-1024|1.0000000000001137|512|1.0000000000001137|512|
|2^-42|2.2737367544323206E-13|0.9999999999997726|-2048|1.0000000000002274|1024|1.0000000000002274|1024|
|2^-41|4.547473508864641E-13|0.9999999999995453|-4096|1.0000000000004547|2048|1.0000000000004547|2048|
|2^-40|9.094947017729282E-13|0.9999999999990905|-8192|1.0000000000009095|4096|1.0000000000009095|4096|
|2^-39|1.8189894035458565E-12|0.999999999998181|-16384|1.000000000001819|8192|1.000000000001819|8192|
|2^-38|3.637978807091713E-12|0.999999999996362|-32768|1.000000000003638|16384|1.000000000003638|16384|
|2^-37|7.275957614183426E-12|0.999999999992724|-65536|1.000000000007276|32768|1.000000000007276|32768|
|2^-36|1.4551915228366852E-11|0.9999999999854481|-131072|1.000000000014552|65536|1.000000000014552|65536|
|2^-35|2.9103830456733704E-11|0.9999999999708962|-262144|1.0000000000291038|131072|1.0000000000291038|131072|
|2^-34|5.820766091346741E-11|0.9999999999417923|-524288|1.0000000000582077|262144|1.0000000000582077|262144|
|2^-33|1.1641532182693481E-10|0.9999999998835847|-1048576|1.0000000001164153|524288|1.0000000001164153|524288|
|2^-32|2.3283064365386963E-10|0.9999999997671694|-2097152|1.0000000002328306|1048576|1.0000000002328306|1048576|
|2^-31|4.6566128730773926E-10|0.9999999995343387|-4194304|1.0000000004656613|2097152|1.0000000004656613|2097152|
|2^-30|9.313225746154785E-10|0.9999999990686774|-8388608|1.0000000009313226|4194304|1.0000000009313226|4194304|
|2^-29|1.862645149230957E-9|0.9999999981373549|-16777216|1.0000000018626451|8388608|1.0000000018626451|8388608|
|2^-28|3.725290298461914E-9|0.9999999962747097|-33554432|1.0000000037252903|16777216|1.0000000037252903|16777216|
|2^-27|7.450580596923828E-9|0.9999999925494194|-67108864|1.0000000074505806|33554432|1.0000000074505806|33554432|
|2^-26|1.4901161193847656E-8|0.9999999850988388|-134217728|1.0000000149011612|67108864|1.0000000149011614|67108865|
|2^-25|2.9802322387695312E-8|0.9999999701976776|-268435456|1.0000000298023224|134217728|1.0000000298023233|134217732|
|2^-24|5.9604644775390625E-8|0.9999999403953552|-536870912|1.0000000596046448|268435456|1.0000000596046483|268435472|
|2^-23|1.1920928955078125E-7|0.9999998807907104|-1073741824|1.0000001192092896|536870912|1.0000001192093038|536870976|
|2^-22|2.384185791015625E-7|0.9999997615814209|-2147483648|1.000000238418579|1073741824|1.000000238418636|1073742080|
|2^-21|4.76837158203125E-7|0.9999995231628418|-4294967296|1.0000004768371582|2147483648|1.0000004768373856|2147484672|
|2^-20|9.5367431640625E-7|0.9999990463256836|-8589934592|1.0000009536743164|4294967296|1.000000953675226|4294971392|
|2^-19|1.9073486328125E-6|0.9999980926513672|-17179869184|1.0000019073486328|8589934592|1.0000019073522708|8589950976|
|2^-18|3.814697265625E-6|0.9999961853027344|-34359738368|1.0000038146972656|17179869184|1.0000038147118175|17179934720|
|2^-17|7.62939453125E-6|0.9999923706054688|-68719476736|1.0000076293945312|34359738368|1.0000076294527394|34360000514|
|2^-16|1.52587890625E-5|0.9999847412109375|-137438953472|1.0000152587890625|68719476736|1.0000152590218967|68720525328|
|2^-15|3.0517578125E-5|0.999969482421875|-274877906944|1.000030517578125|137438953472|1.000030518509476|137443147904|
|2^-14|6.103515625E-5|0.99993896484375|-549755813888|1.00006103515625|274877906944|1.0000610388817677|274894685184|
|2^-13|1.220703125E-4|0.9998779296875|-1099511627776|1.0001220703125|549755813888|1.0001220852154804|549822930945|
|2^-12|2.44140625E-4|0.999755859375|-2199023255552|1.000244140625|1099511627776|1.0002442002442002|1099780128784|
|2^-11|4.8828125E-4|0.99951171875|-4398046511104|1.00048828125|2199023255552|1.0004885197850513|2200097521920|
|2^-10|9.765625E-4|0.9990234375|-8796093022208|1.0009765625|4398046511104|1.0009775171065494|4402345676804|
|2^-9|0.001953125|0.998046875|-17592186044416|1.001953125|8796093022208|1.0019569471624266|8813306511488|
|2^-8|0.00390625|0.99609375|-35184372088832|1.00390625|17592186044416|1.003921568627451|17661175009296|
|2^-7|0.0078125|0.9921875|-70368744177664|1.0078125|35184372088832|1.0078740157480315|35461414388744|
|2^-6|0.015625|0.984375|-140737488355328|1.015625|70368744177664|1.0158730158730158|71485708370960|
|2^-5|0.03125|0.96875|-281474976710656|1.03125|140737488355328|1.032258064516129|145277407334532|
|2^-4|0.0625|0.9375|-562949953421312|1.0625|281474976710656|1.0666666666666667|300239975158033|
|2^-3|0.125|0.875|-1125899906842624|1.125|562949953421312|1.1428571428571428|643371375338642|
|2^-2|0.25|0.75|-2251799813685248|1.25|1125899906842624|1.3333333333333333|1501199875790165|
|2^-1|0.5|0.5|-4503599627370496|1.5|2251799813685248|2.0|4503599627370496|
Note that this sets a similar upper bound when eps is small. When the limit of
eps is 2^-53 the reciprocal method sets an allowed ULP of 1 in both directions.
The old method has 1 and 0 ULP.
When eps is large then the bounds are quite different. However since this is
relative change then the bounds are correct. A relative change of 0.5 would
allow a value to half or double.
The algorithm can be updated to compute the upper and lower bound on delta for
convergence and stop when the relative error is achieved. This allows the
minimum eps of 2^-53 to support 1 ulp changes in both directions.
> 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)
