Changeset 999
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/real.cpp
r998 r999 668 668 } 669 669 670 real atan(real const &x) 671 { 672 /* Computing atan(x): we choose a different Taylor series depending on 673 * the value of x to help with convergence. 674 * 675 * If |x| < 0.5 we evaluate atan(y) near 0: 676 * atan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 ... 677 * 678 * If 0.5 <= |x| < 1.5 we evaluate atan(1+y) near 0: 679 * atan(1+y) = π/4 + y/(1*2^1) - y^2/(2*2^1) + y^3/(3*2^2) 680 * - y^5/(5*2^3) + y^6/(6*2^3) - y^7/(7*2^4) 681 * + y^9/(9*2^5) - y^10/(10*2^5) + y^11/(11*2^6) ... 682 * 683 * If 1.5 <= |x| < 2 we evaluate atan(sqrt(3)+2y) near 0: 684 * atan(sqrt(3)+2y) = π/3 + 1/2 y - sqrt(3)/2 y^2/2 685 * + y^3/3 - sqrt(3)/2 y^4/4 + 1/2 y^5/5 686 * - 1/2 y^7/7 + sqrt(3)/2 y^8/8 687 * - y^9/9 + sqrt(3)/2 y^10/10 - 1/2 y^11/11 688 * + 1/2 y^13/13 - sqrt(3)/2 y^14/14 689 * + y^15/15 - sqrt(3)/2 y^16/16 + 1/2 y^17/17 ... 690 * 691 * If |x| >= 2 we evaluate atan(y) near +∞: 692 * atan(y) = π/2 - y^-1 + y^-3/3 - y^-5/5 + y^-7/7 - y^-9/9 ... 693 */ 694 real absx = fabs(x); 695 696 if (absx < (real::R_1 >> 1)) 697 { 698 real ret = x, xn = x, mx2 = -x * x; 699 for (int i = 3; i < 100; i += 2) 700 { 701 xn *= mx2; 702 ret += xn / (real)i; 703 } 704 return ret; 705 } 706 707 real ret = 0; 708 709 if (absx < (real::R_3 >> 1)) 710 { 711 real y = real::R_1 - absx; 712 real yn = y, my2 = -y * y; 713 for (int i = 0; i < 200; i += 2) 714 { 715 ret += (yn / (real)(2 * i + 1)) >> (i + 1); 716 yn *= y; 717 ret += (yn / (real)(2 * i + 2)) >> (i + 1); 718 yn *= y; 719 ret += (yn / (real)(2 * i + 3)) >> (i + 2); 720 yn *= my2; 721 } 722 ret = real::R_PI_4 - ret; 723 } 724 else if (absx < real::R_2) 725 { 726 real y = (absx - real::R_SQRT3) >> 1; 727 real yn = y, my2 = -y * y; 728 for (int i = 1; i < 200; i += 6) 729 { 730 ret += (yn / (real)i) >> 1; 731 yn *= y; 732 ret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1); 733 yn *= y; 734 ret += yn / (real)(i + 2); 735 yn *= y; 736 ret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3); 737 yn *= y; 738 ret += (yn / (real)(i + 4)) >> 1; 739 yn *= my2; 740 } 741 ret = real::R_PI_3 + ret; 742 } 743 else 744 { 745 real y = re(absx); 746 real yn = y, my2 = -y * y; 747 ret = y; 748 for (int i = 3; i < 120; i += 2) 749 { 750 yn *= my2; 751 ret += yn / (real)i; 752 } 753 ret = real::R_PI_2 - ret; 754 } 755 756 /* Propagate sign */ 757 ret.m_signexp |= (x.m_signexp & 0x80000000u); 758 return ret; 759 } 760 670 761 void real::print(int ndigits) const 671 762 { … … 745 836 real const real::R_1 = (real)1.0; 746 837 real const real::R_2 = (real)2.0; 838 real const real::R_3 = (real)3.0; 747 839 real const real::R_10 = (real)10.0; 748 840 … … 754 846 real const real::R_PI = fast_pi(); 755 847 real const real::R_PI_2 = R_PI >> 1; 848 real const real::R_PI_3 = R_PI / R_3; 756 849 real const real::R_PI_4 = R_PI >> 2; 757 850 real const real::R_1_PI = re(R_PI); … … 759 852 real const real::R_2_SQRTPI = re(sqrt(R_PI)) << 1; 760 853 real const real::R_SQRT2 = sqrt(R_2); 854 real const real::R_SQRT3 = sqrt(R_3); 761 855 real const real::R_SQRT1_2 = R_SQRT2 >> 1; 762 856 -
trunk/src/real.h
r998 r999 71 71 friend real sin(real const &x); 72 72 friend real cos(real const &x); 73 friend real atan(real const &x); 73 74 74 75 void print(int ndigits = 150) const; … … 77 78 static real const R_1; 78 79 static real const R_2; 80 static real const R_3; 79 81 static real const R_10; 80 82 … … 86 88 static real const R_PI; 87 89 static real const R_PI_2; 90 static real const R_PI_3; 88 91 static real const R_PI_4; 89 92 static real const R_1_PI; … … 91 94 static real const R_2_SQRTPI; 92 95 static real const R_SQRT2; 96 static real const R_SQRT3; 93 97 static real const R_SQRT1_2; 94 98
Note: See TracChangeset
for help on using the changeset viewer.