// Copyright Michael Monagan 2019-2020
// Compile with  gcc -O3 -shared -o linalg.so -fPIC linalg.c

#define LONG long long int
#include <stdio.h>
#include <stdlib.h>
#include "int128g.c"


/******************************************************************************************/
/*       Zp arithmetic                                                                    */
/******************************************************************************************/

LONG add64s(LONG a, LONG b, LONG p) { LONG t; t = (a-p)+b; t += (t>>63) & p; return t; }
LONG sub64s(LONG a, LONG b, LONG p) { LONG t; t = a-b; t += (t>>63) & p; return t; }
LONG neg64s(LONG a, LONG p) { return (a==0) ? 0 : p-a; }
LONG mul64s(LONG a, LONG b, LONG p) {
        LONG q, r;
        __asm__ __volatile__(           \
        "       mulq    %%rdx           \n\t" \
        "       divq    %4              \n\t" \
        : "=a"(q), "=d"(r) : "0"(a), "1"(b), "rm"(p));
        return r;
}

/* c^(-1) mod p assuming 0 < c < p < 2^63 */
LONG inv64s( LONG c, LONG p )
{   LONG d,r,q,r1,c1,d1;
    d = p;
    c1 = 1;
    d1 = 0;
    while( d != 0 ) {
        q = c / d;
        r = c - q*d; r1 = c1 - q*d1;
        c = d; c1 = d1;
        d = r; d1 = r1;
    }
    if( c!=1 ) return( 0 );
    if( c1 < 0 ) c1 += p;
    return( c1 );
}

/******************************************************************************************/
/* Linear algebra routines                                                                */
/******************************************************************************************/

LONG * matrix64s( int n )
{   LONG *A; LONG N;
    N = n;
    N = n*N;
    N = sizeof(LONG) * N;
    A = (LONG *) malloc(N);
    return A;
}

/* print an array in form [a0,a1,...,an-1] */
void vecprint64s( LONG *A, int n )
{   int i;
    printf("[");
    for( i=0; i<n; i++ ) { printf("%lld",A[i]); if( i<n-1 ) printf(", "); }
    printf("]");
    return;
}

void matprint64s( LONG *A, int n )
{   int i;
    printf("[");
    for( i=0; i<n-1; i++ ) { vecprint64s(A+i*n,n); printf(",\n"); }
    vecprint64s(A+(n-1)*n,n); 
    printf("];\n");
    return;
}

LONG det64s( LONG * A, int n, LONG p ) {
// M = [[a,b,c],[d,e,f],[1,2,3],[4,5,6]]  M[i][j]
// A = [a,b,c,d,e,f,1,2,3,4,5,6]          A[m*i+j]   m = #columns 

    int i,j,k;
    LONG d,t;
    recint P;
    P = recip1(p);
    d = 1;
    for( k=0; k<n; k++ ) {
        //if( k>0 && k%100 == 0 ) printf("elimination at row %d\n",k);
        for( i=k; i<n && A[n*i+k]==0; i++ ); // look for non-zero pivot
        if( i>=n ) { d = 0; break; }
        if( i!=k ) { // interchange row k with row i
             for( j=k; j<n; j++ ) { t = A[k*n+j]; A[k*n+j] = A[i*n+j]; A[i*n+j] = t; }
             d = neg64s(d,p);
        };
        d = mulrec64(d,A[k*n+k],P);
        if( A[k*n+k]==0 ) { printf("division by 0\n"); return -1; }
        t = inv64s(A[k*n+k],p);
        A[n*k+k] = 1;
        for( j=k+1; j<n; j++ ) { A[n*k+j] = mulrec64(A[n*k+j],t,P); }
        for( i=k+1; i<n; i++ ) { // row i
            if( A[i*n+k]!=0 )
                for( j=k+1; j<n; j++ ) 
                     // A[i,j] = A[i,j] - A[i,k] A[k,j] mod p
                    // A[i*n+j] = sub64s(A[i*n+j],mulrec64(A[i*n+k],A[k*n+j],P), p );
                    if( A[k*n+j]!=0 )  A[i*n+j] = sub64s(A[i*n+j],mulrec64(A[i*n+k],A[k*n+j],P), p );
            A[i*n+k] = 0;
        }
    }
    return d;
}


int minor64s( LONG * A, int n, int m, int *rows, int *cols, LONG p ) {
// Input A[n,m], rows[n], cols[m] 
// Output  r = rank A and rows and cols s.t. rank(B)=r where B=A[rows[0..r-1],cols[0..r-1]]
// Example  A = [[1,1,1],[2,2,2],[1,2,3],[1,0,-1]]
// r = rank(A) = 2 rows = [1,3,-,-], cols = [1,2,-,-] ==> M = [[1,1],[1,2]]
    int i,j,k,r,c,rank;
    LONG t,m64;
//printf("n=%d  m=%d\n",n,m);
//printf("A=\n"); matprint64s( A, n, m );
    recint P;
    P = recip1(p);
    for( i=0; i<n; i++ ) rows[i] = i;
    rank = 0;
    m64 = m; // so the matrix can have more than 2^31 entries
    r = 0; // current row
    c = 0; // current column
    while ( r<n && c<m ) {
        for( i=r; i<n && A[rows[i]*m64+c]==0; i++ ); // look for non-zero pivot
        if( i>=n ) { c++; continue; } // move to next column
        if( i!=r ) { j = rows[i]; rows[i] = rows[r]; rows[r] = j; }; // row interchange
        cols[rank] = c;
        rank ++;
        k = rows[r]; // pivot row
        t = inv64s(A[k*m64+c],p);
        A[k*m64+c] = 1; // make the pivot row have a pivot = 1
        for( j=c+1; j<m; j++ ) A[k*m64+j] = mulrec64(A[k*m64+j],t,P);
        for( i=r+1; i<n; i++ ) {
            if( A[rows[i]*m64+c]==0 ) continue;
            for( j=c+1; j<m; j++ ) {
                if( A[k*m64+j]==0 ) continue;
                t = mulrec64(A[rows[i]*m64+c],A[k*m64+j],P);
                A[rows[i]*m64+j] = sub64s(A[rows[i]*m64+j],t,p);
            }
            A[rows[i]*m64+c] = 0;
        }
        r++; // next row
        c++; // next column
//printf("r=%d  c=%d\n", r, c);
//printf("A=\n"); matprint64s( A, n, m );
    }
    return rank;
}