3 * Grantham's Frobenius primality test
5 * (c) 2018 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
37 /*----- Main code ---------------------------------------------------------*/
39 static int squarep(mp
*n
)
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
);
49 /* --- @pgen_granfrob@ --- *
51 * Arguments: @mp *n@ = an integer to test
52 * @int a, b@ = coefficients; if @a@ is zero then choose
55 * Returns: One of the @PGEN_...@ codes.
57 * Use: Performs a quadratic versoin of Grantham's Frobenius
58 * primality test, which is a simple extension of the standard
61 * If %$a^2 - 4 b$% is a perfect square then the test can't
62 * work; this function returns @PGEN_ABORT@ under these
66 int pgen_granfrob(mp
*n
, int a
, int b
)
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
;
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
);
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
);
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$%.
92 /* Explicit parameters. Just set them and check that they'll work. */
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
; }
97 if (d
>= 0) wd
= d
; else { wd
= -d
; md
.f
|= MP_NEG
; }
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
103 e
= mp_jacobi(&md
, n
);
105 /* If %$\Delta$% is a perfect square then the test can't work. */
106 if (e
== 1 && squarep(&md
)) { rc
= PGEN_ABORT
; goto end
; }
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$%.
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
;
120 if (md
.f
&MP_NEG
) { wb
= (wd
+ 1)/4; d
= -d
; }
121 else { wb
= (wd
- 1)/4; mb
.f
|= MP_NEG
; }
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
; }
131 /* Now we use binary a Lucas chain to evaluate %$V_{n-e}(a, b) \pmod{n}$%.
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
139 * To compute this, we use the handy identities
141 * %$V_{i+j}(a, b) = V_i(a, b) V_j(a, b) - b^i V_{j-i}(a, b)$%
145 * %$U_i(a, b) = (2 V_{i+1}(a, b) - a V_i(a, b))/\Delta$%.
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)$%.
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$%.
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
);
164 for (mpscan_rinitx(&msc
, k
->v
, k
->vl
); mpscan_rstep(&msc
); ) {
165 bit
= mpscan_rbit(&msc
);
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}$%.
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
);
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$%.
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
;
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$%.
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
;
200 /* The Lucas test is that %$U_k(a, b) \equiv 0 \pmod{n}$% if %$n$% is
201 * prime. I'm assured that
203 * %$U_k(a, b) = (2 V_{k+1}(a, b) - a V_k(a, b))/\Delta$%
205 * so this is just a matter of checking that %$2 w - a v \equiv 0$%.
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
; }
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$%.
218 if (MP_ODDP(v
)) v
= mp_add(v
, v
, n
);
220 if (!MP_EQ(v
, e
== +1 ? mm
.r
: bb
)) { rc
= PGEN_FAIL
; goto end
; }
222 /* Looks like we made it. */
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
);
231 /*----- Test rig ----------------------------------------------------------*/
235 #include <mLib/testrig.h>
239 static int verify(dstr
*v
)
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
;
245 rc
= pgen_granfrob(n
, a
, b
);
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
);
257 assert(mparena_count(MPARENA_GLOBAL
) == 0);
261 static test_chunk tests
[] = {
262 { "pgen-granfrob", verify
,
263 { &type_mp
, &type_int
, &type_int
, &type_int
, 0 } },
267 int main(int argc
, char *argv
[])
270 test_run(argc
, argv
, tests
, SRCDIR
"/t/pgen");
276 /*----- That's all, folks -------------------------------------------------*/