Primality proving program based on Pocklington's theorem. GNU MP library and GNU C Compiler are required.
Command line:
echo data data ... 0 | pock
data := n f f ...
n := number
f := prime factor of n-1
Example:
echo 23 2 11 0 0 | pockAugust 17, 2003 ... Version 0.2.1 is available.
pock021.c
//Primality proving program based on Pocklington's theorem
// powered by GMP 4.1.2
// version 0.2.1 by M.Kamada
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "gmp.h"
#define TITLE "Primality proving program based on Pocklington's theorem"
#define VERSION "0.2.1"
#define AUTHOR "M.Kamada"
#define USAGE "usage: echo data data ... 0 | %s\n\
data: n f f ...\n\
n: number\n\
f: prime factor of n-1\n\
example: echo 23 2 11 0 0 | %s\n"
#define MAX_M 10000000
#define PROBAB_P_PARAM 5
#define FLUSH() fflush(stdout)
//primes up to 1000
const unsigned int primes[] = {
2,3,5,7,11,13,17,19,23,29,
31,37,41,43,47,53,59,61,67,71,
73,79,83,89,97,101,103,107,109,113,
127,131,137,139,149,151,157,163,167,173,
179,181,191,193,197,199,211,223,227,229,
233,239,241,251,257,263,269,271,277,281,
283,293,307,311,313,317,331,337,347,349,
353,359,367,373,379,383,389,397,401,409,
419,421,431,433,439,443,449,457,461,463,
467,479,487,491,499,503,509,521,523,541,
547,557,563,569,571,577,587,593,599,601,
607,613,617,619,631,641,643,647,653,659,
661,673,677,683,691,701,709,719,727,733,
739,743,751,757,761,769,773,787,797,809,
811,821,823,827,829,839,853,857,859,863,
877,881,883,887,907,911,919,929,937,941,
947,953,967,971,977,983,991,997,
};
#define FACTORLIST_LENGTH 200
//prime factor list
typedef struct {
int length;
mpz_t f[FACTORLIST_LENGTH];
int e[FACTORLIST_LENGTH];
} factorlist_t;
//main proof
// mpz_t n : number
// factorlist_t f : prime factor list of n-1
// int result : -1=unknown, 0=definitely composite, 2=definitely prime
int pocklington(mpz_t n, factorlist_t *f) {
int result = -1;
unsigned int *a, *a_fermat_tested;
const unsigned int *a_limit =
(unsigned int *)primes + sizeof(primes) / sizeof(unsigned int);
mpz_t n1; //n-1
mpz_t t, u;
int i; //factor #
mpz_init(n1);
mpz_init(t);
mpz_init(u);
mpz_sub_ui(n1, n, 1); //n-1
a = a_fermat_tested = (unsigned int *)primes;
i = 0;
for (;;) {
//check gcd(a,n)==1
if (mpz_cmp_ui(n, *a) == 0) { //n==a
result = 2; //definitely prime
break;
}
if (mpz_divisible_ui_p(n, *a)) { //n is divisible by a
printf("n is divisible by %u\n", *a);
FLUSH();
result = 0; //definitely composite
break;
}
//check a^(n-1)==1 (mod n)
if (a >= a_fermat_tested) {
mpz_set_ui(t, *a); //t=a
mpz_powm(t, t, n1, n); //t=a^n1 mod n
gmp_printf("%u^(n-1)=%Zd (mod n)\n", *a, t);
FLUSH();
if (mpz_cmp_ui(t, 1) != 0) { //a^(n-1)!=1 (mod n)
result = 0; //definitely composite
break;
}
a_fermat_tested++;
}
//check gcd(a^((n-1)/f[i])-1,n)==1
mpz_set_ui(t, *a); //t=a
mpz_div(u, n1, f->f[i]); //u=n1/f[i]
mpz_powm(t, t, u, n); //t=a^(n1/f[i]) mod n, 0<t<n
mpz_sub_ui(t, t, 1); //t=a^(n1/f[i])-1, 0<=t<n-1
mpz_gcd(t, t, n); //t=gcd(a^(n1/f[i])-1,n), 0<t<=n
if (mpz_cmp_ui(t, 1) != 0) {
if (mpz_cmp(t, n) < 0) {
gmp_printf("n is divisible by %Zd\n", t);
FLUSH();
result = 0; //definitely composite
break;
}
printf("gcd(%u^((n-1)/f[%d])-1,n)=n\n", *a, i);
FLUSH();
a++;
if (a >= a_limit) {
result = -1; //unknown
break;
}
} else {
printf("gcd(%u^((n-1)/f[%d])-1,n)=1\n", *a, i);
FLUSH();
i++;
if (i >= f->length) {
result = 2; //definitely prime
break;
}
a = (unsigned int *)primes;
}
}
mpz_clear(n1);
mpz_clear(t);
mpz_clear(u);
return result;
}
//main proof 2
// mpz_t n : N
// mpz_t u : F
// mpz_t t : R
int pocklington2(mpz_t n, mpz_t u, mpz_t t) {
int result = -1;
mpz_t s, r, d, x, y;
unsigned int m, lambda;
mpz_init(s);
mpz_init(r);
mpz_init(d);
mpz_init(x);
mpz_init(y);
do {
mpz_set(x, u); //F
mpz_add(x, x, x); //2F
mpz_fdiv_qr(s, r, t, x); //R=2Fs+r, 1<=r<2F
printf("R=2Fs+r, 1<=r<2F\n");
gmp_printf("s=%Zd\n", s);
gmp_printf("r=%Zd\n", r);
FLUSH();
// (mF+1)(2F^2+(r-m)F+1) >= N+1
// (mF+1)(2F^2+rF-mF+1) >= N+1
// mF(2F^2+rF-mF+1)+2F^2+rF-mF+1 >= N+1
// mF(2F^2+rF-mF+1)+2F^2+rF-mF-N >= 0
// -F^2*m^2 + F^2(2F+r)m + F(2F+r)-N >= 0
// F^2*m^2 - F^2(2F+r)m + N-F(2F+r) <= 0
// D=F^2((F(2F+r))^2-4(N-F(2F+r)))
// ceil((F^2(2F+r)-isqrt(D))/2F^2) <= m <=
// floor((F^2(2F+r)+isqrt(D))/2F^2)
// m=max(1,ceil((F^2(2F+r)-isqrt(D))/2F^2))
mpz_set(x, u); //F
mpz_add(x, x, x); //2F
mpz_add(x, x, r); //2F+r
mpz_mul(x, x, u); //F(2F+r)
mpz_mul(d, x, x); //(F(2F+r))^2
mpz_sub(y, n, x); //N-F(2F+r)
mpz_mul_ui(y, y, 4); //4(N-F(2F+r))
mpz_sub(d, d, y); //(F(2F+r))^2-4(N-F(2F+r))
mpz_mul(y, u, u); //F^2
mpz_mul(d, d, y); //F^2((F(2F+r))^2-4(N-F(2F+r)))
printf("let D be discriminant of (mF+1)(2F^2+(r-m)F+1)>=N+1\n");
gmp_printf("D=F^2((F(2F+r))^2-4(N-F(2F+r)))=%Zd\n", d);
FLUSH();
if (mpz_sgn(d) < 0) {
printf("D is negative. probably there is insufficient factor of n-1\n");
FLUSH();
break;
}
mpz_mul(x, x, u); //F^2(2F+r)
mpz_sqrt(d, d); //isqrt(D)
mpz_sub(x, x, d); //F^2(2F+r)-isqrt(D)
mpz_add(y, y, y); //2F^2
mpz_add(x, x, y); //F^2(2F+r)-isqrt(D)+2F^2
mpz_sub_ui(x, x, 1); //F^2(2F+r)-isqrt(D)+2F^2-1
mpz_div(x, x, y); //ceil((F^2(2F+r)-isqrt(D))/2F^2)
if (mpz_cmp_ui(x, 1) < 0) {
mpz_set_ui(x, 1);
}
gmp_printf("m=max(1,ceil((F^2(2F+r)-isqrt(D))/2F^2))=%Zd\n", x);
FLUSH();
if (mpz_cmp_ui(x, MAX_M) > 0) {
printf("m is too large\n");
break;
}
m = mpz_get_ui(x);
if (m > 1) {
mpz_set(x, u); //F
mpz_add_ui(x, x, 1); //F+1
for (lambda = 1; lambda < m; lambda++) {
if (mpz_divisible_p(n, x)) {
gmp_printf("n is divisible by %Zd\n", x);
FLUSH();
result = 0;
break;
}
mpz_add(x, x, u); //next lambda*F+1
}
if (lambda < m) {
break;
}
printf("satisfied that n is not divisible by lambda*F+1, 1<=lambda<m\n");
FLUSH();
}
mpz_set(x, u); //F
mpz_mul(x, x, x); //F^2
mpz_add(x, x, x); //2F^2
mpz_set(y, r); //r
mpz_sub_ui(y, y, m); //r-m
mpz_mul(y, y, u); //(r-m)F
mpz_add(x, x, y); //2F^2+(r-m)F
mpz_add_ui(x, x, 1); //2F^2+(r-m)F+1
mpz_set(y, u); //F
mpz_mul_ui(y, y, m); //mF
mpz_add_ui(y, y, 1); //mF+1
mpz_mul(x, x, y); //(mF+1)(2F^2+(r-m)F+1)
gmp_printf("(mF+1)(2F^2+(r-m)F+1)=%Zd\n", x);
FLUSH();
if (mpz_cmp(n, x) < 0) {
printf("(mF+1)(2F^2+(r-m)F+1) is greater than n\n");
FLUSH();
} else {
printf("(mF+1)(2F^2+(r-m)F+1) is not greater than n\n");
FLUSH();
break;
}
//final test
if (mpz_sgn(s) == 0) {
printf("s is zero\n");
result = 2; //definitely prime
FLUSH();
break;
}
printf("s is not zero\n");
mpz_set(x, r); //r
mpz_mul(x, x, x); //r^2
mpz_set(y, s); //s
mpz_mul_ui(y, y, 8); //8s
mpz_sub(x, x, y); //r^2-8s
gmp_printf("r^2-8s=%Zd\n", x);
mpz_sqrtrem(x, y, x);
printf("r^2-8s=x^2+y, 0<=y<2x+1\n");
gmp_printf("x=%Zd\n", x);
gmp_printf("y=%Zd\n", y);
if (mpz_sgn(y) != 0) {
printf("r^2-8s is not a square\n");
result = 2; //definitely prime
} else {
printf("r^2-8s is a perfect square\n");
result = 0; //definitely composite
}
FLUSH();
} while (0);
mpz_clear(s);
mpz_clear(r);
mpz_clear(d);
mpz_clear(x);
mpz_clear(y);
return result;
}
const char *pp_string[3] = {
"definitely composite",
"probably prime",
"definitely prime",
};
int main(int argc, char *argv[]) {
clock_t tm0, tm1;
factorlist_t f;
mpz_t n, n1, t, u;
int i;
int pp;
printf("%s\n", TITLE);
printf(" powered by GMP %s\n", gmp_version);
printf(" version %s by %s\n", VERSION, AUTHOR);
FLUSH();
switch (argc) {
case 1:
break;
default:
printf(USAGE, argv[0], argv[0]);
FLUSH();
return 0;
}
mpz_init(n);
mpz_init(n1);
mpz_init(t);
mpz_init(u);
for (i = 0; i < FACTORLIST_LENGTH; i++) {
mpz_init(f.f[i]);
}
{
clock_t t = clock();
while ((tm0 = clock()) == t);
}
for (;;) {
//input n and f[i]
FLUSH();
if (gmp_scanf("%Zd", n) != 1) {
break;
}
if (mpz_sgn(n) == 0) {
break;
}
for (i = 0; i < FACTORLIST_LENGTH; i++) {
if (gmp_scanf("%Zd", f.f[i]) != 1) {
break;
}
if (mpz_sgn(f.f[i]) == 0) {
break;
}
}
if (i == 0) {
printf("error: no prime factors specified\n");
FLUSH();
break;
}
if (i == FACTORLIST_LENGTH) { //last box is reserved
printf("error: too many factors specified\n");
FLUSH();
break;
}
f.length = i;
//output n and f[i]
gmp_printf("n=%Zd\n", n);
FLUSH();
for (i = 0; i < f.length; i++) {
gmp_printf("f[%d]=%Zd\n", i, f.f[i]);
FLUSH();
}
//range check
if (mpz_cmp_ui(n, 2) < 0) {
printf("error: n is out of range\n");
FLUSH();
break;
}
for (i = 0; i < f.length; i++) {
if (mpz_cmp_ui(f.f[i], 2) < 0) {
printf("error: f[%d] is out of range\n", i);
FLUSH();
break;
}
}
if (i < f.length) {
break;
}
//pretest
if (mpz_cmp_ui(n, 2) == 0) {
printf("n is definitely prime\n");
FLUSH();
continue;
}
if (mpz_even_p(n)) {
printf("n is divisible by 2\n");
FLUSH();
continue;
}
mpz_sqrtrem(t, u, n);
if (mpz_sgn(u) == 0) {
printf("n is a perfect square\n");
FLUSH();
continue;
}
//prime factor 2 is required
for (i = 0; i < f.length; i++) {
if (mpz_cmp_ui(f.f[i], 2) == 0) {
break;
}
}
if (i == f.length) {
mpz_set_ui(f.f[i], 2); //perhaps reserved last box
f.length++;
gmp_printf("f[%d]=%Zd\n", i, f.f[i]);
FLUSH();
}
//n-1
mpz_sub_ui(n1, n, 1);
//prime factor check
printf("prime factor check\n");
FLUSH();
mpz_set(t, n1);
for (i = 0; i < f.length; i++) {
f.e[i] = 0;
while (mpz_divisible_p(t, f.f[i])) {
mpz_div(t, t, f.f[i]);
f.e[i]++;
}
if (f.e[i] == 0) {
printf("error: f[%d] is not a divisor of n-1 or duplicated\n", i);
FLUSH();
break;
}
pp = mpz_probab_prime_p(f.f[i], PROBAB_P_PARAM);
if (pp == 0) {
printf("error: f[%d] is definitely composite\n", i);
FLUSH();
break;
}
printf("f[%d] is a %s factor of n-1\n", i, pp_string[pp]);
FLUSH();
}
if (i < f.length) {
break;
}
//F
printf("F");
for (i = 0; i < f.length; i++) {
if (f.e[i] == 1) {
printf("%cf[%d]", (i == 0 ? '=' : '*'), i);
} else {
printf("%cf[%d]^%d", (i == 0 ? '=' : '*'), i, f.e[i]);
}
}
printf("\n");
FLUSH();
//F and R
mpz_div(u, n1, t); //F
printf("n-1=F*R\n");
gmp_printf("F=%Zd\n", u);
gmp_printf("R=%Zd\n", t);
FLUSH();
if (mpz_cmp(u, t) <= 0) {
printf("F is not greater than R\n");
} else {
printf("F is greater than R\n");
}
FLUSH();
//main proof
printf("main proof\n");
FLUSH();
pp = pocklington(n, &f);
if (pp == 2 && mpz_cmp(u, t) <= 0) {
pp = pocklington2(n, u, t);
}
//result
if (pp >= 0) {
printf("n is %s\n", pp_string[pp]);
} else {
printf("primality of n has not proved\n");
}
FLUSH();
}
tm1 = clock();
printf("%.8g sec.\n", ((double)(tm1 - tm0)) / CLK_TCK);
FLUSH();
mpz_clear(n);
mpz_clear(n1);
mpz_clear(t);
mpz_clear(u);
for (i = 0; i < FACTORLIST_LENGTH; i++) {
mpz_clear(f.f[i]);
}
return 0;
}