root / Pi / C / OpenMP / Pi_OpenMP.c @ 31
History  View  Annotate  Download (2.3 kB)
1 
//


2 
// Estimation of Pi using Monte Carlo exploration process

3 
// gcc std=c99 O3 o Pi Pi.c lm

4 
//

5  
6 
#include <math.h> 
7 
#include <stdio.h> 
8 
#include <stdlib.h> 
9 
#include <omp.h> 
10 
#include <limits.h> 
11  
12 
// Marsaglia RNG very simple implementation

13 
#define znew ((z=36969*(z&65535)+(z>>16))<<16) 
14 
#define wnew ((w=18000*(w&65535)+(w>>16))&65535) 
15 
#define MWC (znew+wnew)

16 
#define SHR3 (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5)) 
17 
#define CONG (jcong=69069*jcong+1234567) 
18 
#define KISS ((MWC^CONG)+SHR3)

19  
20 
#define MWCfp MWC * 2.328306435454494e10f 
21 
#define KISSfp KISS * 2.328306435454494e10f 
22  
23 
#define ITERATIONS 1000000000 
24  
25 
#define PROCESS 4 
26  
27 
#ifdef LONG

28 
#define LENGTH long long 
29 
#else

30 
#define LENGTH int 
31 
#endif

32  
33 
LENGTH MainLoopGlobal(LENGTH iterations,unsigned int seed_w,unsigned int seed_z) 
34 
{ 
35 
unsigned int z=seed_z; 
36 
unsigned int w=seed_w; 
37  
38 
LENGTH total=0;

39  
40 
for (LENGTH i=0;i<iterations;i++) { 
41  
42 
float x=MWCfp ;

43 
float y=MWCfp ;

44  
45 
// Matching test

46 
int inside=((x*x+y*y) < 1.0f) ? 1:0; 
47 
total+=inside; 
48 
} 
49  
50 
return(total);

51 
} 
52  
53 
int main(int argc, char *argv[]) { 
54  
55 
unsigned int seed_w=10,seed_z=10,process=PROCESS; 
56 
LENGTH iterations=ITERATIONS; 
57 
LENGTH inside[1024],insides=0; 
58  
59 
if (argc > 1) { 
60 
iterations=(LENGTH)atoll(argv[1]);

61 
process=atoi(argv[2]);

62 
} 
63 
else {

64 
printf("\n\tPi : Estimate Pi with Monte Carlo exploration\n\n");

65 
printf("\t\t#1 : number of iterations (default 1 billion)\n");

66 
printf("\t\t#2 : number of process (default 4)\n\n");

67 
} 
68  
69 
printf ("\n\tInformation about architecture:\n\n");

70  
71 
printf ("\tSizeof int = %lld bytes.\n", (long long)sizeof(int)); 
72 
printf ("\tSizeof long = %lld bytes.\n", (long long)sizeof(long)); 
73 
printf ("\tSizeof long long = %lld bytes.\n\n", (long long)sizeof(long long)); 
74  
75 
printf ("\tMax int = %u\n", INT_MAX);

76 
printf ("\tMax long = %ld\n", LONG_MAX);

77 
printf ("\tMax long long = %lld\n\n", LLONG_MAX);

78  
79 
#pragma omp parallel for 
80 
for (int i=0 ; i<process; i++) { 
81 
inside[i]=MainLoopGlobal(iterations/process,seed_w,seed_z); 
82 
printf("\tFound %lld for process %i\n",(long long)inside[i],i); 
83 
} 
84 
printf("\n");

85  
86 
for (int i=0 ; i<process; i++) { 
87 
insides+=inside[i]; 
88 
} 
89 

90 
float pi=4.*(float)insides/(float)iterations; 
91  
92 
printf("\tPi=%f with error %f and %lld iterations\n\n",pi,

93 
fabs(pi4*atan(1))/pi,(long long)iterations); 
94 

95 
} 