```#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);
}
```