Download code

From LiteratePrograms

Jump to: navigation, search

Back to Miller-Rabin_primality_test_(C)

Download for Windows: zip

Download for UNIX: zip, tar.gz, tar.bz2

integer.c

  1 /* Copyright (c) 2010 the authors listed at the following URL, and/or
  2 the authors of referenced articles or incorporated external code:
  3 http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?action=history&offset=20060620214840
  4 
  5 Permission is hereby granted, free of charge, to any person obtaining
  6 a copy of this software and associated documentation files (the
  7 "Software"), to deal in the Software without restriction, including
  8 without limitation the rights to use, copy, modify, merge, publish,
  9 distribute, sublicense, and/or sell copies of the Software, and to
 10 permit persons to whom the Software is furnished to do so, subject to
 11 the following conditions:
 12 
 13 The above copyright notice and this permission notice shall be
 14 included in all copies or substantial portions of the Software.
 15 
 16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 19 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 20 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 21 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 22 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 23 
 24 Retrieved from: http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?oldid=6349
 25 */
 26 
 27 #include <stdlib.h>   /* for malloc() */

 28 
 29 #include <string.h>   /* for memset(), memmove() */

 30 
 31 #include <math.h>  /* for ceil */

 32 
 33 #include "integer.h"

 34 
 35 integer create_integer(int components) {
 36     integer result;
 37     result.num_components = components;
 38     result.c = (component_t*)malloc(sizeof(component_t)*components);
 39     return result;
 40 }
 41 
 42 
 43 void free_integer(integer i) {
 44     free(i.c);
 45 }
 46 
 47 
 48 void set_zero_integer(integer i) {
 49     memset(i.c, 0, sizeof(component_t)*i.num_components);
 50 }
 51 
 52 
 53 int is_zero_integer(integer x) {
 54     int i;
 55     for(i=0; i < x.num_components; i++) {
 56         if (x.c[i] != 0) return 0;
 57     }
 58     return 1;
 59 }
 60 
 61 
 62 void copy_integer(integer source, integer target) {
 63     memmove(target.c, source.c,
 64             sizeof(component_t)*MIN(source.num_components, target.num_components));
 65     if (target.num_components > source.num_components) {
 66         memset(target.c + source.num_components, 0,
 67                sizeof(component_t)*(target.num_components - source.num_components));
 68     }
 69 }
 70 
 71 
 72 void add_integer(integer left, integer right, integer result) {
 73     double_component_t carry = 0;
 74     int i;
 75     for(i=0; i<left.num_components || i<right.num_components || carry != 0; i++) {
 76         double_component_t partial_sum = carry;
 77         carry = 0;
 78         if (i < left.num_components)  partial_sum += left.c[i];
 79         if (i < right.num_components) partial_sum += right.c[i];
 80         if (partial_sum > MAX_COMPONENT) {
 81             partial_sum &= MAX_COMPONENT;
 82             carry = 1;
 83         }
 84         result.c[i] = (component_t)partial_sum;
 85     }
 86     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
 87 }
 88 
 89 
 90 void subtract_integer(integer left, integer right, integer result) {
 91     int borrow = 0;
 92     int i;
 93     for(i=0; i<left.num_components; i++) {
 94         double_component_t lhs = left.c[i];
 95         double_component_t rhs = (i < right.num_components) ? right.c[i] : 0;
 96         if (borrow) {
 97 	    if (lhs <= rhs) {
 98 	        /* leave borrow set to 1 */
 99 	        lhs += (MAX_COMPONENT + 1) - 1;
100 	    } else {
101 	        borrow = 0;
102 	        lhs--;
103 	    }
104 	}
105 
106         if (lhs < rhs) {
107             borrow = 1;
108             lhs += MAX_COMPONENT + 1;
109         }
110         result.c[i] = lhs - rhs;
111     }
112     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
113 }
114 
115 
116 void multiply_small_integer(integer left, component_t right, integer result) {
117     double_component_t carry = 0;
118     int i;
119     for(i=0; i<left.num_components || carry != 0; i++) {
120         double_component_t partial_sum = carry;
121         carry = 0;
122         if (i < left.num_components)  partial_sum += left.c[i]*right;
123         carry = partial_sum >> COMPONENT_BITS;
124         result.c[i] = (component_t)(partial_sum & MAX_COMPONENT);
125     }
126     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
127 }
128 
129 
130 void multiply_integer(integer left, integer right, integer result) {
131     int i, lidx, ridx;
132     double_component_t carry = 0;
133     int max_size_no_carry;
134     int left_max_component  = left.num_components - 1;
135     int right_max_component = right.num_components - 1;
136     while(left.c[left_max_component]   == 0) left_max_component--;
137     while(right.c[right_max_component] == 0) right_max_component--;
138 
139     max_size_no_carry = left_max_component + right_max_component;
140     for(i=0; i <= max_size_no_carry || carry != 0; i++) {
141         double_component_t partial_sum = carry;
142         carry = 0;
143         lidx = MIN(i, left_max_component);
144         ridx = i - lidx;
145         while(lidx >= 0 && ridx <= right_max_component) {
146             partial_sum += ((double_component_t)left.c[lidx])*right.c[ridx];
147             carry += partial_sum >> COMPONENT_BITS;
148             partial_sum &= MAX_COMPONENT;
149             lidx--; ridx++;
150         }
151         result.c[i] = partial_sum;
152     }
153     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
154 }
155 
156 
157 int compare_integers(integer left, integer right) {
158     int i = MAX(left.num_components - 1, right.num_components - 1);
159     for ( ; i >= 0; i--) {
160         component_t left_comp =
161             (i < left.num_components) ? left.c[i] : 0;
162         component_t right_comp =
163             (i < right.num_components) ? right.c[i] : 0;
164         if (left_comp < right_comp)
165             return -1;
166         else if (left_comp > right_comp)
167             return  1;
168     }
169     return 0;
170 }
171 
172 
173 void shift_left_one_integer(integer arg) {
174     int i;
175     arg.c[arg.num_components - 1] <<= 1;
176     for (i = arg.num_components - 2; i >= 0; i--) {
177         arg.c[i + 1] |= arg.c[i] >> (COMPONENT_BITS - 1);
178         arg.c[i] <<= 1;
179     }
180 }
181 
182 
183 void shift_right_one_integer(integer arg) {
184     int i;
185     arg.c[0] >>= 1;
186     for (i = 1; i < arg.num_components; i++) {
187         arg.c[i - 1] |= (arg.c[i] & 1) << (COMPONENT_BITS - 1);
188         arg.c[i] >>= 1;
189     }
190 }
191 
192 
193 component_t mod_small_integer(integer left, component_t right) {
194     double_component_t mod_two_power = 1;
195     double_component_t result = 0;
196     int i, bit;
197     for(i=0; i<left.num_components; i++) {
198         for(bit=0; bit<COMPONENT_BITS; bit++) {
199             if ((left.c[i] & (1 << bit)) != 0) {
200                 result += mod_two_power;
201                 if (result >= right) {
202                     result -= right;
203                 }
204             }
205             mod_two_power <<= 1;
206             if (mod_two_power >= right) {
207                 mod_two_power -= right;
208             }
209         }
210     }
211     return (component_t)result;
212 }
213 
214 
215 void mod_integer(integer left, integer right, integer result) {
216     integer mod_two_power = create_integer(right.num_components + 1);
217     int i, bit;
218     set_zero_integer(result);
219     set_zero_integer(mod_two_power);
220     mod_two_power.c[0] = 1;
221     for(i=0; i<left.num_components; i++) {
222         for(bit=0; bit<COMPONENT_BITS; bit++) {
223             if ((left.c[i] & (1 << bit)) != 0) {
224                 add_integer(result, mod_two_power, result);
225                 if (compare_integers(result, right) >= 0) {
226                     subtract_integer(result, right, result);
227                 }
228             }
229             shift_left_one_integer(mod_two_power);
230             if (compare_integers(mod_two_power, right) >= 0) {
231                 subtract_integer(mod_two_power, right, mod_two_power);
232             }
233         }
234     }
235     free_integer(mod_two_power);
236 }
237 
238 
239 void divide_small_integer(integer left, component_t right, integer result) {
240     double_component_t dividend = 0;
241     int i;
242     for (i = left.num_components - 1; i >= 0; i--) {
243         dividend |= left.c[i];
244         result.c[i] = dividend/right;
245         dividend = (dividend % right) << COMPONENT_BITS;
246     }
247 }
248 
249 
250 integer string_to_integer(char* s) {
251     integer result = create_integer((int)ceil(LOG_2_10*strlen(s)/COMPONENT_BITS));
252     integer digit = create_integer(1);
253     int i;
254     for (i = 0; s[i] != '\0'; i++) {
255         multiply_small_integer(result, 10, result);
256         digit.c[0] = s[i] - '0';
257         add_integer(result, digit, result);
258     }
259     free_integer(digit);
260     return result;
261 }
262 
263 
264 char* integer_to_string(integer x) {
265     int i, result_len;
266     char* result =
267         (char*)malloc((int)ceil(COMPONENT_BITS*x.num_components/LOG_2_10) + 2);
268     integer ten = create_integer(1);
269     ten.c[0] = 10;
270 
271     if (is_zero_integer(x)) {
272         strcpy(result, "0");
273     } else {
274         for (i = 0; !is_zero_integer(x); i++) {
275             result[i] = (char)mod_small_integer(x, 10) + '0';
276             divide_small_integer(x, 10, x);
277         }
278         result[i] = '\0';
279     }
280 
281     result_len = strlen(result);
282     for(i=0; i < result_len/2; i++) {
283         char temp = result[i];
284         result[i] = result[result_len - i - 1];
285         result[result_len - i - 1] = temp;
286     }
287 
288 
289     free_integer(ten);
290     return result;
291 }
292 
293 
294 
295


integer.h

 1 /* Copyright (c) 2010 the authors listed at the following URL, and/or
 2 the authors of referenced articles or incorporated external code:
 3 http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?action=history&offset=20060620214840
 4 
 5 Permission is hereby granted, free of charge, to any person obtaining
 6 a copy of this software and associated documentation files (the
 7 "Software"), to deal in the Software without restriction, including
 8 without limitation the rights to use, copy, modify, merge, publish,
 9 distribute, sublicense, and/or sell copies of the Software, and to
10 permit persons to whom the Software is furnished to do so, subject to
11 the following conditions:
12 
13 The above copyright notice and this permission notice shall be
14 included in all copies or substantial portions of the Software.
15 
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
19 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
20 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
21 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
22 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 
24 Retrieved from: http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?oldid=6349
25 */
26 
27 #ifndef _INTEGER_H_

28 #define _INTEGER_H_

29 
30 #include <limits.h>  /* for CHAR_BIT */

31 
32 typedef unsigned short component_t;
33 typedef unsigned long double_component_t;
34 
35 #define MAX_COMPONENT  ((component_t)(-1))

36 
37 #define COMPONENT_BITS  (sizeof(component_t)*CHAR_BIT)

38 
39 #define LOG_2_10    3.3219280948873623478703194294894

40 
41 
42 #define MIN(x,y)  ((x)<(y) ? (x) : (y))

43 #define MAX(x,y)  ((x)>(y) ? (x) : (y))

44 
45 
46 typedef struct {
47     component_t* c;    /* least-significant word first */
48     int num_components;
49 } integer;
50 
51 
52 integer create_integer(int components);
53 void free_integer(integer i);
54 void set_zero_integer(integer i);
55 void copy_integer(integer source, integer target);
56 void add_integer(integer left, integer right, integer result);
57 void subtract_integer(integer left, integer right, integer result);
58 void multiply_small_integer(integer left, component_t right, integer result);
59 void multiply_integer(integer left, integer right, integer result);
60 int compare_integers(integer left, integer right);
61 void shift_left_one_integer(integer arg);
62 void shift_right_one_integer(integer arg);
63 component_t mod_small_integer(integer left, component_t right);
64 void mod_integer(integer left, integer right, integer result);
65 void divide_small_integer(integer left, component_t right, integer result);
66 integer string_to_integer(char* s);
67 char* integer_to_string(integer x);
68 int is_zero_integer(integer x);
69 
70 #endif /* #ifndef _INTEGER_H_ */

71 
72


miller-rabin.c

  1 /* Copyright (c) 2010 the authors listed at the following URL, and/or
  2 the authors of referenced articles or incorporated external code:
  3 http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?action=history&offset=20060620214840
  4 
  5 Permission is hereby granted, free of charge, to any person obtaining
  6 a copy of this software and associated documentation files (the
  7 "Software"), to deal in the Software without restriction, including
  8 without limitation the rights to use, copy, modify, merge, publish,
  9 distribute, sublicense, and/or sell copies of the Software, and to
 10 permit persons to whom the Software is furnished to do so, subject to
 11 the following conditions:
 12 
 13 The above copyright notice and this permission notice shall be
 14 included in all copies or substantial portions of the Software.
 15 
 16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 19 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 20 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 21 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 22 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 23 
 24 Retrieved from: http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?oldid=6349
 25 */
 26 
 27 #include <stdio.h>

 28 #include <stdlib.h>

 29 #include <string.h>

 30 #include <time.h>

 31 #include "integer.h"

 32 
 33 #define COMPOSITE 0

 34 #define PRIME     1

 35 
 36 
 37 integer modular_exponent(integer base, integer power, integer modulus) {
 38     int i, bit;
 39     integer result = create_integer(modulus.num_components + 1);
 40     integer temp   = create_integer(modulus.num_components*2 + 1);
 41     set_zero_integer(result);
 42     result.c[0] = 1;
 43 
 44     for(i=power.num_components - 1; i>=0; i--) {
 45         for(bit=COMPONENT_BITS-1; bit>=0; bit--) {
 46             multiply_integer(result, result, temp);
 47             mod_integer(temp, modulus, result);
 48             if ((power.c[i] & (1 << bit)) != 0) {
 49                 multiply_integer(result, base, temp);
 50                 mod_integer(temp, modulus, result);
 51             }
 52         }
 53     }
 54 
 55     free_integer(temp);
 56     return result;
 57 }
 58 
 59 
 60 int miller_rabin_pass(integer a, integer n) {
 61     int i, s, result;
 62     integer a_to_power, d, one, n_minus_one;
 63     integer temp = create_integer(n.num_components*2 + 1);
 64 
 65     one = create_integer(1);
 66     set_zero_integer(one);
 67     one.c[0] = 1;
 68     n_minus_one = create_integer(n.num_components);
 69     subtract_integer(n, one, n_minus_one);
 70 
 71     s = 0;
 72     d = create_integer(n_minus_one.num_components);
 73     copy_integer(n_minus_one, d);
 74     while ((d.c[0] % 2) == 0) {
 75         shift_right_one_integer(d);
 76         s++;
 77     }
 78 
 79     a_to_power = modular_exponent(a, d, n);
 80     if (compare_integers(a_to_power, one) == 0)  { result=PRIME; goto exit; }
 81     for(i=0; i < s-1; i++) {
 82         if (compare_integers(a_to_power, n_minus_one) == 0) { result=PRIME; goto exit; }
 83         multiply_integer(a_to_power, a_to_power, temp);
 84         mod_integer(temp, n, a_to_power);
 85     }
 86     if (compare_integers(a_to_power, n_minus_one) == 0) { result=PRIME; goto exit; }
 87     result = COMPOSITE;
 88 
 89 exit:
 90     free_integer(temp);
 91     free_integer(a_to_power);
 92     free_integer(one);
 93     free_integer(n_minus_one);
 94     return result;
 95 }
 96 
 97 
 98 void random_integer(integer max, integer result) {
 99     int i, most_sig = max.num_components - 1;
100     while (max.c[most_sig] == 0) most_sig--;
101     for (i=0; i<most_sig; i++) {
102         result.c[i] = rand() % (MAX_COMPONENT + 1);
103     }
104     result.c[most_sig] = rand() % max.c[most_sig];
105 }
106 
107 
108 int miller_rabin(integer n) {
109     integer a = create_integer(n.num_components);
110     int repeat;
111     for(repeat=0; repeat<20; repeat++) {
112         do {
113             random_integer(n, a);
114         } while (is_zero_integer(a));
115         if (miller_rabin_pass(a, n) == COMPOSITE) {
116             return COMPOSITE;
117         }
118     }
119     return PRIME;
120 }
121 
122 
123 int main(int argc, char* argv[]) {
124     srand(time(NULL));
125     if (strcmp(argv[1], "test") == 0) {
126         integer n = string_to_integer(argv[2]);
127         puts(miller_rabin(n) == PRIME ? "PRIME" : "COMPOSITE");
128     } else if (strcmp(argv[1], "genprime") == 0) {
129         integer max = create_integer(atoi(argv[2])/COMPONENT_BITS);
130         integer p   = create_integer(max.num_components);
131         set_zero_integer(max);
132         max.c[max.num_components-1] = MAX_COMPONENT;
133         do {
134             random_integer(max, p);
135             if ((p.c[0] % 2) == 0) continue;
136 	    if (mod_small_integer(p, 3) == 0) continue;
137 	    if (mod_small_integer(p, 5) == 0) continue;
138 	    if (mod_small_integer(p, 7) == 0) continue;
139 
140         } while (miller_rabin(p) == COMPOSITE);
141         puts(integer_to_string(p));
142     }
143     return 0;
144 }
145


miller-rabin_16.c

 1 /* Copyright (c) 2010 the authors listed at the following URL, and/or
 2 the authors of referenced articles or incorporated external code:
 3 http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?action=history&offset=20060620214840
 4 
 5 Permission is hereby granted, free of charge, to any person obtaining
 6 a copy of this software and associated documentation files (the
 7 "Software"), to deal in the Software without restriction, including
 8 without limitation the rights to use, copy, modify, merge, publish,
 9 distribute, sublicense, and/or sell copies of the Software, and to
10 permit persons to whom the Software is furnished to do so, subject to
11 the following conditions:
12 
13 The above copyright notice and this permission notice shall be
14 included in all copies or substantial portions of the Software.
15 
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
19 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
20 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
21 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
22 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 
24 Retrieved from: http://en.literateprograms.org/Miller-Rabin_primality_test_(C)?oldid=6349
25 */
26 
27 #include <stdio.h>

28 #include <stdlib.h>

29 
30 typedef unsigned short u16;
31 typedef unsigned long  u32;
32 
33 #define COMPOSITE 0

34 #define PRIME     1

35 
36 
37 u16 modular_exponent_16(u16 base, u16 power, u16 modulus) {
38     u32 result = 1;
39     int i;
40     for (i=15; i>=0; i--) {
41         result = (result*result) % modulus;
42         if (power & (1 << i)) {
43             result = (result*base) % modulus;
44         }
45     }
46     return (u16)result; /* Will not truncate since modulus is a u16 */
47 }
48 
49 
50 int miller_rabin_pass_16(u16 a, u16 n) {
51     u16 a_to_power, s, d, i;
52     s = 0;
53     d = n - 1;
54     while ((d % 2) == 0) {
55         d /= 2;
56         s++;
57     }
58 
59     a_to_power = modular_exponent_16(a, d, n);
60     if (a_to_power == 1) return PRIME;
61     for(i=0; i < s-1; i++) {
62         if (a_to_power == n-1) return PRIME;
63         a_to_power = modular_exponent_16(a_to_power, 2, n);
64     }
65     if (a_to_power == n-1) return PRIME;
66     return COMPOSITE;
67 }
68 
69 
70 int miller_rabin_16(u16 n) {
71     if (n <= 1) return COMPOSITE;
72     if (n == 2) return PRIME;
73     if (miller_rabin_pass_16( 2, n) == PRIME &&
74         (n <= 7  || miller_rabin_pass_16( 7, n) == PRIME) &&
75         (n <= 61 || miller_rabin_pass_16(61, n) == PRIME))
76         return PRIME;
77     else
78         return COMPOSITE;
79 }
80 
81 
82 int main(int argc, char* argv[]) {
83     u16 n = (u16)atoi(argv[1]);
84     puts(miller_rabin_16(n) == PRIME ? "PRIME" : "COMPOSITE");
85     return 0;
86 }
87


Retrieved from "http://en.literateprograms.org/Special:Downloadcode"
Views
Personal tools