1 /* Exercise mpz_probab_prime_p.
2
3 Copyright 2002, 2018-2019, 2022 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include "gmp-impl.h"
23 #include "tests.h"
24
25
26 /* Enhancements:
27
28 - Test some big primes don't come back claimed to be composite.
29 - Test some big composites don't come back claimed to be certainly prime.
30 - Test some big composites with small factors are identified as certainly
31 composite. */
32
33
34 /* return 2 if prime, 0 if composite */
35 int
36 isprime (unsigned long n)
37 {
38 if (n < 4)
39 return (n & 2);
40 if ((n & 1) == 0)
41 return 0;
42
43 for (unsigned long i = 3; i*i <= n; i+=2)
44 if ((n % i) == 0)
45 return 0;
46
47 return 2;
48 }
49
50 void
51 check_one (mpz_srcptr n, int want)
52 {
53 int got;
54
55 got = mpz_probab_prime_p (n, 25);
56
57 /* "definitely prime" (2) is fine if we only wanted "probably prime" (1) */
58 if ((got != want) && (got != want * 2))
59 {
60 printf ("mpz_probab_prime_p\n");
61 mpz_trace (" n ", n);
62 printf (" got =%d", got);
63 printf (" want=%d", want);
64 abort ();
65 }
66 }
67
68 void
69 check_pn (mpz_ptr n, int want)
70 {
71 check_one (n, want);
72 mpz_neg (n, n);
73 check_one (n, want);
74 }
75
76 /* expect certainty for small n */
77 void
78 check_small (void)
79 {
80 mpz_t n;
81 long i;
82
83 mpz_init (n);
84
85 for (i = 0; i < 300; i++)
86 {
87 mpz_set_si (n, i);
88 check_pn (n, isprime (i));
89 }
90
91 mpz_clear (n);
92 }
93
94 void
95 check_composites (int count)
96 {
97 int i;
98 mpz_t a, b, n, bs;
99 unsigned long size_range, size;
100 gmp_randstate_ptr rands = RANDS;
101
102 mpz_init (a);
103 mpz_init (b);
104 mpz_init (n);
105 mpz_init (bs);
106
107 static const char * const composites[] = {
108 "225670644213750121", /* n=61*C16, if D < 61, (n/D) = 1. */
109 "2386342059899637841", /* n=61*C17, if D < 61, (n/D) = 1. */
110 "1194649", /* A square, but strong base-2 pseudoprime, */
111 "12327121", /* another base-2 pseudoprime square. */
112 "18446744066047760377", /* Should trigger Fibonacci's test; */
113 "10323769", /* &3==1, Lucas' test with D=37; */
114 "1397419", /* &3==3, Lucas' test with D=43; */
115 "11708069165918597341", /* &3==1, Lucas' test with large D=107; */
116 "395009109077493751", /* &3==3, Lucas' test with large D=113. */
117 NULL
118 };
119
120 for (i = 0; composites[i]; i++)
121 {
122 mpz_set_str_or_abort (n, composites[i], 0);
123 check_one (n, 0);
124 }
125
126 for (i = 0; i < count; i++)
127 {
128 mpz_urandomb (bs, rands, 32);
129 size_range = mpz_get_ui (bs) % 13 + 1; /* 0..8192 bit operands */
130
131 mpz_urandomb (bs, rands, size_range);
132 size = mpz_get_ui (bs);
133 mpz_rrandomb (a, rands, size);
134
135 mpz_urandomb (bs, rands, 32);
136 size_range = mpz_get_ui (bs) % 13 + 1; /* 0..8192 bit operands */
137 mpz_rrandomb (b, rands, size);
138
139 /* Exclude trivial factors */
140 if (mpz_cmp_ui (a, 1) == 0)
141 mpz_set_ui (a, 2);
142 if (mpz_cmp_ui (b, 1) == 0)
143 mpz_set_ui (b, 2);
144
145 mpz_mul (n, a, b);
146
147 check_pn (n, 0);
148 }
149 mpz_clear (a);
150 mpz_clear (b);
151 mpz_clear (n);
152 mpz_clear (bs);
153 }
154
155 static void
156 check_primes (void)
157 {
158 static const char * const primes[] = {
159 "2", "53", "1234567891",
160 "2055693949", "1125899906842597", "16412292043871650369",
161 "18446744075358702679", /* Lucas' test with large D=107. */
162 /* diffie-hellman-group1-sha1, also "Well known group 2" in RFC
163 2412, 2^1024 - 2^960 - 1 + 2^64 * { [2^894 pi] + 129093 } */
164 "0xFFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD1"
165 "29024E088A67CC74020BBEA63B139B22514A08798E3404DD"
166 "EF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245"
167 "E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7ED"
168 "EE386BFB5A899FA5AE9F24117C4B1FE649286651ECE65381"
169 "FFFFFFFFFFFFFFFF",
170 NULL
171 };
172
173 mpz_t n;
174 int i;
175
176 mpz_init (n);
177
178 for (i = 0; primes[i]; i++)
179 {
180 mpz_set_str_or_abort (n, primes[i], 0);
181 check_one (n, 1);
182 }
183 mpz_clear (n);
184 }
185
186 static void
187 check_fermat_mersenne (int count)
188 {
189 int fermat_exponents [] = {1, 2, 4, 8, 16};
190 int mersenne_exponents [] = {2, 3, 5, 7, 13, 17, 19, 31, 61, 89,
191 107, 127, 521, 607, 1279, 2203, 2281,
192 3217, 4253, 4423, 9689, 9941, 11213,
193 19937, 21701, 23209, 44497, 86243};
194 mpz_t pp;
195 int i, j, want;
196
197 mpz_init (pp);
198 count = MIN (110000, count);
199
200 for (i=1; i<count; ++i)
201 {
202 mpz_set_ui (pp, 1);
203 mpz_setbit (pp, i); /* 2^i + 1 */
204 want = 0;
205 for (j = 0; j < numberof (fermat_exponents); j++)
206 if (fermat_exponents[j] == i)
207 {
208 /* Fermat's primes are small enough for a definite answer. */
209 want = 2;
210 break;
211 }
212 check_one (pp, want);
213
214 mpz_sub_ui (pp, pp, 2); /* 2^i - 1 */
215 want = 0;
216 for (j = 0; j < numberof (mersenne_exponents); j++)
217 if (mersenne_exponents[j] == i)
218 {
219 want = 1 << (i < 50);
220 break;
221 }
222 check_one (pp, want);
223 }
224 mpz_clear (pp);
225 }
226
227 int
228 main (int argc, char **argv)
229 {
230 int count = 1000;
231
232 TESTS_REPS (count, argv, argc);
233
234 tests_start ();
235
236 check_small ();
237 check_fermat_mersenne (count >> 3);
238 check_composites (count);
239 check_primes ();
240
241 tests_end ();
242 exit (0);
243 }