//V35.0

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include <sys/time.h>

#define LONG long long int
#define ULONG unsigned long long int
#define ULONG unsigned long long int
#define UINT32 unsigned int
#define UINT64 unsigned long long

typedef struct {
	UINT64 s;	/* shift */
	UINT64 v;	/* reciprocal of d */
	UINT64 d0;	/* divisor shifted up */
	UINT64 d1;
} recint;

recint recip1(UINT64  p);

ULONG seed, mult;
LONG rand64s(LONG p);


//Custom Functions

//create an array of length n of 64 bit integers
LONG * array( int n );
LONG * array64s( LONG n );

//basic functions
LONG log2N(LONG n);
LONG ceil2(LONG n);

//fills an array with 0s
void cleanArray64s(LONG *A, int len);

//prints univariate and bivariate Polynomials in Maple orientation
void polprint64s( LONG *A, int d );
void polprintbivar64s( LONG *A, int dx, int dy, char * x, char * y);

//Compare solution output by the Hensel lifting functions to the actual answer
void ComparePolyFactorsBivar(LONG *A,LONG *F,int dx,int dz,int n,LONG alpha,LONG p);

//Generate polynomials to test the
void createProblem(LONG *A, LONG *F0, int n, int dx, int dz, int dxA, int dzA, LONG alpha, LONG p);


//This preforms hensel lifting on n factors to factor polynomial A
int HenselLiftCubic(LONG *A, int dx, int dz, LONG *F0, int n, int du, LONG *F, LONG alpha, LONG *TempSpace, LONG p );
int HenselLiftQuartic(LONG *A, int dx, int dz, LONG *F0, int n, int du, LONG *F, LONG alpha, LONG *W, LONG p );

//main
int main () {

    clock_t T1,T2;
    LONG p;

    int dx,dy,n,np,du,i,j,k,t;
    LONG *A,*F0,*F,*W,*disp,alpha,tempArraySize;
    recint P;

   p = 1009;
   //p = pow(2,31)-1;
   P = recip1(p);

        n = 4; //number of factors
        dx = 8; //deg(A,x);
        dy = 8; //deg(A,y);
        du = 2;  //max{degrees of each initial function in x}
        alpha = 3;

        //calculates np, the number of polynomial products
        np = 1; t = n;
        while (t>1){
           np +=  t;
           t = ceil2(t);
        }


        A = array((dx+1)*(dy+1));
        disp = array((dx+1)*(dy+1));
        F0 = array((du+1)*n);
        F = array((dx+1)*(dy+1)*n);
        tempArraySize = (n-1)*n*du/2 + n + (dx+1)*(dy+1) +5*n*(du+1) + 3*(dx+1) + (dx+1)*(dx+1) + (du+1)*(dx/2 + 1) + n*(dx+1) + 4*(dx+1)*(dy+n)*(log2N(n)+1) + 2*np*(dx+1) + 32*(dy+1) + (dx+1)*((dx+1)/2+1);
        W = array(tempArraySize);

        //printf("There are %d factors, each of which has degree(x) = %d and degree(y) = %d\n",n,du,du);
        printf("We wish to factor the polynomial A\n");

        cleanArray64s(A,(dx+1)*(dy+1));
        cleanArray64s(F,(dx+1)*(dy+1)*n);
        cleanArray64s(W,tempArraySize);
        cleanArray64s(F0,(du+1)*n);
        cleanArray64s(disp,(dx+1)*(dy+1));

        //Populate polynomial A

        A[0*(dx+1)+0] = 191; A[0*(dx+1)+1] = 979; A[0*(dx+1)+2] = 97; A[0*(dx+1)+3] = 376; A[0*(dx+1)+4] = 92; A[0*(dx+1)+5] = 370; A[0*(dx+1)+6] = 595; A[0*(dx+1)+7] = 237; A[0*(dx+1)+8] = 1;
        A[1*(dx+1)+0] = 947; A[1*(dx+1)+1] = 48; A[1*(dx+1)+2] = 990; A[1*(dx+1)+3] = 365; A[1*(dx+1)+4] = 373; A[1*(dx+1)+5] = 213; A[1*(dx+1)+6] = 358; A[1*(dx+1)+7] = 46; A[1*(dx+1)+8] = 0;
        A[2*(dx+1)+0] = 338; A[2*(dx+1)+1] = 806; A[2*(dx+1)+2] = 487; A[2*(dx+1)+3] = 957; A[2*(dx+1)+4] = 538; A[2*(dx+1)+5] = 78; A[2*(dx+1)+6] = 109; A[2*(dx+1)+7] = 739; A[2*(dx+1)+8] = 0;
        A[3*(dx+1)+0] = 622; A[3*(dx+1)+1] = 497; A[3*(dx+1)+2] = 517; A[3*(dx+1)+3] = 650; A[3*(dx+1)+4] = 143; A[3*(dx+1)+5] = 282; A[3*(dx+1)+6] = 582; A[3*(dx+1)+7] = 0; A[3*(dx+1)+8] = 0;
        A[4*(dx+1)+0] = 555; A[4*(dx+1)+1] = 26; A[4*(dx+1)+2] = 179; A[4*(dx+1)+3] = 911; A[4*(dx+1)+4] = 479; A[4*(dx+1)+5] = 79; A[4*(dx+1)+6] = 945; A[4*(dx+1)+7] = 0; A[4*(dx+1)+8] = 0;
        A[5*(dx+1)+0] = 142; A[5*(dx+1)+1] = 598; A[5*(dx+1)+2] = 999; A[5*(dx+1)+3] = 386; A[5*(dx+1)+4] = 540; A[5*(dx+1)+5] = 994; A[5*(dx+1)+6] = 0; A[5*(dx+1)+7] = 0; A[5*(dx+1)+8] = 0;
        A[6*(dx+1)+0] = 630; A[6*(dx+1)+1] = 158; A[6*(dx+1)+2] = 966; A[6*(dx+1)+3] = 320; A[6*(dx+1)+4] = 1002; A[6*(dx+1)+5] = 539; A[6*(dx+1)+6] = 0; A[6*(dx+1)+7] = 0; A[6*(dx+1)+8] = 0;
        A[7*(dx+1)+0] = 62; A[7*(dx+1)+1] = 599; A[7*(dx+1)+2] = 464; A[7*(dx+1)+3] = 404; A[7*(dx+1)+4] = 568; A[7*(dx+1)+5] = 0; A[7*(dx+1)+6] = 0; A[7*(dx+1)+7] = 0; A[7*(dx+1)+8] = 0;
        A[8*(dx+1)+0] = 154; A[8*(dx+1)+1] = 61; A[8*(dx+1)+2] = 592; A[8*(dx+1)+3] = 209; A[8*(dx+1)+4] = 16; A[8*(dx+1)+5] = 0; A[8*(dx+1)+6] = 0; A[8*(dx+1)+7] = 0; A[8*(dx+1)+8] = 0;

        //Populate the initial factors F0
        F0[0*(du+1)+0] = 19; F0[0*(du+1)+1] = 289; F0[0*(du+1)+2] = 1;
        F0[1*(du+1)+0] = 950; F0[1*(du+1)+1] = 35; F0[1*(du+1)+2] = 1;
        F0[2*(du+1)+0] = 660; F0[2*(du+1)+1] = 89; F0[2*(du+1)+2] = 1;
        F0[3*(du+1)+0] = 37; F0[3*(du+1)+1] = 559; F0[3*(du+1)+2] = 1;

        printf("A := ");polprintbivar64s(A,dx,dy,"x","y");
        printf("\n\n");
        for(i=0;i<n;i++){
           printf("f%d :=",i+1); polprint64s(&F0[i*(du+1)],du); printf(":");
           printf("\n");
        }
        printf("The total degree of A is degree(x) = %d and degree(z) = %d\n",dx,dy);



            T1 = clock();
            HenselLiftCubic(A, dx, dy, F0, n, du, F, alpha, W, p); //Our cubic Linear Hensel lifting algorithm
            //HenselLiftQuartic(A, dx, dy, F0, n, du, F, alpha, W, p); //Bernardin's Quartic Linear Hensel lifting algorithm
            T2 = clock();

        printf("Time for Cubic Algorithm = %8.2fms\n",(T2-T1)/1000.0/1.0);
        printf("\n");

        printf("After executing our algorithm, we get the factors:\n");

        for(i=0;i<n;i++){
             printf("f%d :=",i+1); polprintbivar64s(&F[i*(dx+1)*(dy+1)],dx,dy,"x","(y-3)");
            printf(";\n");
        }
        printf("\n");

        printf("A - f1*f2*f3*f4 = \n");
        ComparePolyFactorsBivar(A,F,dx,dy,n,alpha,p); //Checks if we factored correctly
        printf("\n");

        free(A);
        free(disp);
        free(F0);
        free(F);
        free(W);

    return 0;

}