```#include <pari/pari.h>

#define swap(x,y)  {GEN  _z=(GEN)x; x=y; y=(long)_z;}
/* Assumes 'a' square Flm. Returns det(a) mod l. IN PLACE (destroys a) */
ulong
Flm_det(GEN a, ulong l)
{
long i, j, k, nbco = lg(a)-1;
ulong p, invp, x = 1UL;

for (i = 1; i < nbco; i++)
{
p = coeff(a,i,i); k=i;
if (!p)
{
do k++; while(k <= nbco && !gcoeff(a,i,k));
if (k > nbco) return 0;
}
if (k != i)
{
swap(a[i], a[k]);
x = l - x;
p = coeff(a,i,i); /* pivot */
}
x = Fl_mul(x, p, l);

invp = Fl_inv(p, l);
for (k = i+1; k <= nbco; k++)
{
ulong m = coeff(a,i,k);
if (m)
{
m = l - Fl_mul(m, invp, l);
for (j = i+1; j <= nbco; j++)
coeff(a,j,k) = Fl_add(coeff(a,j,k), Fl_mul(m, coeff(a,j,i), l), l);
}
}
}
return Fl_mul(x, coeff(a,nbco,nbco), l);
}

static long
center(ulong u, ulong p) { return (long) (u > (p >> 1))? u - p: u; }

/* 27449 = prime(3000) */
byteptr
init_modular(ulong *p) { *p = 27449; return diffptr + 3000; }

/* assume A is a ZM */
GEN
ZM_det(GEN A)
{
pari_sp av = avma;
ulong Dp, p; /* p will loop over prime numbers, Dp = D mod p */
byteptr ptr = init_modular(&p);
GEN D = NULL, qp, q  = gun;

for(;;)
{
pari_sp av2 = avma;

NEXT_PRIME_VIADIFF_CHECK(p, ptr);
Dp = Flm_det(ZM_to_Flm(A, p), p);
avma = av2;
qp = muliu(q, p);
if (!D)
D = stoi( center(Dp, p) );
else
if (Z_incremental_CRT(&D, Dp, q, qp, p)) break;
fprintferr("p = %lu, Dp = %lu, D = %Z\n", p, Dp, D);
q = qp;
}
return gerepileuptoint(av, D);
}

int
main()
{
GEN A;
pari_init(1000000,2);
printf("A = "); A = lisGEN(stdin);
pariputsf( "det(A) is probably %Z\n", ZM_det(A) );
}
```