OLD | NEW |
1 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm), | 1 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm), |
2 // | 2 // |
3 // ==================================================== | 3 // ==================================================== |
4 // Copyright (C) 1993-2004 by Sun Microsystems, Inc. All rights reserved. | 4 // Copyright (C) 1993-2004 by Sun Microsystems, Inc. All rights reserved. |
5 // | 5 // |
6 // Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 // Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 // Permission to use, copy, modify, and distribute this | 7 // Permission to use, copy, modify, and distribute this |
8 // software is freely granted, provided that this notice | 8 // software is freely granted, provided that this notice |
9 // is preserved. | 9 // is preserved. |
10 // ==================================================== | 10 // ==================================================== |
(...skipping 379 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
390 function MathTan(x) { | 390 function MathTan(x) { |
391 x = x * 1; // Convert to number. | 391 x = x * 1; // Convert to number. |
392 if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) { | 392 if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) { |
393 // |x| < pi/4, approximately. No reduction needed. | 393 // |x| < pi/4, approximately. No reduction needed. |
394 return KernelTan(x, 0, 1); | 394 return KernelTan(x, 0, 1); |
395 } | 395 } |
396 REMPIO2(x); | 396 REMPIO2(x); |
397 return KernelTan(y0, y1, (n & 1) ? -1 : 1); | 397 return KernelTan(y0, y1, (n & 1) ? -1 : 1); |
398 } | 398 } |
399 | 399 |
400 // ES6 draft 09-27-13, section 20.2.2.20. | |
401 // Math.log1p | |
402 // | |
403 // Method : | |
404 // 1. Argument Reduction: find k and f such that | |
405 // 1+x = 2^k * (1+f), | |
406 // where sqrt(2)/2 < 1+f < sqrt(2) . | |
407 // | |
408 // Note. If k=0, then f=x is exact. However, if k!=0, then f | |
409 // may not be representable exactly. In that case, a correction | |
410 // term is need. Let u=1+x rounded. Let c = (1+x)-u, then | |
411 // log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u), | |
412 // and add back the correction term c/u. | |
413 // (Note: when x > 2**53, one can simply return log(x)) | |
414 // | |
415 // 2. Approximation of log1p(f). | |
416 // Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) | |
417 // = 2s + 2/3 s**3 + 2/5 s**5 + ....., | |
418 // = 2s + s*R | |
419 // We use a special Reme algorithm on [0,0.1716] to generate | |
420 // a polynomial of degree 14 to approximate R The maximum error | |
421 // of this polynomial approximation is bounded by 2**-58.45. In | |
422 // other words, | |
423 // 2 4 6 8 10 12 14 | |
424 // R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s | |
425 // (the values of Lp1 to Lp7 are listed in the program) | |
426 // and | |
427 // | 2 14 | -58.45 | |
428 // | Lp1*s +...+Lp7*s - R(z) | <= 2 | |
429 // | | | |
430 // Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. | |
431 // In order to guarantee error in log below 1ulp, we compute log | |
432 // by | |
433 // log1p(f) = f - (hfsq - s*(hfsq+R)). | |
434 // | |
435 // 3. Finally, log1p(x) = k*ln2 + log1p(f). | |
436 // = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) | |
437 // Here ln2 is split into two floating point number: | |
438 // ln2_hi + ln2_lo, | |
439 // where n*ln2_hi is always exact for |n| < 2000. | |
440 // | |
441 // Special cases: | |
442 // log1p(x) is NaN with signal if x < -1 (including -INF) ; | |
443 // log1p(+INF) is +INF; log1p(-1) is -INF with signal; | |
444 // log1p(NaN) is that NaN with no signal. | |
445 // | |
446 // Accuracy: | |
447 // according to an error analysis, the error is always less than | |
448 // 1 ulp (unit in the last place). | |
449 // | |
450 // Constants: | |
451 // Constants are found in fdlibm.cc. We assume the C++ compiler to convert | |
452 // from decimal to binary accurately enough to produce the intended values. | |
453 // | |
454 // Note: Assuming log() return accurate answer, the following | |
455 // algorithm can be used to compute log1p(x) to within a few ULP: | |
456 // | |
457 // u = 1+x; | |
458 // if (u==1.0) return x ; else | |
459 // return log(u)*(x/(u-1.0)); | |
460 // | |
461 // See HP-15C Advanced Functions Handbook, p.193. | |
462 // | |
463 define LN2_HI = 6.93147180369123816490e-01; | 400 define LN2_HI = 6.93147180369123816490e-01; |
464 define LN2_LO = 1.90821492927058770002e-10; | 401 define LN2_LO = 1.90821492927058770002e-10; |
465 define TWO_THIRD = 6.666666666666666666e-01; | |
466 define LP1 = 6.666666666666735130e-01; | |
467 define LP2 = 3.999999999940941908e-01; | |
468 define LP3 = 2.857142874366239149e-01; | |
469 define LP4 = 2.222219843214978396e-01; | |
470 define LP5 = 1.818357216161805012e-01; | |
471 define LP6 = 1.531383769920937332e-01; | |
472 define LP7 = 1.479819860511658591e-01; | |
473 | 402 |
474 // 2^54 | 403 // 2^54 |
475 define TWO54 = 18014398509481984; | 404 define TWO54 = 18014398509481984; |
476 | 405 |
477 function MathLog1p(x) { | |
478 x = x * 1; // Convert to number. | |
479 var hx = %_DoubleHi(x); | |
480 var ax = hx & 0x7fffffff; | |
481 var k = 1; | |
482 var f = x; | |
483 var hu = 1; | |
484 var c = 0; | |
485 var u = x; | |
486 | |
487 if (hx < 0x3fda827a) { | |
488 // x < 0.41422 | |
489 if (ax >= 0x3ff00000) { // |x| >= 1 | |
490 if (x === -1) { | |
491 return -INFINITY; // log1p(-1) = -inf | |
492 } else { | |
493 return NaN; // log1p(x<-1) = NaN | |
494 } | |
495 } else if (ax < 0x3c900000) { | |
496 // For |x| < 2^-54 we can return x. | |
497 return x; | |
498 } else if (ax < 0x3e200000) { | |
499 // For |x| < 2^-29 we can use a simple two-term Taylor series. | |
500 return x - x * x * 0.5; | |
501 } | |
502 | |
503 if ((hx > 0) || (hx <= -0x402D413D)) { // (int) 0xbfd2bec3 = -0x402d413d | |
504 // -.2929 < x < 0.41422 | |
505 k = 0; | |
506 } | |
507 } | |
508 | |
509 // Handle Infinity and NaN | |
510 if (hx >= 0x7ff00000) return x; | |
511 | |
512 if (k !== 0) { | |
513 if (hx < 0x43400000) { | |
514 // x < 2^53 | |
515 u = 1 + x; | |
516 hu = %_DoubleHi(u); | |
517 k = (hu >> 20) - 1023; | |
518 c = (k > 0) ? 1 - (u - x) : x - (u - 1); | |
519 c = c / u; | |
520 } else { | |
521 hu = %_DoubleHi(u); | |
522 k = (hu >> 20) - 1023; | |
523 } | |
524 hu = hu & 0xfffff; | |
525 if (hu < 0x6a09e) { | |
526 u = %_ConstructDouble(hu | 0x3ff00000, %_DoubleLo(u)); // Normalize u. | |
527 } else { | |
528 ++k; | |
529 u = %_ConstructDouble(hu | 0x3fe00000, %_DoubleLo(u)); // Normalize u/2. | |
530 hu = (0x00100000 - hu) >> 2; | |
531 } | |
532 f = u - 1; | |
533 } | |
534 | |
535 var hfsq = 0.5 * f * f; | |
536 if (hu === 0) { | |
537 // |f| < 2^-20; | |
538 if (f === 0) { | |
539 if (k === 0) { | |
540 return 0.0; | |
541 } else { | |
542 return k * LN2_HI + (c + k * LN2_LO); | |
543 } | |
544 } | |
545 var R = hfsq * (1 - TWO_THIRD * f); | |
546 if (k === 0) { | |
547 return f - R; | |
548 } else { | |
549 return k * LN2_HI - ((R - (k * LN2_LO + c)) - f); | |
550 } | |
551 } | |
552 | |
553 var s = f / (2 + f); | |
554 var z = s * s; | |
555 var R = z * (LP1 + z * (LP2 + z * (LP3 + z * (LP4 + | |
556 z * (LP5 + z * (LP6 + z * LP7)))))); | |
557 if (k === 0) { | |
558 return f - (hfsq - s * (hfsq + R)); | |
559 } else { | |
560 return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f); | |
561 } | |
562 } | |
563 | |
564 // ES6 draft 09-27-13, section 20.2.2.14. | 406 // ES6 draft 09-27-13, section 20.2.2.14. |
565 // Math.expm1 | 407 // Math.expm1 |
566 // Returns exp(x)-1, the exponential of x minus 1. | 408 // Returns exp(x)-1, the exponential of x minus 1. |
567 // | 409 // |
568 // Method | 410 // Method |
569 // 1. Argument reduction: | 411 // 1. Argument reduction: |
570 // Given x, find r and integer k such that | 412 // Given x, find r and integer k such that |
571 // | 413 // |
572 // x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 | 414 // x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 |
573 // | 415 // |
(...skipping 526 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1100 | 942 |
1101 utils.InstallFunctions(GlobalMath, DONT_ENUM, [ | 943 utils.InstallFunctions(GlobalMath, DONT_ENUM, [ |
1102 "cos", MathCos, | 944 "cos", MathCos, |
1103 "sin", MathSin, | 945 "sin", MathSin, |
1104 "tan", MathTan, | 946 "tan", MathTan, |
1105 "sinh", MathSinh, | 947 "sinh", MathSinh, |
1106 "cosh", MathCosh, | 948 "cosh", MathCosh, |
1107 "tanh", MathTanh, | 949 "tanh", MathTanh, |
1108 "log10", MathLog10, | 950 "log10", MathLog10, |
1109 "log2", MathLog2, | 951 "log2", MathLog2, |
1110 "log1p", MathLog1p, | |
1111 "expm1", MathExpm1 | 952 "expm1", MathExpm1 |
1112 ]); | 953 ]); |
1113 | 954 |
1114 %SetForceInlineFlag(MathSin); | 955 %SetForceInlineFlag(MathSin); |
1115 %SetForceInlineFlag(MathCos); | 956 %SetForceInlineFlag(MathCos); |
1116 | 957 |
1117 }) | 958 }) |
OLD | NEW |