basic.js (4715B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2018 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 'use strict'; 20 21 // MODULES // 22 23 var abs = require( './../../../../base/special/abs' ); 24 var EPS = require( '@stdlib/constants/float64/eps' ); 25 var FLOAT32_SMALLEST_NORMAL = require( '@stdlib/constants/float32/smallest-normal' ); 26 27 28 // VARIABLES // 29 30 var MAX_ITER = 1000000; 31 32 33 // FUNCTIONS // 34 35 /** 36 * Evaluates a continued fraction expansion. 37 * 38 * ```text 39 * a1 40 * --------------- 41 * b1 + a2 42 * ---------- 43 * b2 + a3 44 * ----- 45 * b3 + ... 46 * ``` 47 * 48 * @private 49 * @param {Function} gen - function giving terms of continued fraction expansion 50 * @param {PositiveNumber} factor - further terms are only added as long as factor*result is smaller than the next term 51 * @param {PositiveInteger} maxIter - maximum number of iterations 52 * @returns {number} evaluated expansion 53 */ 54 function continuedFractionA( gen, factor, maxIter ) { 55 var delta; 56 var a0; 57 var C; 58 var D; 59 var f; 60 var v; 61 62 v = gen(); 63 f = v[ 1 ]; 64 a0 = v[ 0 ]; 65 if ( f === 0 ) { 66 f = FLOAT32_SMALLEST_NORMAL; 67 } 68 C = f; 69 D = 0.0; 70 71 do { 72 v = gen(); 73 if ( v ) { 74 D = v[ 1 ] + ( v[ 0 ] * D ); 75 if ( D === 0.0 ) { 76 D = FLOAT32_SMALLEST_NORMAL; 77 } 78 C = v[ 1 ] + ( v[ 0 ] / C ); 79 if ( C === 0.0 ) { 80 C = FLOAT32_SMALLEST_NORMAL; 81 } 82 D = 1.0 / D; 83 delta = C * D; 84 f *= delta; 85 } 86 } while ( v && ( abs( delta - 1.0 ) > factor ) && --maxIter ); // eslint-disable-line no-plusplus 87 88 return a0 / f; 89 } 90 91 /** 92 * Evaluates a continued fraction expansion. 93 * 94 * ```text 95 * b0 + a1 96 * --------------- 97 * b1 + a2 98 * ---------- 99 * b2 + a3 100 * ----- 101 * b3 + ... 102 * ``` 103 * 104 * @private 105 * @param {Function} gen - function giving terms of continued fraction expansion 106 * @param {PositiveNumber} factor - further terms are only added as long as factor*result is smaller than the next term 107 * @param {PositiveInteger} maxIter - maximum number of iterations 108 * @returns {number} evaluated expansion 109 */ 110 function continuedFractionB( gen, factor, maxIter ) { 111 var delta; 112 var C; 113 var D; 114 var f; 115 var v; 116 117 v = gen(); 118 f = v[ 1 ]; 119 if ( f === 0.0 ) { 120 f = FLOAT32_SMALLEST_NORMAL; 121 } 122 C = f; 123 D = 0.0; 124 do { 125 v = gen(); 126 if ( v ) { 127 D = v[ 1 ] + ( v[ 0 ] * D ); 128 if ( D === 0.0 ) { 129 D = FLOAT32_SMALLEST_NORMAL; 130 } 131 C = v[ 1 ] + ( v[ 0 ] / C ); 132 if ( C === 0.0 ) { 133 C = FLOAT32_SMALLEST_NORMAL; 134 } 135 D = 1.0 / D; 136 delta = C * D; 137 f *= delta; 138 } 139 } while ( v && ( abs( delta - 1.0 ) > factor ) && --maxIter ); // eslint-disable-line no-plusplus 140 return f; 141 } 142 143 144 // MAIN // 145 146 /** 147 * Evaluates the continued fraction approximation for the supplied series generator using the modified Lentz algorithm. 148 * 149 * ## References 150 * 151 * - Lentz, William J. 1976. "Generating bessel functions in Mie scattering calculations using continued fractions." _Applied Optics_ 15 (3): 668–71. doi:[10.1364/AO.15.000668](https://doi.org/10.1364/AO.15.000668). 152 * 153 * @param {Function} generator - function returning terms of continued fraction expansion 154 * @param {Object} [options] - function options 155 * @param {PositiveInteger} [options.maxIter=1000000] - maximum number of iterations 156 * @param {PositiveNumber} [options.tolerance=2.22e-16] - further terms are only added as long as the next term is greater than current term times the tolerance 157 * @param {boolean} [options.keep=false] - whether to keep the leading b term 158 * @returns {number} value of continued fraction 159 * 160 * @example 161 * // Continued fraction for (e-1)^(-1): 162 * var gen = generator(); 163 * var out = continuedFraction( gen ); 164 * // returns ~0.582 165 * 166 * function generator() { 167 * var i = 0; 168 * return function() { 169 * i++; 170 * return [ i, i ]; 171 * }; 172 * } 173 */ 174 function continuedFraction( generator, options ) { 175 var maxIter; 176 var opts; 177 var eps; 178 179 opts = {}; 180 if ( arguments.length > 1 ) { 181 opts = options; 182 } 183 eps = opts.tolerance || EPS; 184 maxIter = opts.maxIter || MAX_ITER; 185 186 if ( opts.keep ) { 187 return continuedFractionB( generator, eps, maxIter ); 188 } 189 return continuedFractionA( generator, eps, maxIter ); 190 } 191 192 193 // EXPORTS // 194 195 module.exports = continuedFraction;