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