<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:x="urn:schemas-microsoft-com:office:excel" xmlns:p="urn:schemas-microsoft-com:office:powerpoint" xmlns:a="urn:schemas-microsoft-com:office:access" xmlns:dt="uuid:C2F41010-65B3-11d1-A29F-00AA00C14882" xmlns:s="uuid:BDC6E3F0-6DA3-11d1-A2A3-00AA00C14882" xmlns:rs="urn:schemas-microsoft-com:rowset" xmlns:z="#RowsetSchema" xmlns:b="urn:schemas-microsoft-com:office:publisher" xmlns:ss="urn:schemas-microsoft-com:office:spreadsheet" xmlns:c="urn:schemas-microsoft-com:office:component:spreadsheet" xmlns:odc="urn:schemas-microsoft-com:office:odc" xmlns:oa="urn:schemas-microsoft-com:office:activation" xmlns:html="http://www.w3.org/TR/REC-html40" xmlns:q="http://schemas.xmlsoap.org/soap/envelope/" xmlns:D="DAV:" xmlns:x2="http://schemas.microsoft.com/office/excel/2003/xml" xmlns:ois="http://schemas.microsoft.com/sharepoint/soap/ois/" xmlns:dir="http://schemas.microsoft.com/sharepoint/soap/directory/" xmlns:ds="http://www.w3.org/2000/09/xmldsig#" xmlns:dsp="http://schemas.microsoft.com/sharepoint/dsp" xmlns:udc="http://schemas.microsoft.com/data/udc" xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:sub="http://schemas.microsoft.com/sharepoint/soap/2002/1/alerts/" xmlns:ec="http://www.w3.org/2001/04/xmlenc#" xmlns:sp="http://schemas.microsoft.com/sharepoint/" xmlns:sps="http://schemas.microsoft.com/sharepoint/soap/" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:udcxf="http://schemas.microsoft.com/data/udc/xmlfile" xmlns:wf="http://schemas.microsoft.com/sharepoint/soap/workflow/" xmlns:mver="http://schemas.openxmlformats.org/markup-compatibility/2006" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns:mrels="http://schemas.openxmlformats.org/package/2006/relationships" xmlns:ex12t="http://schemas.microsoft.com/exchange/services/2006/types" xmlns:ex12m="http://schemas.microsoft.com/exchange/services/2006/messages" xmlns:Z="urn:schemas-microsoft-com:" xmlns:st="" xmlns="http://www.w3.org/TR/REC-html40"> <head> <meta http-equiv=Content-Type content="text/html; charset=us-ascii"> <meta name=Generator content="Microsoft Word 12 (filtered medium)"> <style>  </style>  </head> <body lang=EN-US link=blue vlink=purple> <div class=Section1> Hi!<o:p></o:p> Since GNU GCC 4.3 now supports OpenMP and most new AMD and Intel processors have more than one core, I thought a parallel version of the gmp-chudnovsky program would be an interesting exercise.  This simple parallel version of gmp-chudnovsky has the factorization performance enhancement removed, so running with one core will not be as fast as the original version.  I hope to post a parallel version with the factorization performance enhancement in the near future.<o:p></o:p> <o:p> </o:p> Later,<o:p></o:p> David Carver<o:p></o:p>  <o:p></o:p> <pre>/* Pi computation using Chudnovsky's algortithm.<o:p></o:p></pre><pre><o:p> </o:p></pre><pre> * Copyright 2002, 2005 Hanhong Xue (macroxue at yahoo dot com)<o:p></o:p></pre><pre><o:p> </o:p></pre><pre> * Slightly modified 2005 by Torbjorn Granlund (tege at swox dot com) to allow<o:p></o:p></pre><pre>   more than 2G digits to be computed.<o:p></o:p></pre><pre><o:p> </o:p></pre><pre> * Modifed 2008 by David Carver (dcarver at tacc dot utexas dot edu) to enable<o:p></o:p></pre><pre>   multi-threading using the algorithm from "Computation of High-Precision <o:p></o:p></pre><pre>   Mathematical Constants in a Combined Cluster and Grid Environment" by <o:p></o:p></pre><pre>   Daisuke Takahashi, Mitsuhisa Sato, and Taisuke Boku.  <o:p></o:p></pre><pre><o:p> </o:p></pre><pre>   For GNU gcc 4.3<o:p></o:p></pre><pre>   gcc -fopenmp -Wall -O2 -o pgmp-chudnovsky pchudnovsky.c -lgmp -lm<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>   For Intel 10.1 compiler<o:p></o:p></pre><pre>   icc -openmp -O2 -o pgmp-chudnovsky pgmp-chudnovsky.c -lgmp -lm<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>   For AIX xlc<o:p></o:p></pre><pre>   xlc_r -qsmp=omp -O2 -o pgmp-chudnovsky pchudnovsky.c  -lgmp -lm<o:p></o:p></pre><pre><o:p> </o:p></pre><pre> * Redistribution and use in source and binary forms, with or without<o:p></o:p></pre><pre> * modification, are permitted provided that the following conditions are met:<o:p></o:p></pre><pre> * 1. Redistributions of source code must retain the above copyright notice,<o:p></o:p></pre><pre> * this list of conditions and the following disclaimer.<o:p></o:p></pre><pre> * 2. Redistributions in binary form must reproduce the above copyright notice,<o:p></o:p></pre><pre> * this list of conditions and the following disclaimer in the documentation<o:p></o:p></pre><pre> * and/or other materials provided with the distribution.<o:p></o:p></pre><pre> *<o:p></o:p></pre><pre> * THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR<o:p></o:p></pre><pre> * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF<o:p></o:p></pre><pre> * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO<o:p></o:p></pre><pre> * EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,<o:p></o:p></pre><pre> * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,<o:p></o:p></pre><pre> * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;<o:p></o:p></pre><pre> * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,<o:p></o:p></pre><pre> * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR<o:p></o:p></pre><pre> * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF<o:p></o:p></pre><pre> * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.<o:p></o:p></pre><pre> */<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>#include <assert.h><o:p></o:p></pre><pre>#include <math.h><o:p></o:p></pre><pre>#include <stdio.h><o:p></o:p></pre><pre>#include <stdlib.h><o:p></o:p></pre><pre>#include <string.h><o:p></o:p></pre><pre>#include <time.h><o:p></o:p></pre><pre>#include "gmp.h"<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>#define A   13591409<o:p></o:p></pre><pre>#define B   545140134<o:p></o:p></pre><pre>#define C   640320<o:p></o:p></pre><pre>#define D   12<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>#define BITS_PER_DIGIT   3.32192809488736234787<o:p></o:p></pre><pre>#define DIGITS_PER_ITER  14.1816474627254776555<o:p></o:p></pre><pre>#define DOUBLE_PREC      53<o:p></o:p></pre><pre><o:p> </o:p></pre><pre><o:p> </o:p></pre><pre>char *prog_name;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>#if CHECK_MEMUSAGE<o:p></o:p></pre><pre>#undef CHECK_MEMUSAGE<o:p></o:p></pre><pre>#define CHECK_MEMUSAGE                                               \<o:p></o:p></pre><pre>  do {                                                               \<o:p></o:p></pre><pre>    char buf[100];                                                   \<o:p></o:p></pre><pre>    snprintf (buf, 100,                                                     \<o:p></o:p></pre><pre>              "ps aguxw | grep '[%c]%s'", prog_name[0], prog_name+1);        \<o:p></o:p></pre><pre>    system (buf);                                                    \<o:p></o:p></pre><pre>  } while (0)<o:p></o:p></pre><pre>#else<o:p></o:p></pre><pre>#undef CHECK_MEMUSAGE<o:p></o:p></pre><pre>#define CHECK_MEMUSAGE<o:p></o:p></pre><pre>#endif<o:p></o:p></pre><pre><o:p> </o:p></pre><pre><o:p> </o:p></pre><pre>////////////////////////////////////////////////////////////////////////////<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>mpf_t t1, t2;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>// r = sqrt(x)<o:p></o:p></pre><pre>void<o:p></o:p></pre><pre>my_sqrt_ui(mpf_t r, unsigned long x)<o:p></o:p></pre><pre>{<o:p></o:p></pre><pre>  unsigned long prec, bits, prec0;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  prec0 = mpf_get_prec(r);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  if (prec0<=DOUBLE_PREC) {<o:p></o:p></pre><pre>    mpf_set_d(r, sqrt(x));<o:p></o:p></pre><pre>    return;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  bits = 0;<o:p></o:p></pre><pre>  for (prec=prec0; prec>DOUBLE_PREC;) {<o:p></o:p></pre><pre>    int bit = prec&1;<o:p></o:p></pre><pre>    prec = (prec+bit)/2;<o:p></o:p></pre><pre>    bits = bits*2+bit;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mpf_set_prec_raw(t1, DOUBLE_PREC);<o:p></o:p></pre><pre>  mpf_set_d(t1, 1/sqrt(x));<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  while (prec<prec0) {<o:p></o:p></pre><pre>    prec *=2;<o:p></o:p></pre><pre>    if (prec<prec0) {<o:p></o:p></pre><pre>      /* t1 = t1+t1*(1-x*t1*t1)/2; */<o:p></o:p></pre><pre>      mpf_set_prec_raw(t2, prec);<o:p></o:p></pre><pre>      mpf_mul(t2, t1, t1);         // half x half -> full<o:p></o:p></pre><pre>      mpf_mul_ui(t2, t2, x);<o:p></o:p></pre><pre>      mpf_ui_sub(t2, 1, t2);<o:p></o:p></pre><pre>      mpf_set_prec_raw(t2, prec/2);<o:p></o:p></pre><pre>      mpf_div_2exp(t2, t2, 1);<o:p></o:p></pre><pre>      mpf_mul(t2, t2, t1);         // half x half -> half<o:p></o:p></pre><pre>      mpf_set_prec_raw(t1, prec);<o:p></o:p></pre><pre>      mpf_add(t1, t1, t2);<o:p></o:p></pre><pre>    } else {<o:p></o:p></pre><pre>      prec = prec0;<o:p></o:p></pre><pre>      /* t2=x*t1, t1 = t2+t1*(x-t2*t2)/2; */<o:p></o:p></pre><pre>      mpf_set_prec_raw(t2, prec/2);<o:p></o:p></pre><pre>      mpf_mul_ui(t2, t1, x);<o:p></o:p></pre><pre>      mpf_mul(r, t2, t2);          // half x half -> full<o:p></o:p></pre><pre>      mpf_ui_sub(r, x, r);<o:p></o:p></pre><pre>      mpf_mul(t1, t1, r);          // half x half -> half<o:p></o:p></pre><pre>      mpf_div_2exp(t1, t1, 1);<o:p></o:p></pre><pre>      mpf_add(r, t1, t2);<o:p></o:p></pre><pre>      break;<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre>    prec -= (bits&1);<o:p></o:p></pre><pre>    bits /=2;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>}<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>// r = y/x   WARNING: r cannot be the same as y.<o:p></o:p></pre><pre>void<o:p></o:p></pre><pre>my_div(mpf_t r, mpf_t y, mpf_t x)<o:p></o:p></pre><pre>{<o:p></o:p></pre><pre>  unsigned long prec, bits, prec0;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  prec0 = mpf_get_prec(r);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  if (prec0<=DOUBLE_PREC) {<o:p></o:p></pre><pre>    mpf_set_d(r, mpf_get_d(y)/mpf_get_d(x));<o:p></o:p></pre><pre>    return;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  bits = 0;<o:p></o:p></pre><pre>  for (prec=prec0; prec>DOUBLE_PREC;) {<o:p></o:p></pre><pre>    int bit = prec&1;<o:p></o:p></pre><pre>    prec = (prec+bit)/2;<o:p></o:p></pre><pre>    bits = bits*2+bit;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mpf_set_prec_raw(t1, DOUBLE_PREC);<o:p></o:p></pre><pre>  mpf_ui_div(t1, 1, x);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  while (prec<prec0) {<o:p></o:p></pre><pre>    prec *=2;<o:p></o:p></pre><pre>    if (prec<prec0) {<o:p></o:p></pre><pre>      /* t1 = t1+t1*(1-x*t1); */<o:p></o:p></pre><pre>      mpf_set_prec_raw(t2, prec);<o:p></o:p></pre><pre>      mpf_mul(t2, x, t1);          // full x half -> full<o:p></o:p></pre><pre>      mpf_ui_sub(t2, 1, t2);<o:p></o:p></pre><pre>      mpf_set_prec_raw(t2, prec/2);<o:p></o:p></pre><pre>      mpf_mul(t2, t2, t1);         // half x half -> half<o:p></o:p></pre><pre>      mpf_set_prec_raw(t1, prec);<o:p></o:p></pre><pre>      mpf_add(t1, t1, t2);<o:p></o:p></pre><pre>    } else {<o:p></o:p></pre><pre>      prec = prec0;<o:p></o:p></pre><pre>      /* t2=y*t1, t1 = t2+t1*(y-x*t2); */<o:p></o:p></pre><pre>      mpf_set_prec_raw(t2, prec/2);<o:p></o:p></pre><pre>      mpf_mul(t2, t1, y);          // half x half -> half<o:p></o:p></pre><pre>      mpf_mul(r, x, t2);           // full x half -> full<o:p></o:p></pre><pre>      mpf_sub(r, y, r);<o:p></o:p></pre><pre>      mpf_mul(t1, t1, r);          // half x half -> half<o:p></o:p></pre><pre>      mpf_add(r, t1, t2);<o:p></o:p></pre><pre>      break;<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre>    prec -= (bits&1);<o:p></o:p></pre><pre>    bits /=2;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>}<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>////////////////////////////////////////////////////////////////////////////<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>int      out=0;<o:p></o:p></pre><pre>mpz_t   **pstack, **qstack, **gstack;<o:p></o:p></pre><pre>long int cores=1, depth, cores_depth;<o:p></o:p></pre><pre>double   progress=0, percent;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>// binary splitting<o:p></o:p></pre><pre>void<o:p></o:p></pre><pre>sum(unsigned long i, unsigned long j, unsigned long gflag)<o:p></o:p></pre><pre>{<o:p></o:p></pre><pre>  mpz_mul(pstack[i][0], pstack[i][0], pstack[j][0]);<o:p></o:p></pre><pre>  mpz_mul(qstack[i][0], qstack[i][0], pstack[j][0]);<o:p></o:p></pre><pre>  mpz_mul(qstack[j][0], qstack[j][0], gstack[i][0]);<o:p></o:p></pre><pre>  mpz_add(qstack[i][0], qstack[i][0], qstack[j][0]);<o:p></o:p></pre><pre>  if (gflag) {<o:p></o:p></pre><pre>     mpz_mul(gstack[i][0], gstack[i][0], gstack[j][0]);<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>}<o:p></o:p></pre><pre>void<o:p></o:p></pre><pre>bs(unsigned long a, unsigned long b, unsigned long gflag, unsigned long level, unsigned long index, unsigned long top)<o:p></o:p></pre><pre>{<o:p></o:p></pre><pre>  unsigned long mid;<o:p></o:p></pre><pre>  int ccc;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  if (out&2) {<o:p></o:p></pre><pre>    fprintf(stderr,"bs: a = %ld b = %ld gflag = %ld index = %ld level = %ld top = %ld \n", a,b,gflag,index,level,top);<o:p></o:p></pre><pre>    fflush(stderr);<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  if ((b > a) && (b-a==1)) {<o:p></o:p></pre><pre>    /*<o:p></o:p></pre><pre>      g(b-1,b) = (6b-5)(2b-1)(6b-1)<o:p></o:p></pre><pre>      p(b-1,b) = b^3 * C^3 / 24<o:p></o:p></pre><pre>      q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).<o:p></o:p></pre><pre>    */<o:p></o:p></pre><pre>    mpz_set_ui(pstack[index][top], b);<o:p></o:p></pre><pre>    mpz_mul_ui(pstack[index][top], pstack[index][top], b);<o:p></o:p></pre><pre>    mpz_mul_ui(pstack[index][top], pstack[index][top], b);<o:p></o:p></pre><pre>    mpz_mul_ui(pstack[index][top], pstack[index][top], (C/24)*(C/24));<o:p></o:p></pre><pre>    mpz_mul_ui(pstack[index][top], pstack[index][top], C*24);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    mpz_set_ui(gstack[index][top], 2*b-1);<o:p></o:p></pre><pre>    mpz_mul_ui(gstack[index][top], gstack[index][top], 6*b-1);<o:p></o:p></pre><pre>    mpz_mul_ui(gstack[index][top], gstack[index][top], 6*b-5);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    mpz_set_ui(qstack[index][top], b);<o:p></o:p></pre><pre>    mpz_mul_ui(qstack[index][top], qstack[index][top], B);<o:p></o:p></pre><pre>    mpz_add_ui(qstack[index][top], qstack[index][top], A);<o:p></o:p></pre><pre>    mpz_mul   (qstack[index][top], qstack[index][top], gstack[index][top]);<o:p></o:p></pre><pre>    if (b%2)<o:p></o:p></pre><pre>      mpz_neg(qstack[index][top], qstack[index][top]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    if (b>(int)(progress)) {<o:p></o:p></pre><pre>      fprintf(stderr,"."); fflush(stderr);<o:p></o:p></pre><pre>      progress += percent*2;<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  } else {<o:p></o:p></pre><pre>    /*<o:p></o:p></pre><pre>      p(a,b) = p(a,m) * p(m,b)<o:p></o:p></pre><pre>      g(a,b) = g(a,m) * g(m,b)<o:p></o:p></pre><pre>      q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)<o:p></o:p></pre><pre>    */<o:p></o:p></pre><pre>    mid = a+(b-a)*0.5224;     // tuning parameter<o:p></o:p></pre><pre>    bs(a, mid, 1, level+1, index, top);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    bs(mid, b, gflag, level+1, index, top+1);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    ccc = level == 0;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    if (ccc) CHECK_MEMUSAGE;<o:p></o:p></pre><pre>    mpz_mul(pstack[index][top], pstack[index][top], pstack[index][top+1]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    if (ccc) CHECK_MEMUSAGE;<o:p></o:p></pre><pre>    mpz_mul(qstack[index][top], qstack[index][top], pstack[index][top+1]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    if (ccc) CHECK_MEMUSAGE;<o:p></o:p></pre><pre>    mpz_mul(qstack[index][top+1], qstack[index][top+1], gstack[index][top]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    if (ccc) CHECK_MEMUSAGE;<o:p></o:p></pre><pre>    mpz_add(qstack[index][top], qstack[index][top], qstack[index][top+1]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    if (gflag) {<o:p></o:p></pre><pre>      mpz_mul(gstack[index][top], gstack[index][top], gstack[index][top+1]);<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>}<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>int<o:p></o:p></pre><pre>main(int argc, char *argv[])<o:p></o:p></pre><pre>{<o:p></o:p></pre><pre>  mpf_t  pi, qi;<o:p></o:p></pre><pre>  long int d=100, terms, i, j, k, cores_size;<o:p></o:p></pre><pre>  unsigned long psize, qsize, mid;<o:p></o:p></pre><pre>  clock_t begin, mid0, mid1, mid2, mid3, mid4, end;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  prog_name = argv[0];<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  if (argc==1) {<o:p></o:p></pre><pre>    fprintf(stderr,"\nSyntax: %s <digits> <option> <cores>\n",prog_name);<o:p></o:p></pre><pre>    fprintf(stderr,"      <digits> digits of pi to output\n");<o:p></o:p></pre><pre>    fprintf(stderr,"      <option> 0 - just run (default)\n");<o:p></o:p></pre><pre>    fprintf(stderr,"               1 - output digits\n");<o:p></o:p></pre><pre>    fprintf(stderr,"               2 - debug\n");<o:p></o:p></pre><pre>    fprintf(stderr,"      <cores> number of cores (default 1)\n");<o:p></o:p></pre><pre>    exit(1);<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>  if (argc>1)<o:p></o:p></pre><pre>    d = strtoul(argv[1], 0, 0);<o:p></o:p></pre><pre>  if (argc>2)<o:p></o:p></pre><pre>    out = atoi(argv[2]);<o:p></o:p></pre><pre>  if (argc>3)<o:p></o:p></pre><pre>    cores = atoi(argv[3]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  terms = d/DIGITS_PER_ITER;<o:p></o:p></pre><pre>  depth = 0;<o:p></o:p></pre><pre>  while ((1L<<depth)<terms)<o:p></o:p></pre><pre>    depth++;<o:p></o:p></pre><pre>  depth++;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  if (cores < 1) {<o:p></o:p></pre><pre>        fprintf(stderr,"Number of cores reset from %ld to 1\n",cores); <o:p></o:p></pre><pre>        fflush(stderr);<o:p></o:p></pre><pre>        cores = 1;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>  if ((terms > 0) && (terms < cores)) {<o:p></o:p></pre><pre>        fprintf(stderr,"Number of cores reset from %ld to %ld\n",cores,terms); <o:p></o:p></pre><pre>        fflush(stderr);<o:p></o:p></pre><pre>        cores = terms;<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre>  cores_depth = 0;<o:p></o:p></pre><pre>  while ((1L<<cores_depth)<cores)<o:p></o:p></pre><pre>    cores_depth++;<o:p></o:p></pre><pre>  cores_size=pow(2,cores_depth);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  percent = terms/100.0;<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  fprintf(stderr,"#terms=%ld, depth=%ld, cores=%ld\n", terms, depth, cores);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  begin = mid0 = clock();<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  pstack = malloc(sizeof(mpz_t)*cores);<o:p></o:p></pre><pre>  qstack = malloc(sizeof(mpz_t)*cores);<o:p></o:p></pre><pre>  gstack = malloc(sizeof(mpz_t)*cores);<o:p></o:p></pre><pre>  /* allocate stacks */<o:p></o:p></pre><pre>  for (j = 0; j < cores_size; j++) {<o:p></o:p></pre><pre>    pstack[j] = malloc(sizeof(mpz_t)*depth);<o:p></o:p></pre><pre>    qstack[j] = malloc(sizeof(mpz_t)*depth);<o:p></o:p></pre><pre>    gstack[j] = malloc(sizeof(mpz_t)*depth);<o:p></o:p></pre><pre>    for (i = 0; i < depth; i++) {<o:p></o:p></pre><pre>      mpz_init(pstack[j][i]);<o:p></o:p></pre><pre>      mpz_init(qstack[j][i]);<o:p></o:p></pre><pre>      mpz_init(gstack[j][i]);<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre> <o:p></o:p></pre><pre>  /* begin binary splitting process */<o:p></o:p></pre><pre>  if (terms<=0) {<o:p></o:p></pre><pre>    mpz_set_ui(pstack[0][0],1);<o:p></o:p></pre><pre>    mpz_set_ui(qstack[0][0],0);<o:p></o:p></pre><pre>    mpz_set_ui(gstack[0][0],1);<o:p></o:p></pre><pre>  } else {<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    mid0 = clock();<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    mid = terms / cores; <o:p></o:p></pre><pre><o:p> </o:p></pre><pre>#ifdef _OPENMP<o:p></o:p></pre><pre>#pragma omp parallel for<o:p></o:p></pre><pre>#endif<o:p></o:p></pre><pre>    for (i = 0; i < cores; i++) {<o:p></o:p></pre><pre>      if (i < (cores-1))<o:p></o:p></pre><pre>         bs(i*mid, (i+1)*mid, 1, 0, i, 0);<o:p></o:p></pre><pre>      else<o:p></o:p></pre><pre>         bs(i*mid, terms, 1, 0, i, 0);<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    for (k = 1; k < cores_size; k*=2) {<o:p></o:p></pre><pre>#ifdef _OPENMP<o:p></o:p></pre><pre>#pragma omp parallel for<o:p></o:p></pre><pre>#endif<o:p></o:p></pre><pre>      for (i = 0; i < cores; i=i+2*k) {<o:p></o:p></pre><pre>        if (i+k < cores) {<o:p></o:p></pre><pre>          sum( i, i+k, 1);<o:p></o:p></pre><pre>        }<o:p></o:p></pre><pre>      }<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>    for (j=0; j<cores_size; j++) {<o:p></o:p></pre><pre>      free(gstack[j]);<o:p></o:p></pre><pre>    }<o:p></o:p></pre><pre>    free(gstack);<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mid1 = clock();<o:p></o:p></pre><pre>  fprintf(stderr,"\nbs      time = %6.3f\n", (double)(mid1-mid0)/CLOCKS_PER_SEC);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  /* prepare to convert integers to floats */<o:p></o:p></pre><pre>  mpf_set_default_prec((long int)(d*BITS_PER_DIGIT+16));<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  /*<o:p></o:p></pre><pre>          p*(C/D)*sqrt(C)<o:p></o:p></pre><pre>    pi = -----------------<o:p></o:p></pre><pre>             (q+A*p)<o:p></o:p></pre><pre>  */<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  psize = mpz_sizeinbase(pstack[0][0],10);<o:p></o:p></pre><pre>  qsize = mpz_sizeinbase(qstack[0][0],10);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mpz_addmul_ui(qstack[0][0], pstack[0][0], A);<o:p></o:p></pre><pre>  mpz_mul_ui(pstack[0][0], pstack[0][0], C/D);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mpf_init(pi);<o:p></o:p></pre><pre>  mpf_set_z(pi, pstack[0][0]);<o:p></o:p></pre><pre>  mpz_clear(pstack[0][0]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mpf_init(qi);<o:p></o:p></pre><pre>  mpf_set_z(qi, qstack[0][0]);<o:p></o:p></pre><pre>  mpz_clear(qstack[0][0]);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  free(pstack[0]);<o:p></o:p></pre><pre>  free(qstack[0]);<o:p></o:p></pre><pre>  free(pstack);<o:p></o:p></pre><pre>  free(qstack);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mid2 = clock();<o:p></o:p></pre><pre>  //fprintf(stderr,"time = %6.3f\n", (double)(mid2-mid1)/CLOCKS_PER_SEC);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  /* initialize temp float variables for sqrt & div */<o:p></o:p></pre><pre>  mpf_init(t1);<o:p></o:p></pre><pre>  mpf_init(t2);<o:p></o:p></pre><pre>  //mpf_set_prec_raw(t1, mpf_get_prec(pi));<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  /* final step */<o:p></o:p></pre><pre>  fprintf(stderr,"div     ");  fflush(stderr);<o:p></o:p></pre><pre>  my_div(qi, pi, qi);<o:p></o:p></pre><pre>  mid3 = clock();<o:p></o:p></pre><pre>  fprintf(stderr,"time = %6.3f\n", (double)(mid3-mid2)/CLOCKS_PER_SEC);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  fprintf(stderr,"sqrt    ");  fflush(stderr);<o:p></o:p></pre><pre>  my_sqrt_ui(pi, C);<o:p></o:p></pre><pre>  mid4 = clock();<o:p></o:p></pre><pre>  fprintf(stderr,"time = %6.3f\n", (double)(mid4-mid3)/CLOCKS_PER_SEC);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  fprintf(stderr,"mul     ");  fflush(stderr);<o:p></o:p></pre><pre>  mpf_mul(qi, qi, pi);<o:p></o:p></pre><pre>  end = clock();<o:p></o:p></pre><pre>  fprintf(stderr,"time = %6.3f\n", (double)(end-mid4)/CLOCKS_PER_SEC);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  fprintf(stderr,"total   time = %6.3f\n", (double)(end-begin)/CLOCKS_PER_SEC);<o:p></o:p></pre><pre>  fflush(stderr);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  fprintf(stderr,"   P size=%ld digits (%f)\n"<o:p></o:p></pre><pre>        "   Q size=%ld digits (%f)\n",<o:p></o:p></pre><pre>        psize, (double)psize/d, qsize, (double)qsize/d);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  /* output Pi and timing statistics */<o:p></o:p></pre><pre>  if (out&1)  {<o:p></o:p></pre><pre>    fprintf(stdout,"pi(0,%ld)=\n", terms);<o:p></o:p></pre><pre>    mpf_out_str(stdout, 10, d+2, qi);<o:p></o:p></pre><pre>    fprintf(stdout,"\n");<o:p></o:p></pre><pre>  }<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  /* free float resources */<o:p></o:p></pre><pre>  mpf_clear(pi);<o:p></o:p></pre><pre>  mpf_clear(qi);<o:p></o:p></pre><pre><o:p> </o:p></pre><pre>  mpf_clear(t1);<o:p></o:p></pre><pre>  mpf_clear(t2);<o:p></o:p></pre><pre>  exit (0);<o:p></o:p></pre><pre>}<o:p></o:p></pre> <o:p> </o:p> </div> </body> </html>