main.js (7017B)
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 * ## Notice 20 * 21 * The following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/s_cbrt.c?view=markup}. The implementation follows the original, but has been modified for JavaScript. 22 * 23 * ```text 24 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 25 * 26 * Developed at SunPro, a Sun Microsystems, Inc. business. 27 * Permission to use, copy, modify, and distribute this 28 * software is freely granted, provided that this notice 29 * is preserved. 30 * 31 * Optimized by Bruce D. Evans. 32 * ``` 33 */ 34 35 'use strict'; 36 37 // MODULES // 38 39 var FLOAT64_SMALLEST_NORMAL = require( '@stdlib/constants/float64/smallest-normal' ); 40 var getHighWord = require( '@stdlib/number/float64/base/get-high-word' ); 41 var setHighWord = require( '@stdlib/number/float64/base/set-high-word' ); 42 var isinfinite = require( './../../../../base/assert/is-infinite' ); 43 var fromWords = require( '@stdlib/number/float64/base/from-words' ); 44 var toWords = require( '@stdlib/number/float64/base/to-words' ); 45 var isnan = require( './../../../../base/assert/is-nan' ); 46 var polyval = require( './polyval_p.js' ); 47 48 49 // VARIABLES // 50 51 // 0x80000000 = 2147483648 => 1 00000000000 00000000000000000000 52 var SIGN_MASK = 0x80000000>>>0; // asm type annotation 53 54 // 0x7fffffff = 2147483647 => 0 11111111111 11111111111111111111 55 var ABS_MASK = 0x7fffffff>>>0; // asm type annotation 56 57 // 2**32 - 1 = 4294967295 => 11111111111111111111111111111111 58 var HIGH_WORD_MASK = 4294967295>>>0; // asm type annotation 59 60 // 2**31 + 2**30 = 3221225472 => 11000000000000000000000000000000 61 var LOW_WORD_MASK = 3221225472>>>0; // asm type annotation 62 63 // 2**54 64 var TWO_54 = 18014398509481984.0; 65 66 // 2**31 = 0x80000000 = 2147483648 => 1 00000000000 00000000000000000000 67 var TWO_31 = 0x80000000>>>0; // asm type annotation 68 69 // 0x00000001 = 1 => 0 00000000000 00000000000000000001 70 var ONE = 0x00000001>>>0; // asm type annotation 71 72 // B1 = (1023-1023/3-0.03306235651)*2**20 73 var B1 = 715094163>>>0; // asm type annotation 74 75 // B2 = (1023-1023/3-54/3-0.03306235651)*2**20 76 var B2 = 696219795>>>0; // asm type annotation 77 78 // 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 79 var FLOAT64_SMALLEST_NORMAL_HIGH_WORD = getHighWord( FLOAT64_SMALLEST_NORMAL ); // eslint-disable-line id-length 80 81 // Words workspace: 82 var WORDS = [ 0>>>0, 0>>>0 ]; // asm type annotations 83 84 85 // MAIN // 86 87 /** 88 * Computes the cube root of a double-precision floating-point number. 89 * 90 * ## Method 91 * 92 * 1. Rough cube root to \\( 5 \\) bits: 93 * 94 * ```tex 95 * \sqrt\[3\]{2^e (1+m)} \approx 2^(e/3) \biggl(1 + \frac{(e \mathrm{mod}\ 3) + m}{3}\biggr) 96 * ``` 97 * 98 * where \\( e \\) is a nonnegative integer, \\( m \\) is real and in \\( \[0, 1) \\), and \\( / \\) and \\( \mathrm{mod} \\) are integer division and modulus with rounding toward \\( -\infty \\). 99 * 100 * The RHS is always greater than or equal to the LHS and has a maximum relative error of about \\( 1 \\) in \\( 16 \\). 101 * 102 * Adding a bias of \\( -0.03306235651 \\) to the \\( (e \mathrm{mod} 3+ m )/ 3 \\) term reduces the error to about \\( 1 \\) in \\( 32 \\). 103 * 104 * With the IEEE floating point representation, for finite positive normal values, ordinary integer division of the value in bits magically gives almost exactly the RHS of the above provided we first subtract the exponent bias (\\( 1023 \\) for doubles) and later add it back. 105 * 106 * We do the subtraction virtually to keep \\( e \geq 0 \\) so that ordinary integer division rounds toward \\( -\infty \\); this is also efficient. 107 * 108 * 2. New cube root to \\( 23 \\) bits: 109 * 110 * ```tex 111 * \sqrt[3]{x} = t \cdot \sqrt\[3\]{x/t^3} \approx t \mathrm{P}(t^3/x) 112 * ``` 113 * 114 * where \\( \mathrm{P}(r) \\) is a polynomial of degree \\( 4 \\) that approximates \\( 1 / \sqrt\[3\]{r} \\) to within \\( 2^{-23.5} \\) when \\( |r - 1| < 1/10 \\). 115 * 116 * The rough approximation has produced \\( t \\) such than \\( |t/sqrt\[3\]{x} - 1| \lesssim 1/32 \\), and cubing this gives us bounds for \\( r = t^3/x \\). 117 * 118 * 3. Round \\( t \\) away from \\( 0 \\) to \\( 23 \\) bits (sloppily except for ensuring that the result is larger in magnitude than \\( \sqrt\[3\]{x} \\) but not much more than \\( 2 \\) 23-bit ulps larger). 119 * 120 * With rounding toward zero, the error bound would be \\( \approx 5/6 \\) instead of \\( \approx 4/6 \\). 121 * 122 * With a maximum error of \\( 2 \\) 23-bit ulps in the rounded \\( t \\), the infinite-precision error in the Newton approximation barely affects the third digit in the final error \\( 0.667 \\); the error in the rounded \\( t \\) can be up to about \\( 3 \\) 23-bit ulps before the final error is larger than \\( 0.667 \\) ulps. 123 * 124 * 4. Perform one step of a Newton iteration to get \\( 53 \\) bits with an error of \\( < 0.667 \\) ulps. 125 * 126 * 127 * @param {number} x - input value 128 * @returns {number} cube root 129 * 130 * @example 131 * var v = cbrt( 64.0 ); 132 * // returns 4.0 133 * 134 * @example 135 * var v = cbrt( 27.0 ); 136 * // returns 3.0 137 * 138 * @example 139 * var v = cbrt( 0.0 ); 140 * // returns 0.0 141 * 142 * @example 143 * var v = cbrt( -9.0 ); 144 * // returns ~-2.08 145 * 146 * @example 147 * var v = cbrt( NaN ); 148 * // returns NaN 149 */ 150 function cbrt( x ) { 151 var sgn; 152 var hx; 153 var hw; 154 var r; 155 var s; 156 var t; 157 var w; 158 if ( 159 x === 0.0 || // handles +-0 160 isnan( x ) || 161 isinfinite( x ) 162 ) { 163 return x; 164 } 165 hx = getHighWord( x )>>>0; 166 sgn = (hx & SIGN_MASK)>>>0; 167 hx &= ABS_MASK; 168 169 // Rough cbrt... 170 if ( hx < FLOAT64_SMALLEST_NORMAL_HIGH_WORD ) { 171 t = TWO_54 * x; 172 hw = ( getHighWord( t )&ABS_MASK )>>>0; 173 hw = ( ( (hw/3)>>>0 ) + B2 )>>>0; 174 t = fromWords( sgn|hw, 0 ); 175 } else { 176 t = 0.0; 177 hw = ( ( (hx/3)>>>0 ) + B1 )>>>0; 178 t = setHighWord( t, sgn|hw ); 179 } 180 // New cbrt... 181 r = ( t*t ) * ( t/x ); 182 t *= polyval( r ); 183 184 // Round `t` away from `0` to `23` bits... 185 toWords( WORDS, t ); 186 if ( WORDS[ 1 ]&TWO_31 ) { 187 // Perform manual addition, since we are split across two words... 188 WORDS[ 0 ] += ONE; // carry the one 189 WORDS[ 1 ] &= ~TWO_31; // clear the bit 190 } else { 191 WORDS[ 1 ] |= TWO_31; 192 } 193 t = fromWords( WORDS[0]&HIGH_WORD_MASK, WORDS[1]&LOW_WORD_MASK ); 194 195 // Newton iteration... 196 s = t * t; // `t*t` is exact 197 r = x / s; // error `<= 0.5` ulps; `|r| < |t|` 198 w = t + t; // `t+t` is exact 199 r = ( r-t ) / ( w+r ); // `r-t` is exact; `w+r ~= 3*t` 200 t += t * r; // error `<= 0.5 + 0.5/3 + eps` 201 202 return t; 203 } 204 205 206 // EXPORTS // 207 208 module.exports = cbrt;