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