/*
  Author:  Pate Williams (c) 1997

  Lucas-Lehmer primality test.
  See "Implementation of a New Primality Test",
  Mathematics of Computation, Volume 48, Number 177,
  January 1987, pages 103 - 121 by H. Cohen and
  A. K. Lenstra.
*/

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"

#define BOUND 1000000l
#define t 25200l

int trial_division(long bound, verylong zn,
                   long *lm, long *lmSize,
                   long *lp, long *lpSize,
                   verylong *zfm,
                   verylong *zfp,
                   verylong *znp,
                   verylong *znm,
                   verylong *zrm,
                   verylong *zrp)
/* returns 0 if number is composite 1 if prime
   2 if proved prime by trial division */
{
  int factored = 0, flag = 0;
  long b, p, r;
  verylong za = 0, zs = 0, zn1 = 0;

  *lmSize = *lpSize = 0;
  zsadd(zn, - 1l, znm);
  zsadd(zn, + 1l, znp);
  zcopy(*znm, zrm);
  zcopy(*znp, zrp);
  zcopy(*znp, &zn1);
  zpstart2();
  p = zpnext();
  zsqrt(zn, &zs, &za);
  if (zscompare(zs, bound) > 0)
    b = bound;
  else {
    flag = 1;
    b = ztoint(zs);
  }
  while (!factored && p <= b) {
    r = zsmod(*znp, p);
    factored = r == 1;
    if (!factored) {
      if (r == 0) {
        do {
          zsdiv(*zrp, p, &za);
          zcopy(za, zrp);
        } while (zsmod(*zrp, p) == 0);
        lp[*lpSize] = p;
        *lpSize = *lpSize + 1;
      }
      if (r == 2) {
        do {
          zsdiv(*zrm, p, &za);
          zcopy(za, zrm);
        } while (zsmod(*zrm, p) == 0);
        lm[*lmSize] = p;
        *lmSize = *lmSize + 1;
      }
    }
    p = zpnext();
  }
  zdiv(*znm, *zrm, zfm, &za);
  zdiv(*znp, *zrp, zfp, &za);
  zfree(&za);
  zfree(&zs);
  zfree(&zn1);
  if (!factored && flag && p > b) factored = 2;
  else factored = !factored;
  return factored;
}

void ring_mul(long a, long r, long u,
              verylong zn,
              verylong *zx, verylong *zy,
              verylong *zz)
{
  verylong za = 0, zb = 0, zc = 0;
  verylong zp0 = 0, zp1 = 0, zs0 = 0, zs1 = 0;

  zmulmod(zx[0], zy[0], zn, &zp0);
  zmulmod(zx[1], zy[1], zn, &zp1);
  zaddmod(zx[0], zx[1], zn, &zs0);
  zaddmod(zy[0], zy[1], zn, &zs1);
  if (r == 1) {
    zsmulmod(zp1, a, zn, &za);
    zaddmod(zp0, za, zn, &zz[0]);
    zmulmod(zs0, zs1, zn, &za);
    zsubmod(za, zp0, zn, &zb);
    zsubmod(zb, zp1, zn, &zz[1]);
  }
  else if (r == 3) {
    zaddmod(zp0, zp1, zn, &zz[0]);
    zsmulmod(zp1, u - 1, zn, &za);
    zmulmod(zs0, zs1, zn, &zb);
    zaddmod(zb, za, zn, &zc);
    zsubmod(zc, zp0, zn, &zz[1]);
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zp0);
  zfree(&zp1);
  zfree(&zs0);
  zfree(&zs1);
}

void ring_sqr(long u, verylong zn,
              verylong *zx, verylong *zz)
{
  verylong za = 0, zb = 0, zs = 0;

  zsmulmod(zx[1], u, zn, &za);
  zsmulmod(zx[0], 2l, zn, &zb);
  zaddmod(za, zb, zn, &zs);
  zmulmod(zx[0], zs, zn, &za);
  zsadd(za, - 1l, &zb);
  zmod(zb, zn, &zz[0]);
  zmulmod(zx[1], zs, zn, &zz[1]);
  zfree(&za);
  zfree(&zb);
  zfree(&zs);
}

void ring_pow(long a, long u,
              verylong ze,
              verylong zn,
              verylong *zx,
              verylong *zz)
/* calculates z = x ^ e in the ring
   (Z / nZ)[T]/(T ^ 2 - u * T - 1) */
{
  long i, r = zsmod(zn, 4l);
  verylong za = 0, zb = 0;
  verylong zA[2], zS[2], zT[2];

  zA[0] = zA[1] = zS[0] = zS[1] = zT[0] = zT[1] = 0;
  zcopy(ze, &za);
  zzero(&zA[1]);
  zone(&zA[0]);
  zcopy(zx[0], &zS[0]);
  zcopy(zx[1], &zS[1]);
  while (zscompare(za, 0l) != 0) {
    if (zodd(za)) {
      ring_mul(a, r, u, zn, zA, zS, zz);
      zcopy(zz[0], &zA[0]);
      zcopy(zz[1], &zA[1]);
    }
    zrshift(za, 1l, &zb);
    zcopy(zb, &za);
    if (zscompare(za, 0l) != 0) {
      ring_sqr(u, zn, zS, zT);
      zcopy(zT[0], &zS[0]);
      zcopy(zT[1], &zS[1]);
    }
  }
  zfree(&za);
  zfree(&zb);
  for (i = 0; i < 2; i++) {
    zfree(&zA[i]);
    zfree(&zS[i]);
    zfree(&zT[i]);
  }
}

int test_for_nm(int *flag,
                verylong zn,
                verylong znm,
                verylong zp,
                verylong *zprod,
                verylong **zbeta)
/* returns - 1 if test fails, 0 if n is composite,
   1 otherwise */
{
  int found = 0, value;
  long i, l, p, q, x;
  verylong za = 0, zb = 0, zo = 0, zx = 0;

  zpstart2();
  for (i = 0; !found && i < 50; i++) {
    x = zpnext();
    zdiv(znm, zp, &za, &zb);
    zintoz(x, &zx);
    zexpmod(zx, za, zn, &zb);
    found = zscompare(zb, 1l) != 0;
  }
  if (found) {
    zexpmod(zx, znm, zn, &za);
    value = zscompare(za, 1l) == 0;
    if (value) {
      zcopy(*zprod, &zo);
      zsadd(zb, - 1l, &za);
      zmulmod(*zprod, za, zn, &zb);
      zcopy(zb, zprod);
      value = zscompare(*zprod, 0l) == 0;
      if (value) {
        zgcd(zo, zn, &za);
        printf("factor = ");  zwriteln(za);
        value = 0;
      }
      else if (zprobprime(zp, 5l)) {
        zpstart();
        l = 1;
        p = ztoint(zp);
        do {
          q = pow(p, l);
          if (t % q == 0) {
            flag[q] = 1;
            for (i = 0; i < q - 1; i++) {
              zsmul(znm, i, &za);
              zsdiv(za, q, &zb);
              zexpmod(zx, zb, zn, &zbeta[q][i]);
            }
          }
          l++;
        } while (q < t);
        value = 1;
      }
      else value = 1;
    }
  }
  else value = - 1;
  if (!value) printf("composite by n - 1 test\n");
  zfree(&za);
  zfree(&zb);
  zfree(&zo);
  zfree(&zx);
  return value;
}

int test_for_np(long a, long u,
                verylong zn,
                verylong znp,
                verylong zp,
                verylong *zprod)
/* returns
  - 1 if test fails, 0 if n is composite,
  + 1 otherwise */
{
  int factored = 0, found = 0, value;
  long d, m;
  verylong za = 0, zd = 0, ze = 0, zi = 0;
  verylong zx[2], zz[2];

  zx[0] = zx[1] = zz[0] = zz[1] = 0;
  zdiv(znp, zp, &ze, &za);
  for (m = 1; !factored && !found && m <= 50; m++) {
    d = m * (m + u) - a;
    zintoz(d, &zd);
    zinvmod(zd, zn, &zi);
    factored = zscompare(zi, 0l) == 0;
    if (!factored) {
      zintoz(m * m + a, &za);
      zmulmod(za, zi, zn, &zx[0]);
      zintoz((2 * m + u) * a, &za);
      zmulmod(za, zi, zn, &zx[1]);
      ring_pow(a, u, ze, zn, zx, zz);
      found = zscompare(zz[0], 1l) != 0 ||
              zscompare(zz[1], 0l) != 0;
    }
  }
  if (!factored && found) {
    ring_pow(a, u, znp, zn, zx, zz);
    factored = zscompare(zz[0], 1l) == 1 &&
               zscompare(zz[1], 0l) == 0;
    if (!factored) {
      if (zscompare(zx[0], 0l) != 0)
        zmulmod(*zprod, zx[0], zn, &za);
      else
        zmulmod(*zprod, zx[1], zn, &za);
      factored = zscompare(za, 0l) == 0;
      if (factored) {
        zgcd(*zprod, zn, &za);
        printf("factor: ");
        zwriteln(za);
      }
      else zcopy(za, zprod);
      if (factored) value = 0; else value = 1;
    }
  }
  else if (factored) value = 0;
  else value = - 1;
  zfree(&za);
  zfree(&zd);
  zfree(&ze);
  zfree(&zi);
  zfree(&zx[0]);
  zfree(&zx[1]);
  zfree(&zz[0]);
  zfree(&zz[1]);
  if (!value) printf("composite by n + 1 test\n");
  return value;
}

void norm(long a, long u, verylong zn,
          verylong *znorm, verylong *zx)
{
  verylong za = 0, zb = 0, zc = 0;

  zmulmod(zx[0], zx[0], zn, &za);
  zmulmod(zx[0], zx[1], zn, &zb);
  zsmulmod(zb, u, zn, &zc);
  zaddmod(za, zc, zn, &zb);
  zmulmod(zx[1], zx[1], zn, &za);
  zsmulmod(za, - a, zn, &zc);
  zaddmod(zb, zc, zn, znorm);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
}

int Lucas_Lehmer(int *flag, long bound, long lmSize,
                 long *lm, long lpSize, long *lp,
                 verylong zn,
                 verylong zfm, verylong zfp,
                 verylong znm, verylong znp,
                 verylong zrm, verylong zrp,
                 verylong **zbeta)
/* returns
  - 4 could not find suitable a,
  - 3 could not find suitable u,
  - 2 if n - 1 test fails,
  - 1 n + 1 test fails, 0 if n is composite,
  + 1 if passed  */
{
  int f4_1, factored = 0, found = 0, test = 1;
  long a, d, i, l, m, p, r, u;
  verylong za = 0, zd = 0, ze = 0, zi = 0, zm = 0;
  verylong zp = 0, zu = 0, zprod = 0, zx[2], zz[2];

  zx[0] = zx[1] = zz[0] = zz[1] = 0;
  zone(&zprod);
  for (i = 0; test && i < lmSize; i++) {
    p = lm[i];
    zintoz(p, &zp);
    test = test_for_nm(flag, zn, znm, zp, &zprod, zbeta);
  }
  if (test == 1) {
    if (zcompare(zfm, zfp) > 0) zcopy(zfm, &za);
    else zcopy(zfp, &za);
    zmul(za, zfm, &zd);
    zmul(zd, zfp, &ze);
    zintoz(bound, &za);
    zsexp(za, 3l, &zm);
    zmul(ze, zm, &za);
    f4_1 = zcompare(zn, za) < 0;
    if (f4_1 && zscompare(zrm, 1l) != 0)
      test = test_for_nm(flag, zn, znm, zrm, &zprod, zbeta);
  }
  if (test == 1) {
    r = zsmod(zn, 4l);
    if (r == 1) {
      zrshift(znm, 1l, &ze);
      zpstart2();
      for (i = 1; !found && i <= 50; i++) {
        a = zpnext();
        zintoz(a, &za);
        zexpmod(za, ze, zn, &zm);
        found = zcompare(zm, znm) == 0;
      }
      if (!found) test = - 4;
      else {
        l = 1;
        do {
          p = pow(2, l);
          zintoz(a, &za);
          if (t % p == 0 && flag[p]) {
            for (i = 0; i < p; i++) {
              zsmul(znm, i, &zd);
              zsdiv(zd, p, &ze);
              zexpmod(za, ze, zn, &zbeta[p][i]);
            }
          }
          l++;
        } while (p < sqrt(t));
      }
    }
    else {
      a = 1;
      for (i = 1; !found && i <= 50; i++) {
        zintoz(i * i + 4l, &zu);
        found = zjacobi(zu, zn) == - 1;
        if (found) u = i;
      }
      if (!found) test = - 3;
      else {
        found = 0;
        for (m = 1; !factored && !found && m <= 50; m++) {
          d = m * (m + u) - a;
          zintoz(d, &zd);
          zinvmod(zd, zn, &zi);
          factored = zscompare(zi, 0l) == 0;
          if (!factored) {
            found = 1;
            zintoz(m * m + a, &za);
            zmulmod(za, zi, zn, &zx[0]);
            zintoz((2 * m + u) * a, &za);
            zmulmod(za, zi, zn, &zx[1]);
          }
        }
        if (!factored && !found) test = - 2;
        else if (found) {
          zone(&zx[0]);
          zintoz(u, &zx[1]);
          zrshift(znp, 1l, &ze);
          ring_pow(a, u, ze, zn, zx, zz);
          test = zcompare(zz[0], znm) == 0 &&
                 zscompare(zz[1], 0l) == 0;
        }
        else test = 0;
      }
    }
  }
  else if (test == - 1) test = - 2;
  if (test == 1) {
    for (i = 0; test && i < lpSize; i++) {
      zintoz(lp[i], &zp);
      test = test_for_np(a, u, zn, znp, zp, &zprod);
    }
    if (test == 1) {
      if (f4_1) {
        if (zscompare(zrp, 1l) == 0) test = 2;
        else  {
          test = test_for_np(a, u, zn, znp, zrp, &zprod);
          if (test) test = 2;
        }
      }
    }
  }
  if (test == 1) {
    zgcd(zprod, zn, &zd);
    test = zscompare(zd, 1l) == 0;
  }
  zfree(&za);
  zfree(&zd);
  zfree(&ze);
  zfree(&zi);
  zfree(&zm);
  zfree(&zp);
  zfree(&zu);
  zfree(&zprod);
  zfree(&zx[0]);
  zfree(&zx[1]);
  zfree(&zz[0]);
  zfree(&zz[1]);
  return test;
}

int main(void)
{
  char answer[256];
  int flag[256], td;
  long cols, rows;
  long lm[256], lp[256], lmSize, lpSize;
  long i, j, l, n = 0, p, q, r[32], tt;
  verylong zn = 0;
  verylong zfm = 0, zfp = 0, znm = 0, znp = 0;
  verylong zrm = 0, zrp = 0;
  verylong **zbeta;

  zpstart2();
  do {
    p = zpnext();
    if (t % p == 0) {
      r[n++] = p;
      l = 2;
      do {
        q = pow(p, l);
        if (t % q == 0) r[n++] = q;
        l++;
      } while (q < sqrt(t));
    }
  } while (p < sqrt(t));
  for (i = 0; i < n - 1; i++)
    for (j = i + 1; j < n; j++)
      if (r[i] > r[j])
        tt = r[i], r[i] = r[j], r[j] = tt;
  cols = rows = r[n - 1] + 1;
  zbeta = calloc(rows, sizeof(verylong *));
  if (!zbeta) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    exit(1);
  }
  for (i = 0; i < rows; i++) {
    zbeta[i] = calloc(cols, sizeof(verylong));
    if (!zbeta[i]) {
      fprintf(stderr, "fatal error\ninsufficient memory\n");
      exit(1);
    }
  }
  do {
    printf("enter the number to be tested below:\n");
    zread(&zn);
    zwriteln(zn);
    if (zprobprime(zn, 5l)) {
      td = trial_division(BOUND, zn, lm, &lmSize,
                         lp, &lpSize,  &zfm, &zfp, &znp,
                         &znm, &zrm, &zrp);
      if (td == 1) {
        for (i = 0; i < 256; i++) flag[i] = 0;
        switch (Lucas_Lehmer(flag, BOUND, lmSize,
                             lm, lpSize, lp, zn,
                             zfm, zfp, znm, znp,
                             zrm, zrp, zbeta)) {
          case - 4 :
            printf("failed could not find suitable a\n");
            break;
          case - 3 :
            printf("failed could not find suitable u\n");
            break;
          case - 2 :
            printf("Lucas-Lehmer n - 1 test failed\n");
            break;
          case - 1 :
            printf("Lucas-Lehmer n + 1 test failed\n");
            break;
          case + 0 :
            break;
          case + 1 :
            printf("number passed Lucas-Lehmer test\n");
            break;
          case + 2 :
            printf("number proven prime by Lucas-Lehmer test\n");
            break;
        }
      }
      else if (!td)
        printf("number composite by trial division\n");
      else
        printf("number proven prime by trial division\n");
    }
    else printf("number composite by Miller-Rabin\n");
    printf("another number (n or y) ? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zn);
  zfree(&zfm);
  zfree(&zfp);
  zfree(&znm);
  zfree(&znp);
  zfree(&zrm);
  zfree(&zrp);
  for (i = 0; i < rows; i++) {
    for (j = 0; j < cols; j++)
      zfree(&zbeta[i][j]);
    free(zbeta[i]);
  }
  free(zbeta);
  return 0;
}