#include <pari/pari.h>

/* A in Z[X], return kernel of (Frob - Id) over Fp[X] / A */
GEN
Berlekamp_Ker(GEN A, GEN p)
{
  long j, N = degpol(A);
  GEN v, w, Q = zeromat(N, N);

  w = v = FpXQ_pow(polx[varn(A)], p, A, p); /* x^p mod (A,p) */
  for (j = 2; j <= N; j++)
  {
    gel(Q, j) = RX_to_RV(w, N);
    gcoeff(Q, j, j) = addis(gcoeff(Q, j, j), -1);
    if (j < N)
    {
      pari_sp av = avma; /*FpXQ_mul is not stack clean*/
      w = gerepileupto(av, FpXQ_mul(w, v, A, p)); /* w *= v (mod A,p)*/
    }
  }
  return FpM_ker(Q, p);
}

long
nbfact(GEN A, GEN p)
{
  pari_sp av = avma;
  GEN K = Berlekamp_Ker(A, p);
  long dim = glength(K);
  avma = av; return dim;
}

int
main()
{
  GEN p, A;
  pari_init(1000000, 2);
  printf("p = "); p = lisGEN(stdin);
  printf("A = "); A = lisGEN(stdin);
  if (!p || !A) err(talker, "read failed");

  if (!BSW_psp(p)) err(talker, "%Z is not prime", p);

  if (typ(A) != t_POL || !FpX_is_squarefree(A, p))
    err(talker, "%Z is not a squarefree polynomial", A);

  pariputsf("... has %ld factor(s) mod %Z.\n", nbfact(A, p), p);
}