c_det.c

#include "c_det.h"

typedef __int128  int128_t;

// -----------------------------------------------------------------------------

// we assume a and b are already mod p
inline int64_t sub_modp( int64_t p , int64_t a , int64_t b )
{
  if (b <= a) 
    { return (a - b); }
  else
    { return (a + p - b); }
}

inline int64_t mul_modp( int64_t p0 , int64_t a0 , int64_t b0 )
{
  int128_t p = p0;
  int128_t a = a0;
  int128_t b = b0;
  int128_t c = a*b;
  c = c % p;
  return ((int64_t)c);
}

// -----------------------------------------------------------------------------

int64_t euclid( int64_t p , int64_t x1_ , int64_t x2_ , int64_t u_ , int64_t v_ )
{
  int64_t halfp1 = (p + 1) >> 1;

  int64_t x1 = x1_; 
  int64_t x2 = x2_;
  int64_t u  = u_;
  int64_t v  = v_;

  while( (u!=1) && (v!=1) )
  {
    while (!(u & 1))
    { // u even
      u = u >> 1;
      if (x1 & 1) { /* x1 odd */  x1 = (x1 >> 1) + halfp1; } else { x1 = x1 >> 1; }
    }
      
    while (!(v & 1))
    { // v even
      v = v >> 1;
      if (x2 & 1) { /* x2 odd */  x2 = (x2 >> 1) + halfp1; } else { x2 = x2 >> 1; }
    }

    if (u >= v)
    {
      u  = u - v;
      if ( x1 >= x2 ) { x1 = (x1 - x2); } else { x1 = (x1 + p - x2); } 
    }
    else 
    {
      v  = v - u;
      if ( x2 >= x1)  { x2 = (x2 - x1); } else { x2 = (x2 + p - x1); } 
    }

  }

  if (u==1) { return x1; }
  if (v==1) { return x2; }
  return 0;                       // shouldn't happen
}

// -----------------------------------------------------------------------------

inline int64_t div_modp( int64_t p , int64_t a , int64_t b )
{
  // return mul_modp( p , a , inv_modp( p , b ) );
  return euclid( p , a , 0 , b , p );}

// mod p inverse using the binary Euclidean algorithm 
int64_t inv_modp( int64_t p , int64_t a )
{
  return euclid( p , 1 , 0 , a , p );
}

// -----------------------------------------------------------------------------

// determinant mod p (64 bit), using Gauss elimination
int64_t det_modp(int64_t p, int n, int64_t *mat)
{
  // safety first
  for (int i=0;i<n*n;i++) { if ((mat[i] >= p) || (mat[i]<0)) { mat[i] = mat[i] % p; } }

  int negative = 0;

  for (int i=0;i<n-1;i++)
  {
    int64_t *row = mat + i*n;
   
    // find pivot element
    int j; 
    for (j=i;j<n;j++) { if (row[j] != 0) break; }
    if ( (j >= n) || (row[j] == 0) ) { return 0; }
     
    if (j > i)
    { // exchange columns
      int64_t *q = row;      
      for (int k=i;k<n;k++)
      { 
        int64_t x;
        x    = q[i];
        q[i] = q[j];
        q[j] = x;
        q += n;     
      }
      negative = negative ^ 1;     // track the sign changes
    }

    // zero out the i-th column
    int64_t *q = row + n;
    for (int k=i+1;k<n;k++)
    { 
      int64_t m = div_modp( p , q[i] , row[i] );
      q[i] = 0;  
      for (int l=i+1;l<n;l++) 
      { 
        q[l] = sub_modp( p , q[l] , mul_modp( p , m , row[l] ) );
      }
      q += n;
    }
  }

  int64_t det = mat[0];
  for (int i=1;i<n;i++) { det = mul_modp( p , det , mat[i*(n+1)] ); }

  if ((negative) && (det!=0)) { return (p-det); } else { return det; }
}