Changeset 1013
 Timestamp:
 Oct 10, 2011, 3:33:09 AM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/real.cpp
r1006 r1013 488 488 ret.m_signexp = (exponent  1) & 0x7fffffffu; 489 489 490 /* Five steps of NewtonRaphson seems enough for 32bigit reals. */ 491 real two = 2; 492 ret = ret * (two  ret * x); 493 ret = ret * (two  ret * x); 494 ret = ret * (two  ret * x); 495 ret = ret * (two  ret * x); 496 ret = ret * (two  ret * x); 490 /* FIXME: log2(BIGITS) steps of NewtonRaphson seems to be enough for 491 * convergence, but this hasn't been checked seriously. */ 492 for (int i = 1; i < real::BIGITS; i *= 2) 493 ret = ret * (real::R_2  ret * x); 497 494 498 495 return ret; … … 533 530 534 531 uint32_t exponent = (x.m_signexp & 0x7fffffffu); 535 exponent = ((1 << 30) + (1 << 29)  1)  (exponent + 1) / 2;532 exponent = ((1 << 30) + (1 << 29)  1)  (exponent + 1) / 2; 536 533 exponent = exponent + (v.x >> 23)  (u.x >> 23); 537 534 ret.m_signexp = exponent & 0x7fffffffu; 538 535 539 /* Five steps of NewtonRaphson seems enough for 32bigit reals. */ 540 real three = 3; 541 ret = ret * (three  ret * ret * x); 542 ret.m_signexp; 543 ret = ret * (three  ret * ret * x); 544 ret.m_signexp; 545 ret = ret * (three  ret * ret * x); 546 ret.m_signexp; 547 ret = ret * (three  ret * ret * x); 548 ret.m_signexp; 549 ret = ret * (three  ret * ret * x); 550 ret.m_signexp; 536 /* FIXME: log2(BIGITS) steps of NewtonRaphson seems to be enough for 537 * convergence, but this hasn't been checked seriously. */ 538 for (int i = 1; i < real::BIGITS; i *= 2) 539 { 540 ret = ret * (real::R_3  ret * ret * x); 541 ret.m_signexp; 542 } 551 543 552 544 return ret * x; … … 560 552 } 561 553 562 static real fast log(real const &x)554 static real fast_log(real const &x) 563 555 { 564 556 /* This fast log method is tuned to work on the [1..2] range and … … 578 570 * sqrt() call. */ 579 571 real y = sqrt(x); 580 real z = (y  (real)1) / (y + (real)1), z2 = z * z, zn = z2; 581 real sum = 1.0; 582 583 for (int i = 3; i < 200; i += 2) 584 { 585 sum += zn / (real)i; 572 real z = (y  real::R_1) / (y + real::R_1), z2 = z * z, zn = z2; 573 real sum = real::R_1; 574 575 for (int i = 3; ; i += 2) 576 { 577 real newsum = sum + zn / (real)i; 578 if (newsum == sum) 579 break; 580 sum = newsum; 586 581 zn *= z2; 587 582 } … … 590 585 } 591 586 592 static real LOG_2 = fastlog((real)2);593 594 587 real log(real const &x) 595 588 { 596 589 /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M), 597 * with the property that M is in [1..2[, so fast log() applies here. */590 * with the property that M is in [1..2[, so fast_log() applies here. */ 598 591 real tmp = x; 599 592 if (x.m_signexp >> 31  x.m_signexp == 0) … … 604 597 } 605 598 tmp.m_signexp = (1 << 30)  1; 606 return (real)(x.m_signexp  (1 << 30) + 1) * LOG_2 + fastlog(tmp);599 return (real)(x.m_signexp  (1 << 30) + 1) * real::R_LN2 + fast_log(tmp); 607 600 } 608 601 … … 625 618 * return x1 * 2^E0 626 619 */ 627 int e0 = x / LOG_2;628 real x0 = x  (real)e0 * LOG_2;629 real x1 = 1.0, fact = 1.0, xn = x0;630 631 for (int i = 1; i < 100; i++)620 int e0 = x / real::R_LN2; 621 real x0 = x  (real)e0 * real::R_LN2; 622 real x1 = real::R_1, fact = real::R_1, xn = x0; 623 624 for (int i = 1; ; i++) 632 625 { 633 626 fact *= (real)i; 634 x1 += xn / fact; 627 real newx1 = x1 + xn / fact; 628 if (newx1 == x1) 629 break; 630 x1 = newx1; 635 631 xn *= x0; 636 632 } … … 720 716 absx = real::R_PI  absx; 721 717 722 real ret = 0.0, fact = 1.0, xn = absx, x2 = absx * absx;718 real ret = real::R_0, fact = real::R_1, xn = absx, x2 = absx * absx; 723 719 for (int i = 1; ; i += 2) 724 720 { 725 721 real newret = ret + xn / fact; 726 if ( ret == newret)722 if (newret == ret) 727 723 break; 728 724 ret = newret; … … 756 752 757 753 real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1; 758 for (int i = 1; i < 280; i++)754 for (int i = 1; ; i++) 759 755 { 760 756 xn *= x2; 761 ret += (fact1 * xn / ((real)(2 * i + 1) * fact2)) >> (i * 2); 757 real mul = (real)(2 * i + 1); 758 real newret = ret + ((fact1 * xn / (mul * fact2)) >> (i * 2)); 759 if (newret == ret) 760 break; 761 ret = newret; 762 762 fact1 *= (real)((2 * i + 1) * (2 * i + 2)); 763 763 fact2 *= (real)((i + 1) * (i + 1)); … … 820 820 { 821 821 real ret = x, xn = x, mx2 = x * x; 822 for (int i = 3; i < 100; i += 2)822 for (int i = 3; ; i += 2) 823 823 { 824 824 xn *= mx2; 825 ret += xn / (real)i; 825 real newret = ret + xn / (real)i; 826 if (newret == ret) 827 break; 828 ret = newret; 826 829 } 827 830 return ret; … … 834 837 real y = real::R_1  absx; 835 838 real yn = y, my2 = y * y; 836 for (int i = 0; i < 200; i += 2)839 for (int i = 0; ; i += 2) 837 840 { 838 re t += (yn / (real)(2 * i + 1)) >> (i + 1);841 real newret = ret + ((yn / (real)(2 * i + 1)) >> (i + 1)); 839 842 yn *= y; 840 ret += (yn / (real)(2 * i + 2)) >> (i + 1);843 newret += (yn / (real)(2 * i + 2)) >> (i + 1); 841 844 yn *= y; 842 ret += (yn / (real)(2 * i + 3)) >> (i + 2); 845 newret += (yn / (real)(2 * i + 3)) >> (i + 2); 846 if (newret == ret) 847 break; 848 ret = newret; 843 849 yn *= my2; 844 850 } … … 849 855 real y = (absx  real::R_SQRT3) >> 1; 850 856 real yn = y, my2 = y * y; 851 for (int i = 1; i < 200; i += 6)857 for (int i = 1; ; i += 6) 852 858 { 853 re t += (yn / (real)i) >> 1;859 real newret = ret + ((yn / (real)i) >> 1); 854 860 yn *= y; 855 ret = (real::R_SQRT3 >> 1) * yn / (real)(i + 1);861 newret = (real::R_SQRT3 >> 1) * yn / (real)(i + 1); 856 862 yn *= y; 857 ret += yn / (real)(i + 2);863 newret += yn / (real)(i + 2); 858 864 yn *= y; 859 ret = (real::R_SQRT3 >> 1) * yn / (real)(i + 3);865 newret = (real::R_SQRT3 >> 1) * yn / (real)(i + 3); 860 866 yn *= y; 861 ret += (yn / (real)(i + 4)) >> 1; 867 newret += (yn / (real)(i + 4)) >> 1; 868 if (newret == ret) 869 break; 870 ret = newret; 862 871 yn *= my2; 863 872 } … … 869 878 real yn = y, my2 = y * y; 870 879 ret = y; 871 for (int i = 3; i < 120; i += 2)880 for (int i = 3; ; i += 2) 872 881 { 873 882 yn *= my2; 874 ret += yn / (real)i; 883 real newret = ret + yn / (real)i; 884 if (newret == ret) 885 break; 886 ret = newret; 875 887 } 876 888 ret = real::R_PI_2  ret; … … 945 957 real const m0 = x0 * x0, m1 = x1 * x1, r16 = 16.0, r4 = 4.0; 946 958 947 /* Degree 240 is required for 512bit mantissa precision */ 948 for (int i = 1; i < 240; i += 2) 949 { 950 ret += r16 / (x0 * (real)i)  r4 / (x1 * (real)i); 959 for (int i = 1; ; i += 2) 960 { 961 real newret = ret + r16 / (x0 * (real)i)  r4 / (x1 * (real)i); 962 if (newret == ret) 963 break; 964 ret = newret; 951 965 x0 *= m0; 952 966 x1 *= m1; … … 962 976 real const real::R_10 = (real)10.0; 963 977 964 real const real::R_E = exp(R_1); 965 real const real::R_LN2 = log(R_2); 978 real const real::R_LN2 = fast_log(R_2); 966 979 real const real::R_LN10 = log(R_10); 967 980 real const real::R_LOG2E = re(R_LN2); 968 981 real const real::R_LOG10E = re(R_LN10); 982 real const real::R_E = exp(R_1); 969 983 real const real::R_PI = fast_pi(); 970 984 real const real::R_PI_2 = R_PI >> 1;
Note: See TracChangeset
for help on using the changeset viewer.