1.     #include <stdio.h>
2.     #include <stdlib.h>
3.     #include <string.h>
4.     #include <math.h>
5.     #include <omp.h>
6.     
7.     // fill p[] with primes up to max plus one more
8.     int
9.     primearray(int max, int p[]) {
10.        p[0] = 3;
11.        int top = 0, n, i;
12.        for (n = 5;; n += 2)
13.            for (i = 0;; i++)
14.                if (n % p[i] == 0)
15.                    break;
16.                else
17.                    if (p[i]*p[i] > n) {
18.                        p[++top] = n;
19.        printf("prime: %d top: %d %d %ld  \n",
20.                 n, top, p[top], (long)p[top]*p[top]);
21.                        if (n > max)
22.                            return top; // last n
23.                        break;
24.                    }
25.    }
26.    // use primes in p[] to check up to max
27.    // 5 is p[1] is the 3rd prime, so start count from 2
28.    int
29.    primecount_omp(long max, int p[]) {
30.        int n, i, sum = 2;
31.        #pragma omp parallel for private(i) schedule(guided) 
      reduction(+:sum) 
32.        for (n = 5; n <= max; n += 2)
33.            for (i = 0;; i++)
34.                if (n % p[i] == 0)
35.                    break;
36.                else
37.                    if (p[i]*p[i] > n) {
38.                        sum++;
39.                        break;
40.                    }
41.        return sum;
42.    }
43.    // call primearray() with SQR and then primecount_omp() with MAX
44.    // PNT (x/logx) to find approx. size of array p[] psize
45.    int
46.    main(int argc, char** argv) {
47.        if (argc < 2) {
48.            printf("You need to supply a limit\n");
49.            return 1;
50.        }
51.        long MAX = strtoul(argv[1], 0, 10);
52.        int SQR = sqrt(MAX);
53.        int psize = 50 + 1.2 * SQR / log(SQR);
54.        int *p = malloc(psize * sizeof*p);
55.        int top = primearray(SQR, p);
56.        printf("psize: %d top: %d %d %ld\n", psize, top, p[top], 
      (long)p[top]*p[top]); 
57.        int sum = primecount_omp(MAX, p);
58.        printf("Count: %d (%.2f%%)\n", sum, (double) 100*sum / MAX) ;
59.        double lnx = MAX / log(MAX);
60.        printf("LnX:   %.2f (%f)\n",  lnx, sum / lnx);
61.        return 0;
62.    }