## 1. Abstract

Primality proving program based on Pocklington's theorem. GNU MP library and GNU C Compiler are required.

## 3. Usage

```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 | pock```

## 4. History

August 17, 2003 ... Version 0.2.1 is available.

## 5. Program

pock021.c

```//Primality proving program based on Pocklington's theorem

#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 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_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_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_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
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;
}
}
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_set(y, r);  //r
mpz_sub_ui(y, y, m);  //r-m
mpz_mul(y, y, u);  //(r-m)F
mpz_set(y, u);  //F
mpz_mul_ui(y, y, m);  //mF
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("  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;
}```

## 6. References

• J. Brillhart, D. H. Lehmer and J. L. Selfridge, "New primality criteria and factorizations of 2m±1," Math. Comp., 29 (1975) 620--647.