chiark - git - mdw - catacomb/blob - math/pgen-granfrob.c

   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 -------------------------------------------------*/