lambdaeta.js (2609B)
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 exp = require( './../../../../base/special/exp' ); 25 var ln = require( './../../../../base/special/ln' ); 26 var evalpoly = require( './../../../../base/tools/evalpoly' ); 27 var polyvalAK1 = require( './polyval_ak1.js' ); 28 var polyvalAK2 = require( './polyval_ak2.js' ); 29 30 31 // VARIABLES // 32 33 var THRESHOLD = 1.0e-8; 34 var ONEO12 = 0.0833333333333333333333333333333; 35 var ONEO120 = 0.00833333333333333333333333333333; 36 37 // Polynomial coefficient workspace: 38 var AK = [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]; // WARNING: not thread safe 39 40 41 // MAIN // 42 43 /** 44 * Returns the positive number satisfying \\( \eta^2/2=\lambda-1-\ln(\lambda) \\) with \\( \operatorname{sign}(\lambda-1)=\operatorname{sign}(\eta) \\). 45 * 46 * @private 47 * @param {number} eta - eta value 48 * @returns {number} value satisfying equation 49 */ 50 function lambdaeta( eta ) { 51 var L2; 52 var L3; 53 var L4; 54 var L5; 55 var la; 56 var L; 57 var q; 58 var r; 59 var s; 60 61 s = eta * eta * 0.5; 62 if ( eta === 0.0 ) { 63 la = 0.0; 64 } 65 else if ( eta < -1.0 ) { 66 r = exp( -1.0 - s ); 67 la = polyvalAK1( r ); 68 } 69 else if ( eta < 1.0 ) { 70 r = eta; 71 la = polyvalAK2( r ); 72 } 73 else { 74 r = 11.0 + s; 75 L = ln( r ); 76 la = r + L; 77 r = 1.0 / r; 78 L2 = L * L; 79 L3 = L2 * L; 80 L4 = L3 * L; 81 L5 = L4 * L; 82 AK[ 1 ] = ( 2.0-L ) * 0.5; 83 AK[ 2 ] = ( ( -9.0*L ) + 6.0 + ( 2.0*L2 ) ) / 6.0; 84 AK[ 3 ] = -( (3*L3)+ (36*L) - (22*L2) - 12 ) * ONEO12; 85 AK[ 4 ] = ( 60.0 + (350.0*L2) - (300.0*L) - (125.0*L3) + (12.0*L4) ) / 60.0; // eslint-disable-line max-len 86 AK[ 5 ] = -(-120 - (274*L4) + (900*L) - (1700*L2) + (1125*L3) + (20*L5)) * ONEO120; // eslint-disable-line max-len 87 la += ( L * r * evalpoly( AK, r ) ); 88 } 89 r = 1.0; 90 if ( 91 ( eta > -3.5 && eta < -0.03 ) || 92 ( eta > 0.03 && eta < 40.0 ) 93 ) { 94 r = 1.0; 95 q = la; 96 do { 97 la = q * ( s+ln(q) ) / ( q-1.0 ); 98 r = abs( ( q/la ) - 1.0 ); 99 q = la; 100 } while ( r > THRESHOLD ); 101 } 102 return la; 103 } 104 105 106 // EXPORTS // 107 108 module.exports = lambdaeta;