001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.numbers.fraction; 018 019import java.util.function.Supplier; 020 021/** 022 * Provides a means to evaluate 023 * <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">generalized continued fractions</a>. 024 * 025 * <p>The continued fraction uses the following form for the numerator ({@code a}) and 026 * denominator ({@code b}) coefficients: 027 * <pre> 028 * a1 029 * b0 + ------------------ 030 * b1 + a2 031 * ------------- 032 * b2 + a3 033 * -------- 034 * b3 + ... 035 * </pre> 036 * 037 * <p>A generator of the coefficients must be provided to evaluate the continued fraction. 038 * 039 * <p>The implementation of the fraction evaluation is based on the modified Lentz algorithm 040 * as described on page 508 in: 041 * 042 * <ul> 043 * <li> 044 * I. J. Thompson, A. R. Barnett (1986). 045 * "Coulomb and Bessel Functions of Complex Arguments and Order." 046 * Journal of Computational Physics 64, 490-509. 047 * <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf"> 048 * https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a> 049 * </li> 050 * </ul> 051 * 052 * @see <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">Wikipedia: Generalized continued fraction</a> 053 * @see <a href="https://en.wikipedia.org/wiki/Generalized_continued_fraction">MathWorld: Generalized continued fraction</a> 054 * @since 1.1 055 */ 056public final class GeneralizedContinuedFraction { 057 /** 058 * The value for any number close to zero. 059 * 060 * <p>"The parameter small should be some non-zero number less than typical values of 061 * eps * |b_n|, e.g., 1e-50". 062 */ 063 static final double SMALL = 1e-50; 064 /** Default maximum number of iterations. */ 065 static final int DEFAULT_ITERATIONS = Integer.MAX_VALUE; 066 /** 067 * Minimum relative error epsilon. Equal to 1 - Math.nextDown(1.0), or 2^-53. 068 * 069 * <p>The epsilon is used to compare the change in the magnitude of the fraction 070 * convergent to 1.0. In theory eps can be 2^-53 reflecting the smallest reduction in 071 * magnitude possible i.e. {@code next = previous * Math.nextDown(1.0)}, or zero 072 * reflecting exact convergence. 073 * 074 * <p>If set to zero then the algorithm requires exact convergence which may not be possible 075 * due to floating point error in the algorithm. For example the golden ratio will not 076 * converge. 077 * 078 * <p>The minimum value will stop the recursive evaluation at the smallest possible 079 * increase or decrease in the convergent. 080 */ 081 private static final double MIN_EPSILON = 0x1.0p-53; 082 /** Maximum relative error epsilon. This is configured to prevent incorrect usage. Values 083 * higher than 1.0 invalidate the relative error lower bound of {@code (1 - eps) / 1}. 084 * Set to 0.5 which is a very weak relative error tolerance. */ 085 private static final double MAX_EPSILON = 0.5; 086 /** Default low threshold for change in magnitude. Precomputed using MIN_EPSILON. 087 * Equal to 1 - 2^-53. */ 088 private static final double DEFAULT_LOW = 1 - MIN_EPSILON; 089 /** Default absolute difference threshold for change in magnitude. Precomputed using MIN_EPSILON. 090 * Equal to {@code 1 / (1 - 2^-53) = 2^-52}. */ 091 private static final double DEFAULT_EPS = 0x1.0p-52; 092 093 /** 094 * Defines the <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html"> 095 * {@code n}-th "a" and "b" coefficients</a> of the continued fraction. 096 * 097 * @since 1.1 098 */ 099 public static final class Coefficient { 100 /** "a" coefficient. */ 101 private final double a; 102 /** "b" coefficient. */ 103 private final double b; 104 105 /** 106 * @param a "a" coefficient 107 * @param b "b" coefficient 108 */ 109 private Coefficient(double a, double b) { 110 this.a = a; 111 this.b = b; 112 } 113 114 /** 115 * Returns the {@code n}-th "a" coefficient of the continued fraction. 116 * 117 * @return the coefficient <code>a<sub>n</sub></code>. 118 */ 119 public double getA() { 120 return a; 121 } 122 123 /** 124 * Returns the {@code n}-th "b" coefficient of the continued fraction. 125 * 126 * @return the coefficient <code>b<sub>n</sub></code>. 127 */ 128 public double getB() { 129 return b; 130 } 131 132 /** 133 * Create a new coefficient. 134 * 135 * @param a "a" coefficient 136 * @param b "b" coefficient 137 * @return the coefficient 138 */ 139 public static Coefficient of(double a, double b) { 140 return new Coefficient(a, b); 141 } 142 } 143 144 /** No instances. */ 145 private GeneralizedContinuedFraction() {} 146 147 /** 148 * Evaluates the continued fraction. 149 * 150 * <p>Note: The first generated partial numerator a<sub>0</sub> is discarded. 151 * 152 * @param gen Generator of coefficients. 153 * @return the value of the continued fraction. 154 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number of 155 * iterations is reached before the expected convergence is achieved. 156 * @see #value(Supplier,double,int) 157 */ 158 public static double value(Supplier<Coefficient> gen) { 159 return value(gen, MIN_EPSILON, DEFAULT_ITERATIONS); 160 } 161 162 /** 163 * Evaluates the continued fraction. 164 * 165 * <p>Note: The first generated partial numerator a<sub>0</sub> is discarded. 166 * 167 * @param gen Generator of coefficients. 168 * @param epsilon Maximum relative error allowed. 169 * @return the value of the continued fraction. 170 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number of 171 * iterations is reached before the expected convergence is achieved. 172 * @see #value(Supplier,double,int) 173 */ 174 public static double value(Supplier<Coefficient> gen, double epsilon) { 175 return value(gen, epsilon, DEFAULT_ITERATIONS); 176 } 177 178 /** 179 * Evaluates the continued fraction. 180 * <pre> 181 * a1 182 * b0 + ------------------ 183 * b1 + a2 184 * ------------- 185 * b2 + a3 186 * -------- 187 * b3 + ... 188 * </pre> 189 * 190 * <p>Setting coefficient a<sub>n</sub> to zero will signal the end of the recursive evaluation. 191 * 192 * <p>Note: The first generated partial numerator a<sub>0</sub> is discarded. 193 * 194 * <p><b>Usage Note</b> 195 * 196 * <p>This method is not functionally identical to calling 197 * {@link #value(double, Supplier, double, int)} with the generator configured to 198 * provide coefficients from n=1 and supplying b<sub>0</sub> separately. In some cases 199 * the computed result from the two variations may be different by more than the 200 * provided epsilon. The other method should be used if b<sub>0</sub> is zero or very 201 * small. See the corresponding javadoc for details. 202 * 203 * @param gen Generator of coefficients. 204 * @param epsilon Maximum relative error allowed. 205 * @param maxIterations Maximum number of iterations. 206 * @return the value of the continued fraction. 207 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number of 208 * iterations is reached before the expected convergence is achieved. 209 * @see #value(double, Supplier, double, int) 210 */ 211 public static double value(Supplier<Coefficient> gen, double epsilon, int maxIterations) { 212 // Use the first b coefficient to seed the evaluation of the fraction. 213 // Coefficient a is discarded. 214 final Coefficient c = gen.get(); 215 return evaluate(c.getB(), gen, epsilon, maxIterations); 216 } 217 218 /** 219 * Evaluates the continued fraction. 220 * 221 * <p>Note: The initial term b<sub>0</sub> is supplied as an argument. 222 * Both of the first generated terms a and b are used. This fraction evaluation 223 * can be used when: 224 * <ul> 225 * <li>b<sub>0</sub> is not part of a regular series 226 * <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component 227 * <li>b<sub>0</sub> is very small and the result is expected to approach zero 228 * </ul> 229 * 230 * @param b0 Coefficient b<sub>0</sub>. 231 * @param gen Generator of coefficients. 232 * @return the value of the continued fraction. 233 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number 234 * of iterations is reached before the expected convergence is achieved. 235 * @see #value(double,Supplier,double,int) 236 */ 237 public static double value(double b0, Supplier<Coefficient> gen) { 238 return value(b0, gen, MIN_EPSILON, DEFAULT_ITERATIONS); 239 } 240 241 /** 242 * Evaluates the continued fraction. 243 * 244 * <p>Note: The initial term b<sub>0</sub> is supplied as an argument. 245 * Both of the first generated terms a and b are used. This fraction evaluation 246 * can be used when: 247 * <ul> 248 * <li>b<sub>0</sub> is not part of a regular series 249 * <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component 250 * <li>b<sub>0</sub> is very small and the result is expected to approach zero 251 * </ul> 252 * 253 * @param b0 Coefficient b<sub>0</sub>. 254 * @param gen Generator of coefficients. 255 * @param epsilon Maximum relative error allowed. 256 * @return the value of the continued fraction. 257 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number 258 * of iterations is reached before the expected convergence is achieved. 259 * @see #value(double,Supplier,double,int) 260 */ 261 public static double value(double b0, Supplier<Coefficient> gen, double epsilon) { 262 return value(b0, gen, epsilon, DEFAULT_ITERATIONS); 263 } 264 265 /** 266 * Evaluates the continued fraction. 267 * <pre> 268 * a1 269 * b0 + ------------------ 270 * b1 + a2 271 * ------------- 272 * b2 + a3 273 * -------- 274 * b3 + ... 275 * </pre> 276 * 277 * <p>Setting coefficient a<sub>n</sub> to zero will signal the end of the recursive evaluation. 278 * 279 * <p>Note: The initial term b<sub>0</sub> is supplied as an argument. 280 * Both of the first generated terms a and b are used. This fraction evaluation 281 * can be used when: 282 * <ul> 283 * <li>b<sub>0</sub> is not part of a regular series 284 * <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component 285 * <li>b<sub>0</sub> is very small and the result is expected to approach zero 286 * </ul> 287 * 288 * <p><b>Usage Note</b> 289 * 290 * <p>This method is not functionally identical to calling 291 * {@link #value(Supplier, double, int)} with the generator configured to provide term 292 * "b<sub>0</sub>" in the first coefficient. In some cases the computed result from 293 * the two variations may be different by more than the provided epsilon. The 294 * convergence of the continued fraction algorithm relies on computing an update 295 * multiplier applied to the current value. Convergence is faster if the initial value 296 * is close to the final value. The {@link #value(Supplier, double, int)} method will 297 * initialise the current value using b<sub>0</sub> and evaluate the continued 298 * fraction using updates computed from the generated coefficients. This method 299 * initialises the algorithm using b1 to evaluate part of the continued fraction and 300 * computes the result as: 301 * 302 * <pre> 303 * a1 304 * b0 + ------ 305 * part 306 * </pre> 307 * 308 * <p>This is preferred if b<sub>0</sub> is smaller in magnitude than the continued 309 * fraction component. In particular the evaluation algorithm sets a bound on the 310 * minimum initial value as {@code 1e-50}. If b<sub>0</sub> is smaller than this value 311 * then using this method is the preferred evaluation. 312 * 313 * @param b0 Coefficient b<sub>0</sub>. 314 * @param gen Generator of coefficients. 315 * @param epsilon Maximum relative error allowed. 316 * @param maxIterations Maximum number of iterations. 317 * @return the value of the continued fraction. 318 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number 319 * of iterations is reached before the expected convergence is achieved. 320 * @see #value(Supplier,double,int) 321 */ 322 public static double value(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) { 323 // Use the first b coefficient to seed the evaluation of the fraction. 324 // Coefficient a is used to compute the final result as the numerator term a1. 325 // The supplied b0 is added to the result. 326 final Coefficient c = gen.get(); 327 return b0 + c.getA() / evaluate(c.getB(), gen, epsilon, maxIterations); 328 } 329 330 /** 331 * Evaluates the continued fraction using the modified Lentz algorithm described in 332 * Thompson and Barnett (1986) Journal of Computational Physics 64, 490-509. 333 * <pre> 334 * a1 335 * b0 + ------------------ 336 * b1 + a2 337 * ------------- 338 * b2 + a3 339 * -------- 340 * b3 + ... 341 * </pre> 342 * 343 * <p>Note: The initial term b<sub>0</sub> is supplied as an argument. 344 * Both of the first generated terms a and b are used. 345 * 346 * <p><b>Implementation Note</b> 347 * 348 * <p>This method is private and functionally different from 349 * {@link #value(double, Supplier, double, int)}. The convergence of the algorithm relies on 350 * computing an update multiplier applied to the current value, initialised as b0. Accuracy 351 * of the evaluation can be effected if the magnitude of b0 is very different from later 352 * terms. In particular if initialised as 0 the algorithm will not function and so must 353 * set b0 to a small non-zero number. The public methods with the leading b0 term 354 * provide evaluation of the fraction if the term b0 is zero. 355 * 356 * @param b0 Coefficient b<sub>0</sub>. 357 * @param gen Generator of coefficients. 358 * @param epsilon Maximum relative error allowed. 359 * @param maxIterations Maximum number of iterations. 360 * @return the value of the continued fraction. 361 * @throws ArithmeticException if the algorithm fails to converge or if the maximal number 362 * of iterations is reached before the expected convergence is achieved. 363 */ 364 static double evaluate(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) { 365 // Relative error epsilon should not be zero to prevent drift in the event 366 // that the update ratio never achieves 1.0. 367 368 // Epsilon is the relative change allowed from 1. Configure the absolute limits so 369 // convergence requires: low <= deltaN <= high 370 // low = 1 - eps 371 // high = 1 / (1 - eps) 372 // High is always further from 1 than low in absolute distance. Do not store high 373 // but store the maximum absolute deviation from 1 for convergence = high - 1. 374 // If this is achieved a second check is made against low. 375 double low; 376 double eps; 377 if (epsilon > MIN_EPSILON && epsilon <= MAX_EPSILON) { 378 low = 1 - epsilon; 379 eps = 1 / low - 1; 380 } else { 381 // Precomputed defaults. Used when epsilon <= MIN_EPSILON 382 low = DEFAULT_LOW; 383 eps = DEFAULT_EPS; 384 } 385 386 double hPrev = updateIfCloseToZero(b0); 387 388 // Notes from Thompson and Barnett: 389 // 390 // Fraction convergent: hn = An / Bn 391 // A(-1) = 1, A0 = b0, B(-1) = 0, B0 = 1 392 393 // Compute the ratios: 394 // Dn = B(n-1) / Bn = 1 / (an * D(n-1) + bn) 395 // Cn = An / A(n-1) = an / C(n-1) + bn 396 // 397 // Ratio of successive convergents: 398 // delta n = hn / h(n-1) 399 // = Cn / Dn 400 401 // Avoid divisors being zero (less than machine precision) by shifting them to e.g. 1e-50. 402 403 double dPrev = 0.0; 404 double cPrev = hPrev; 405 406 for (int n = maxIterations; n > 0; n--) { 407 final Coefficient c = gen.get(); 408 final double a = c.getA(); 409 final double b = c.getB(); 410 411 double dN = updateIfCloseToZero(b + a * dPrev); 412 final double cN = updateIfCloseToZero(b + a / cPrev); 413 414 dN = 1 / dN; 415 final double deltaN = cN * dN; 416 final double hN = hPrev * deltaN; 417 418 // If the fraction is convergent then deltaN -> 1. 419 // Computation of deltaN = 0 or deltaN = big will result in zero or overflow. 420 // Directly check for overflow on hN (this ensures the result is finite). 421 422 if (!Double.isFinite(hN)) { 423 throw new FractionException("Continued fraction diverged to " + hN); 424 } 425 426 // Check for underflow on deltaN. This allows fractions to compute zero 427 // if this is the convergent limit. 428 // Note: deltaN is only zero if dN > 1e-50 / min_value, or 2.02e273. 429 // Since dN is the ratio of convergent denominators this magnitude of 430 // ratio is a presumed to be an error. 431 if (deltaN == 0) { 432 throw new FractionException("Ratio of successive convergents is zero"); 433 } 434 435 // Update from Thompson and Barnett to use <= eps in place of < eps. 436 // eps = high - 1 437 // A second check is made to ensure: 438 // low <= deltaN <= high 439 if (Math.abs(deltaN - 1) <= eps && deltaN >= low) { 440 return hN; 441 } 442 443 dPrev = dN; 444 cPrev = cN; 445 hPrev = hN; 446 } 447 448 throw new FractionException("Maximum iterations (%d) exceeded", maxIterations); 449 } 450 451 /** 452 * Returns the value, or if close to zero returns a small epsilon of the same sign. 453 * 454 * <p>This method is used in Thompson & Barnett to monitor both the numerator and denominator 455 * ratios for approaches to zero. 456 * 457 * @param value the value 458 * @return the value (or small epsilon) 459 */ 460 private static double updateIfCloseToZero(double value) { 461 return Math.abs(value) < SMALL ? Math.copySign(SMALL, value) : value; 462 } 463}