Changeset 1019
 Timestamp:
 Oct 12, 2011, 12:01:51 AM (9 years ago)
 Location:
 trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/real.cpp
r1018 r1019 552 552 553 553 return ret * x; 554 } 555 556 real cbrt(real const &x) 557 { 558 /* if zero, return x */ 559 if (!(x.m_signexp << 1)) 560 return x; 561 562 /* Use the system's float inversion to approximate cbrt(x). First 563 * we construct a float in the [1..8[ range that has roughly the same 564 * mantissa as our real. Its exponent is 0, 1 or 2, depending on the 565 * value of x. The final exponent is 0 or 1 (special case). We use 566 * the final exponent and final mantissa to prefill the result. */ 567 union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f }; 568 v.x += ((x.m_signexp % 3) << 23); 569 v.x = x.m_mantissa[0] >> 9; 570 v.f = powf(v.f, 0.33333333333333333f); 571 572 real ret; 573 ret.m_mantissa[0] = v.x << 9; 574 575 uint32_t sign = x.m_signexp & 0x80000000u; 576 ret.m_signexp = sign; 577 578 int exponent = (x.m_signexp & 0x7fffffffu)  (1 << 30) + 1; 579 exponent = exponent / 3 + (v.x >> 23)  (u.x >> 23); 580 ret.m_signexp = (exponent + (1 << 30)  1) & 0x7fffffffu; 581 582 /* FIXME: 1+log2(BIGITS) steps of NewtonRaphson seems to be enough for 583 * convergence, but this hasn't been checked seriously. */ 584 for (int i = 1; i <= real::BIGITS; i *= 2) 585 { 586 static real third = re(real::R_3); 587 ret = third * (x / (ret * ret) + (ret << 1)); 588 } 589 590 return ret; 554 591 } 555 592 … … 609 646 } 610 647 648 real log2(real const &x) 649 { 650 /* Strategy for log2(x): see log(x). */ 651 real tmp = x; 652 if (x.m_signexp >> 31  x.m_signexp == 0) 653 { 654 tmp.m_signexp = 0xffffffffu; 655 tmp.m_mantissa[0] = 0xffffffffu; 656 return tmp; 657 } 658 tmp.m_signexp = (1 << 30)  1; 659 return (real)(x.m_signexp  (1 << 30) + 1) + fast_log(tmp) * real::R_LOG2E; 660 } 661 662 static real fast_exp(real const &x) 663 { 664 /* This fast exp method is tuned to work on the [0..1] range and 665 * no effort whatsoever was made to improve convergence outside this 666 * domain of validity. */ 667 real ret = real::R_1, fact = real::R_1, xn = x; 668 669 for (int i = 1; ; i++) 670 { 671 fact *= (real)i; 672 real newret = ret + xn / fact; 673 if (newret == ret) 674 break; 675 ret = newret; 676 xn *= x; 677 } 678 679 return ret; 680 } 681 611 682 real exp(real const &x) 612 683 { … … 629 700 int e0 = x / real::R_LN2; 630 701 real x0 = x  (real)e0 * real::R_LN2; 631 real x1 = real::R_1, fact = real::R_1, xn = x0; 632 633 for (int i = 1; ; i++) 634 { 635 fact *= (real)i; 636 real newx1 = x1 + xn / fact; 637 if (newx1 == x1) 638 break; 639 x1 = newx1; 640 xn *= x0; 641 } 642 702 real x1 = fast_exp(x0); 643 703 x1.m_signexp += e0; 644 704 return x1; 705 } 706 707 real exp2(real const &x) 708 { 709 /* Strategy for exp2(x): see strategy in exp(). */ 710 int e0 = x; 711 real x0 = x  (real)e0; 712 real x1 = fast_exp(x0 * real::R_LN2); 713 x1.m_signexp += e0; 714 return x1; 715 } 716 717 real copysign(real const &x, real const &y) 718 { 719 real ret = x; 720 ret.m_signexp &= 0x7fffffffu; 721 ret.m_signexp = y.m_signexp & 0x80000000u; 722 return ret; 645 723 } 646 724 … … 929 1007 void real::print(int ndigits) const 930 1008 { 931 real const r1 = 1, r10 = 10;932 1009 real x = *this; 933 1010 … … 942 1019 if (x.m_signexp) 943 1020 { 944 for (real div = r1, newdiv; true; div = newdiv)1021 for (real div = R_1, newdiv; true; div = newdiv) 945 1022 { 946 newdiv = div * r10;1023 newdiv = div * R_10; 947 1024 if (x < newdiv) 948 1025 { … … 952 1029 exponent++; 953 1030 } 954 for (real mul = 1, newx; true; mul *= r10)1031 for (real mul = 1, newx; true; mul *= R_10) 955 1032 { 956 1033 newx = x * mul; 957 if (newx >= r1)1034 if (newx >= R_1) 958 1035 { 959 1036 x = newx; … … 972 1049 printf("."); 973 1050 x = real(digit); 974 x *= r10;1051 x *= R_10; 975 1052 } 976 1053 
trunk/src/real.h
r1018 r1019 63 63 operator bool() const; 64 64 65 friend real fabs(real const &x); 66 67 friend real re(real const &x); 68 friend real sqrt(real const &x); 69 friend real log(real const &x); 70 friend real exp(real const &x); 71 72 friend real floor(real const &x); 73 friend real ceil(real const &x); 74 friend real round(real const &x); 75 friend real fmod(real const &x, real const &y); 76 65 /* Trigonometric functions */ 77 66 friend real sin(real const &x); 78 67 friend real cos(real const &x); … … 81 70 friend real acos(real const &x); 82 71 friend real atan(real const &x); 72 /* FIXME: atan2 */ 73 74 /* Hyperbolic functions */ 75 /* FIXME: sinh cosh tanh */ 76 77 /* Exponential and logarithmic functions */ 78 friend real exp(real const &x); 79 friend real exp2(real const &x); 80 friend real log(real const &x); 81 friend real log2(real const &x); 82 /* FIXME: frexp ldexp log10 modf */ 83 84 /* Power functions */ 85 friend real re(real const &x); 86 friend real sqrt(real const &x); 87 friend real cbrt(real const &x); 88 /* FIXME: pow */ 89 90 /* Rounding, absolute value, remainder etc. */ 91 friend real ceil(real const &x); 92 friend real copysign(real const &x, real const &y); 93 friend real floor(real const &x); 94 friend real fabs(real const &x); 95 friend real round(real const &x); 96 friend real fmod(real const &x, real const &y); 83 97 84 98 void print(int ndigits = 150) const; 
trunk/test/unit/real.cpp
r1017 r1019 35 35 LOLUNIT_ASSERT_EQUAL(a10, 10.0); 36 36 37 double b0 = log(real::R_E); 38 LOLUNIT_ASSERT_EQUAL(b0, 1.0); 39 40 double b1 = exp(re(real::R_LOG2E)); 41 LOLUNIT_ASSERT_EQUAL(b1, 2.0); 42 43 double b2 = exp(re(real::R_LOG10E)); 44 LOLUNIT_ASSERT_EQUAL(b2, 10.0); 45 46 double b3 = exp(real::R_LN2); 47 LOLUNIT_ASSERT_EQUAL(b3, 2.0); 48 49 double b4 = exp(real::R_LN10); 50 LOLUNIT_ASSERT_EQUAL(b4, 10.0); 51 52 double b5 = sin(real::R_PI); 53 double b6 = cos(real::R_PI); 54 LOLUNIT_ASSERT_DOUBLES_EQUAL(b5, 0.0, 1e100); 55 LOLUNIT_ASSERT_EQUAL(b6, 1.0); 56 57 double b7 = sin(real::R_PI_2); 58 double b8 = cos(real::R_PI_2); 59 LOLUNIT_ASSERT_EQUAL(b7, 1.0); 60 LOLUNIT_ASSERT_DOUBLES_EQUAL(b8, 0.0, 1e100); 61 62 double b9 = sin(real::R_PI_4) * sin(real::R_PI_4); 63 double b10 = cos(real::R_PI_4) * cos(real::R_PI_4); 64 LOLUNIT_ASSERT_EQUAL(b9, 0.5); 65 LOLUNIT_ASSERT_EQUAL(b10, 0.5); 37 double b1 = log(real::R_E); 38 double b2 = log2(real::R_2); 39 LOLUNIT_ASSERT_EQUAL(b1, 1.0); 40 LOLUNIT_ASSERT_EQUAL(b2, 1.0); 41 42 double c1 = exp(re(real::R_LOG2E)); 43 double c2 = log(exp2(real::R_LOG2E)); 44 LOLUNIT_ASSERT_EQUAL(c1, 2.0); 45 LOLUNIT_ASSERT_EQUAL(c2, 1.0); 46 47 double d1 = exp(re(real::R_LOG10E)); 48 LOLUNIT_ASSERT_EQUAL(d1, 10.0); 49 50 double e1 = exp(real::R_LN2); 51 LOLUNIT_ASSERT_EQUAL(e1, 2.0); 52 53 double f1 = exp(real::R_LN10); 54 LOLUNIT_ASSERT_EQUAL(f1, 10.0); 55 56 double g1 = sin(real::R_PI); 57 double g2 = cos(real::R_PI); 58 LOLUNIT_ASSERT_DOUBLES_EQUAL(g1, 0.0, 1e100); 59 LOLUNIT_ASSERT_EQUAL(g2, 1.0); 60 61 double h1 = sin(real::R_PI_2); 62 double h2 = cos(real::R_PI_2); 63 LOLUNIT_ASSERT_EQUAL(h1, 1.0); 64 LOLUNIT_ASSERT_DOUBLES_EQUAL(h2, 0.0, 1e100); 65 66 double i1 = sin(real::R_PI_4) * sin(real::R_PI_4); 67 double i2 = cos(real::R_PI_4) * cos(real::R_PI_4); 68 LOLUNIT_ASSERT_EQUAL(i1, 0.5); 69 LOLUNIT_ASSERT_EQUAL(i2, 0.5); 66 70 } 67 71
Note: See TracChangeset
for help on using the changeset viewer.