Part of Slepp's Projects — Pastebin — TURL — Imagebin — Filebin
Feedback -- English French German Japanese
Create Upload Newest Tools Donate
Sign In | Create Account

Paste Actions

Download Raw
Embed
Post to Twitter
View Tree
Remove Paste

Paste Details

C Source format, 2 years 7 months old, and 124 hits.

Advertising

Paste Description for mp-chudnovsky.c

morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 1000000000 0 # 1 000 000 000
#terms=70513669, depth=28
sieve Segmentation Fault (core dumped)
morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 100000000 0 # 100 000 000
#terms=7051366, depth=24
sieve time = 1.179
..................................................

bs time = 512.977
gcd time = 0.000
div time = 134.902
sqrt time = 25.185
mul time = 17.535
total time = 691.889
P size=145605885 digits (1.456059)
Q size=145605879 digits (1.456059)
morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 500000000 0 # 500 000 000
#terms=35256834, depth=27
sieve time = 6.331
...................................................
GNU MP: Cannot allocate memory (size=738213888)
Abort (core dumped)
morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 500000000 0 # 400 000 000
#terms=35256834, depth=27
sieve time = 6.523
...................................................
GNU MP: Cannot allocate memory (size=738213888)



^CAbort (core dumped)
morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 400000000 0 # 400 000 000
#terms=28205467, depth=26
sieve time = 5.128
...................................................

bs time = 2906.343
gcd time = 0.000
div time = 679.893
sqrt time = 118.614
mul time = 80.787
total time = 3791.225
morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 400000000 0 # 400 000 000
#terms=28205467, depth=26
sieve time = 5.307
........GNU MP: Cannot allocate memory (size=37765120)
^CAbort (core dumped)
morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out 400000000 0 # 400 000 000
#terms=28205467, depth=26
sieve time = 5.147
...................................................

bs time = 2904.718
gcd time = 0.000
div time = 679.525
sqrt time = 118.587
mul time = 81.760
total time = 3790.189
P size=582427117 digits (1.456068)
Q size=582427111 digits (1.456068)
morten@opensolaris.2009.06:~$

mp-chudnovsky.c
Wednesday, October 27th, 2010 at 4:19:08pm MDT 

  1. /* Pi computation using Chudnovsky's algortithm.
  2.  
  3. * Copyright 2002, 2005 Hanhong Xue (macroxue at yahoo dot com)
  4.  
  5. * Slightly modified 2005 by Torbjorn Granlund to allow more than 2G
  6.    digits to be computed.
  7.  
  8. * Redistribution and use in source and binary forms, with or without
  9. * modification, are permitted provided that the following conditions are met:
  10. * 1. Redistributions of source code must retain the above copyright notice,
  11. * this list of conditions and the following disclaimer.
  12. * 2. Redistributions in binary form must reproduce the above copyright notice,
  13. * this list of conditions and the following disclaimer in the documentation
  14. * and/or other materials provided with the distribution.
  15. *
  16. * THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
  17. * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  18. * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO
  19. * EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  20. * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  21. * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
  22. * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
  23. * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
  24. * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  25. * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26. */
  27.  
  28. #include <assert.h>
  29. #include <math.h>
  30. #include <stdio.h>
  31. #include <stdlib.h>
  32. #include <string.h>
  33. #include <time.h>
  34. #include "gmp.h"
  35.  
  36. #define A   13591409
  37. #define B   545140134
  38. #define C   640320
  39. #define D   12
  40.  
  41. #define BITS_PER_DIGIT   3.32192809488736234787
  42. #define DIGITS_PER_ITER  14.1816474627254776555
  43. #define DOUBLE_PREC      53
  44.  
  45. #ifdef __GNUC__
  46. #define inline __inline__
  47. #endif
  48.  
  49. char *prog_name;
  50.  
  51. #if CHECK_MEMUSAGE
  52. #undef CHECK_MEMUSAGE
  53. #define CHECK_MEMUSAGE              \
  54.   do {                  \
  55.     char buf[100];                                          \
  56.     snprintf (buf, 100,       \
  57.               "ps aguxw | grep '[%c]%s'", prog_name[0], prog_name+1);   \
  58.     system (buf);                                                 \
  59.   } while (0)
  60. #else
  61. #undef CHECK_MEMUSAGE
  62. #define CHECK_MEMUSAGE
  63. #endif
  64.  
  65.  
  66. /* Return user CPU time measured in milliseconds.  */
  67.  
  68. #if !defined (__sun) \
  69.     && (defined (USG) || defined (__SVR4) || defined (_UNICOS) \
  70.         || defined (__hpux))
  71. int
  72. cputime ()
  73. {
  74.   return (int) ((double) clock () * 1000 / CLOCKS_PER_SEC);
  75. }
  76. #else
  77. #include <sys/types.h>
  78. #include <sys/time.h>
  79. #include <sys/resource.h>
  80.  
  81. int
  82. cputime ()
  83. {
  84.   struct rusage rus;
  85.  
  86.   getrusage (0, &rus);
  87.   return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
  88. }
  89. #endif
  90.  
  91. /*///////////////////////////////////////////////////////////////////////////*/
  92.  
  93. mpf_t t1, t2;
  94.  
  95. /* r = sqrt(x) */
  96. void
  97. my_sqrt_ui(mpf_t r, unsigned long x)
  98. {
  99.   unsigned long prec, bits, prec0;
  100.  
  101.   prec0 = mpf_get_prec(r);
  102.  
  103.   if (prec0<=DOUBLE_PREC) {
  104.     mpf_set_d(r, sqrt(x));
  105.     return;
  106.   }
  107.  
  108.   bits = 0;
  109.   for (prec=prec0; prec>DOUBLE_PREC;)
  110.     {
  111.       int bit = prec&1;
  112.       prec = (prec+bit)/2;
  113.       bits = bits*2+bit;
  114.     }
  115.  
  116.   mpf_set_prec_raw(t1, DOUBLE_PREC);
  117.   mpf_set_d(t1, 1/sqrt(x));
  118.  
  119.   while (prec<prec0)
  120.     {
  121.       prec *=2;
  122.       if (prec<prec0)
  123.         {
  124.           /* t1 = t1+t1*(1-x*t1*t1)/2; */
  125.           mpf_set_prec_raw(t2, prec);
  126.           mpf_mul(t2, t1, t1);         /* half x half -> full */
  127.           mpf_mul_ui(t2, t2, x);
  128.           mpf_ui_sub(t2, 1, t2);
  129.           mpf_set_prec_raw(t2, prec/2);
  130.           mpf_div_2exp(t2, t2, 1);
  131.           mpf_mul(t2, t2, t1);         /* half x half -> half */
  132.           mpf_set_prec_raw(t1, prec);
  133.           mpf_add(t1, t1, t2);
  134.         }
  135.       else
  136.         {
  137.           break;
  138.         }
  139.       prec -= (bits&1);
  140.       bits /=2;
  141.     }
  142.   /* t2=x*t1, t1 = t2+t1*(x-t2*t2)/2; */
  143.   mpf_set_prec_raw(t2, prec0/2);
  144.   mpf_mul_ui(t2, t1, x);
  145.   mpf_mul(r, t2, t2);          /* half x half -> full */
  146.   mpf_ui_sub(r, x, r);
  147.   mpf_mul(t1, t1, r);          /* half x half -> half */
  148.   mpf_div_2exp(t1, t1, 1);
  149.   mpf_add(r, t1, t2);
  150. }
  151.  
  152. /* r = y/x   WARNING: r cannot be the same as y. */
  153. #if __GMP_MP_RELEASE >= 50001
  154. #define my_div mpf_div
  155. #else
  156. void
  157. my_div(mpf_t r, mpf_t y, mpf_t x)
  158. {
  159.   unsigned long prec, bits, prec0;
  160.  
  161.   prec0 = mpf_get_prec(r);
  162.  
  163.   if (prec0<=DOUBLE_PREC) {
  164.     mpf_set_d(r, mpf_get_d(y)/mpf_get_d(x));
  165.     return;
  166.   }
  167.  
  168.   bits = 0;
  169.   for (prec=prec0; prec>DOUBLE_PREC;) {
  170.     int bit = prec&1;
  171.     prec = (prec+bit)/2;
  172.     bits = bits*2+bit;
  173.   }
  174.  
  175.   mpf_set_prec_raw(t1, DOUBLE_PREC);
  176.   mpf_ui_div(t1, 1, x);
  177.  
  178.   while (prec<prec0) {
  179.     prec *=2;
  180.     if (prec<prec0) {
  181.       /* t1 = t1+t1*(1-x*t1); */
  182.       mpf_set_prec_raw(t2, prec);
  183.       mpf_mul(t2, x, t1);          /* full x half -> full */
  184.       mpf_ui_sub(t2, 1, t2);
  185.       mpf_set_prec_raw(t2, prec/2);
  186.       mpf_mul(t2, t2, t1);         /* half x half -> half */
  187.       mpf_set_prec_raw(t1, prec);
  188.       mpf_add(t1, t1, t2);
  189.     } else {
  190.       prec = prec0;
  191.       /* t2=y*t1, t1 = t2+t1*(y-x*t2); */
  192.       mpf_set_prec_raw(t2, prec/2);
  193.       mpf_mul(t2, t1, y);          /* half x half -> half */
  194.       mpf_mul(r, x, t2);           /* full x half -> full */
  195.       mpf_sub(r, y, r);
  196.       mpf_mul(t1, t1, r);          /* half x half -> half */
  197.       mpf_add(r, t1, t2);
  198.       break;
  199.     }
  200.     prec -= (bits&1);
  201.     bits /=2;
  202.   }
  203. }
  204. #endif
  205.  
  206. /*///////////////////////////////////////////////////////////////////////////*/
  207.  
  208. #define min(x,y) ((x)<(y)?(x):(y))
  209. #define max(x,y) ((x)>(y)?(x):(y))
  210.  
  211. typedef struct {
  212.   unsigned long max_facs;
  213.   unsigned long num_facs;
  214.   unsigned long *fac;
  215.   unsigned long *pow;
  216. } fac_t[1];
  217.  
  218. typedef struct {
  219.   long int fac;
  220.   long int pow;
  221.   long int nxt;
  222. } sieve_t;
  223.  
  224. sieve_t *sieve;
  225. long int sieve_size;
  226. fac_t   ftmp, fmul;
  227.  
  228. #define INIT_FACS 32
  229.  
  230. void
  231. fac_show(fac_t f)
  232. {
  233.   long int i;
  234.   for (i=0; i<f[0].num_facs; i++)
  235.     if (f[0].pow[i]==1)
  236.       printf("%ld ", f[0].fac[i]);
  237.     else
  238.       printf("%ld^%ld ", f[0].fac[i], f[0].pow[i]);
  239.   printf("\n");
  240. }
  241.  
  242. inline void
  243. fac_reset(fac_t f)
  244. {
  245.   f[0].num_facs = 0;
  246. }
  247.  
  248. inline void
  249. fac_init_size(fac_t f, long int s)
  250. {
  251.   if (s<INIT_FACS)
  252.     s=INIT_FACS;
  253.  
  254.   f[0].fac  = malloc(s*sizeof(unsigned long)*2);
  255.   f[0].pow  = f[0].fac + s;
  256.   f[0].max_facs = s;
  257.  
  258.   fac_reset(f);
  259. }
  260.  
  261. inline void
  262. fac_init(fac_t f)
  263. {
  264.   fac_init_size(f, INIT_FACS);
  265. }
  266.  
  267. inline void
  268. fac_clear(fac_t f)
  269. {
  270.   free(f[0].fac);
  271. }
  272.  
  273. inline void
  274. fac_resize(fac_t f, long int s)
  275. {
  276.   if (f[0].max_facs < s) {
  277.     fac_clear(f);
  278.     fac_init_size(f, s);
  279.   }
  280. }
  281.  
  282. /* f = base^pow */
  283. inline void
  284. fac_set_bp(fac_t f, unsigned long base, long int pow)
  285. {
  286.   long int i;
  287.   assert(base<sieve_size);
  288.   for (i=0, base/=2; base>0; i++, base = sieve[base].nxt) {
  289.     f[0].fac[i] = sieve[base].fac;
  290.     f[0].pow[i] = sieve[base].pow*pow;
  291.   }
  292.   f[0].num_facs = i;
  293.   assert(i<=f[0].max_facs);
  294. }
  295.  
  296. /* r = f*g */
  297. inline void
  298. fac_mul2(fac_t r, fac_t f, fac_t g)
  299. {
  300.   long int i, j, k;
  301.  
  302.   for (i=j=k=0; i<f[0].num_facs && j<g[0].num_facs; k++) {
  303.     if (f[0].fac[i] == g[0].fac[j]) {
  304.       r[0].fac[k] = f[0].fac[i];
  305.       r[0].pow[k] = f[0].pow[i] + g[0].pow[j];
  306.       i++; j++;
  307.     } else if (f[0].fac[i] < g[0].fac[j]) {
  308.       r[0].fac[k] = f[0].fac[i];
  309.       r[0].pow[k] = f[0].pow[i];
  310.       i++;
  311.     } else {
  312.       r[0].fac[k] = g[0].fac[j];
  313.       r[0].pow[k] = g[0].pow[j];
  314.       j++;
  315.     }
  316.   }
  317.   for (; i<f[0].num_facs; i++, k++) {
  318.     r[0].fac[k] = f[0].fac[i];
  319.     r[0].pow[k] = f[0].pow[i];
  320.   }
  321.   for (; j<g[0].num_facs; j++, k++) {
  322.     r[0].fac[k] = g[0].fac[j];
  323.     r[0].pow[k] = g[0].pow[j];
  324.   }
  325.   r[0].num_facs = k;
  326.   assert(k<=r[0].max_facs);
  327. }
  328.  
  329. /* f *= g */
  330. inline void
  331. fac_mul(fac_t f, fac_t g)
  332. {
  333.   fac_t tmp;
  334.   fac_resize(fmul, f[0].num_facs + g[0].num_facs);
  335.   fac_mul2(fmul, f, g);
  336.   tmp[0]  = f[0];
  337.   f[0]    = fmul[0];
  338.   fmul[0] = tmp[0];
  339. }
  340.  
  341. /* f *= base^pow */
  342. inline void
  343. fac_mul_bp(fac_t f, unsigned long base, unsigned long pow)
  344. {
  345.   fac_set_bp(ftmp, base, pow);
  346.   fac_mul(f, ftmp);
  347. }
  348.  
  349. /* remove factors of power 0 */
  350. inline void
  351. fac_compact(fac_t f)
  352. {
  353.   long int i, j;
  354.   for (i=0, j=0; i<f[0].num_facs; i++) {
  355.     if (f[0].pow[i]>0) {
  356.       if (j<i) {
  357.               f[0].fac[j] = f[0].fac[i];
  358.         f[0].pow[j] = f[0].pow[i];
  359.       }
  360.       j++;
  361.     }
  362.   }
  363.   f[0].num_facs = j;
  364. }
  365.  
  366. /* convert factorized form to number */
  367. void
  368. bs_mul(mpz_t r, long int a, long int b)
  369. {
  370.   long int i, j;
  371.   if (b-a<=32) {
  372.     mpz_set_ui(r, 1);
  373.     for (i=a; i<b; i++)
  374.       for (j=0; j<fmul[0].pow[i]; j++)
  375.         mpz_mul_ui(r, r, fmul[0].fac[i]);
  376.   } else {
  377.     mpz_t r2;
  378.     mpz_init(r2);
  379.     bs_mul(r2, a, (a+b)/2);
  380.     bs_mul(r, (a+b)/2, b);
  381.     mpz_mul(r, r, r2);
  382.     mpz_clear(r2);
  383.   }
  384. }
  385.  
  386. mpz_t    gcd, mgcd;
  387.  
  388. #if HAVE_DIVEXACT_PREINV
  389. void mpz_invert_mod_2exp (mpz_ptr, mpz_srcptr);
  390. void mpz_divexact_pre (mpz_ptr, mpz_srcptr, mpz_srcptr, mpz_srcptr);
  391.  
  392. #endif
  393.  
  394. /* f /= gcd(f,g), g /= gcd(f,g) */
  395. void
  396. fac_remove_gcd(mpz_t p, fac_t fp, mpz_t g, fac_t fg)
  397. {
  398.   long int i, j, k, c;
  399.   fac_resize(fmul, min(fp->num_facs, fg->num_facs));
  400.   for (i=j=k=0; i<fp->num_facs && j<fg->num_facs; ) {
  401.     if (fp->fac[i] == fg->fac[j]) {
  402.       c = min(fp->pow[i], fg->pow[j]);
  403.       fp->pow[i] -= c;
  404.       fg->pow[j] -= c;
  405.       fmul->fac[k] = fp->fac[i];
  406.       fmul->pow[k] = c;
  407.       i++; j++; k++;
  408.     } else if (fp->fac[i] < fg->fac[j]) {
  409.       i++;
  410.     } else {
  411.       j++;
  412.     }
  413.   }
  414.   fmul->num_facs = k;
  415.   assert(k <= fmul->max_facs);
  416.  
  417.   if (fmul->num_facs) {
  418.     bs_mul(gcd, 0, fmul->num_facs);
  419. #if HAVE_DIVEXACT_PREINV
  420.     mpz_invert_mod_2exp (mgcd, gcd);
  421.     mpz_divexact_pre (p, p, gcd, mgcd);
  422.     mpz_divexact_pre (g, g, gcd, mgcd);
  423. #else
  424. #define SIZ(x) x->_mp_size
  425.     mpz_divexact(p, p, gcd);
  426.     mpz_divexact(g, g, gcd);
  427. #endif
  428.     fac_compact(fp);
  429.     fac_compact(fg);
  430.   }
  431. }
  432.  
  433. /*///////////////////////////////////////////////////////////////////////////*/
  434.  
  435. int      out=0;
  436. mpz_t   *pstack, *qstack, *gstack;
  437. fac_t  *fpstack, *fgstack;
  438. long int      top = 0;
  439. double   progress=0, percent;
  440.  
  441. #define p1 (pstack[top])
  442. #define q1 (qstack[top])
  443. #define g1 (gstack[top])
  444. #define fp1 (fpstack[top])
  445. #define fg1 (fgstack[top])
  446.  
  447. #define p2 (pstack[top+1])
  448. #define q2 (qstack[top+1])
  449. #define g2 (gstack[top+1])
  450. #define fp2 (fpstack[top+1])
  451. #define fg2 (fgstack[top+1])
  452.  
  453. long gcd_time = 0;
  454.  
  455. /* binary splitting */
  456. void
  457. bs(unsigned long a, unsigned long b, unsigned gflag, long int level)
  458. {
  459.   unsigned long i, mid;
  460.   int ccc;
  461.  
  462.   if (b-a==1) {
  463.     /*
  464.       g(b-1,b) = (6b-5)(2b-1)(6b-1)
  465.       p(b-1,b) = b^3 * C^3 / 24
  466.       q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
  467.     */
  468.     mpz_set_ui(p1, b);
  469.     mpz_mul_ui(p1, p1, b);
  470.     mpz_mul_ui(p1, p1, b);
  471.     mpz_mul_ui(p1, p1, (C/24)*(C/24));
  472.     mpz_mul_ui(p1, p1, C*24);
  473.  
  474.     mpz_set_ui(g1, 2*b-1);
  475.     mpz_mul_ui(g1, g1, 6*b-1);
  476.     mpz_mul_ui(g1, g1, 6*b-5);
  477.  
  478.     mpz_set_ui(q1, b);
  479.     mpz_mul_ui(q1, q1, B);
  480.     mpz_add_ui(q1, q1, A);
  481.     mpz_mul   (q1, q1, g1);
  482.     if (b%2)
  483.       mpz_neg(q1, q1);
  484.  
  485.     i=b;
  486.     while ((i&1)==0) i>>=1;
  487.     fac_set_bp(fp1, i, 3);      /*  b^3 */
  488.     fac_mul_bp(fp1, 3*5*23*29, 3);
  489.     fp1[0].pow[0]--;
  490.  
  491.     fac_set_bp(fg1, 2*b-1, 1);  /* 2b-1 */
  492.     fac_mul_bp(fg1, 6*b-1, 1);  /* 6b-1 */
  493.     fac_mul_bp(fg1, 6*b-5, 1);  /* 6b-5 */
  494.  
  495.     if (b>(int)(progress)) {
  496.       printf("."); fflush(stdout);
  497.       progress += percent*2;
  498.     }
  499.  
  500.   } else {
  501.     /*
  502.       p(a,b) = p(a,m) * p(m,b)
  503.       g(a,b) = g(a,m) * g(m,b)
  504.       q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
  505.     */
  506.     mid = a+(b-a)*0.5224;     /* tuning parameter */
  507.     bs(a, mid, 1, level+1);
  508.  
  509.     top++;
  510.     bs(mid, b, gflag, level+1);
  511.     top--;
  512.  
  513.     if (level == 0)
  514.       puts ("");
  515.  
  516.     ccc = level == 0;
  517.  
  518.     if (ccc) CHECK_MEMUSAGE;
  519.  
  520.     if (level>=4) {           /* tuning parameter */
  521. #if 0
  522.       long t = cputime();
  523. #endif
  524.       fac_remove_gcd(p2, fp2, g1, fg1);
  525. #if 0
  526.       gcd_time += cputime()-t;
  527. #endif
  528.     }
  529.  
  530.     if (ccc) CHECK_MEMUSAGE;
  531.     mpz_mul(p1, p1, p2);
  532.  
  533.     if (ccc) CHECK_MEMUSAGE;
  534.     mpz_mul(q1, q1, p2);
  535.  
  536.     if (ccc) CHECK_MEMUSAGE;
  537.     mpz_mul(q2, q2, g1);
  538.  
  539.     if (ccc) CHECK_MEMUSAGE;
  540.     mpz_add(q1, q1, q2);
  541.  
  542.     if (ccc) CHECK_MEMUSAGE;
  543.     fac_mul(fp1, fp2);
  544.  
  545.     if (gflag) {
  546.       mpz_mul(g1, g1, g2);
  547.       fac_mul(fg1, fg2);
  548.     }
  549.   }
  550.  
  551.   if (out&2) {
  552.     printf("p(%ld,%ld)=",a,b); fac_show(fp1);
  553.     if (gflag)
  554.       printf("g(%ld,%ld)=",a,b); fac_show(fg1);
  555.   }
  556. }
  557.  
  558. void
  559. build_sieve(long int n, sieve_t *s)
  560. {
  561.   long int m, i, j, k;
  562.  
  563.   sieve_size = n;
  564.   m = (long int)sqrt(n);
  565.   memset(s, 0, sizeof(sieve_t)*n/2);
  566.  
  567.   s[1/2].fac = 1;
  568.   s[1/2].pow = 1;
  569.  
  570.   for (i=3; i<=n; i+=2) {
  571.     if (s[i/2].fac == 0) {
  572.       s[i/2].fac = i;
  573.       s[i/2].pow = 1;
  574.       if (i<=m) {
  575.         for (j=i*i, k=i/2; j<=n; j+=i+i, k++) {
  576.           if (s[j/2].fac==0) {
  577.             s[j/2].fac = i;
  578.             if (s[k].fac == i) {
  579.               s[j/2].pow = s[k].pow + 1;
  580.               s[j/2].nxt = s[k].nxt;
  581.             } else {
  582.               s[j/2].pow = 1;
  583.               s[j/2].nxt = k;
  584.             }
  585.           }
  586.         }
  587.       }
  588.     }
  589.   }
  590. }
  591.  
  592. int
  593. main(int argc, char *argv[])
  594. {
  595.   mpf_t  pi, qi;
  596.   long int d=100, i, depth=1, terms;
  597.   unsigned long psize, qsize;
  598.   long begin, mid0, mid1, mid2, mid3, mid4, end;
  599.  
  600.   prog_name = argv[0];
  601.  
  602.   if (argc>1)
  603.     d = strtoul(argv[1], 0, 0);
  604.   if (argc>2)
  605.     out = atoi(argv[2]);
  606.  
  607.   terms = d/DIGITS_PER_ITER;
  608.   while ((1L<<depth)<terms)
  609.     depth++;
  610.   depth++;
  611.   percent = terms/100.0;
  612.   printf("#terms=%ld, depth=%ld\n", terms, depth);
  613.  
  614.   begin = cputime();
  615.   printf("sieve   "); fflush(stdout);
  616.  
  617.   sieve_size = max(3*5*23*29+1, terms*6);
  618.   sieve = (sieve_t *)malloc(sizeof(sieve_t)*sieve_size/2);
  619.   build_sieve(sieve_size, sieve);
  620.  
  621.   mid0 = cputime();
  622.   printf("time = %6.3f\n", (double)(mid0-begin)/1000);
  623.  
  624.   /* allocate stacks */
  625.   pstack = malloc(sizeof(mpz_t)*depth);
  626.   qstack = malloc(sizeof(mpz_t)*depth);
  627.   gstack = malloc(sizeof(mpz_t)*depth);
  628.   fpstack = malloc(sizeof(fac_t)*depth);
  629.   fgstack = malloc(sizeof(fac_t)*depth);
  630.   for (i=0; i<depth; i++) {
  631.     mpz_init(pstack[i]);
  632.     mpz_init(qstack[i]);
  633.     mpz_init(gstack[i]);
  634.     fac_init(fpstack[i]);
  635.     fac_init(fgstack[i]);
  636.   }
  637.   mpz_init(gcd);
  638. #if HAVE_DIVEXACT_PREINV
  639.   mpz_init(mgcd);
  640. #endif
  641.   fac_init(ftmp);
  642.   fac_init(fmul);
  643.  
  644.   /* begin binary splitting process */
  645.   if (terms<=0) {
  646.     mpz_set_ui(p2,1);
  647.     mpz_set_ui(q2,0);
  648.     mpz_set_ui(g2,1);
  649.   } else {
  650.     bs(0,terms,0,0);
  651.   }
  652.  
  653.   mid1 = cputime();
  654.   printf("\nbs      time = %6.3f\n", (double)(mid1-mid0)/1000);
  655.   printf("   gcd  time = %6.3f\n", (double)(gcd_time)/1000);
  656.  
  657.   /* printf("misc    "); fflush(stdout); */
  658.  
  659.   /* free some resources */
  660.   free(sieve);
  661.  
  662. #if HAVE_DIVEXACT_PREINV
  663.   mpz_clear(mgcd);
  664. #endif
  665.   mpz_clear(gcd);
  666.   fac_clear(ftmp);
  667.   fac_clear(fmul);
  668.  
  669.   for (i=1; i<depth; i++) {
  670.     mpz_clear(pstack[i]);
  671.     mpz_clear(qstack[i]);
  672.     mpz_clear(gstack[i]);
  673.     fac_clear(fpstack[i]);
  674.     fac_clear(fgstack[i]);
  675.   }
  676.  
  677.   mpz_clear(gstack[0]);
  678.   fac_clear(fpstack[0]);
  679.   fac_clear(fgstack[0]);
  680.  
  681.   free(gstack);
  682.   free(fpstack);
  683.   free(fgstack);
  684.  
  685.   /* prepare to convert integers to floats */
  686.   mpf_set_default_prec((long int)(d*BITS_PER_DIGIT+16));
  687.  
  688.   /*
  689.           p*(C/D)*sqrt(C)
  690.     pi = -----------------
  691.              (q+A*p)
  692.   */
  693.  
  694.   psize = mpz_sizeinbase(p1,10);
  695.   qsize = mpz_sizeinbase(q1,10);
  696.  
  697.   mpz_addmul_ui(q1, p1, A);
  698.   mpz_mul_ui(p1, p1, C/D);
  699.  
  700.   mpf_init(pi);
  701.   mpf_set_z(pi, p1);
  702.   mpz_clear(p1);
  703.  
  704.   mpf_init(qi);
  705.   mpf_set_z(qi, q1);
  706.   mpz_clear(q1);
  707.  
  708.   free(pstack);
  709.   free(qstack);
  710.  
  711.   mid2 = cputime();
  712.   /* printf("time = %6.3f\n", (double)(mid2-mid1)/1000); */
  713.  
  714.   /* initialize temp float variables for sqrt & div */
  715.   mpf_init(t1);
  716.   mpf_init(t2);
  717.   /* mpf_set_prec_raw(t1, mpf_get_prec(pi)); */
  718.  
  719.   /* final step */
  720.   printf("div     ");  fflush(stdout);
  721.   my_div(qi, pi, qi);
  722.   mid3 = cputime();
  723.   printf("time = %6.3f\n", (double)(mid3-mid2)/1000);
  724.  
  725.   printf("sqrt    ");  fflush(stdout);
  726.   my_sqrt_ui(pi, C);
  727.   mid4 = cputime();
  728.   printf("time = %6.3f\n", (double)(mid4-mid3)/1000);
  729.  
  730.   printf("mul     ");  fflush(stdout);
  731.   mpf_mul(qi, qi, pi);
  732.   end = cputime();
  733.   printf("time = %6.3f\n", (double)(end-mid4)/1000);
  734.  
  735.   printf("total   time = %6.3f\n", (double)(end-begin)/1000);
  736.   fflush(stdout);
  737.  
  738.   printf("   P size=%ld digits (%f)\n"
  739.          "   Q size=%ld digits (%f)\n",
  740.          psize, (double)psize/d, qsize, (double)qsize/d);
  741.  
  742.   /* output Pi and timing statistics */
  743.   if (out&1)  {
  744.     printf("pi(0,%ld)=\n", terms);
  745.     mpf_out_str(stdout, 10, d+2, qi);
  746.     printf("\n");
  747.   }
  748.  
  749.   /* free float resources */
  750.   mpf_clear(pi);
  751.   mpf_clear(qi);
  752.  
  753.   mpf_clear(t1);
  754.   mpf_clear(t2);
  755.   exit (0);
  756. }
  757.  
  758.  
  759.  
  760. /*
  761.  
  762.  
  763.  
  764. gcc gmp-chudnovsky.c -I /opt/csw/include  -lm  -lgmp -o  gmp-chudnovsky.out
  765.  
  766. ./gmp-chudnovsky.out 50  1
  767. #terms=3, depth=3
  768. sieve   time =  0.001
  769. ...
  770.  
  771. bs      time =  0.000
  772.    gcd  time =  0.000
  773. div     time =  0.000
  774. sqrt    time =  0.000
  775. mul     time =  0.000
  776. total   time =  0.001
  777.    P size=51 digits (1.020000)
  778.    Q size=44 digits (0.880000)
  779. pi(0,3)=
  780. 0.3141592653589793238462643383279502884197169399375106e1
  781.  
  782. */
  783.  
  784.  
  785.  
  786. /*
  787.  
  788. ./gmp-chudnovsky.out 100000  0 | grep total
  789. total   time =  0.138
  790.  
  791. ./gmp-chudnovsky.out 1000000  0 | grep total
  792. total   time =  2.221
  793.  
  794. I expect        I got
  795.  
  796.    0.03          0.140
  797.    0.48          2.169
  798.    8.2          39.323
  799. 134           691.670
  800. 2097          nothing,
  801.  
  802. */
  803.  
  804.  
  805. /*
  806. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  100000 0 | grep total
  807. total   time =  0.140
  808. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  1000000 0 | grep total
  809. total   time =  2.169
  810. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  10000000 0 | grep total
  811. total   time = 39.323
  812. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  100000000 0 | grep total
  813. total   time = 691.670
  814. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  1000000000 0 | grep total
  815. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  1000000000 0 | grep total
  816. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  100000000 0 | grep total
  817. total   time = 691.197
  818. morten@opensolaris.2009.06:~$ ./gmp-chudnovsky.out  1000000000 0 | grep total
  819. morten@opensolaris.2009.06:~$
  820.  
  821.  
  822. please confer http://gmplib.org/pi-with-gmp.html
  823.  
  824.  
  825. */

advertising

Update the Post

Either update this post and resubmit it with changes, or make a new post.

You may also comment on this post.

update paste below
details of the post (optional)

Note: Only the paste content is required, though the following information can be useful to others.

Save name / title?

(space separated, optional)



Please note that information posted here will expire by default in one month. If you do not want it to expire, please set the expiry time above. If it is set to expire, web search engines will not be allowed to index it prior to it expiring. Items that are not marked to expire will be indexable by search engines. Be careful with your passwords. All illegal activities will be reported and any information will be handed over to the authorities, so be good.

worth-right
fantasy-obligation fantasy-obligation