2 * Copyright (C) 1999-2000 Harri Porten (porten@kde.org)
3 * Copyright (C) 2007, 2008, 2013 Apple Inc. All Rights Reserved.
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2 of the License, or (at your option) any later version.
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with this library; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 #include "MathObject.h"
25 #include "ObjectPrototype.h"
26 #include "JSCInlines.h"
28 #include <wtf/Assertions.h>
29 #include <wtf/MathExtras.h>
30 #include <wtf/RandomNumber.h>
31 #include <wtf/RandomNumberSeed.h>
32 #include <wtf/Vector.h>
36 STATIC_ASSERT_IS_TRIVIALLY_DESTRUCTIBLE(MathObject
);
38 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncAbs(ExecState
*);
39 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncACos(ExecState
*);
40 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncACosh(ExecState
*);
41 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncASin(ExecState
*);
42 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncASinh(ExecState
*);
43 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan(ExecState
*);
44 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncATanh(ExecState
*);
45 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan2(ExecState
*);
46 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncCbrt(ExecState
*);
47 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncCeil(ExecState
*);
48 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncCos(ExecState
*);
49 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncCosh(ExecState
*);
50 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncExp(ExecState
*);
51 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncExpm1(ExecState
*);
52 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncFloor(ExecState
*);
53 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncFround(ExecState
*);
54 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncHypot(ExecState
*);
55 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog(ExecState
*);
56 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog1p(ExecState
*);
57 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog10(ExecState
*);
58 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog2(ExecState
*);
59 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncMax(ExecState
*);
60 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncMin(ExecState
*);
61 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncPow(ExecState
*);
62 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncRandom(ExecState
*);
63 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncRound(ExecState
*);
64 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncSin(ExecState
*);
65 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncSinh(ExecState
*);
66 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncSqrt(ExecState
*);
67 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncTan(ExecState
*);
68 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncTanh(ExecState
*);
69 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncTrunc(ExecState
*);
70 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncIMul(ExecState
*);
76 const ClassInfo
MathObject::s_info
= { "Math", &Base::s_info
, 0, 0, CREATE_METHOD_TABLE(MathObject
) };
78 MathObject::MathObject(VM
& vm
, Structure
* structure
)
79 : JSNonFinalObject(vm
, structure
)
83 void MathObject::finishCreation(VM
& vm
, JSGlobalObject
* globalObject
)
85 Base::finishCreation(vm
);
86 ASSERT(inherits(info()));
88 putDirectWithoutTransition(vm
, Identifier(&vm
, "E"), jsNumber(exp(1.0)), DontDelete
| DontEnum
| ReadOnly
);
89 putDirectWithoutTransition(vm
, Identifier(&vm
, "LN2"), jsNumber(log(2.0)), DontDelete
| DontEnum
| ReadOnly
);
90 putDirectWithoutTransition(vm
, Identifier(&vm
, "LN10"), jsNumber(log(10.0)), DontDelete
| DontEnum
| ReadOnly
);
91 putDirectWithoutTransition(vm
, Identifier(&vm
, "LOG2E"), jsNumber(1.0 / log(2.0)), DontDelete
| DontEnum
| ReadOnly
);
92 putDirectWithoutTransition(vm
, Identifier(&vm
, "LOG10E"), jsNumber(0.4342944819032518), DontDelete
| DontEnum
| ReadOnly
);
93 putDirectWithoutTransition(vm
, Identifier(&vm
, "PI"), jsNumber(piDouble
), DontDelete
| DontEnum
| ReadOnly
);
94 putDirectWithoutTransition(vm
, Identifier(&vm
, "SQRT1_2"), jsNumber(sqrt(0.5)), DontDelete
| DontEnum
| ReadOnly
);
95 putDirectWithoutTransition(vm
, Identifier(&vm
, "SQRT2"), jsNumber(sqrt(2.0)), DontDelete
| DontEnum
| ReadOnly
);
97 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "abs"), 1, mathProtoFuncAbs
, AbsIntrinsic
, DontEnum
| Function
);
98 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "acos"), 1, mathProtoFuncACos
, NoIntrinsic
, DontEnum
| Function
);
99 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "asin"), 1, mathProtoFuncASin
, NoIntrinsic
, DontEnum
| Function
);
100 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "atan"), 1, mathProtoFuncATan
, NoIntrinsic
, DontEnum
| Function
);
101 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "acosh"), 1, mathProtoFuncACosh
, NoIntrinsic
, DontEnum
| Function
);
102 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "asinh"), 1, mathProtoFuncASinh
, NoIntrinsic
, DontEnum
| Function
);
103 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "atanh"), 1, mathProtoFuncATanh
, NoIntrinsic
, DontEnum
| Function
);
104 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "atan2"), 2, mathProtoFuncATan2
, NoIntrinsic
, DontEnum
| Function
);
105 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "cbrt"), 1, mathProtoFuncCbrt
, NoIntrinsic
, DontEnum
| Function
);
106 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "ceil"), 1, mathProtoFuncCeil
, CeilIntrinsic
, DontEnum
| Function
);
107 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "cos"), 1, mathProtoFuncCos
, CosIntrinsic
, DontEnum
| Function
);
108 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "cosh"), 1, mathProtoFuncCosh
, NoIntrinsic
, DontEnum
| Function
);
109 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "exp"), 1, mathProtoFuncExp
, ExpIntrinsic
, DontEnum
| Function
);
110 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "expm1"), 1, mathProtoFuncExpm1
, NoIntrinsic
, DontEnum
| Function
);
111 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "floor"), 1, mathProtoFuncFloor
, FloorIntrinsic
, DontEnum
| Function
);
112 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "fround"), 1, mathProtoFuncFround
, FRoundIntrinsic
, DontEnum
| Function
);
113 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "hypot"), 2, mathProtoFuncHypot
, NoIntrinsic
, DontEnum
| Function
);
114 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "log"), 1, mathProtoFuncLog
, LogIntrinsic
, DontEnum
| Function
);
115 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "log10"), 1, mathProtoFuncLog10
, NoIntrinsic
, DontEnum
| Function
);
116 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "log1p"), 1, mathProtoFuncLog1p
, NoIntrinsic
, DontEnum
| Function
);
117 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "log2"), 1, mathProtoFuncLog2
, NoIntrinsic
, DontEnum
| Function
);
118 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "max"), 2, mathProtoFuncMax
, MaxIntrinsic
, DontEnum
| Function
);
119 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "min"), 2, mathProtoFuncMin
, MinIntrinsic
, DontEnum
| Function
);
120 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "pow"), 2, mathProtoFuncPow
, PowIntrinsic
, DontEnum
| Function
);
121 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "random"), 0, mathProtoFuncRandom
, NoIntrinsic
, DontEnum
| Function
);
122 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "round"), 1, mathProtoFuncRound
, RoundIntrinsic
, DontEnum
| Function
);
123 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "sin"), 1, mathProtoFuncSin
, SinIntrinsic
, DontEnum
| Function
);
124 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "sinh"), 1, mathProtoFuncSinh
, NoIntrinsic
, DontEnum
| Function
);
125 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "sqrt"), 1, mathProtoFuncSqrt
, SqrtIntrinsic
, DontEnum
| Function
);
126 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "tan"), 1, mathProtoFuncTan
, NoIntrinsic
, DontEnum
| Function
);
127 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "tanh"), 1, mathProtoFuncTanh
, NoIntrinsic
, DontEnum
| Function
);
128 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "trunc"), 1, mathProtoFuncTrunc
, NoIntrinsic
, DontEnum
| Function
);
129 putDirectNativeFunctionWithoutTransition(vm
, globalObject
, Identifier(&vm
, "imul"), 1, mathProtoFuncIMul
, IMulIntrinsic
, DontEnum
| Function
);
132 // ------------------------------ Functions --------------------------------
134 EncodedJSValue JSC_HOST_CALL
mathProtoFuncAbs(ExecState
* exec
)
136 return JSValue::encode(jsNumber(fabs(exec
->argument(0).toNumber(exec
))));
139 EncodedJSValue JSC_HOST_CALL
mathProtoFuncACos(ExecState
* exec
)
141 return JSValue::encode(jsDoubleNumber(acos(exec
->argument(0).toNumber(exec
))));
144 EncodedJSValue JSC_HOST_CALL
mathProtoFuncASin(ExecState
* exec
)
146 return JSValue::encode(jsDoubleNumber(asin(exec
->argument(0).toNumber(exec
))));
149 EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan(ExecState
* exec
)
151 return JSValue::encode(jsDoubleNumber(atan(exec
->argument(0).toNumber(exec
))));
154 EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan2(ExecState
* exec
)
156 double arg0
= exec
->argument(0).toNumber(exec
);
157 double arg1
= exec
->argument(1).toNumber(exec
);
158 return JSValue::encode(jsDoubleNumber(atan2(arg0
, arg1
)));
161 EncodedJSValue JSC_HOST_CALL
mathProtoFuncCeil(ExecState
* exec
)
163 return JSValue::encode(jsNumber(ceil(exec
->argument(0).toNumber(exec
))));
166 EncodedJSValue JSC_HOST_CALL
mathProtoFuncCos(ExecState
* exec
)
168 return JSValue::encode(jsDoubleNumber(cos(exec
->argument(0).toNumber(exec
))));
171 EncodedJSValue JSC_HOST_CALL
mathProtoFuncExp(ExecState
* exec
)
173 return JSValue::encode(jsDoubleNumber(exp(exec
->argument(0).toNumber(exec
))));
176 EncodedJSValue JSC_HOST_CALL
mathProtoFuncFloor(ExecState
* exec
)
178 return JSValue::encode(jsNumber(floor(exec
->argument(0).toNumber(exec
))));
181 EncodedJSValue JSC_HOST_CALL
mathProtoFuncHypot(ExecState
* exec
)
183 unsigned argsCount
= exec
->argumentCount();
185 Vector
<double, 8> args
;
186 args
.reserveInitialCapacity(argsCount
);
187 for (unsigned i
= 0; i
< argsCount
; ++i
) {
188 args
.uncheckedAppend(exec
->uncheckedArgument(i
).toNumber(exec
));
189 if (exec
->hadException())
190 return JSValue::encode(jsNull());
191 if (std::isinf(args
[i
]))
192 return JSValue::encode(jsDoubleNumber(+std::numeric_limits
<double>::infinity()));
193 max
= std::max(fabs(args
[i
]), max
);
197 // Kahan summation algorithm significantly reduces the numerical error in the total obtained.
199 double compensation
= 0;
200 for (double argument
: args
) {
201 double scaledArgument
= argument
/ max
;
202 double summand
= scaledArgument
* scaledArgument
- compensation
;
203 double preliminary
= sum
+ summand
;
204 compensation
= (preliminary
- sum
) - summand
;
207 return JSValue::encode(jsDoubleNumber(sqrt(sum
) * max
));
210 EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog(ExecState
* exec
)
212 return JSValue::encode(jsDoubleNumber(log(exec
->argument(0).toNumber(exec
))));
215 EncodedJSValue JSC_HOST_CALL
mathProtoFuncMax(ExecState
* exec
)
217 unsigned argsCount
= exec
->argumentCount();
218 double result
= -std::numeric_limits
<double>::infinity();
219 for (unsigned k
= 0; k
< argsCount
; ++k
) {
220 double val
= exec
->uncheckedArgument(k
).toNumber(exec
);
221 if (std::isnan(val
)) {
223 } else if (val
> result
|| (!val
&& !result
&& !std::signbit(val
)))
226 return JSValue::encode(jsNumber(result
));
229 EncodedJSValue JSC_HOST_CALL
mathProtoFuncMin(ExecState
* exec
)
231 unsigned argsCount
= exec
->argumentCount();
232 double result
= +std::numeric_limits
<double>::infinity();
233 for (unsigned k
= 0; k
< argsCount
; ++k
) {
234 double val
= exec
->uncheckedArgument(k
).toNumber(exec
);
235 if (std::isnan(val
)) {
237 } else if (val
< result
|| (!val
&& !result
&& std::signbit(val
)))
240 return JSValue::encode(jsNumber(result
));
243 #if PLATFORM(IOS) && CPU(ARM_THUMB2)
245 static double fdlibmPow(double x
, double y
);
247 static ALWAYS_INLINE
bool isDenormal(double x
)
249 static const uint64_t signbit
= 0x8000000000000000ULL
;
250 static const uint64_t minNormal
= 0x0001000000000000ULL
;
251 return (bitwise_cast
<uint64_t>(x
) & ~signbit
) - 1 < minNormal
- 1;
254 static ALWAYS_INLINE
bool isEdgeCase(double x
)
256 static const uint64_t signbit
= 0x8000000000000000ULL
;
257 static const uint64_t infinity
= 0x7fffffffffffffffULL
;
258 return (bitwise_cast
<uint64_t>(x
) & ~signbit
) - 1 >= infinity
- 1;
261 static ALWAYS_INLINE
double mathPow(double x
, double y
)
263 if (!isDenormal(x
) && !isDenormal(y
)) {
264 double libmResult
= pow(x
,y
);
265 if (libmResult
|| isEdgeCase(x
) || isEdgeCase(y
))
268 return fdlibmPow(x
,y
);
273 ALWAYS_INLINE
double mathPow(double x
, double y
)
280 EncodedJSValue JSC_HOST_CALL
mathProtoFuncPow(ExecState
* exec
)
284 double arg
= exec
->argument(0).toNumber(exec
);
285 double arg2
= exec
->argument(1).toNumber(exec
);
287 if (std::isnan(arg2
))
288 return JSValue::encode(jsNaN());
289 if (std::isinf(arg2
) && fabs(arg
) == 1)
290 return JSValue::encode(jsNaN());
291 return JSValue::encode(jsNumber(mathPow(arg
, arg2
)));
294 EncodedJSValue JSC_HOST_CALL
mathProtoFuncRandom(ExecState
* exec
)
296 return JSValue::encode(jsDoubleNumber(exec
->lexicalGlobalObject()->weakRandomNumber()));
299 EncodedJSValue JSC_HOST_CALL
mathProtoFuncRound(ExecState
* exec
)
301 double arg
= exec
->argument(0).toNumber(exec
);
302 double integer
= ceil(arg
);
303 return JSValue::encode(jsNumber(integer
- (integer
- arg
> 0.5)));
306 EncodedJSValue JSC_HOST_CALL
mathProtoFuncSin(ExecState
* exec
)
308 return JSValue::encode(jsDoubleNumber(sin(exec
->argument(0).toNumber(exec
))));
311 EncodedJSValue JSC_HOST_CALL
mathProtoFuncSqrt(ExecState
* exec
)
313 return JSValue::encode(jsDoubleNumber(sqrt(exec
->argument(0).toNumber(exec
))));
316 EncodedJSValue JSC_HOST_CALL
mathProtoFuncTan(ExecState
* exec
)
318 return JSValue::encode(jsDoubleNumber(tan(exec
->argument(0).toNumber(exec
))));
321 EncodedJSValue JSC_HOST_CALL
mathProtoFuncIMul(ExecState
* exec
)
323 int32_t left
= exec
->argument(0).toInt32(exec
);
324 if (exec
->hadException())
325 return JSValue::encode(jsNull());
326 int32_t right
= exec
->argument(1).toInt32(exec
);
327 return JSValue::encode(jsNumber(left
* right
));
330 EncodedJSValue JSC_HOST_CALL
mathProtoFuncACosh(ExecState
* exec
)
332 return JSValue::encode(jsDoubleNumber(acosh(exec
->argument(0).toNumber(exec
))));
335 EncodedJSValue JSC_HOST_CALL
mathProtoFuncASinh(ExecState
* exec
)
337 return JSValue::encode(jsDoubleNumber(asinh(exec
->argument(0).toNumber(exec
))));
340 EncodedJSValue JSC_HOST_CALL
mathProtoFuncATanh(ExecState
* exec
)
342 return JSValue::encode(jsDoubleNumber(atanh(exec
->argument(0).toNumber(exec
))));
345 EncodedJSValue JSC_HOST_CALL
mathProtoFuncCbrt(ExecState
* exec
)
347 return JSValue::encode(jsDoubleNumber(cbrt(exec
->argument(0).toNumber(exec
))));
350 EncodedJSValue JSC_HOST_CALL
mathProtoFuncCosh(ExecState
* exec
)
352 return JSValue::encode(jsDoubleNumber(cosh(exec
->argument(0).toNumber(exec
))));
355 EncodedJSValue JSC_HOST_CALL
mathProtoFuncExpm1(ExecState
* exec
)
357 return JSValue::encode(jsDoubleNumber(expm1(exec
->argument(0).toNumber(exec
))));
360 EncodedJSValue JSC_HOST_CALL
mathProtoFuncFround(ExecState
* exec
)
362 return JSValue::encode(jsDoubleNumber(static_cast<float>(exec
->argument(0).toNumber(exec
))));
365 EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog1p(ExecState
* exec
)
367 double value
= exec
->argument(0).toNumber(exec
);
369 return JSValue::encode(jsDoubleNumber(value
));
370 return JSValue::encode(jsDoubleNumber(log1p(value
)));
373 EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog10(ExecState
* exec
)
375 return JSValue::encode(jsDoubleNumber(log10(exec
->argument(0).toNumber(exec
))));
378 EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog2(ExecState
* exec
)
380 return JSValue::encode(jsDoubleNumber(log2(exec
->argument(0).toNumber(exec
))));
383 EncodedJSValue JSC_HOST_CALL
mathProtoFuncSinh(ExecState
* exec
)
385 return JSValue::encode(jsDoubleNumber(sinh(exec
->argument(0).toNumber(exec
))));
388 EncodedJSValue JSC_HOST_CALL
mathProtoFuncTanh(ExecState
* exec
)
390 return JSValue::encode(jsDoubleNumber(tanh(exec
->argument(0).toNumber(exec
))));
393 EncodedJSValue JSC_HOST_CALL
mathProtoFuncTrunc(ExecState
*exec
)
395 return JSValue::encode(jsNumber(exec
->argument(0).toIntegerPreserveNaN(exec
)));
399 #if PLATFORM(IOS) && CPU(ARM_THUMB2)
401 // The following code is taken from netlib.org:
402 // http://www.netlib.org/fdlibm/fdlibm.h
403 // http://www.netlib.org/fdlibm/e_pow.c
404 // http://www.netlib.org/fdlibm/s_scalbn.c
406 // And was originally distributed under the following license:
409 * ====================================================
410 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
412 * Developed at SunSoft, a Sun Microsystems, Inc. business.
413 * Permission to use, copy, modify, and distribute this
414 * software is freely granted, provided that this notice
416 * ====================================================
419 * ====================================================
420 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
422 * Permission to use, copy, modify, and distribute this
423 * software is freely granted, provided that this notice
425 * ====================================================
428 /* __ieee754_pow(x,y) return x**y
431 * Method: Let x = 2 * (1+f)
432 * 1. Compute and return log2(x) in two pieces:
434 * where w1 has 53-24 = 29 bit trailing zeros.
435 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
436 * arithmetic, where |y'|<=0.5.
437 * 3. Return x**y = 2**n*exp(y'*log2)
440 * 1. (anything) ** 0 is 1
441 * 2. (anything) ** 1 is itself
442 * 3. (anything) ** NAN is NAN
443 * 4. NAN ** (anything except 0) is NAN
444 * 5. +-(|x| > 1) ** +INF is +INF
445 * 6. +-(|x| > 1) ** -INF is +0
446 * 7. +-(|x| < 1) ** +INF is +0
447 * 8. +-(|x| < 1) ** -INF is +INF
448 * 9. +-1 ** +-INF is NAN
449 * 10. +0 ** (+anything except 0, NAN) is +0
450 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
451 * 12. +0 ** (-anything except 0, NAN) is +INF
452 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
453 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
454 * 15. +INF ** (+anything except 0,NAN) is +INF
455 * 16. +INF ** (-anything except 0,NAN) is +0
456 * 17. -INF ** (anything) = -0 ** (-anything)
457 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
458 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
461 * pow(x,y) returns x**y nearly rounded. In particular
462 * pow(integer,integer)
463 * always returns the correct integer provided it is
467 * The hexadecimal values are the intended ones for the following
468 * constants. The decimal values may be used, provided that the
469 * compiler will convert from decimal to binary accurately enough
470 * to produce the hexadecimal values shown.
473 #define __HI(x) *(1+(int*)&x)
474 #define __LO(x) *(int*)&x
478 dp_h
[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
479 dp_l
[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
483 two53
= 9007199254740992.0, /* 0x43400000, 0x00000000 */
487 two54
= 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
488 twom54
= 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
489 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
490 L1
= 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
491 L2
= 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
492 L3
= 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
493 L4
= 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
494 L5
= 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
495 L6
= 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
496 P1
= 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
497 P2
= -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
498 P3
= 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
499 P4
= -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
500 P5
= 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
501 lg2
= 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
502 lg2_h
= 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
503 lg2_l
= -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
504 ovt
= 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
505 cp
= 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
506 cp_h
= 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
507 cp_l
= -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
508 ivln2
= 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
509 ivln2_h
= 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
510 ivln2_l
= 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
512 inline double fdlibmScalbn (double x
, int n
)
517 k
= (hx
&0x7ff00000)>>20; /* extract exponent */
518 if (k
==0) { /* 0 or subnormal x */
519 if ((lx
|(hx
&0x7fffffff))==0) return x
; /* +-0 */
522 k
= ((hx
&0x7ff00000)>>20) - 54;
523 if (n
< -50000) return tiny
*x
; /*underflow*/
525 if (k
==0x7ff) return x
+x
; /* NaN or Inf */
527 if (k
> 0x7fe) return huge
*copysign(huge
,x
); /* overflow */
528 if (k
> 0) /* normal result */
529 {__HI(x
) = (hx
&0x800fffff)|(k
<<20); return x
;}
531 if (n
> 50000) /* in case integer overflow in n+k */
532 return huge
*copysign(huge
,x
); /*overflow*/
533 else return tiny
*copysign(tiny
,x
); /*underflow*/
535 k
+= 54; /* subnormal result */
536 __HI(x
) = (hx
&0x800fffff)|(k
<<20);
540 double fdlibmPow(double x
, double y
)
542 double z
,ax
,z_h
,z_l
,p_h
,p_l
;
543 double y1
,t1
,t2
,r
,s
,t
,u
,v
,w
;
544 int i0
,i1
,i
,j
,k
,yisint
,n
;
548 i0
= ((*(int*)&one
)>>29)^1; i1
=1-i0
;
549 hx
= __HI(x
); lx
= __LO(x
);
550 hy
= __HI(y
); ly
= __LO(y
);
551 ix
= hx
&0x7fffffff; iy
= hy
&0x7fffffff;
553 /* y==zero: x**0 = 1 */
554 if((iy
|ly
)==0) return one
;
556 /* +-NaN return x+y */
557 if(ix
> 0x7ff00000 || ((ix
==0x7ff00000)&&(lx
!=0)) ||
558 iy
> 0x7ff00000 || ((iy
==0x7ff00000)&&(ly
!=0)))
561 /* determine if y is an odd int when x < 0
562 * yisint = 0 ... y is not an integer
563 * yisint = 1 ... y is an odd int
564 * yisint = 2 ... y is an even int
568 if(iy
>=0x43400000) yisint
= 2; /* even integer y */
569 else if(iy
>=0x3ff00000) {
570 k
= (iy
>>20)-0x3ff; /* exponent */
573 if(static_cast<unsigned>(j
<<(52-k
))==ly
) yisint
= 2-(j
&1);
576 if((j
<<(20-k
))==iy
) yisint
= 2-(j
&1);
581 /* special value of y */
583 if (iy
==0x7ff00000) { /* y is +-inf */
584 if(((ix
-0x3ff00000)|lx
)==0)
585 return y
- y
; /* inf**+-1 is NaN */
586 else if (ix
>= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
587 return (hy
>=0)? y
: zero
;
588 else /* (|x|<1)**-,+inf = inf,0 */
589 return (hy
<0)?-y
: zero
;
591 if(iy
==0x3ff00000) { /* y is +-1 */
592 if(hy
<0) return one
/x
; else return x
;
594 if(hy
==0x40000000) return x
*x
; /* y is 2 */
595 if(hy
==0x3fe00000) { /* y is 0.5 */
596 if(hx
>=0) /* x >= +0 */
602 /* special value of x */
604 if(ix
==0x7ff00000||ix
==0||ix
==0x3ff00000){
605 z
= ax
; /*x is +-0,+-inf,+-1*/
606 if(hy
<0) z
= one
/z
; /* z = (1/|x|) */
608 if(((ix
-0x3ff00000)|yisint
)==0) {
609 z
= (z
-z
)/(z
-z
); /* (-1)**non-int is NaN */
611 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
619 /* (x<0)**(non-int) is NaN */
620 if((n
|yisint
)==0) return (x
-x
)/(x
-x
);
622 s
= one
; /* s (sign of result -ve**odd) = -1 else = 1 */
623 if((n
|(yisint
-1))==0) s
= -one
;/* (-ve)**(odd int) */
626 if(iy
>0x41e00000) { /* if |y| > 2**31 */
627 if(iy
>0x43f00000){ /* if |y| > 2**64, must o/uflow */
628 if(ix
<=0x3fefffff) return (hy
<0)? huge
*huge
:tiny
*tiny
;
629 if(ix
>=0x3ff00000) return (hy
>0)? huge
*huge
:tiny
*tiny
;
631 /* over/underflow if x is not close to one */
632 if(ix
<0x3fefffff) return (hy
<0)? s
*huge
*huge
:s
*tiny
*tiny
;
633 if(ix
>0x3ff00000) return (hy
>0)? s
*huge
*huge
:s
*tiny
*tiny
;
634 /* now |1-x| is tiny <= 2**-20, suffice to compute
635 log(x) by x-x^2/2+x^3/3-x^4/4 */
636 t
= ax
-one
; /* t has 20 trailing zeros */
637 w
= (t
*t
)*(0.5-t
*(0.3333333333333333333333-t
*0.25));
638 u
= ivln2_h
*t
; /* ivln2_h has 21 sig. bits */
639 v
= t
*ivln2_l
-w
*ivln2
;
644 double ss
,s2
,s_h
,s_l
,t_h
,t_l
;
646 /* take care subnormal number */
648 {ax
*= two53
; n
-= 53; ix
= __HI(ax
); }
649 n
+= ((ix
)>>20)-0x3ff;
651 /* determine interval */
652 ix
= j
|0x3ff00000; /* normalize ix */
653 if(j
<=0x3988E) k
=0; /* |x|<sqrt(3/2) */
654 else if(j
<0xBB67A) k
=1; /* |x|<sqrt(3) */
655 else {k
=0;n
+=1;ix
-= 0x00100000;}
658 /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
659 u
= ax
-bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
664 /* t_h=ax+bp[k] High */
666 __HI(t_h
)=((ix
>>1)|0x20000000)+0x00080000+(k
<<18);
667 t_l
= ax
- (t_h
-bp
[k
]);
668 s_l
= v
*((u
-s_h
*t_h
)-s_h
*t_l
);
669 /* compute log(ax) */
671 r
= s2
*s2
*(L1
+s2
*(L2
+s2
*(L3
+s2
*(L4
+s2
*(L5
+s2
*L6
)))));
676 t_l
= r
-((t_h
-3.0)-s2
);
677 /* u+v = ss*(1+...) */
680 /* 2/(3log2)*(ss+...) */
684 z_h
= cp_h
*p_h
; /* cp_h+cp_l = 2/(3*log2) */
685 z_l
= cp_l
*p_h
+p_l
*cp
+dp_l
[k
];
686 /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
688 t1
= (((z_h
+z_l
)+dp_h
[k
])+t
);
690 t2
= z_l
-(((t1
-t
)-dp_h
[k
])-z_h
);
693 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
696 p_l
= (y
-y1
)*t1
+y
*t2
;
701 if (j
>=0x40900000) { /* z >= 1024 */
702 if(((j
-0x40900000)|i
)!=0) /* if z > 1024 */
703 return s
*huge
*huge
; /* overflow */
705 if(p_l
+ovt
>z
-p_h
) return s
*huge
*huge
; /* overflow */
707 } else if((j
&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
708 if(((j
-0xc090cc00)|i
)!=0) /* z < -1075 */
709 return s
*tiny
*tiny
; /* underflow */
711 if(p_l
<=z
-p_h
) return s
*tiny
*tiny
; /* underflow */
715 * compute 2**(p_h+p_l)
720 if(i
>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
721 n
= j
+(0x00100000>>(k
+1));
722 k
= ((n
&0x7fffffff)>>20)-0x3ff; /* new k for n */
724 __HI(t
) = (n
&~(0x000fffff>>k
));
725 n
= ((n
&0x000fffff)|0x00100000)>>(20-k
);
732 v
= (p_l
-(t
-p_h
))*lg2
+t
*lg2_l
;
736 t1
= z
- t
*(P1
+t
*(P2
+t
*(P3
+t
*(P4
+t
*P5
))));
737 r
= (z
*t1
)/(t1
-two
)-(w
+z
*w
);
741 if((j
>>20)<=0) z
= fdlibmScalbn(z
,n
); /* subnormal output */
742 else __HI(z
) += (n
<<20);