1. /* sdcrpi.cpp by K.Tsuru */
  2. // function ID 3508 DRADIX, constant pi since version 2.30
  3. /**************************************************************
  4. pi by Ramanujan's formula used by W.Gosper
  5. 2sqrt(2) (4n)!(1103+26390n)
  6. 1/pi = --------- sum_{n=0}^{\inf} ------------------
  7. 99^2 (4^n 99^n n!)^4
  8. sqrt(8) B(0)...B(n-1)
  9. =-------- sum_{n=0}^{\inf} --------------- A(n)
  10. 99^2 C(0)...C(n)
  11. Adding one term the presision rises by eight digits.
  12. See "sdcchpi.cpp" for detail.
  13. Binary Splitting method Nov.13,2016
  14. 28.4(sec) for 4190029 digits
  15. 69.5(sec) 8000029
  16. ***************************************************************/
  17. #ifndef SN_H
  18. #include "sn.h"
  19. #endif
  20. static const SLong D=1103L;
  21. static const SLong Ec=26390L; // To distinguish from E().
  22. static const SLong F=4L*99L;
  23. static const SLong G=99L*99L;
  24. static double upPerTerm = 7.98;
  25. #define BS_RJ_PI_SLONG 1 // SLong version is faster
  26. #if BS_RJ_PI_SLONG
  27. ////// SLong version /////// 2.36(sec) for 100000 digits
  28. class BSRamanujanPi : public BinarySplitting<SLong> {
  29. public:
  30. /* Constructer : 'L' is upToTerm, 'f' is precision*/
  31. BSRamanujanPi(long L, long f) : BinarySplitting<SLong>(L, f){}
  32. void setABC(long k, SLong& a, SLong& b, SLong& c) {
  33. // A(k)
  34. a = D + Ec * k;
  35. // B(k)= 8(k+1)(4k+3)(2k+1)(4k+1)
  36. if(k+1L < LONG_MAX/8L){
  37. long d =8*(k+1L), p = 4L*k+3L, q = 2L*k+1L, r = p-2L;
  38. b= d; b *= p; b *= q; b *= r;
  39. }else{
  40. SLong d =k+1L, p = 4L*k+3L, q = 2L*k+1L, r = p-2L;
  41. b= 8; b *= d; b *= p; b *= q; b *= r;
  42. }
  43. // C(k)
  44. if(k){ c = k*F; c *= c; c*=c; } // c = Lpow(k*F, 4);
  45. // if(k) c = Lpow(k*F, 4);
  46. else c = 1;
  47. }
  48. SDouble getValue() {
  49. putTogether();
  50. SDouble s = BinarySplitting<SLong>::getValue(true), pi = G*RecSqrt(8.0)*s;
  51. return pi;
  52. }
  53. };
  54. SDouble RamanujanPi() {
  55. SDouble C;
  56. long f = long(C.EffFig() + C.Hidden()+ 1u)*DFIGURES;
  57. long L = (long)((double)f/upPerTerm);
  58. BSRamanujanPi pi(L, f);
  59. return pi.getValue();
  60. }
  61. #else
  62. ////// SDouble version /////// 3.75(sec) for 100000 digits
  63. class BSRamanujanPi : public BinarySplitting<SDouble> {
  64. public:
  65. /* Constructer : 'L' is upToTerm, 'f' is precision*/
  66. BSRamanujanPi(long L, long f) : BinarySplitting<SDouble>(L, f){}
  67. void setABC(long k, SDouble& a, SDouble& b, SDouble& c) {
  68. // A(k)
  69. a = D + Ec * k;
  70. // B(k)= 8(k+1)(4k+3)(2k+1)(4k+1)
  71. SLong d =8L*(k+1L), p = 4L*k+3L, q = 2L*k+1L, r = p-2L;
  72. b = d; b *= p; b *= q; b *= r;
  73. // C(k)
  74. if(k){ c = k*F; c *= c; c*=c; } // c = Lpow(k*F, 4);
  75. else c = 1;
  76. }
  77. SDouble getValue() {
  78. putTogether();
  79. SDouble s = BinarySplitting<SDouble>::getValue(true), pi = G*RecSqrt(8.0)*s;
  80. return pi;
  81. }
  82. };
  83. SDouble RamanujanPi() {
  84. SDouble C;
  85. long f = long( C.EffFig() + C.Hidden() ), L = long(upRate*double(f));
  86. BSRamanujanPi pi(L, f);
  87. return pi.getValue();
  88. }
  89. #endif // BS_RJ_PI_SLONG

sdcRjnpi.cpp : last modifiled at 2016/11/20 10:05:50(2,983 bytes)
created at 2017/10/07 10:21:15
The creation time of this html file is 2017/10/07 10:30:03 (Sat Oct 07 10:30:03 2017).