From: Mark Wooding Date: Sat, 17 Nov 2018 19:21:43 +0000 (+0000) Subject: math/: Implement Grantham's Frobenius (primality) test. X-Git-Tag: 2.5.0~20 X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/commitdiff_plain/6001a9ffafa1e77b2d192938d79e6da80febdc43 math/: Implement Grantham's Frobenius (primality) test. This is a rather heavyweight test which is effective when checking possibly adversarial numbers. There are no known composites which pass both this test and the Miller--Rabin test with witness 2 (although infinitely many are conjectured to exist); the combination is called the `Baillie--PSW' test (after Baillie, Pomerance, Selfridge, and Wagstaff). Modify `pgen_primep' to use Baillie--PSW. Since Baillie--PSW is somewhat faster than the many rounds of Miller-- Rabin which `pgen_primep' used to use, celebrate by raising the `keen' threshold in the `dh-param.c' test. This work was prompted by the paper `Prime and Prejudice', by Martin R. Albrecht, Jake Massimo, Kenneth G. Paterson, and Juraj Somorovsky; though, since Catacomb already used 32 iterations of Miller--Rabin with random witnesses, I can confidently state that the previous implementation was inefficient but secure when used with a good randomness source. --- diff --git a/math/Makefile.am b/math/Makefile.am index ac654d91..dd48057a 100644 --- a/math/Makefile.am +++ b/math/Makefile.am @@ -245,9 +245,10 @@ libmath_la_SOURCES += pfilt.c pkginclude_HEADERS += pgen.h libmath_la_SOURCES += pgen.c libmath_la_SOURCES += pgen-gcd.c +libmath_la_SOURCES += pgen-granfrob.c libmath_la_SOURCES += pgen-simul.c libmath_la_SOURCES += pgen-stdev.c -TESTS += pgen.t$(EXEEXT) +TESTS += pgen.t$(EXEEXT) pgen-granfrob.t$(EXEEXT) EXTRA_DIST += t/pgen ## Finding primitive elements in finite fields. diff --git a/math/pgen-granfrob.c b/math/pgen-granfrob.c new file mode 100644 index 00000000..8e49c731 --- /dev/null +++ b/math/pgen-granfrob.c @@ -0,0 +1,276 @@ +/* -*-c-*- + * + * Grantham's Frobenius primality test + * + * (c) 2018 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software: you can redistribute it and/or modify it + * under the terms of the GNU Library General Public License as published + * by the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb. If not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, + * USA. + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" +#include "mpmont.h" +#include "mpscan.h" +#include "pgen.h" + +#include "mptext.h" + +/*----- Main code ---------------------------------------------------------*/ + +static int squarep(mp *n) +{ + mp *t = MP_NEW; + int rc; + + if (MP_NEGP(n)) return (0); + t = mp_sqrt(t, n); t = mp_sqr(t, t); + rc = MP_EQ(t, n); mp_drop(t); return (rc); +} + +/* --- @pgen_granfrob@ --- * + * + * Arguments: @mp *n@ = an integer to test + * @int a, b@ = coefficients; if @a@ is zero then choose + * automatically + * + * Returns: One of the @PGEN_...@ codes. + * + * Use: Performs a quadratic versoin of Grantham's Frobenius + * primality test, which is a simple extension of the standard + * Lucas test. + * + * If %$a^2 - 4 b$% is a perfect square then the test can't + * work; this function returns @PGEN_ABORT@ under these + * circumstances. + */ + +int pgen_granfrob(mp *n, int a, int b) +{ + mp *v = MP_NEW, *w = MP_NEW, *aa = MP_NEW, *bb = MP_NEW, *bi = MP_NEW, + *k = MP_NEW, *x = MP_NEW, *y = MP_NEW, *z = MP_NEW, *t, *u; + mp ma; mpw wa; + mp mb; mpw wb; + mp md; mpw wd; int d; + mpmont mm; + mpscan msc; + int e, bit, rc; + + /* Maybe this is a no-hoper. */ + if (MP_NEGP(n)) return (PGEN_FAIL); + if (MP_EQ(n, MP_TWO)) return (PGEN_DONE); + if (!MP_ODDP(n)) return (PGEN_FAIL); + + /* First, build the parameters as large integers. */ + mp_build(&ma, &wa, &wa + 1); mp_build(&mb, &wb, &wb + 1); + mp_build(&md, &wd, &wd + 1); + mpmont_create(&mm, n); + + /* Prepare the Lucas sequence parameters. Here, %$\Delta$% is the + * disciminant of the polynomial %$p(x) = x^2 - a x + b$%, i.e., + * %$\Delta = a^2 - 4 b$%. + */ + if (a) { + /* Explicit parameters. Just set them and check that they'll work. */ + + if (a >= 0) wa = a; else { wa = -a; ma.f |= MP_NEG; } + if (b >= 0) wb = b; else { wb = -b; mb.f |= MP_NEG; } + d = a*a - 4*b; + if (d >= 0) wd = d; else { wd = -d; md.f |= MP_NEG; } + + /* Determine the quadratic character of %$\Delta$%. If %$(\Delta | n)$% + * is zero then we'll have a problem, but we'll catch that case with the + * GCD check below. + */ + e = mp_jacobi(&md, n); + + /* If %$\Delta$% is a perfect square then the test can't work. */ + if (e == 1 && squarep(&md)) { rc = PGEN_ABORT; goto end; } + } else { + /* Determine parameters. Use Selfridge's `Method A': choose the first + * %$\Delta$% from the sequence %$5$%, %$-7$%, %%\dots%%, such that + * %$(\Delta | n) = -1$%. + */ + + wa = 1; wd = 5; + for (;;) { + e = mp_jacobi(&md, n); if (e != +1) break; + if (wd == 25 && squarep(n)) { rc = PGEN_FAIL; goto end; } + wd += 2; md.f ^= MP_NEG; + } + a = 1; d = wd; + if (md.f&MP_NEG) { wb = (wd + 1)/4; d = -d; } + else { wb = (wd - 1)/4; mb.f |= MP_NEG; } + b = (1 - d)/4; + } + + /* The test won't work if %$\gcd(2 a b \Delta, n) \ne 1$%. */ + x = mp_lsl(x, &ma, 1); x = mp_mul(x, x, &mb); x = mp_mul(x, x, &md); + mp_gcd(&y, 0, 0, x, n); + if (!MP_EQ(y, MP_ONE)) + { rc = MP_EQ(y, n) ? PGEN_ABORT : PGEN_FAIL; goto end; } + + /* Now we use binary a Lucas chain to evaluate %$V_{n-e}(a, b) \pmod{n}$%. + * Here, + * + * * %$U_{i+1}(a, b) = a U_i(a, b) - b U_{i-1}(a, b)$%, and + * * %$V_{i+1}(a, b) = a V_i(a, b) - b V_{i-1}(a, b)$%; with + * * %$U_0(a, b) = 0$%, $%U_1(a, b) = 1$%, %$V_0(a, b) = 2$%, and + * %$V_1(a, b) = a$%. + * + * To compute this, we use the handy identities + * + * %$V_{i+j}(a, b) = V_i(a, b) V_j(a, b) - b^i V_{j-i}(a, b)$% + * + * and + * + * %$U_i(a, b) = (2 V_{i+1}(a, b) - a V_i(a, b))/\Delta$%. + * + * Let %$k = n - e$%. Given %$V_i(a, b)$% and %$V_{i+1}(a, b)$%, we can + * determine either %$V_{2i}(a, b)$% and %$V_{2i+1}(a, b)$%, or + * %$V_{2i+1}(a, b)$% and %$V_{2i+2}(a, b)$%. + * + * To do this, suppose that %$n < 2^\ell$% and %$0 \le i \le \ell%; we'll + * start with %$i = 0$%. Divide %$n = n_i 2^{\ell-i} + n'_i$% with + * %$n'_i < 2^{\ell-i}$%. To do this, we maintain %$v_i = V_{n_i}(a, b)$%, + * %$w_i = V_{n_i+1}(a, b)$%, and %$b_i = b^{n_i}$%, all modulo %$n$%. If + * %$n'_i < 2^{\ell-i-1}$% then we have %$n'_{i+1} = n'_i$% and + * %$n_{i+i} = 2 n_i$%; otherwise %$n'_{i+1} = n'_i - 2^{\ell-i-1}$% and + * %$n_{i+i} = 2 n_i + 1$%. + */ + k = mp_add(k, n, e > 0 ? MP_MONE : MP_ONE); + aa = mpmont_mul(&mm, aa, &ma, mm.r2); + bb = mpmont_mul(&mm, bb, &mb, mm.r2); bi = MP_COPY(mm.r); + v = mpmont_mul(&mm, v, MP_TWO, mm.r2); w = MP_COPY(aa); + + for (mpscan_rinitx(&msc, k->v, k->vl); mpscan_rstep(&msc); ) { + bit = mpscan_rbit(&msc); + + /* We will want %$x = V_{n_i+1}(a, b) = V_{n_i} V_{n_i+1} - a b^{n_i}$%, + * but we don't yet know whether this is %$v_{i+1}$% or %$w_{i+1}$%. + */ + x = mpmont_mul(&mm, x, v, w); + if (a == 1) x = mp_sub(x, x, bi); + else { y = mpmont_mul(&mm, y, aa, bi); x = mp_sub(x, x, y); } + if (MP_NEGP(x)) x = mp_add(x, x, n); + + if (!bit) { + /* We're in the former case: %$n_{i+i} = 2 n_i$%. So %$w_{i+1} = x$%, + * %$v_{i+1} = (v_i^2 - 2 b_i$%, and %$b_{i+1} = b_i^2$%. + */ + + y = mp_sqr(y, v); y = mpmont_reduce(&mm, y, y); + y = mp_sub(y, y, bi); if (MP_NEGP(y)) y = mp_add(y, y, n); + y = mp_sub(y, y, bi); if (MP_NEGP(y)) y = mp_add(y, y, n); + bi = mp_sqr(bi, bi); bi = mpmont_reduce(&mm, bi, bi); + t = v; u = w; v = y; w = x; x = t; y = u; + } else { + /* We're in the former case: %$n_{i+i} = 2 n_i + 1$%. So + * %$v_{i+1} = x$%, %$w_{i+1} = w_i^2 - 2 b b^i$%$%, and + * %$b_{i+1} = b b_i^2$%. + */ + + y = mp_sqr(y, w); y = mpmont_reduce(&mm, y, y); + z = mpmont_mul(&mm, z, bi, bb); + y = mp_sub(y, y, z); if (MP_NEGP(y)) y = mp_add(y, y, n); + y = mp_sub(y, y, z); if (MP_NEGP(y)) y = mp_add(y, y, n); + bi = mpmont_mul(&mm, bi, bi, z); + t = v; u = w; v = x; w = y; x = t; y = u; + } + } + + /* The Lucas test is that %$U_k(a, b) \equiv 0 \pmod{n}$% if %$n$% is + * prime. I'm assured that + * + * %$U_k(a, b) = (2 V_{k+1}(a, b) - a V_k(a, b))/\Delta$% + * + * so this is just a matter of checking that %$2 w - a v \equiv 0$%. + */ + x = mp_add(x, w, w); y = mp_sub(y, x, n); + if (!MP_NEGP(y)) { t = x; x = y; y = t; } + if (a == 1) x = mp_sub(x, x, v); + else { y = mpmont_mul(&mm, y, v, aa); x = mp_sub(x, x, y); } + if (MP_NEGP(x)) x = mp_add(x, x, n); + if (!MP_ZEROP(x)) { rc = PGEN_FAIL; goto end; } + + /* Grantham's Frobenius test is that, also, %$V_k(a, b) v = \equiv 2 b$% + * if %$n$% is prime and %$(\Delta | n) = -1$%, or %$v \equiv 2$% if + * %$(\Delta | n) = +1$%. + */ + if (MP_ODDP(v)) v = mp_add(v, v, n); + v = mp_lsr(v, v, 1); + if (!MP_EQ(v, e == +1 ? mm.r : bb)) { rc = PGEN_FAIL; goto end; } + + /* Looks like we made it. */ + rc = PGEN_PASS; +end: + mp_drop(v); mp_drop(w); mp_drop(aa); mp_drop(bb); mp_drop(bi); + mp_drop(k); mp_drop(x); mp_drop(y); mp_drop(z); + mpmont_destroy(&mm); + return (rc); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +#include "mptext.h" + +static int verify(dstr *v) +{ + mp *n = *(mp **)v[0].buf; + int a = *(int *)v[1].buf, b = *(int *)v[2].buf, xrc = *(int *)v[3].buf, rc; + int ok = 1; + + rc = pgen_granfrob(n, a, b); + if (rc != xrc) { + fputs("\n*** pgen_granfrob failed", stdout); + fputs("\nn = ", stdout); mp_writefile(n, stdout, 10); + printf("\na = %d", a); + printf("\nb = %d", a); + printf("\nexp rc = %d", xrc); + printf("\ncalc rc = %d\n", rc); + ok = 0; + } + + mp_drop(n); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static test_chunk tests[] = { + { "pgen-granfrob", verify, + { &type_mp, &type_int, &type_int, &type_int, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/t/pgen"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/math/pgen.c b/math/pgen.c index 9a822f57..f10d585f 100644 --- a/math/pgen.c +++ b/math/pgen.c @@ -319,28 +319,33 @@ mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx, * @grand *gr@ = a random number source * * Returns: Nonzero if @p@ is really prime. + * + * Use: Checks the primality of @p@. If @p@ is prime, then this + * function returns nonzero; if @p@ is really composite then it + * %%\emph{probably}%% returns zero, but might not. + * + * Currently, this function uses the Baillie--PSW test, which + * combines a single Miller--Rabin test with witness 2 with a + * single Frobenius test with parameters chosen using + * Selfridge's `Method A'. No composites are known which pass + * this test, though it's conjectured that infinitely many + * exist. */ int pgen_primep(mp *p, grand *gr) { - int i; rabin r; - mp *x = MP_NEW; + int rc; if (MP_NEGP(p)) return (0); switch (pfilt_smallfactor(p)) { case PGEN_DONE: return (1); case PGEN_FAIL: return (0); } - rabin_create(&r, p); - for (i = 32; i; i--) { - x = mprand_range(x, p, gr, 0); - if (rabin_rtest(&r, x) == PGEN_FAIL) - break; - } - MP_DROP(x); - rabin_destroy(&r); - return (!i); + rabin_create(&r, p); rc = rabin_test(&r, MP_TWO); rabin_destroy(&r); + if (rc == PGEN_FAIL) return (0); + rc = pgen_granfrob(p, 0, 0); if (rc == PGEN_FAIL) return (0); + return (1); } /*----- Test rig ----------------------------------------------------------*/ diff --git a/math/pgen.h b/math/pgen.h index b103c30e..219b261b 100644 --- a/math/pgen.h +++ b/math/pgen.h @@ -270,6 +270,29 @@ extern mp *pgen(const char */*name*/, mp */*d*/, mp */*m*/, unsigned /*steps*/, pgen_proc */*step*/, void */*sctx*/, unsigned /*tests*/, pgen_proc */*test*/, void */*tctx*/); +/* --- @pgen_granfrob@ --- * + * + * Arguments: @mp *n@ = an integer to test + * @int a, b@ = coefficients; if @a@ is zero then choose + * automatically + * + * Returns: One of the @PGEN_...@ codes. + * + * Use: Performs a quadratic versoin of Grantham's Frobenius + * primality test, which is a simple extension of the standard + * Lucas test. + * + * If %$a^2 - 4 b$% is a perfect square then the test can't + * work; this function returns @PGEN_ABORT@ under these + * circumstances. + * + * If @a@ is zero on entry, then the function will choose + * suitable parameters deterministically -- i.e., it always + * chooses the same parameters for a given %$n$%. + */ + +extern int pgen_granfrob(mp */*n*/, int /*a*/, int /*b*/); + /* --- @pgen_primep@ --- * * * Arguments: @mp *p@ = a number to check diff --git a/math/t/pgen b/math/t/pgen index 8d52d197..be26c183 100644 --- a/math/t/pgen +++ b/math/t/pgen @@ -14,6 +14,29 @@ pgen { 166359567317705838255275971708060308423814413741683015010175247351623188739655446196925981468626681882384215574706593049022467680136399439302347043107836749816290369600677730213469006507173065402294688841278559283358390567733443050775707749725690534182003442070447739085348456478911335969765393755383551520173 166359567317705838255275971708060308423814413741683015010175247351623188739655446196925981468626681882384215574706593049022467680136399439302347043107836749816290369600677730213469006507173065402294688841278559283358390567733443050775707749725690534182003442070447739085348456478911335969765393755383551520257; } +pgen-granfrob { + 5 0 0 -1; + 7 0 0 4; + 15 0 0 3; + 5777 1 -1 4; # pseudoprime + 40301809 0 0 4; + 86059163416987297647409667483582114939806237974424324409828198660056356336227 1 5 4; + 102508420970861015999300753620309481186457893679971500520427161277511389396803 1 5 4; + 72291866454056552194087337607224612505157525245486245416393486917859196707519 1 5 4; + 72291866454056552194087337607224612505157525245486265416393486917859196707519 1 5 3; + + ## A large Frobenius pseudoprime: call the first number p_1; then p_2 = 31 + ## (p_1 + 1) - 1 and p_3 = 43 (p_1 + 1) - 1. These three are all prime. + ## Their product is a strong Lucas, and Frobenius, pseudoprime. + ## + ## See `Prime and Prejudice' by Martin R. Albrecht, Jake Massimo, Kenneth + ## G. Paterson, and Juraj Somorovsky. + 3690125385954346893658786222051913500627130245213169388019826598097107079718295481926241398412699320815932808015860263240282855670239765686869973444864115322609857375876438922226372746215468824202413623127 0 0 4; + 114393886964584753703422372883609318519441037601608251028614624541010319471267159939713483350793678945293917048491668160448768525777432736292969176790787575000905578652169606589017555132679533550274822316967 0 0 4; + 158675391596036916427327807548232280526966600544166283684852543718175604427886705722828380131746070795085110744681991319332162793820309924535408858129156958872223867162686873655734028087265159440703785794503 0 0 4; + 66981291792500223036804182765508448534715465524671325885174850970812009004775815201151227900130153990294748113034471984909912807896550069799856170439734910206802409847773026240559371480115711600866989845251707737806461503879250232804362190067578216069266197879151809743235261582813331022213587929425243163096486125825510076936556242805690400001899138503900919499414951069309064408305196756524628693684938044145785145327821174180933033293089394794328963673467918652042794300291355500468079109432376296868174257674548727592142782202898031102246775544402811199608266683925072825828225074019194302318324623049819212337927 0 0 4; +} + primep { -5 0; -1 0; @@ -24,6 +47,30 @@ primep { 4 0; 40301809 1; 40301811 0; + + ## A small Lucas pseudoprime: 5777 = 53*109. + 5777 0; + + ## A large strong pseudoprime: this is the product of + ## + ## p_1 = 142445387161415482404826365418175962266689133006163 + ## p_2 = 5840260873618034778597880982145214452934254453252643 + ## p_3 = 14386984103302963722887462907235772188935602433622363 + ## + ## See `Prime and Prejudice' by Martin R. Albrecht, Jake Massimo, Kenneth + ## G. Paterson, and Juraj Somorovsky. + 142445387161415482404826365418175962266689133006163 1; + 5840260873618034778597880982145214452934254453252643 1; + 14386984103302963722887462907235772188935602433622363 1; + 11968794224604718293549908104759518204343930652759288592987578098131927050572705181539873293848476235393230314654912729920657864630317971562727057595285667 0; + + ## A large Lucas pseudoprime: call the first number p_1; then p_2 = 31 (p_1 + ## + 1) - 1 and p_3 = 43 (p_1 + 1) - 1. These three are all prime. Their + ## product is a strong Lucas pseudoprime. + 3690125385954346893658786222051913500627130245213169388019826598097107079718295481926241398412699320815932808015860263240282855670239765686869973444864115322609857375876438922226372746215468824202413623127 1; + 114393886964584753703422372883609318519441037601608251028614624541010319471267159939713483350793678945293917048491668160448768525777432736292969176790787575000905578652169606589017555132679533550274822316967 1; + 158675391596036916427327807548232280526966600544166283684852543718175604427886705722828380131746070795085110744681991319332162793820309924535408858129156958872223867162686873655734028087265159440703785794503 1; + 66981291792500223036804182765508448534715465524671325885174850970812009004775815201151227900130153990294748113034471984909912807896550069799856170439734910206802409847773026240559371480115711600866989845251707737806461503879250232804362190067578216069266197879151809743235261582813331022213587929425243163096486125825510076936556242805690400001899138503900919499414951069309064408305196756524628693684938044145785145327821174180933033293089394794328963673467918652042794300291355500468079109432376296868174257674548727592142782202898031102246775544402811199608266683925072825828225074019194302318324623049819212337927 0; } primeiter { diff --git a/pub/dh-param.c b/pub/dh-param.c index b593fb8c..69e24376 100644 --- a/pub/dh-param.c +++ b/pub/dh-param.c @@ -126,7 +126,7 @@ int main(int argc, char *argv[]) group *g; dh_infofromdata(&dp, pe->data); g = group_prime(&dp); - if (mp_bits(dp.p) > 2048 && + if (mp_bits(dp.p) > 3072 && (!argv[1] || strcmp(argv[1], "keen") != 0)) { printf(" [%s skipped]", pe->name); fflush(stdout);