1818*
1919* ## Notice
2020*
21- * The following copyright and license 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_pow .c}. The implementation follows the original, but has been modified for JavaScript.
21+ * The following copyright and license were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/e_powf .c}. The implementation follows the original, but has been modified for JavaScript.
2222*
2323* ```text
2424* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
@@ -39,57 +39,50 @@ var isInfinitef = require( '@stdlib/math/base/assert/is-infinitef' );
3939var isIntegerf = require ( '@stdlib/math/base/assert/is-integerf' ) ;
4040var sqrtf = require ( '@stdlib/math/base/special/sqrtf' ) ;
4141var absf = require ( '@stdlib/math/base/special/absf' ) ;
42- var toWords = require ( '@stdlib/number/float32/base/to-words ' ) ;
43- var setLowWord = require ( '@stdlib/number/float32/base/set-low -word' ) ;
44- var uint32ToInt32 = require ( '@stdlib/number/uint32 /base/to-int32 ' ) ;
42+ var fromWordf = require ( '@stdlib/number/float32/base/from-word ' ) ;
43+ var toWordf = require ( '@stdlib/number/float32/base/to -word' ) ;
44+ var float64ToFloat32 = require ( '@stdlib/number/float64 /base/to-float32 ' ) ;
4545var NINF = require ( '@stdlib/constants/float32/ninf' ) ;
4646var PINF = require ( '@stdlib/constants/float32/pinf' ) ;
47- var ABS_MASK = require ( '@stdlib/constants/float32/high-word-abs-mask' ) ;
4847var FLOAT32_ABS_MASK = require ( '@stdlib/constants/float32/abs-mask' ) ;
4948var xIsZero = require ( './x_is_zero.js' ) ;
50- var yIsHuge = require ( './y_is_huge.js' ) ;
5149var yIsInfinite = require ( './y_is_infinite.js' ) ;
52- var log2axf = require ( './log2axf .js' ) ;
53- var logxf = require ( './logxf .js' ) ;
54- var pow2f = require ( './pow2f .js' ) ;
50+ var log2ax = require ( './log2ax .js' ) ;
51+ var logx = require ( './logx .js' ) ;
52+ var pow2 = require ( './pow2 .js' ) ;
5553
5654
5755// VARIABLES //
5856
59- // 0x3fefffff = 1072693247 => 0 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1
60- var HIGH_MAX_NEAR_UNITY = 0x3fefffff | 0 ; // asm type annotation
57+ // 0x3f7ffff8 = 1065353216 => 0 01111111011 11111111111111111000 => biased exponent: 126 = -1+127 => 2^-1 * (1.999999)
58+ var HIGH_MAX_NEAR_UNITY = 0x3f7ffff8 ; // asm type annotation
6159
62- // 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31
63- var HIGH_BIASED_EXP_31 = 0x41e00000 | 0 ; // asm type annotation
60+ // 0x4d000000 = 1291845632 => 0 10011010000 00000000000000000000 => biased exponent: 208 = 81+127 => 2^81
61+ var HIGH_BIASED_EXP_27 = 0x4d000000 ; // asm type annotation
6462
65- // 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64
66- var HIGH_BIASED_EXP_64 = 0x43f00000 | 0 ; // asm type annotation
63+ // 0x43000000 = 1124073472 => 0 10000110000 00000000000000000000 => biased exponent: 96 = -31+127 => 2^16
64+ var SMALL_BIASED_EXP_10 = 0x43000000 ; // asm type annotation
6765
68- // 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024
69- var HIGH_BIASED_EXP_10 = 0x40900000 | 0 ; // asm type annotation
66+ // 0x3f800007 = 1065353223 => 0 01111111100 00000000000000000111 => biased exponent: 127 = 0+127 => 2^0 * (1.0000000007)
67+ var BIASED_EXP_0 = 0x3f800007 ; // asm type annotation
7068
71- // 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1
72- var HIGH_BIASED_EXP_0 = 0x3ff00000 | 0 ; // asm type annotation
69+ // 0xfffff000 = 4294963200 => 1 11111111111 11111111111100000000
70+ var HIGH_BIASED_EXP_NEG_64 = 0xfffff000 ; // asm type annotation
7371
74- // 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000
75- var HIGH_1075 = 0x4090cc00 | 0 ; // asm type annotation
72+ // 0x43160000 = 1128275968 => 0 10000110011 00000000000000000000 => biased exponent: 99 = -28+127 => 2^22
73+ var SMALL_60 = 0x43160000 ; // asm type annotation
7674
77- // 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000
78- var HIGH_NEG_1075 = 0xc090cc00 >>> 0 ; // asm type annotation
79-
80- var HIGH_NUM_NONSIGN_BITS = 31 | 0 ; // asm type annotation
75+ // 0xc3160000 = 3273056256 => 1 10000110011 00000000000000000000 => biased exponent: 99 = -28+127 => -2^22
76+ var SMALL_NEG_60 = 0xc3160000 >>> 0 ; // asm type annotation
8177
8278var HUGE = 1.0e30 ;
8379var TINY = 1.0e-30 ;
8480
8581// -(1024-log2(ovfl+.5ulp))
8682var OVT = 4.2995665694e-08 ;
8783
88- // High/low words workspace:
89- var WORDS = [ 0 | 0 , 0 | 0 ] ;
90-
9184// Log workspace:
92- var LOG_WORKSPACE = [ 0.0 , 0.0 ] ;
85+ var LOG_WORKSPACE = 0.0 ;
9386
9487
9588// MAIN //
@@ -194,64 +187,55 @@ var LOG_WORKSPACE = [ 0.0, 0.0 ];
194187* // returns NaN
195188*/
196189function powf ( x , y ) {
197- var ahx ; // absolute value high word `x`
198- var ahy ; // absolute value high word `y`
199- var ax ; // absolute value `x`
200- var hx ; // high word `x`
201- var lx ; // low word `x`
202- var hy ; // high word `y`
203- var ly ; // low word `y`
204- var sx ; // sign `x`
205- var sy ; // sign `y`
190+ var ahs ;
191+ var ahx ;
192+ var ahy ;
193+ var ax ;
194+ var hx ;
195+ var hy ;
196+ var hp ;
197+ var lp ;
198+ var sn ;
206199 var y1 ;
200+ var t1 ;
201+ var t2 ;
207202 var hp ;
208203 var lp ;
209204 var z ; // y prime
210205 var j ;
211- var i ;
212- if ( isnanf ( x ) || isnanf ( y ) ) {
213- return NaN ;
214- }
215- // Split `y` into high and low words:
216- toWords . assign ( y , WORDS , 1 , 0 ) ;
217- hy = WORDS [ 0 ] ;
218- ly = WORDS [ 1 ] ;
206+
207+ hx = fromWordf ( x ) ;
208+ hy = fromWordf ( y ) ;
219209
220210 // Special cases `y`...
221- if ( ly === 0 ) {
222- if ( y === 0.0 ) {
223- return 1.0 ;
224- }
225- if ( y === 1.0 ) {
226- return x ;
227- }
228- if ( y === - 1.0 ) {
229- return 1.0 / x ;
230- }
231- if ( y === 0.5 ) {
232- return sqrtf ( x ) ;
233- }
234- if ( y === - 0.5 ) {
235- return 1.0 / sqrtf ( x ) ;
236- }
237- if ( y === 2.0 ) {
238- return x * x ;
239- }
240- if ( y === 3.0 ) {
241- return x * x * x ;
242- }
243- if ( y === 4.0 ) {
244- x *= x ;
245- return x * x ;
246- }
247- if ( isInfinitef ( y ) ) {
248- return yIsInfinite ( x , y ) ;
249- }
211+ if ( y === 0.0 ) {
212+ return 1.0 ;
213+ }
214+ if ( y === 1.0 ) {
215+ return x ;
216+ }
217+ if ( y === - 1.0 ) {
218+ return 1.0 / x ;
219+ }
220+ if ( y === 0.5 ) {
221+ return sqrtf ( x ) ;
222+ }
223+ if ( y === - 0.5 ) {
224+ return 1.0 / sqrtf ( x ) ;
225+ }
226+ if ( y === 2.0 ) {
227+ return x * x ;
228+ }
229+ if ( y === 3.0 ) {
230+ return x * x * x ;
231+ }
232+ if ( y === 4.0 ) {
233+ x *= x ;
234+ return x * x ;
235+ }
236+ if ( isInfinitef ( y ) ) {
237+ return yIsInfinite ( x , y ) ;
250238 }
251- // Split `x` into high and low words:
252- toWords . assign ( x , WORDS , 1 , 0 ) ;
253- hx = WORDS [ 0 ] ;
254- lx = WORDS [ 1 ] ;
255239
256240 // Special cases `x`...
257241 if ( lx === 0 ) {
@@ -291,68 +275,60 @@ function powf( x, y ) {
291275 ahx = ( hx & FLOAT32_ABS_MASK ) ; // asm type annotation
292276 ahy = ( hy & FLOAT32_ABS_MASK ) ; // asm type annotation
293277
294- // Determine the sign of the result...
295- if ( sx && isOddf ( y ) ) {
296- sx = - 1.0 ;
278+ if ( sn && isOddf ( y ) ) {
279+ sn = - 1.0 ;
297280 } else {
298- sx = 1.0 ;
281+ sn = 1.0 ;
299282 }
300- // Case 1: `|y|` is huge...
301-
302- // |y| > 2^31
303- if ( ahy > HIGH_BIASED_EXP_31 ) {
304- // `|y| > 2^64`, then must over- or underflow...
305- if ( ahy > HIGH_BIASED_EXP_64 ) {
306- return yIsHuge ( x , y ) ;
307- }
308- // Over- or underflow if `x` is not close to unity...
283+ // |y| is huge
309284
285+ // if |y| > 2^27
286+ if ( ahy > HIGH_BIASED_EXP_27 ) {
287+ // Over- or underflow if `x` is not close to unity...
310288 if ( ahx < HIGH_MAX_NEAR_UNITY ) {
311289 // y < 0
312- if ( sy === 1 ) {
290+ if ( sn === 1 ) {
313291 // Signal overflow...
314- return sx * HUGE * HUGE ;
292+ return sn * HUGE * HUGE ;
315293 }
316294 // Signal underflow...
317- return sx * TINY * TINY ;
295+ return sn * TINY * TINY ;
318296 }
319- if ( ahx > HIGH_BIASED_EXP_0 ) {
297+ if ( ahx > BIASED_EXP_0 ) {
320298 // y > 0
321- if ( sy === 0 ) {
299+ if ( sn === 0 ) {
322300 // Signal overflow...
323- return sx * HUGE * HUGE ;
301+ return sn * HUGE * HUGE ;
324302 }
325303 // Signal underflow...
326- return sx * TINY * TINY ;
304+ return sn * TINY * TINY ;
327305 }
328306 // At this point, `|1-x|` is tiny (`<= 2^-20`). Suffice to compute `log(x)` by `x - x^2/2 + x^3/3 - x^4/4`.
329- t = logxf ( LOG_WORKSPACE , ax ) ;
307+ t = logx ( LOG_WORKSPACE , ax ) ;
330308 }
331309 // Case 2: `|y|` is not huge...
332310 else {
333- t = log2axf ( LOG_WORKSPACE , ax , ahx ) ;
311+ t = log2ax ( LOG_WORKSPACE , ax , ahx ) ;
334312 }
335313 // Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`...
336- y1 = setLowWord ( y , 0 ) ;
337- lp = ( ( y - y1 ) * t [ 0 ] ) + ( y * t [ 1 ] ) ;
338- hp = y1 * t [ 0 ] ;
339- z = lp + hp ;
340-
341- // Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner.
342- toWords . assign ( z , WORDS , 1 , 0 ) ;
343- j = uint32ToInt32 ( WORDS [ 0 ] ) ;
344- i = uint32ToInt32 ( WORDS [ 1 ] ) ;
345-
346- // z >= 1024
347- if ( j >= HIGH_BIASED_EXP_10 ) {
348- // z > 1024
349- if ( ( ( j - HIGH_BIASED_EXP_10 ) | i ) !== 0 ) {
350- // Signal overflow...
351- return sx * HUGE * HUGE ;
352- }
314+ ahs = fromWordf ( y ) ;
315+ ahs = float64ToFloat32 ( ahs & HIGH_BIASED_EXP_NEG_64 ) ;
316+ y1 = toWordf ( ahs ) ;
317+ lp = float64ToFloat32 ( ( float64ToFloat32 ( y - y1 ) * t1 ) + float64ToFloat32 ( y * t2 ) ) ; // eslint-disable-line max-len
318+ hp = float64ToFloat32 ( y1 * t1 ) ;
319+ z = float64ToFloat32 ( lp + hp ) ;
320+ j = fromWordf ( z ) ;
321+
322+ // z >= 128
323+ if ( j > SMALL_BIASED_EXP_10 ) {
324+ // Signal overflow
325+ return sn * HUGE * HUGE ;
326+ }
327+ // z == 128
328+ else if ( j === SMALL_BIASED_EXP_10 ) {
353329 if ( ( lp + OVT ) > ( z - hp ) ) {
354330 // Signal overflow...
355- return sx * HUGE * HUGE ;
331+ return sn * HUGE * HUGE ;
356332 }
357333 }
358334 // z <= -150
@@ -367,7 +343,7 @@ function powf( x, y ) {
367343 }
368344 }
369345 // Compute `2^(hp+lp)`...
370- z = pow2f ( j , hp , lp ) ;
346+ z = pow2 ( j , hp , lp ) ;
371347
372348 return sx * z ;
373349}
0 commit comments