prime.c 4.67 KB
Newer Older
atuncer @kelenken's avatar
atuncer @kelenken committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

int test_n_against_a(mpz_t n, mpz_t n_1, mpz_t d, mp_bitcnt_t s, mpz_t a, mpz_t x);


int main(int argc, char *argv[]) {
    int world_size;
    int world_rank;
    mpz_t n;
    mpz_t n_1;
    mp_bitcnt_t s;
    mpz_t d;
    mpz_t a;
    mpz_t a_max;
    mpz_t x;
    size_t logn;
    int answer_local; // 0: inconclusive, 1: prime, 2: composite
    int answer_global; // 0: inconclusive, 1: prime, 2: composite

    MPI_Request answer_found_req;
    int answer_found_flag;
    int is_composite;


    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);
    MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

    // print an error if the starting position is not specified
    if (argc != 2) {
        printf("USAGE: %s <start>\n", argv[0]);
        printf("  start: an integer as a strating point (arbitrary precision)\n", argv[0]);
        MPI_Finalize();
        exit(EXIT_FAILURE);
    }


    // print an error if working with a single process
    if (world_size < 2) {
        printf("ERROR: please create at least 2 processes!\n");
        MPI_Finalize();
        exit(EXIT_FAILURE);
    }


    mpz_init_set_str(n, argv[1], 10);

    mpz_init(n_1);
    mpz_init(d);
    mpz_init(a);
    mpz_init(a_max);
    mpz_init(x);


    while(1) {
        // PREPARE VARIABLES

        // n_1 = n - 1
        mpz_sub_ui(n_1, n, 1);

        // s = [ x WHERE n-1 = 2^x * d ]
        s = mpz_scan1(n_1, 0);

        // d = n_1 / 2^s
        mpz_tdiv_q_2exp(d, n_1, s);

        // logn = log(n)
        logn = mpz_sizeinbase(n, 2);

        // a_max = ( (logn^2) / [number of workers] )
        mpz_set_ui(a_max, logn);
        mpz_mul_ui(a_max, a_max, logn);
        // TODO: range starts from 2, we scan a bit extra
        mpz_cdiv_q_ui(a_max, a_max, world_size);

        // a = [start of scanning range for this worker]
        // a_max = [end of scanning range for this worker]
        mpz_mul_ui(a, a_max, world_rank);
        mpz_mul_ui(a_max, a_max, world_rank + 1);
        // range should start from 2
        mpz_add_ui(a, a, 3);
        mpz_add_ui(a_max, a_max, 3);

        answer_local = 0;

#ifdef DEBUG
        if (world_rank == 0) {
            gmp_printf("n     : %Zd\n", n);
            gmp_printf("n-1   : %Zd\n", n_1);
            gmp_printf("s     : %lu\n", s);
            gmp_printf("d     : %Zd\n", d);
            gmp_printf("logn  : %d\n", logn);
        }
        gmp_printf("p%02d   : %Zd - %Zd\n", world_rank, a, a_max);
#endif

        // scan the assigned range
        do {

            // can it be proven that n is composite?
            if (test_n_against_a(n, n_1, d, s, a, x)) {
                // n is composite
                answer_local = 2;
            } 

            // is {the local answer still inconclusive} and {test range is exhausted}?
            if ( (answer_local == 0) && ((mpz_cmp(a, a_max) >= 0))) {
                // n might be a prime
                answer_local = 1;
            }

#ifdef DEBUG
            printf("composition: %d found %d\n", world_rank, answer_local);
#endif

            // reduce local answers. more definite answers override others
            // {2: definitely composite} > {1: might be prime} > {0: inconclusive}
            MPI_Allreduce(&answer_local, &answer_global, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);

            // prepare to check against the next a in range
            mpz_add_ui(a, a, 1);

            // answer_global == {0 :inconclusive} : inconclusive, restart to check against the next a in range
            // answer_global == {1: might be prime} : inconclusive, but no more iterations available, exit
            // answer_global == {2: definitely composite} : conclusive, exit
        } while (answer_global == 0); 

        // print out the agreed answer (only on the first worker)
        if (world_rank == 0) {
#ifdef DEBUG
            gmp_printf(" %d : %Zd\n", answer_global, n);
#else
            if (answer_global == 1) {
                gmp_printf("%Zd\n", n);
            }
#endif
        }

        // increment n
        mpz_add_ui(n, n, 1);

    }

    MPI_Finalize();
}




int test_n_against_a(mpz_t n, mpz_t n_1, mpz_t d, mp_bitcnt_t s, mpz_t a, mpz_t x) {
    // 0: inconclusive
    // 1: number is composite

    mp_bitcnt_t r;

    // x = a^d (mod n)
    mpz_powm(x, a, d, n);

    // if (x == 1) return inconclusive
    if ( ! mpz_cmp_ui(x, 1) ) return 0;

    // if (x == n-1) return inconclusive
    if ( ! mpz_cmp(x, n_1) ) return 0;

    for (r=1; r<s; r++) {

        // x = x^2 (mod n)
        mpz_powm_ui(x, x, 2, n);

        // if (x == n-1) return inconclusive
        if ( ! mpz_cmp(x, n_1) ) return 0;

    }

    // return is_composite
    return 1;

}