/*
  Author:  Pate Williams (c) 1997

  "Algorithm 3.4.6 (Cantor-Zassenhaus Split).
  Given A contained in Fp[X], p > 2, this
  algorithm outputs its irreducible factors
  (which are all of degree d)." -Henri Cohen-
  See "A Course in Computational Algebraic
  Number Theory" by Henri Cohen page 128.
*/

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

#define DEBUG
#define POLY_SIZE 256

void error(char *message)
{
  fprintf(stderr, "%s", message);
  exit(1);
}

verylong **allocate_very_matrix(long m, long n)
{
  long i;
  verylong **matrix = calloc(m, sizeof(verylong *));

  if (!matrix)
    error("Failure in allocate_very_matrix.");
  for (i = 0; i < m; i++) {
    matrix[i] = calloc(n, sizeof(verylong));
    if (!matrix[i])
      error("Failure in allocate_very_matrix.");
  }
  return matrix;
}

void free_very_matrix(verylong **zm, long m, long n)
{
  long i, j;

  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      zfree(&zm[i][j]);
    free(zm[i]);
  }
  free(zm);
}

void cont(long dA, verylong *zA, verylong *zc)
/* content of a polynomial */
{
  long i;
  verylong zd = 0;

  zgcd(zA[0], zA[1], zc);
  for (i = 2; i <= dA; i++) {
    zgcd(*zc, zA[i], &zd);
    zcopy(zd, zc);
  }
  zfree(&zd);
}

void pp(long dA, verylong *zA)
/* primitive part of a polynomial */
{
  long i;
  verylong za = 0, zc = 0, zr = 0;

  cont(dA, zA, &zc);
  for (i = 0; i <= dA; i++) {
    zdiv(zA[i], zc, &za, &zr);
    zcopy(za, &zA[i]);
  }
  zfree(&za);
  zfree(&zc);
  zfree(&zr);
}

void zpoly_mul(long m, long n, verylong *za, verylong *zb,
               verylong *zc, long *p)
{
  long i, j, k;
  verylong zd = 0, zai = 0, zbk = 0, zsum = 0, zterm = 0;

  *p = m + n;
  for (k = 0; k <= *p; k++) {
    zzero(&zsum);
    for (i = 0; i <= k; i++) {
      j = k - i;
      if (i > m) zzero(&zai);
      else zcopy(za[i], &zai);
      if (j > n) zzero(&zbk);
      else zcopy(zb[j], &zbk);
      zmul(zai, zbk, &zterm);
      zcopy(zsum, &zd);
      zadd(zterm, zd, &zsum);
    }
    zcopy(zsum, &zc[k]);
  }
  zfree(&zd);
  zfree(&zai);
  zfree(&zbk);
  zfree(&zsum);
  zfree(&zterm);
}

void zspow(verylong za, long n, verylong *zb)
{
  int c = zscompare(za, 0l);
  verylong zc = 0;

  if (c > 0)
    zsexp(za, n, zb);
  else if (c < 0) {
    zcopy(za, &zc);
    zabs(&zc);
    zsexp(zc, n, zb);
    if (n & 1) znegate(zb);
  }
  else
    zzero(zb);
  zfree(&zc);
}

void zpoly_div(long m, long n, verylong *zu, verylong *zv,
               verylong *zq, verylong *zr, long *p, long *s)
{
  long j, jk, k, nk;
  verylong za = 0, zb = 0, zvn = 0;

  zcopy(zv[n], &zvn);
  for (j = 0; j <= m; j++)
    zcopy(zu[j], &zr[j]);
  if (m < n) {
    *p = 0, *s = m;
    zzero(&zq[0]);
  }
  else {
    *p = m - n, *s = n - 1;
    for (k = *p; k >= 0; k--) {
      nk = n + k;
      zspow(zvn, k, &za);
      zmul(zr[nk], za, &zq[k]);
      for (j = nk - 1; j >= 0; j--) {
        jk = j - k;
        if (jk >= 0) {
          zmul(zvn, zr[j], &za);
          zmul(zr[nk], zv[jk], &zb);
          zsub(za, zb, &zr[j]);
        }
        else {
          zcopy(zr[j], &za);
          zmul(zvn, za, &zr[j]);
        }
      }
    }
    while (*p > 0 && zscompare(zq[*p], 0l) == 0) *p = *p - 1;
    while (*s > 0 && zscompare(zr[*s], 0l) == 0) *s = *s - 1;
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zvn);
}

void zpoly_pow(long degreeA, long degreem, verylong zn,
               verylong *zA, verylong *zm, verylong *zs,
               long *ds)
{
  long dp, dq, dx = degreeA, i;
  verylong za = 0, zb = 0, zp[POLY_SIZE], zq[POLY_SIZE],
           zx[POLY_SIZE], zy[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zp[i] = zq[i] = zx[i] = zy[i] = 0;
  *ds = 0;
  zcopy(zn, &za);
  zone(&zs[0]);
  for (i = 0; i <= dx; i++) zcopy(zA[i], &zx[i]);
  while (zscompare(za, 0l) > 0) {
    if (zodd(za)) {
      /* s = (s * x) % m; */
      zpoly_mul(*ds, dx, zs, zx, zp, &dp);
      zpoly_div(dp, degreem, zp, zm, zq, zs, &dq, ds);
    }
    zcopy(za, &zb);
    zrshift(zb, 1l, &za);
    if (zscompare(za, 0l) > 0) {
      /* x = (x * x) % m; */
      for (i = 0; i <= dx; i++) zcopy(zx[i], &zy[i]);
      zpoly_mul(dx, dx, zx, zy, zp, &dp);
      zpoly_div(dp, degreem, zp, zm, zq, zx, &dq, &dx);
    }
  }
  zfree(&za);
  zfree(&zb);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zp[i]);
    zfree(&zq[i]);
    zfree(&zx[i]);
    zfree(&zy[i]);
  }
}

void zpoly_sub(long da, long db, verylong *za, verylong *zb,
               verylong *zc, long *dc)
/* c = a - b */
{
  long i;
  verylong zz = 0;

  zzero(&zz);
  if (da >= db) {
    for (i = 0; i <= db; i++)
      zsub(za[i], zb[i], &zc[i]);
    for (i = db + 1; i <= da; i++)
      zcopy(za[i], &zc[i]);
    *dc = da;
  }
  else {
    for (i = 0; i <= da; i++)
      zsub(za[i], zb[i], &zc[i]);
    for (i = da + 1; i <= db; i++)
      zsub(zz, zb[i], &zc[i]);
    *dc = db;
  }
  zfree(&zz);
}

void zpoly_gcd(long degreeA, long degreeB, long p,
               verylong *zA, verylong *zB, verylong *za,
               long *da)
{
  int nonzero = 0, zero;
  long db, dq, dr, i;
  verylong zc = 0, zp = 0;
  verylong zb[POLY_SIZE], zq[POLY_SIZE], zr[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zb[i] = zq[i] = zr[i] = 0;
  if (degreeA > degreeB) {
    *da = degreeA;
    db = degreeB;
    for (i = 0; i <= *da; i++) zcopy(zA[i], &za[i]);
    for (i = 0; i <= db; i++) zcopy(zB[i], &zb[i]);
  }
  else {
    *da = degreeB;
    db = degreeA;
    for (i = 0; i <= *da; i++) zcopy(zB[i], &za[i]);
    for (i = 0; i <= db; i++) zcopy(zA[i], &zb[i]);
  }
  for (i = 0; i <= db && !nonzero; i++)
    nonzero = zscompare(zb[i], 0l) != 0;
  while (nonzero) {
    zpoly_div(*da, db, za, zb, zq, zr, &dq, &dr);
    zintoz(p, &zp);
    for (i = 0; i <= dr; i++) {
      zcopy(zr[i], &zc);
      zmod(zc, zp, &zr[i]);
    }
    zero = 1;
    for (i = dr; i >= 0 && zero; i--) {
      zero = zscompare(zr[i], 0l) == 0;
      if (zero && dr > 0) dr--;
    }
    for (i = 0; i <= db; i++) zcopy(zb[i], &za[i]);
    *da = db;
    for (i = 0; i <= dr; i++) zcopy(zr[i], &zb[i]);
    db = dr;
    nonzero = 0;
    for (i = 0; i <= db && !nonzero; i++)
      nonzero = zscompare(zb[i], 0l) != 0;
  }
  zfree(&zc);
  zfree(&zp);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zb[i]);
    zfree(&zq[i]);
    zfree(&zr[i]);
  }
}

void zpoly_write(char *label, long dA, verylong *zA)
{
  long i;

  printf("%s", label);
  for (i = dA; i >= 0; i--) {
    zwrite(zA[i]);
    printf(" ");
  }
  printf("\n");
}

void Cantor_Zassenhaus(long d, long dA, long p,
                       verylong *zA, verylong *zX,
                       verylong **zB, long *dX,
                       long *dB, long *count)
{
  int zero;
  long dC, dP, dR, dT = 2 * d - 1, i, k = dA / d;
  static verylong za = 0, zb = 0, ze = 0, zp = 0;
  static verylong zC[POLY_SIZE] = {0};
  static verylong zP[POLY_SIZE] = {0};
  static verylong zR[POLY_SIZE] = {0};
  static verylong zT[POLY_SIZE] = {0};

  zintoz(p, &zp);
  zsexp(zp, d, &za);
  zsadd(za, - 1, &zb);
  zlshift(zb, 1l, &ze);
  if (k <= 1) {
    pp(dA, zA);
    for (i = 0; i <= dA; i++)
     zcopy(zA[i], &zB[*count][i]);
    dB[*count] = dA;
    *count = *count + 1;
    zfree(&za);
    zfree(&zb);
    zfree(&ze);
    zfree(&zp);
    for (i = 0; i < POLY_SIZE; i++) {
      zfree(&zC[i]);
      zfree(&zP[i]);
      zfree(&zR[i]);
      zfree(&zT[i]);
    }
    return;
  }
  L2:
  zone(&zT[dT]);
  for (i = 0; i < dT; i++) zrandomb(zp, &zT[i]);
  zpoly_pow(dT, dA, ze, zT, zA, zP, &dP);
  for (i = 0; i <= dP; i++)
    zintoz(zsmod(zP[i], p), &zP[i]);
  pp(dP, zP);
  zero = 1;
  for (i = dP; zero && i >= 0; i--) {
    zero = zscompare(zP[i], 0l) == 0;
    if (zero) dP--;
  }
  zsadd(zP[0], - 1l, &za);
  zcopy(za, &zP[0]);
  zpoly_gcd(dA, dP, p, zA, zP, zC, &dC);
  zsadd(zP[0], + 1l, &za);
  zcopy(za, &zP[0]);
  #ifdef DEBUG
  zpoly_write("T = ", dT, zT);
  zpoly_write("B = ", dC, zC);
  #endif
  if (dC == 0 || dC == dA) goto L2;
  pp(dC, zC);
  #ifdef DEBUG
  zpoly_write("A = ", dA, zA);
  zpoly_write("P = ", dP, zP);
  zpoly_write("T = ", dT, zT);
  zpoly_write("B = ", dC, zC);
  #endif
  for (i = 0; i <= dC; i++)
    zintoz(zsmod(zC[i], p), &zC[i]);
  zpoly_div(dA, dC, zA, zC, zP, zR, &dP, &dR);
  for (i = 0; i <= dC; i++)
    zcopy(zC[i], &zA[i]);
  for (i = 0; i <= dP; i++)
    zintoz(zsmod(zP[i], p), &zX[i]);
  *dX = dP;
  #ifdef DEBUG
  zpoly_write("P = ", dP, zP);
  zpoly_write("X = ", *dX, zX);
  #endif
  Cantor_Zassenhaus(d, dC, p, zC, zX, zB, dX, dB, count);
  for (i = 0; i <= *dX; i++)
    zcopy(zX[i], &zA[i]);
  Cantor_Zassenhaus(d, dP, p, zA, zX, zB, dX, dB, count);
}

int main(void)
{
  char answer[256];
  long count, dA, dB[32], dX, i, m = 32, p;
  verylong zA[POLY_SIZE] = {0};
  verylong zX[POLY_SIZE] = {0};
  verylong **zB = allocate_very_matrix(m, POLY_SIZE);

  do {
    printf("degree of polynomial = ");
    scanf("%d", &dA);
    printf("modulus = ");
    scanf("%d", &p);
    for (i = 0; i <= dA; i++) {
      printf("A[%ld] = ", i);
      zread(&zA[i]);
    }
    count = 0;
    Cantor_Zassenhaus(dA / 2, dA, p, zA, zX, zB, &dX, dB, &count);
    for (i = 0; i < count; i++) {
      zpoly_write("factor = ", dB[i], zB[i]);
    }
    printf("another polynomial (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  for (i = 0; i <= dA; i++) zfree(&zA[i]);
  free_very_matrix(zB, m, POLY_SIZE);
  return 0;
}