Commit | Line | Data |
---|---|---|
6001a9ff MW |
1 | /* -*-c-*- |
2 | * | |
3 | * Grantham's Frobenius primality test | |
4 | * | |
5 | * (c) 2018 Straylight/Edgeware | |
6 | */ | |
7 | ||
8 | /*----- Licensing notice --------------------------------------------------* | |
9 | * | |
10 | * This file is part of Catacomb. | |
11 | * | |
12 | * Catacomb is free software: you can redistribute it and/or modify it | |
13 | * under the terms of the GNU Library General Public License as published | |
14 | * by the Free Software Foundation; either version 2 of the License, or | |
15 | * (at your option) any later version. | |
16 | * | |
17 | * Catacomb is distributed in the hope that it will be useful, but | |
18 | * WITHOUT ANY WARRANTY; without even the implied warranty of | |
19 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
20 | * Library General Public License for more details. | |
21 | * | |
22 | * You should have received a copy of the GNU Library General Public | |
23 | * License along with Catacomb. If not, write to the Free Software | |
24 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, | |
25 | * USA. | |
26 | */ | |
27 | ||
28 | /*----- Header files ------------------------------------------------------*/ | |
29 | ||
30 | #include "mp.h" | |
31 | #include "mpmont.h" | |
32 | #include "mpscan.h" | |
33 | #include "pgen.h" | |
34 | ||
35 | #include "mptext.h" | |
36 | ||
37 | /*----- Main code ---------------------------------------------------------*/ | |
38 | ||
6001a9ff MW |
39 | /* --- @pgen_granfrob@ --- * |
40 | * | |
41 | * Arguments: @mp *n@ = an integer to test | |
42 | * @int a, b@ = coefficients; if @a@ is zero then choose | |
43 | * automatically | |
44 | * | |
45 | * Returns: One of the @PGEN_...@ codes. | |
46 | * | |
9bfacb4f | 47 | * Use: Performs a quadratic version of Grantham's Frobenius |
6001a9ff MW |
48 | * primality test, which is a simple extension of the standard |
49 | * Lucas test. | |
50 | * | |
51 | * If %$a^2 - 4 b$% is a perfect square then the test can't | |
52 | * work; this function returns @PGEN_ABORT@ under these | |
53 | * circumstances. | |
54 | */ | |
55 | ||
56 | int pgen_granfrob(mp *n, int a, int b) | |
57 | { | |
58 | mp *v = MP_NEW, *w = MP_NEW, *aa = MP_NEW, *bb = MP_NEW, *bi = MP_NEW, | |
59 | *k = MP_NEW, *x = MP_NEW, *y = MP_NEW, *z = MP_NEW, *t, *u; | |
60 | mp ma; mpw wa; | |
61 | mp mb; mpw wb; | |
62 | mp md; mpw wd; int d; | |
63 | mpmont mm; | |
64 | mpscan msc; | |
65 | int e, bit, rc; | |
66 | ||
67 | /* Maybe this is a no-hoper. */ | |
68 | if (MP_NEGP(n)) return (PGEN_FAIL); | |
69 | if (MP_EQ(n, MP_TWO)) return (PGEN_DONE); | |
70 | if (!MP_ODDP(n)) return (PGEN_FAIL); | |
71 | ||
72 | /* First, build the parameters as large integers. */ | |
73 | mp_build(&ma, &wa, &wa + 1); mp_build(&mb, &wb, &wb + 1); | |
74 | mp_build(&md, &wd, &wd + 1); | |
75 | mpmont_create(&mm, n); | |
76 | ||
77 | /* Prepare the Lucas sequence parameters. Here, %$\Delta$% is the | |
78 | * disciminant of the polynomial %$p(x) = x^2 - a x + b$%, i.e., | |
79 | * %$\Delta = a^2 - 4 b$%. | |
80 | */ | |
81 | if (a) { | |
82 | /* Explicit parameters. Just set them and check that they'll work. */ | |
83 | ||
84 | if (a >= 0) wa = a; else { wa = -a; ma.f |= MP_NEG; } | |
85 | if (b >= 0) wb = b; else { wb = -b; mb.f |= MP_NEG; } | |
86 | d = a*a - 4*b; | |
87 | if (d >= 0) wd = d; else { wd = -d; md.f |= MP_NEG; } | |
88 | ||
89 | /* Determine the quadratic character of %$\Delta$%. If %$(\Delta | n)$% | |
90 | * is zero then we'll have a problem, but we'll catch that case with the | |
91 | * GCD check below. | |
92 | */ | |
93 | e = mp_jacobi(&md, n); | |
94 | ||
95 | /* If %$\Delta$% is a perfect square then the test can't work. */ | |
0e1199ac | 96 | if (e == 1 && mp_squarep(&md)) { rc = PGEN_ABORT; goto end; } |
6001a9ff MW |
97 | } else { |
98 | /* Determine parameters. Use Selfridge's `Method A': choose the first | |
99 | * %$\Delta$% from the sequence %$5$%, %$-7$%, %%\dots%%, such that | |
100 | * %$(\Delta | n) = -1$%. | |
101 | */ | |
102 | ||
103 | wa = 1; wd = 5; | |
104 | for (;;) { | |
105 | e = mp_jacobi(&md, n); if (e != +1) break; | |
0e1199ac | 106 | if (wd == 25 && mp_squarep(n)) { rc = PGEN_FAIL; goto end; } |
6001a9ff MW |
107 | wd += 2; md.f ^= MP_NEG; |
108 | } | |
109 | a = 1; d = wd; | |
110 | if (md.f&MP_NEG) { wb = (wd + 1)/4; d = -d; } | |
111 | else { wb = (wd - 1)/4; mb.f |= MP_NEG; } | |
112 | b = (1 - d)/4; | |
113 | } | |
114 | ||
115 | /* The test won't work if %$\gcd(2 a b \Delta, n) \ne 1$%. */ | |
116 | x = mp_lsl(x, &ma, 1); x = mp_mul(x, x, &mb); x = mp_mul(x, x, &md); | |
117 | mp_gcd(&y, 0, 0, x, n); | |
118 | if (!MP_EQ(y, MP_ONE)) | |
119 | { rc = MP_EQ(y, n) ? PGEN_ABORT : PGEN_FAIL; goto end; } | |
120 | ||
121 | /* Now we use binary a Lucas chain to evaluate %$V_{n-e}(a, b) \pmod{n}$%. | |
122 | * Here, | |
123 | * | |
124 | * * %$U_{i+1}(a, b) = a U_i(a, b) - b U_{i-1}(a, b)$%, and | |
125 | * * %$V_{i+1}(a, b) = a V_i(a, b) - b V_{i-1}(a, b)$%; with | |
126 | * * %$U_0(a, b) = 0$%, $%U_1(a, b) = 1$%, %$V_0(a, b) = 2$%, and | |
127 | * %$V_1(a, b) = a$%. | |
128 | * | |
129 | * To compute this, we use the handy identities | |
130 | * | |
131 | * %$V_{i+j}(a, b) = V_i(a, b) V_j(a, b) - b^i V_{j-i}(a, b)$% | |
132 | * | |
133 | * and | |
134 | * | |
135 | * %$U_i(a, b) = (2 V_{i+1}(a, b) - a V_i(a, b))/\Delta$%. | |
136 | * | |
137 | * Let %$k = n - e$%. Given %$V_i(a, b)$% and %$V_{i+1}(a, b)$%, we can | |
138 | * determine either %$V_{2i}(a, b)$% and %$V_{2i+1}(a, b)$%, or | |
139 | * %$V_{2i+1}(a, b)$% and %$V_{2i+2}(a, b)$%. | |
140 | * | |
141 | * To do this, suppose that %$n < 2^\ell$% and %$0 \le i \le \ell%; we'll | |
142 | * start with %$i = 0$%. Divide %$n = n_i 2^{\ell-i} + n'_i$% with | |
143 | * %$n'_i < 2^{\ell-i}$%. To do this, we maintain %$v_i = V_{n_i}(a, b)$%, | |
144 | * %$w_i = V_{n_i+1}(a, b)$%, and %$b_i = b^{n_i}$%, all modulo %$n$%. If | |
145 | * %$n'_i < 2^{\ell-i-1}$% then we have %$n'_{i+1} = n'_i$% and | |
146 | * %$n_{i+i} = 2 n_i$%; otherwise %$n'_{i+1} = n'_i - 2^{\ell-i-1}$% and | |
147 | * %$n_{i+i} = 2 n_i + 1$%. | |
148 | */ | |
149 | k = mp_add(k, n, e > 0 ? MP_MONE : MP_ONE); | |
150 | aa = mpmont_mul(&mm, aa, &ma, mm.r2); | |
151 | bb = mpmont_mul(&mm, bb, &mb, mm.r2); bi = MP_COPY(mm.r); | |
152 | v = mpmont_mul(&mm, v, MP_TWO, mm.r2); w = MP_COPY(aa); | |
153 | ||
154 | for (mpscan_rinitx(&msc, k->v, k->vl); mpscan_rstep(&msc); ) { | |
155 | bit = mpscan_rbit(&msc); | |
156 | ||
157 | /* We will want %$x = V_{n_i+1}(a, b) = V_{n_i} V_{n_i+1} - a b^{n_i}$%, | |
158 | * but we don't yet know whether this is %$v_{i+1}$% or %$w_{i+1}$%. | |
159 | */ | |
160 | x = mpmont_mul(&mm, x, v, w); | |
161 | if (a == 1) x = mp_sub(x, x, bi); | |
162 | else { y = mpmont_mul(&mm, y, aa, bi); x = mp_sub(x, x, y); } | |
163 | if (MP_NEGP(x)) x = mp_add(x, x, n); | |
164 | ||
165 | if (!bit) { | |
166 | /* We're in the former case: %$n_{i+i} = 2 n_i$%. So %$w_{i+1} = x$%, | |
167 | * %$v_{i+1} = (v_i^2 - 2 b_i$%, and %$b_{i+1} = b_i^2$%. | |
168 | */ | |
169 | ||
170 | y = mp_sqr(y, v); y = mpmont_reduce(&mm, y, y); | |
171 | y = mp_sub(y, y, bi); if (MP_NEGP(y)) y = mp_add(y, y, n); | |
172 | y = mp_sub(y, y, bi); if (MP_NEGP(y)) y = mp_add(y, y, n); | |
173 | bi = mp_sqr(bi, bi); bi = mpmont_reduce(&mm, bi, bi); | |
174 | t = v; u = w; v = y; w = x; x = t; y = u; | |
175 | } else { | |
176 | /* We're in the former case: %$n_{i+i} = 2 n_i + 1$%. So | |
177 | * %$v_{i+1} = x$%, %$w_{i+1} = w_i^2 - 2 b b^i$%$%, and | |
178 | * %$b_{i+1} = b b_i^2$%. | |
179 | */ | |
180 | ||
181 | y = mp_sqr(y, w); y = mpmont_reduce(&mm, y, y); | |
182 | z = mpmont_mul(&mm, z, bi, bb); | |
183 | y = mp_sub(y, y, z); if (MP_NEGP(y)) y = mp_add(y, y, n); | |
184 | y = mp_sub(y, y, z); if (MP_NEGP(y)) y = mp_add(y, y, n); | |
185 | bi = mpmont_mul(&mm, bi, bi, z); | |
186 | t = v; u = w; v = x; w = y; x = t; y = u; | |
187 | } | |
188 | } | |
189 | ||
190 | /* The Lucas test is that %$U_k(a, b) \equiv 0 \pmod{n}$% if %$n$% is | |
191 | * prime. I'm assured that | |
192 | * | |
193 | * %$U_k(a, b) = (2 V_{k+1}(a, b) - a V_k(a, b))/\Delta$% | |
194 | * | |
195 | * so this is just a matter of checking that %$2 w - a v \equiv 0$%. | |
196 | */ | |
197 | x = mp_add(x, w, w); y = mp_sub(y, x, n); | |
198 | if (!MP_NEGP(y)) { t = x; x = y; y = t; } | |
199 | if (a == 1) x = mp_sub(x, x, v); | |
200 | else { y = mpmont_mul(&mm, y, v, aa); x = mp_sub(x, x, y); } | |
201 | if (MP_NEGP(x)) x = mp_add(x, x, n); | |
202 | if (!MP_ZEROP(x)) { rc = PGEN_FAIL; goto end; } | |
203 | ||
204 | /* Grantham's Frobenius test is that, also, %$V_k(a, b) v = \equiv 2 b$% | |
205 | * if %$n$% is prime and %$(\Delta | n) = -1$%, or %$v \equiv 2$% if | |
206 | * %$(\Delta | n) = +1$%. | |
207 | */ | |
208 | if (MP_ODDP(v)) v = mp_add(v, v, n); | |
209 | v = mp_lsr(v, v, 1); | |
210 | if (!MP_EQ(v, e == +1 ? mm.r : bb)) { rc = PGEN_FAIL; goto end; } | |
211 | ||
212 | /* Looks like we made it. */ | |
213 | rc = PGEN_PASS; | |
214 | end: | |
215 | mp_drop(v); mp_drop(w); mp_drop(aa); mp_drop(bb); mp_drop(bi); | |
216 | mp_drop(k); mp_drop(x); mp_drop(y); mp_drop(z); | |
217 | mpmont_destroy(&mm); | |
218 | return (rc); | |
219 | } | |
220 | ||
221 | /*----- Test rig ----------------------------------------------------------*/ | |
222 | ||
223 | #ifdef TEST_RIG | |
224 | ||
225 | #include <mLib/testrig.h> | |
226 | ||
227 | #include "mptext.h" | |
228 | ||
229 | static int verify(dstr *v) | |
230 | { | |
231 | mp *n = *(mp **)v[0].buf; | |
232 | int a = *(int *)v[1].buf, b = *(int *)v[2].buf, xrc = *(int *)v[3].buf, rc; | |
233 | int ok = 1; | |
234 | ||
235 | rc = pgen_granfrob(n, a, b); | |
236 | if (rc != xrc) { | |
237 | fputs("\n*** pgen_granfrob failed", stdout); | |
238 | fputs("\nn = ", stdout); mp_writefile(n, stdout, 10); | |
239 | printf("\na = %d", a); | |
240 | printf("\nb = %d", a); | |
241 | printf("\nexp rc = %d", xrc); | |
242 | printf("\ncalc rc = %d\n", rc); | |
243 | ok = 0; | |
244 | } | |
245 | ||
246 | mp_drop(n); | |
247 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
248 | return (ok); | |
249 | } | |
250 | ||
251 | static test_chunk tests[] = { | |
252 | { "pgen-granfrob", verify, | |
253 | { &type_mp, &type_int, &type_int, &type_int, 0 } }, | |
254 | { 0, 0, { 0 } } | |
255 | }; | |
256 | ||
257 | int main(int argc, char *argv[]) | |
258 | { | |
259 | sub_init(); | |
260 | test_run(argc, argv, tests, SRCDIR "/t/pgen"); | |
261 | return (0); | |
262 | } | |
263 | ||
264 | #endif | |
265 | ||
266 | /*----- That's all, folks -------------------------------------------------*/ |