/* N integer, believed to be prime, is g a generator for F_N^* ? * -- g an integer (t_INT) * -- L = vector of the (N-1)/p, where N-1 = prod p^{e_p}. Passed as argument * to avoid factoring N-1 at each test */ testgen(g, N, L) = { g = Mod(g, N); \\ g^k % N is very inefficient if (!g, return (0)); if (g^(N-1) != 1, error(N, " n'est pas premier")); for (i=1, #L, if (g^L[i] == 1, return (0)); ); return (1); } \\ Return a generator for (Z/NZ)*. Assume that N is prime gen(N) = { my (P = factor(N-1)[,1]); \\ N-1 = # (Z/NZ)^* L = vector(#P, i, (N-1)/P[i]); while (1, \\ oo loop g = 1 + random(N-1); \\ in [1,N-1] if (testgen(g,N,L), return ([g, P])); ); } /************* "p-1" Primality Proof (naïve) ************/ TRIVIAL = "TRIVIAL"; \\ global constant \\ Return a primality certificate for N certificate(N)= { my(G, g, P); if (N == 2, return ([2, TRIVIAL])); \\ trivial certificate for 2 G = gen(N); g = G[1]; P = G[2];; [N, g, P, vector(#P, i, certificate(P[i]))]; } \\ Check a primality certificate C for integer N. \\ If N is omitted, check whether C is consistent check(C, N = 0) = { my(g = C[2], P, V); if (!N, N = C[1]); \\ N omitted, check consistency if (N != C[1], return (0)); \\ not the right N ! if (g == TRIVIAL, return (1)); P = C[3]; V = C[4]; for (i=1, #V, if (!check(V[i], P[i]), return (0))); \\ We now know that the P[i] are primes if (N-1 != prod(i=1, #P, P[i]^valuation(N-1, P[i])), return (0)); \\ We now know that the P[i] are the prime divisors of N-1 \\ N.B: Fails if the P[i] are not distinct ! if (!testgen(g, N, vector(#P,i,(N-1)/P[i])), return (0)); \\ We now know that g has order N-1 in (Z/NZ)^* return (1); \\ Hence N is prime }