18
18
*
19
19
* ## Notice
20
20
*
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.
22
22
*
23
23
* ```text
24
24
* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
@@ -39,57 +39,50 @@ var isInfinitef = require( '@stdlib/math/base/assert/is-infinitef' );
39
39
var isIntegerf = require ( '@stdlib/math/base/assert/is-integerf' ) ;
40
40
var sqrtf = require ( '@stdlib/math/base/special/sqrtf' ) ;
41
41
var 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 ' ) ;
45
45
var NINF = require ( '@stdlib/constants/float32/ninf' ) ;
46
46
var PINF = require ( '@stdlib/constants/float32/pinf' ) ;
47
- var ABS_MASK = require ( '@stdlib/constants/float32/high-word-abs-mask' ) ;
48
47
var FLOAT32_ABS_MASK = require ( '@stdlib/constants/float32/abs-mask' ) ;
49
48
var xIsZero = require ( './x_is_zero.js' ) ;
50
- var yIsHuge = require ( './y_is_huge.js' ) ;
51
49
var 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' ) ;
55
53
56
54
57
55
// VARIABLES //
58
56
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
61
59
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
64
62
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
67
65
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
70
68
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
73
71
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
76
74
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
81
77
82
78
var HUGE = 1.0e30 ;
83
79
var TINY = 1.0e-30 ;
84
80
85
81
// -(1024-log2(ovfl+.5ulp))
86
82
var OVT = 4.2995665694e-08 ;
87
83
88
- // High/low words workspace:
89
- var WORDS = [ 0 | 0 , 0 | 0 ] ;
90
-
91
84
// Log workspace:
92
- var LOG_WORKSPACE = [ 0.0 , 0.0 ] ;
85
+ var LOG_WORKSPACE = 0.0 ;
93
86
94
87
95
88
// MAIN //
@@ -194,64 +187,55 @@ var LOG_WORKSPACE = [ 0.0, 0.0 ];
194
187
* // returns NaN
195
188
*/
196
189
function 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 ;
206
199
var y1 ;
200
+ var t1 ;
201
+ var t2 ;
207
202
var hp ;
208
203
var lp ;
209
204
var z ; // y prime
210
205
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 ) ;
219
209
220
210
// 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 ) ;
250
238
}
251
- // Split `x` into high and low words:
252
- toWords . assign ( x , WORDS , 1 , 0 ) ;
253
- hx = WORDS [ 0 ] ;
254
- lx = WORDS [ 1 ] ;
255
239
256
240
// Special cases `x`...
257
241
if ( lx === 0 ) {
@@ -291,68 +275,60 @@ function powf( x, y ) {
291
275
ahx = ( hx & FLOAT32_ABS_MASK ) ; // asm type annotation
292
276
ahy = ( hy & FLOAT32_ABS_MASK ) ; // asm type annotation
293
277
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 ;
297
280
} else {
298
- sx = 1.0 ;
281
+ sn = 1.0 ;
299
282
}
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
309
284
285
+ // if |y| > 2^27
286
+ if ( ahy > HIGH_BIASED_EXP_27 ) {
287
+ // Over- or underflow if `x` is not close to unity...
310
288
if ( ahx < HIGH_MAX_NEAR_UNITY ) {
311
289
// y < 0
312
- if ( sy === 1 ) {
290
+ if ( sn === 1 ) {
313
291
// Signal overflow...
314
- return sx * HUGE * HUGE ;
292
+ return sn * HUGE * HUGE ;
315
293
}
316
294
// Signal underflow...
317
- return sx * TINY * TINY ;
295
+ return sn * TINY * TINY ;
318
296
}
319
- if ( ahx > HIGH_BIASED_EXP_0 ) {
297
+ if ( ahx > BIASED_EXP_0 ) {
320
298
// y > 0
321
- if ( sy === 0 ) {
299
+ if ( sn === 0 ) {
322
300
// Signal overflow...
323
- return sx * HUGE * HUGE ;
301
+ return sn * HUGE * HUGE ;
324
302
}
325
303
// Signal underflow...
326
- return sx * TINY * TINY ;
304
+ return sn * TINY * TINY ;
327
305
}
328
306
// 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 ) ;
330
308
}
331
309
// Case 2: `|y|` is not huge...
332
310
else {
333
- t = log2axf ( LOG_WORKSPACE , ax , ahx ) ;
311
+ t = log2ax ( LOG_WORKSPACE , ax , ahx ) ;
334
312
}
335
313
// 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 ) {
353
329
if ( ( lp + OVT ) > ( z - hp ) ) {
354
330
// Signal overflow...
355
- return sx * HUGE * HUGE ;
331
+ return sn * HUGE * HUGE ;
356
332
}
357
333
}
358
334
// z <= -150
@@ -367,7 +343,7 @@ function powf( x, y ) {
367
343
}
368
344
}
369
345
// Compute `2^(hp+lp)`...
370
- z = pow2f ( j , hp , lp ) ;
346
+ z = pow2 ( j , hp , lp ) ;
371
347
372
348
return sx * z ;
373
349
}
0 commit comments