- 注册时间
- 2009-7-21
- 最后登录
- 1970-1-1
- 威望
- 星
- 金币
- 枚
- 贡献
- 分
- 经验
- 点
- 鲜花
- 朵
- 魅力
- 点
- 上传
- 次
- 下载
- 次
- 积分
- 4489
- 在线时间
- 小时
|
楼主 |
发表于 2010-9-7 20:47:29
|
显示全部楼层
another code
- /*
- Hi again!
- Here is a simple recursive parallel gmp-chudnovsky program using Cilk or Cilk++.
- The Cilk/Cilk++ compiler is based on GNU C, with a few extra keywords to control
- parallelism and works with the normal GNU GMP library. Using Cilk/Cilk++
- instead of C with Pthreads or OpenMP makes recursive divide-and-conquer code
- easier to implement, debug (see cilkscreen), and maybe faster.
- More information about Cilk is at http://en.wikipedia.org/wiki/Cilk.
- Also, this simple version of gmp-chudnovsky has the factorization performance
- enhancement removed for readability, so running with one core will not be as
- fast as the Hanhong Xue's original version.
- Enjoy,
- David Carver
- */
- /* Pi computation using Chudnovsky's algortithm.
- * Copyright 2002,2005 Hanhong Xue (macroxue at yahoo dot com)
- * Slightly modified 2005 by Torbjorn Granlund to allow more than 2G
- digits to be computed.
- * Modified 2010 by David Carver (dcarver at tacc dot utexas dot edu) to
- demonstrate a parallel and fully recursive version of the gmp-chudnovsky
- program, using cilk/cilk++.
- To compile with GMP:
- cilk++ -Wall -O2 -o pgmp-chudnovsky pgmp-chudnovsky.cilk -lgmp -lm -lmiser
- To run:
- ./pgmp-chudnovsky 1000 1
- or add the Cilk++ option "-cilk_set_worker_count=" to specify the number
- of workers (hardware threads used).
- ./pgmp-chudnovsky 1000 1 -cilk_set_worker_count=1
- To get help run the program with no options:
- ./pgmp-chudnovsky
- Syntax: ./pgmp-chudnovsky.cilk.bin <digits> <option>
- <digits> digits of pi to output
- <option> 0 - just run (default)
- 1 - output decimal digits to stdout
- * Redistribution and use in source and binary forms,with or without
- * modification,are permitted provided that the following conditions are met:
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- * this list of conditions and the following disclaimer in the documentation
- * and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES,INCLUDING,BUT NOT LIMITED TO,THE IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
- * EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT,INDIRECT,INCIDENTAL,
- * SPECIAL,EXEMPLARY,OR CONSEQUENTIAL DAMAGES (INCLUDING,BUT NOT LIMITED TO,
- * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,DATA,OR PROFITS;
- * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
- * WHETHER IN CONTRACT,STRICT LIABILITY,OR TORT (INCLUDING NEGLIGENCE OR
- * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,EVEN IF
- * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
- #include <assert.h>
- #include <math.h>
- #include <stdio.h>
- #include <stdlib.h>
- #include <string.h>
- #include <time.h>
- #include <sys/time.h>
- #include <gmp.h>
- #include <cilk.h>
- #define A 13591409
- #define B 545140134
- #define C 640320
- #define D 12
- #define BITS_PER_DIGIT 3.32192809488736234787
- #define DIGITS_PER_ITER 14.1816474627254776555
- #define DOUBLE_PREC 53
- ////////////////////////////////////////////////////////////////////////////
- double wall_clock()
- {
- struct timeval timeval;
- (void) gettimeofday (&timeval,NULL);
- return (timeval.tv_sec + (timeval.tv_usec / 1000000.0));
- }
- ////////////////////////////////////////////////////////////////////////////
- /* binary splitting */
- void
- bs(unsigned long a,unsigned long b,unsigned gflag,mpz_t pstack1,mpz_t qstack1,
- mpz_t gstack1)
- {
- unsigned long mid;
- mpz_t pstack2,qstack2,gstack2;
- if (b-a==1) {
- /*
- g(b-1,b) = (6b-5)(2b-1)(6b-1)
- p(b-1,b) = b^3 * C^3 / 24
- q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
- */
- mpz_set_ui(pstack1,b);
- mpz_mul_ui(pstack1,pstack1,b);
- mpz_mul_ui(pstack1,pstack1,b);
- mpz_mul_ui(pstack1,pstack1,(C/24)*(C/24));
- mpz_mul_ui(pstack1,pstack1,C*24);
- mpz_set_ui(gstack1,2*b-1);
- mpz_mul_ui(gstack1,gstack1,6*b-1);
- mpz_mul_ui(gstack1,gstack1,6*b-5);
- mpz_set_ui(qstack1,b);
- mpz_mul_ui(qstack1,qstack1,B);
- mpz_add_ui(qstack1,qstack1,A);
- mpz_mul (qstack1,qstack1,gstack1);
- if (b%2)
- mpz_neg(qstack1,qstack1);
- } else {
- /*
- p(a,b) = p(a,m) * p(m,b)
- g(a,b) = g(a,m) * g(m,b)
- q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
- */
- mid = a+(b-a)*0.5224; /* tuning parameter */
- mpz_init(pstack2);
- mpz_init(qstack2);
- mpz_init(gstack2);
- cilk_spawn bs(mid,b,gflag,pstack2,qstack2,gstack2);
- bs(a,mid,1,pstack1,qstack1,gstack1);
- cilk_sync;
- mpz_mul(pstack1,pstack1,pstack2);
- mpz_mul(qstack1,qstack1,pstack2);
- mpz_mul(qstack2,qstack2,gstack1);
- mpz_add(qstack1,qstack1,qstack2);
- if (gflag) {
- mpz_mul(gstack1,gstack1,gstack2);
- }
- mpz_clear(pstack2);
- mpz_clear(qstack2);
- mpz_clear(gstack2);
- }
- }
- int
- cilk_main(int argc,char *argv[])
- {
- mpf_t pi,qi;
- mpz_t pstack,qstack,gstack;
- long int d=100,out=0,terms;
- unsigned long psize,qsize;
- char *prog_name;
- clock_t begin,end;
- double wbegin,wend;
- prog_name = argv[0];
- if (argc==1) {
- fprintf(stderr,"\nSyntax: %s <digits> <option>\n",prog_name);
- fprintf(stderr," <digits> digits of pi to output\n");
- fprintf(stderr," <option> 0 - just run (default)\n");
- fprintf(stderr," 1 - output decimal digits to stdout\n");
- exit(1);
- }
- if (argc>1)
- d = strtoul(argv[1],0,0);
- if (argc>2)
- out = atoi(argv[2]);
- terms = d/DIGITS_PER_ITER;
- fprintf(stderr,"#terms=%ld workers=%d\n",terms,cilk::current_worker_count());
- begin = clock();
- wbegin = wall_clock();
- mpz_init(pstack);
- mpz_init(qstack);
- mpz_init(gstack);
- /* begin binary splitting process */
- if (terms<=0) {
- mpz_set_ui(pstack,1);
- mpz_set_ui(qstack,0);
- mpz_set_ui(gstack,1);
- } else {
- bs(0,terms,0,pstack,qstack,gstack);
- }
- end = clock();
- wend = wall_clock();
- fprintf(stderr,"bs cputime = %6.3f wallclock = %6.3f\n",
- (double)(end-begin)/CLOCKS_PER_SEC,(wend-wbegin));
- fflush(stderr);
- mpz_clear(gstack);
- /* prepare to convert integers to floats */
- mpf_set_default_prec((long int)(d*BITS_PER_DIGIT+16));
- /*
- p*(C/D)*sqrt(C)
- pi = -----------------
- (q+A*p)
- */
- psize = mpz_sizeinbase(pstack,10);
- qsize = mpz_sizeinbase(qstack,10);
- mpz_addmul_ui(qstack,pstack,A);
- mpz_mul_ui(pstack,pstack,C/D);
- mpf_init(pi);
- mpf_set_z(pi,pstack);
- mpz_clear(pstack);
- mpf_init(qi);
- mpf_set_z(qi,qstack);
- mpz_clear(qstack);
- /* final step */
- mpf_div(qi,pi,qi);
- mpf_sqrt_ui(pi,C);
- mpf_mul(qi,qi,pi);
- end = clock();
- wend = wall_clock();
- /* output Pi and timing statistics */
- fprintf(stderr,"total cputime = %6.3f wallclock = %6.3f\n",
- (double)(end-begin)/CLOCKS_PER_SEC,(wend-wbegin));
- fflush(stderr);
- printf(" P size=%ld digits (%f)\n"
- " Q size=%ld digits (%f)\n",
- psize,(double)psize/d,qsize,(double)qsize/d);
- if (out&1) {
- printf("pi(0,%ld)=\n",terms);
- mpf_out_str(stdout,10,d+2,qi);
- printf("\n");
- }
- /* free float resources */
- mpf_clear(pi);
- mpf_clear(qi);
- exit (0);
- }
|
|
IP卡
狗仔卡
发表于 2010-9-7 19:14:54
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
千斤顶
显身卡
发表于 2010-9-8 12:12:08
