// Compile with
// gcc -O3 -shared -o evalGP1.so -fPIC evalGP1.c
// Copyright Michael Monagan 2025

#include <stdio.h>
#include <stdlib.h>

#define LONG long long int
#define ULONG unsigned long long int


/******************************************************************************************/
/*  Zp utilities                                                                          */
/******************************************************************************************/

#define UINT64 unsigned long long

#include "int128g.c"

//typedef struct {
//        UINT64 s;       /* shift */
//        UINT64 v;       /* reciprocal of d */
//        UINT64 d0;      /* divisor shifted up */
//        UINT64 d1;
//} recint;
//
//recint recip1(UINT64  p);
//UINT64 mulrec64(UINT64  a, UINT64  b, recint  v);


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


LONG * array(LONG n) {
    LONG *A;
    A = (LONG *) malloc( n*sizeof(LONG) );
    if( A==0 ) { printf("out of memory\n"); exit(1); }
    return A;
}

int * array32s( LONG n ) {
    int *A;
    A = (int *) malloc( n*sizeof(int) );
    if( A==0 ) { printf("out of memory\n"); exit(1); }
    return A;
}


/*********************************************************************/
/*  These routines work with the Maple DAG directly                  */
/*********************************************************************/

#define NAME 8
#define TABLEREF 10
#define PROD 14
#define SUM 16
#define POLY 17
#define EXPSEQ 31
#define LIST 32
#define SET 38
#define ARRAY 40
#define INTNEG 1
#define INTPOS 2


int ID( LONG *a ) {
    int id;
    id = (a[0] >> 32) & 63;
    return id;
}

LONG LENGTH( LONG *a ) {
    LONG len;
    len = a[0] & ((1LL << 32)-1);
    return len;
}

int ISPOLY( LONG *a ) {
    return( ID(a)==POLY );
}

int findindex( LONG x, LONG *A, int n ) {
    int i;
    for( i=0; i<n; i++ ) if( x==A[i] ) return i;
    return n;
}


int getsumdegs64s( LONG *f, LONG *x, int *D, int n )  {
int i,j,k,t;
LONG c,*m,d,MAXDEG;
    // f is a SUM dag
    // extract degrees of f in all n variables simulateneously
    if( ID(f)!=SUM ) return 0; t = (LENGTH(f)-1)/2; f++;
    for( i=0; i<n; i++ ) D[i] = 0;
    MAXDEG = (1LL<<31)-1;
    for( i=0; i<t; i++ ) {
        c = f[2*i];
        if( c&1 ) continue;
        m = (LONG *) c;
        if( ID(m)==NAME || ID(m)==TABLEREF ) {
            j = findindex((LONG) m,x,n);
            if( j==n ) return 0;
            if( D[j]==0 ) D[j] = 1;
        } else if( ID(m)==PROD ) {
            for( k=1; k<LENGTH(m); k+=2 ) {
                j = findindex(m[k],x,n);
                if( j==n ) return 0;
                d = m[k+1];
                if( d&1 == 0 ) return 0; // d is not an immediate integer
                d = d/2;
                if( d<0 || d>MAXDEG ) return -1;
                if( d>D[j] ) D[j] = d;
            }
        } else return 0;
    }
    return 1;
}

int evalsumall( LONG *f, LONG *x, int n, LONG **POW, LONG p )  {
    int i,j,k,t;
    LONG c,y,z,*m,d;
    recint P;
    if( ID(f)!=SUM ) return -1;
    t = (LENGTH(f)-1)/2;
    f++;
    P = recip1(p);
    y = 0;
    for( i=0; i<t; i++ ) {
        if( f[2*i] & 1 ) {
            c = f[2*i];
            if( c&1 ) z = c/2; else return -2;
            y = add64s(y,z,p);
            continue;
        }
        c = f[2*i+1];
        if( c&1 ) z = c/2; else return -2;
        m = (LONG *) f[2*i];
        if( ID(m)==NAME || ID(m)==TABLEREF ) {
            j = findindex((LONG) m,x,n);
            if( j==n ) return -1;
            z = mulrec64(z,POW[j][1],P);
        } else if( ID(m)==PROD ) {
            for( k=1; k<LENGTH(m); k+=2 ) {
                j = findindex(m[k],x,n);
                if( j==n ) return -1;
                d = m[k+1];
                if( d&1 == 0 ) return -1; // d is not an immediate integer
                d = d/2;
                if( d<0 ) return -1;
                z = mulrec64(z,POW[j][d],P);
            }
        } else return -3;
        y = add64s(y,z,p);
    }
    return y;
}


int getpolydegs64s( LONG *f, int *D, int n )  {
// [  m1, c1, m2, c2 ]   m1 = [ deg(m1) | deg(m1,x1) | deg(m1,x1) ], c1, ...]
int i,j,t,numbits;
LONG mask,z,d,MAXDEG;
    if( ID(f)!=POLY ) return 0;
    t = (LENGTH(f)-2)/2;
    f += 2; // skip over header and variables
    // extract degrees of f in all n variables simulateneously
    numbits = 64/(n+1);
    mask = ( 1LL << numbits ) - 1;
    for( i=0; i<n; i++ ) D[i] = 0;
    MAXDEG = (1LL<<31)-1;
    for( i=0; i<t; i++ ) {
        z = f[2*i];
        for( j=n-1; j>=0; j-- ) {
            d = z & mask;
            if( d>MAXDEG ) return -1;
            if( d>D[j] ) D[j] = d;
            z = z >> numbits;
        }
    }
    return 1;
}

LONG evalpolyall(LONG *f, int n, LONG **POWERS, LONG p) {
// Let f = sum( a_i M_i(x1,...,xn), i=1..t ) stored in the POLY representation.
// Let POWERS[i][d] = alpha^d
// It is assumed f is in the POLY representation, n>1 and 0 <= a_i < p
    LONG i, j, t, x, y, z, numbits, mask, d, h;
    recint P;
    if( ID(f)!=POLY ) return -1;
    t = (LENGTH(f)-2)/2;
    f += 2;
    P = recip1(p);
    if( n==1 ) numbits=64; else numbits = 64/(n+1);
    mask = ((1LL) << numbits) - 1;
    y = 0;
    for( i=0; i<t; i++ ) {
        x = f[2*i]; // monomial
        z = f[2*i+1]; // coefficient
//printf("n=%d  x=%lld  z=%lld\n",n,x,z);
        if( z&1 ) { // Maple immediate integer
            z = z/2;
            if( z<=-p || z>=p ) z = z%p;
            if( z<0 ) z = z+p;
        } else {
            printf("coefficients must be immediate\n");
            exit(1);
        }
        for( j=n-1; j>=0; j-- ) {
            d = x & mask;
            if( d ) z = mulrec64(z,POWERS[j][d],P);
            x = x >> numbits;
        }
//printf("i=%d  z=%lld\n",i,z);
        y = add64s(y,z,p);
    }
    return y;
}


LONG evalpolytop(LONG *a, LONG *alpha, int n, LONG p) {
    int i,j,k,m,s,D[32];
    LONG y,*POW[32],*Y,*B,*W;
    recint P;

    if( ID(a)!=POLY ) { printf("a is not a POLY\n"); exit(1); }
    Y = (LONG *) a[1];
    if( ID(Y)!=EXPSEQ ) { printf("a is a bad POLY\n"); exit(1); }
    m = LENGTH(Y);
    if( m!=n+1 ) { printf("evalpolyall : variable missmatch\n"); exit(1); }

    getpolydegs64s(a,D,n);
    for( s=i=0; i<n; i++ ) s = s+D[i]+1;
    B = W = array(s);

    P = recip1(p);
    for( i=0; i<n; i++ ) {
        POW[i] = B;
        m = D[i]; //printf("deg(a,x[i])=%d\n",m);
        B[0] = 1;
        for( j=1; j<=m; j++ ) B[j] = mulrec64(alpha[i],B[j-1],P);
        B += m+1;
    }

    y = evalpolyall(a,n,POW,p);
    free(W);
    return y;
}

LONG evalsumtop(LONG *a, LONG *x, LONG *alpha, int n, LONG p) {
    int i,j,k,m,*D,s;
    LONG y,**POW,*B,*W;
    recint P;

    if( ID(a)!=SUM ) { printf("evalsumall: SUM expected\n"); exit(1); }
    D = array32s(n);
    m = getsumdegs64s(a,x,D,n);
    if( m==0 || m==-1 ) { free(D); return m; }

    POW = (LONG **) array(n);
    P = recip1(p);
    s = 0;
    for( i=0; i<n; i++ ) s = s + D[i] + 1;
    B = W = array(s);
    for( i=0; i<n; i++ ) {
        POW[i] = B;
        m = D[i];
        B[0] = 1;
        for( j=1; j<=m; j++ ) B[j] = mulrec64(alpha[i],B[j-1],P);
        B = B+m+1;
    }

    y = evalsumall(a,x,n,POW,p);
    free(POW); free(D); free(W);
    return y;
}

LONG eval64s(LONG *a, LONG *x, LONG *alpha, LONG p) {
// evaluate a at x[i] = alpha[i] mod p where x includes indets(a)
    int i,j,n;
    LONG z,b,*y;

    if( ID(x)!=LIST) { printf("x must be a list\n"); exit(1); }
    if( ID(alpha)!=LIST) { printf("alpha must be a list\n"); exit(1); }
    x = (LONG *) x[1];
    alpha = (LONG *) alpha[1];
    //printf("check %d %d\n",LENGTH(x),LENGTH(alpha));
    //if( LENGTH(x)!=LENGTH(alpha) ) { printf("x and alpha inconsisent\n"); exit(1); }
    // reduce alpha mod p from immediate
    n = LENGTH(x)-1;
    x++; alpha++;
    //if( ID(alpha)==LIST){}
    //if(id==LIST) {alpha++;}
    //printf("a=%d %d %d %d %d\n",alpha[0],alpha[1],alpha[2], alpha[3], alpha[4]);
    if( ID(a)==POLY ) { LONG beta[32];
        //printf("a=%lld %lld %lld %lld %lld\n",alpha[0],alpha[1],alpha[2], alpha[3], alpha[4]);
        y = (LONG *) a[1];
        if( ID(y)!=EXPSEQ ){ return -1;}
        //printf("%d\n",LENGTH(y));
        if( LENGTH(y)!=n+1 ) { return -1;}
        y++;
        for( i=0; i<n; i++ ) {
            j = findindex(x[i],y,n);
            if( j==n ) { printf("variable mismatch\n"); exit(1); }
            b = alpha[i];
            if( !(b&1) ) { printf("alpha[i] is not immediate\n"); exit(1); }
            b = (b/2) % p;
            if( b<0 ) b += p;
            beta[j] = b;
        }
        z = evalpolytop(a,beta,n,p);
        return z;
    }
    if( ID(a)==SUM ) { LONG *beta;
        beta = array(n);
        for( i=0; i<n; i++ ) {
            b = alpha[i];
            if( !(b&1) ) { printf("alpha[i] is not immediate\n"); exit(1); }
            b = (b/2) % p;
            if( b<0 ) b += p;
            beta[i] = b;
        }
        z = evalsumtop(a,x,beta,n,p);
        free(beta);
        return z;
    }
    if( ID(a)==NAME || ID(a)==TABLEREF ) {
        for( i=0; i<n; i++ )
            if( a==(LONG *)x[i] ) {
                b = alpha[i];
                if( b&1 ) {
                    b = (b/2) % p;
                    if( b<0 ) b += p;
                    return b;
                }
                printf("alpha is not immediate\n"); exit(1);
            }
    }
    printf("invalid type for f\n");
    exit(1);
}