Statistics
| Revision:

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.328306435454494e-10f
21
#define KISSfp KISS * 2.328306435454494e-10f
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(pi-4*atan(1))/pi,(long long)iterations);
94
  
95
}