1/* 2 * Copyright (C) 1999-2000 Harri Porten (porten@kde.org) 3 * Copyright (C) 2007, 2008, 2013 Apple Inc. All Rights Reserved. 4 * 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. 9 * 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. 14 * 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 18 * 19 */ 20 21#include "config.h" 22#include "MathObject.h" 23 24#include "Lookup.h" 25#include "ObjectPrototype.h" 26#include "JSCInlines.h" 27#include <time.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> 33 34namespace JSC { 35 36STATIC_ASSERT_IS_TRIVIALLY_DESTRUCTIBLE(MathObject); 37 38static EncodedJSValue JSC_HOST_CALL mathProtoFuncAbs(ExecState*); 39static EncodedJSValue JSC_HOST_CALL mathProtoFuncACos(ExecState*); 40static EncodedJSValue JSC_HOST_CALL mathProtoFuncACosh(ExecState*); 41static EncodedJSValue JSC_HOST_CALL mathProtoFuncASin(ExecState*); 42static EncodedJSValue JSC_HOST_CALL mathProtoFuncASinh(ExecState*); 43static EncodedJSValue JSC_HOST_CALL mathProtoFuncATan(ExecState*); 44static EncodedJSValue JSC_HOST_CALL mathProtoFuncATanh(ExecState*); 45static EncodedJSValue JSC_HOST_CALL mathProtoFuncATan2(ExecState*); 46static EncodedJSValue JSC_HOST_CALL mathProtoFuncCbrt(ExecState*); 47static EncodedJSValue JSC_HOST_CALL mathProtoFuncCeil(ExecState*); 48static EncodedJSValue JSC_HOST_CALL mathProtoFuncCos(ExecState*); 49static EncodedJSValue JSC_HOST_CALL mathProtoFuncCosh(ExecState*); 50static EncodedJSValue JSC_HOST_CALL mathProtoFuncExp(ExecState*); 51static EncodedJSValue JSC_HOST_CALL mathProtoFuncExpm1(ExecState*); 52static EncodedJSValue JSC_HOST_CALL mathProtoFuncFloor(ExecState*); 53static EncodedJSValue JSC_HOST_CALL mathProtoFuncFround(ExecState*); 54static EncodedJSValue JSC_HOST_CALL mathProtoFuncHypot(ExecState*); 55static EncodedJSValue JSC_HOST_CALL mathProtoFuncLog(ExecState*); 56static EncodedJSValue JSC_HOST_CALL mathProtoFuncLog1p(ExecState*); 57static EncodedJSValue JSC_HOST_CALL mathProtoFuncLog10(ExecState*); 58static EncodedJSValue JSC_HOST_CALL mathProtoFuncLog2(ExecState*); 59static EncodedJSValue JSC_HOST_CALL mathProtoFuncMax(ExecState*); 60static EncodedJSValue JSC_HOST_CALL mathProtoFuncMin(ExecState*); 61static EncodedJSValue JSC_HOST_CALL mathProtoFuncPow(ExecState*); 62static EncodedJSValue JSC_HOST_CALL mathProtoFuncRandom(ExecState*); 63static EncodedJSValue JSC_HOST_CALL mathProtoFuncRound(ExecState*); 64static EncodedJSValue JSC_HOST_CALL mathProtoFuncSin(ExecState*); 65static EncodedJSValue JSC_HOST_CALL mathProtoFuncSinh(ExecState*); 66static EncodedJSValue JSC_HOST_CALL mathProtoFuncSqrt(ExecState*); 67static EncodedJSValue JSC_HOST_CALL mathProtoFuncTan(ExecState*); 68static EncodedJSValue JSC_HOST_CALL mathProtoFuncTanh(ExecState*); 69static EncodedJSValue JSC_HOST_CALL mathProtoFuncTrunc(ExecState*); 70static EncodedJSValue JSC_HOST_CALL mathProtoFuncIMul(ExecState*); 71 72} 73 74namespace JSC { 75 76const ClassInfo MathObject::s_info = { "Math", &Base::s_info, 0, 0, CREATE_METHOD_TABLE(MathObject) }; 77 78MathObject::MathObject(VM& vm, Structure* structure) 79 : JSNonFinalObject(vm, structure) 80{ 81} 82 83void MathObject::finishCreation(VM& vm, JSGlobalObject* globalObject) 84{ 85 Base::finishCreation(vm); 86 ASSERT(inherits(info())); 87 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); 96 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); 130} 131 132// ------------------------------ Functions -------------------------------- 133 134EncodedJSValue JSC_HOST_CALL mathProtoFuncAbs(ExecState* exec) 135{ 136 return JSValue::encode(jsNumber(fabs(exec->argument(0).toNumber(exec)))); 137} 138 139EncodedJSValue JSC_HOST_CALL mathProtoFuncACos(ExecState* exec) 140{ 141 return JSValue::encode(jsDoubleNumber(acos(exec->argument(0).toNumber(exec)))); 142} 143 144EncodedJSValue JSC_HOST_CALL mathProtoFuncASin(ExecState* exec) 145{ 146 return JSValue::encode(jsDoubleNumber(asin(exec->argument(0).toNumber(exec)))); 147} 148 149EncodedJSValue JSC_HOST_CALL mathProtoFuncATan(ExecState* exec) 150{ 151 return JSValue::encode(jsDoubleNumber(atan(exec->argument(0).toNumber(exec)))); 152} 153 154EncodedJSValue JSC_HOST_CALL mathProtoFuncATan2(ExecState* exec) 155{ 156 double arg0 = exec->argument(0).toNumber(exec); 157 double arg1 = exec->argument(1).toNumber(exec); 158 return JSValue::encode(jsDoubleNumber(atan2(arg0, arg1))); 159} 160 161EncodedJSValue JSC_HOST_CALL mathProtoFuncCeil(ExecState* exec) 162{ 163 return JSValue::encode(jsNumber(ceil(exec->argument(0).toNumber(exec)))); 164} 165 166EncodedJSValue JSC_HOST_CALL mathProtoFuncCos(ExecState* exec) 167{ 168 return JSValue::encode(jsDoubleNumber(cos(exec->argument(0).toNumber(exec)))); 169} 170 171EncodedJSValue JSC_HOST_CALL mathProtoFuncExp(ExecState* exec) 172{ 173 return JSValue::encode(jsDoubleNumber(exp(exec->argument(0).toNumber(exec)))); 174} 175 176EncodedJSValue JSC_HOST_CALL mathProtoFuncFloor(ExecState* exec) 177{ 178 return JSValue::encode(jsNumber(floor(exec->argument(0).toNumber(exec)))); 179} 180 181EncodedJSValue JSC_HOST_CALL mathProtoFuncHypot(ExecState* exec) 182{ 183 unsigned argsCount = exec->argumentCount(); 184 double max = 0; 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); 194 } 195 if (!max) 196 max = 1; 197 // Kahan summation algorithm significantly reduces the numerical error in the total obtained. 198 double sum = 0; 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; 205 sum = preliminary; 206 } 207 return JSValue::encode(jsDoubleNumber(sqrt(sum) * max)); 208} 209 210EncodedJSValue JSC_HOST_CALL mathProtoFuncLog(ExecState* exec) 211{ 212 return JSValue::encode(jsDoubleNumber(log(exec->argument(0).toNumber(exec)))); 213} 214 215EncodedJSValue JSC_HOST_CALL mathProtoFuncMax(ExecState* exec) 216{ 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)) { 222 result = PNaN; 223 } else if (val > result || (!val && !result && !std::signbit(val))) 224 result = val; 225 } 226 return JSValue::encode(jsNumber(result)); 227} 228 229EncodedJSValue JSC_HOST_CALL mathProtoFuncMin(ExecState* exec) 230{ 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)) { 236 result = PNaN; 237 } else if (val < result || (!val && !result && std::signbit(val))) 238 result = val; 239 } 240 return JSValue::encode(jsNumber(result)); 241} 242 243#if PLATFORM(IOS) && CPU(ARM_THUMB2) 244 245static double fdlibmPow(double x, double y); 246 247static ALWAYS_INLINE bool isDenormal(double x) 248{ 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; 252} 253 254static ALWAYS_INLINE bool isEdgeCase(double x) 255{ 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; 259} 260 261static ALWAYS_INLINE double mathPow(double x, double y) 262{ 263 if (!isDenormal(x) && !isDenormal(y)) { 264 double libmResult = pow(x,y); 265 if (libmResult || isEdgeCase(x) || isEdgeCase(y)) 266 return libmResult; 267 } 268 return fdlibmPow(x,y); 269} 270 271#else 272 273ALWAYS_INLINE double mathPow(double x, double y) 274{ 275 return pow(x, y); 276} 277 278#endif 279 280EncodedJSValue JSC_HOST_CALL mathProtoFuncPow(ExecState* exec) 281{ 282 // ECMA 15.8.2.1.13 283 284 double arg = exec->argument(0).toNumber(exec); 285 double arg2 = exec->argument(1).toNumber(exec); 286 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))); 292} 293 294EncodedJSValue JSC_HOST_CALL mathProtoFuncRandom(ExecState* exec) 295{ 296 return JSValue::encode(jsDoubleNumber(exec->lexicalGlobalObject()->weakRandomNumber())); 297} 298 299EncodedJSValue JSC_HOST_CALL mathProtoFuncRound(ExecState* exec) 300{ 301 double arg = exec->argument(0).toNumber(exec); 302 double integer = ceil(arg); 303 return JSValue::encode(jsNumber(integer - (integer - arg > 0.5))); 304} 305 306EncodedJSValue JSC_HOST_CALL mathProtoFuncSin(ExecState* exec) 307{ 308 return JSValue::encode(jsDoubleNumber(sin(exec->argument(0).toNumber(exec)))); 309} 310 311EncodedJSValue JSC_HOST_CALL mathProtoFuncSqrt(ExecState* exec) 312{ 313 return JSValue::encode(jsDoubleNumber(sqrt(exec->argument(0).toNumber(exec)))); 314} 315 316EncodedJSValue JSC_HOST_CALL mathProtoFuncTan(ExecState* exec) 317{ 318 return JSValue::encode(jsDoubleNumber(tan(exec->argument(0).toNumber(exec)))); 319} 320 321EncodedJSValue JSC_HOST_CALL mathProtoFuncIMul(ExecState* exec) 322{ 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)); 328} 329 330EncodedJSValue JSC_HOST_CALL mathProtoFuncACosh(ExecState* exec) 331{ 332 return JSValue::encode(jsDoubleNumber(acosh(exec->argument(0).toNumber(exec)))); 333} 334 335EncodedJSValue JSC_HOST_CALL mathProtoFuncASinh(ExecState* exec) 336{ 337 return JSValue::encode(jsDoubleNumber(asinh(exec->argument(0).toNumber(exec)))); 338} 339 340EncodedJSValue JSC_HOST_CALL mathProtoFuncATanh(ExecState* exec) 341{ 342 return JSValue::encode(jsDoubleNumber(atanh(exec->argument(0).toNumber(exec)))); 343} 344 345EncodedJSValue JSC_HOST_CALL mathProtoFuncCbrt(ExecState* exec) 346{ 347 return JSValue::encode(jsDoubleNumber(cbrt(exec->argument(0).toNumber(exec)))); 348} 349 350EncodedJSValue JSC_HOST_CALL mathProtoFuncCosh(ExecState* exec) 351{ 352 return JSValue::encode(jsDoubleNumber(cosh(exec->argument(0).toNumber(exec)))); 353} 354 355EncodedJSValue JSC_HOST_CALL mathProtoFuncExpm1(ExecState* exec) 356{ 357 return JSValue::encode(jsDoubleNumber(expm1(exec->argument(0).toNumber(exec)))); 358} 359 360EncodedJSValue JSC_HOST_CALL mathProtoFuncFround(ExecState* exec) 361{ 362 return JSValue::encode(jsDoubleNumber(static_cast<float>(exec->argument(0).toNumber(exec)))); 363} 364 365EncodedJSValue JSC_HOST_CALL mathProtoFuncLog1p(ExecState* exec) 366{ 367 double value = exec->argument(0).toNumber(exec); 368 if (value == 0) 369 return JSValue::encode(jsDoubleNumber(value)); 370 return JSValue::encode(jsDoubleNumber(log1p(value))); 371} 372 373EncodedJSValue JSC_HOST_CALL mathProtoFuncLog10(ExecState* exec) 374{ 375 return JSValue::encode(jsDoubleNumber(log10(exec->argument(0).toNumber(exec)))); 376} 377 378EncodedJSValue JSC_HOST_CALL mathProtoFuncLog2(ExecState* exec) 379{ 380 return JSValue::encode(jsDoubleNumber(log2(exec->argument(0).toNumber(exec)))); 381} 382 383EncodedJSValue JSC_HOST_CALL mathProtoFuncSinh(ExecState* exec) 384{ 385 return JSValue::encode(jsDoubleNumber(sinh(exec->argument(0).toNumber(exec)))); 386} 387 388EncodedJSValue JSC_HOST_CALL mathProtoFuncTanh(ExecState* exec) 389{ 390 return JSValue::encode(jsDoubleNumber(tanh(exec->argument(0).toNumber(exec)))); 391} 392 393EncodedJSValue JSC_HOST_CALL mathProtoFuncTrunc(ExecState*exec) 394{ 395 return JSValue::encode(jsNumber(exec->argument(0).toIntegerPreserveNaN(exec))); 396} 397 398 399#if PLATFORM(IOS) && CPU(ARM_THUMB2) 400 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 405// 406// And was originally distributed under the following license: 407 408/* 409 * ==================================================== 410 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 411 * 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 415 * is preserved. 416 * ==================================================== 417 */ 418/* 419 * ==================================================== 420 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. 421 * 422 * Permission to use, copy, modify, and distribute this 423 * software is freely granted, provided that this notice 424 * is preserved. 425 * ==================================================== 426 */ 427 428/* __ieee754_pow(x,y) return x**y 429 * 430 * n 431 * Method: Let x = 2 * (1+f) 432 * 1. Compute and return log2(x) in two pieces: 433 * log2(x) = w1 + w2, 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) 438 * 439 * Special cases: 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 459 * 460 * Accuracy: 461 * pow(x,y) returns x**y nearly rounded. In particular 462 * pow(integer,integer) 463 * always returns the correct integer provided it is 464 * representable. 465 * 466 * Constants : 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. 471 */ 472 473#define __HI(x) *(1+(int*)&x) 474#define __LO(x) *(int*)&x 475 476static const double 477bp[] = {1.0, 1.5,}, 478dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ 479dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ 480zero = 0.0, 481one = 1.0, 482two = 2.0, 483two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ 484huge = 1.0e300, 485tiny = 1.0e-300, 486 /* for scalbn */ 487two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 488twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ 489 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ 490L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ 491L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ 492L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ 493L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ 494L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ 495L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ 496P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 497P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 498P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 499P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 500P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ 501lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ 502lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ 503lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ 504ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ 505cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ 506cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ 507cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ 508ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ 509ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ 510ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ 511 512inline double fdlibmScalbn (double x, int n) 513{ 514 int k,hx,lx; 515 hx = __HI(x); 516 lx = __LO(x); 517 k = (hx&0x7ff00000)>>20; /* extract exponent */ 518 if (k==0) { /* 0 or subnormal x */ 519 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ 520 x *= two54; 521 hx = __HI(x); 522 k = ((hx&0x7ff00000)>>20) - 54; 523 if (n< -50000) return tiny*x; /*underflow*/ 524 } 525 if (k==0x7ff) return x+x; /* NaN or Inf */ 526 k = k+n; 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;} 530 if (k <= -54) { 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*/ 534 } 535 k += 54; /* subnormal result */ 536 __HI(x) = (hx&0x800fffff)|(k<<20); 537 return x*twom54; 538} 539 540double fdlibmPow(double x, double y) 541{ 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; 545 int hx,hy,ix,iy; 546 unsigned lx,ly; 547 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; 552 553 /* y==zero: x**0 = 1 */ 554 if((iy|ly)==0) return one; 555 556 /* +-NaN return x+y */ 557 if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || 558 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 559 return x+y; 560 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 565 */ 566 yisint = 0; 567 if(hx<0) { 568 if(iy>=0x43400000) yisint = 2; /* even integer y */ 569 else if(iy>=0x3ff00000) { 570 k = (iy>>20)-0x3ff; /* exponent */ 571 if(k>20) { 572 j = ly>>(52-k); 573 if(static_cast<unsigned>(j<<(52-k))==ly) yisint = 2-(j&1); 574 } else if(ly==0) { 575 j = iy>>(20-k); 576 if((j<<(20-k))==iy) yisint = 2-(j&1); 577 } 578 } 579 } 580 581 /* special value of y */ 582 if(ly==0) { 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; 590 } 591 if(iy==0x3ff00000) { /* y is +-1 */ 592 if(hy<0) return one/x; else return x; 593 } 594 if(hy==0x40000000) return x*x; /* y is 2 */ 595 if(hy==0x3fe00000) { /* y is 0.5 */ 596 if(hx>=0) /* x >= +0 */ 597 return sqrt(x); 598 } 599 } 600 601 ax = fabs(x); 602 /* special value of x */ 603 if(lx==0) { 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|) */ 607 if(hx<0) { 608 if(((ix-0x3ff00000)|yisint)==0) { 609 z = (z-z)/(z-z); /* (-1)**non-int is NaN */ 610 } else if(yisint==1) 611 z = -z; /* (x<0)**odd = -(|x|**odd) */ 612 } 613 return z; 614 } 615 } 616 617 n = (hx>>31)+1; 618 619 /* (x<0)**(non-int) is NaN */ 620 if((n|yisint)==0) return (x-x)/(x-x); 621 622 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ 623 if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */ 624 625 /* |y| is huge */ 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; 630 } 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; 640 t1 = u+v; 641 __LO(t1) = 0; 642 t2 = v-(t1-u); 643 } else { 644 double ss,s2,s_h,s_l,t_h,t_l; 645 n = 0; 646 /* take care subnormal number */ 647 if(ix<0x00100000) 648 {ax *= two53; n -= 53; ix = __HI(ax); } 649 n += ((ix)>>20)-0x3ff; 650 j = ix&0x000fffff; 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;} 656 __HI(ax) = ix; 657 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 */ 660 v = one/(ax+bp[k]); 661 ss = u*v; 662 s_h = ss; 663 __LO(s_h) = 0; 664 /* t_h=ax+bp[k] High */ 665 t_h = zero; 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) */ 670 s2 = ss*ss; 671 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); 672 r += s_l*(s_h+ss); 673 s2 = s_h*s_h; 674 t_h = 3.0+s2+r; 675 __LO(t_h) = 0; 676 t_l = r-((t_h-3.0)-s2); 677 /* u+v = ss*(1+...) */ 678 u = s_h*t_h; 679 v = s_l*t_h+t_l*ss; 680 /* 2/(3log2)*(ss+...) */ 681 p_h = u+v; 682 __LO(p_h) = 0; 683 p_l = v-(p_h-u); 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 */ 687 t = (double)n; 688 t1 = (((z_h+z_l)+dp_h[k])+t); 689 __LO(t1) = 0; 690 t2 = z_l-(((t1-t)-dp_h[k])-z_h); 691 } 692 693 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ 694 y1 = y; 695 __LO(y1) = 0; 696 p_l = (y-y1)*t1+y*t2; 697 p_h = y1*t1; 698 z = p_l+p_h; 699 j = __HI(z); 700 i = __LO(z); 701 if (j>=0x40900000) { /* z >= 1024 */ 702 if(((j-0x40900000)|i)!=0) /* if z > 1024 */ 703 return s*huge*huge; /* overflow */ 704 else { 705 if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ 706 } 707 } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */ 708 if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ 709 return s*tiny*tiny; /* underflow */ 710 else { 711 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ 712 } 713 } 714 /* 715 * compute 2**(p_h+p_l) 716 */ 717 i = j&0x7fffffff; 718 k = (i>>20)-0x3ff; 719 n = 0; 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 */ 723 t = zero; 724 __HI(t) = (n&~(0x000fffff>>k)); 725 n = ((n&0x000fffff)|0x00100000)>>(20-k); 726 if(j<0) n = -n; 727 p_h -= t; 728 } 729 t = p_l+p_h; 730 __LO(t) = 0; 731 u = t*lg2_h; 732 v = (p_l-(t-p_h))*lg2+t*lg2_l; 733 z = u+v; 734 w = v-(z-u); 735 t = z*z; 736 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 737 r = (z*t1)/(t1-two)-(w+z*w); 738 z = one-(r-z); 739 j = __HI(z); 740 j += (n<<20); 741 if((j>>20)<=0) z = fdlibmScalbn(z,n); /* subnormal output */ 742 else __HI(z) += (n<<20); 743 return s*z; 744} 745 746#endif 747 748} // namespace JSC 749