Table of contents

  1. Abstract
  2. Results
  3. Usage
  4. History
  5. Program
  6. References
  7. Links to related websites

1. Abstract

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

2. Results

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
//  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;
}

6. References

7. Links to related websites