#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <time.h>

/////////////////////////////////////////////////////
//poly.c
/////////////////////////////////////////////////////
//Version 1.0 (3/7/2011)
/////////////////////////////////////////////////////

/////////////////////////////////////////////////////
//To compile: gcc -lm -O3 poly.c
//(-lm to include math.h and complex.h)
//(-O3 for maximum optimization)
//
//Them, Run the built file with arguments "n d"
//
//The output will be a file in the root directory
//of the form ndpq_DD.MM.YY.txt
////////////////////////////////////////////////////

//Groups of polynomials will be referred to as:
//(n, d, P, Q), where
//n is the number of variables,
//d is the maximum degree
//P is the polymod
//Q is the funcmod.

/* This code will work correctly for the following values of n, d, P, Q:

(note: d must be a positive integer)

P = 2; Q = 3                  P = 3; Q = 2
-------------                 -------------
n = 2:  d <= 2                n = 2: d <= 2
n = 3:  d <= 3                n = 3: d <= 3
n = 4:  d <= 4                n = 4: d <= 4
n = 5:  d <= 5                n = 5: d <= 5
n = 6:  d <= 5                n = 6: d <= 2
n = 7:  d <= 2                n = 7: d <= 2
n = 8:  d <= 2                n = 8: d <= 2
n = 9:  d <= 2
n = 10: d <= 2

For the above, the number of terms (monomials) will not exceed 63,
and the number of polynomials will not exceed 2^63.

*/

/////////////Constants/////////////

//polymod p
#define P 2

//funcmod q
#define Q 3

//parameter for the evalTable, that will precompute
//evaluation of polynomials in chunks.
#define T 9

//The evalTable will be of size TABLE_X * TABLE_Y:
#define TABLE_X (int)pow(P, T)
#define TABLE_Y (2 << T)

//Pointer to output file.
FILE *fp;

//full name of file (with date)
char name[20];

//Array used for computing order the order of terms.
//Generates all terms in the extended input.
unsigned int allTerms[65][6]; //65: upper bound on # of terms.
                              //6: upper bound on max degree
//Index for allTerms array.
unsigned int C = 0;

//////////////////////////////////

/* x^n (for long longs)*/
unsigned long long ipow(unsigned long long x, unsigned int n){
   unsigned long long ret = 1;
    while (n > 0){
        if (n % 2){
            ret *= x;
            n -= 1;
        }
        x *= x;
        n /= 2;
    }
    return ret;
}

/* Count the number of set bits in the binary representation
 * of i. Algorithm source:
 * http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
 */
unsigned int countSetBits(unsigned int i){
    i = i - ((i >> 1) & 0x55555555);
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
    return ((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) >> 24;
}

/* Binomial Coefficient: (n choose m)*/
unsigned int binom(unsigned int n, unsigned int m){
    if ((m > n) || (m < 0))
        return 0;

    unsigned int i;
    unsigned int retval = 1;
    for(i = 0; i < m; i++){
        retval *= (n - i);
        retval /= (i + 1);
    }
    return retval;
}

/* Returns the number of possible terms in a
 * polynomial of arity n and max degree d */
unsigned int countTerms(unsigned int n, unsigned int d){
    unsigned int retval = 0;
    unsigned int i;
    for (i = 0; i <= d; i++){
        retval += binom(n, i);
    }
    return retval;
}

/* Returns the number of possible polynomials
 * under the given modulus, with the given
 * arity and max degree */
unsigned long long countPolys(unsigned int termCount){
    return ipow(P, termCount);
}

/* Given (integer representations of) a polynomial
 * and a term, returns the coefficient for that term
 * in the given modulus */
unsigned int termCoeff(unsigned long long poly, unsigned int term){
    return (poly / ipow(P, term)) % P;
}

/* Evaluates the given term on the given input.
 * For {0, 1}:
 * Returns 1 when there is no 1 bit in the term
 * that is in the same position as a 0 bit in the input.
 * In other words: returns 1 if the bitwise logical
 * implication [term => inpt] is all 1's.
 * For {-1, 1}
 * Keeps track of parity by reversing it for every 1 bit
 * in the input that matches the corresponding bit in
 * the term.*/
 unsigned int evalTerm(unsigned long long inpt, unsigned int term){
    //{0,1}: AND
    #if P == 2
        return (term & ~inpt) ? 0 : 1;
    //{-1, 1}: Parity
    #else
        unsigned int r = 0;
        while (0 < countSetBits(term)){
             if ((term & 1) & (inpt & 1)){
                r = !r;
                }
            term >>= 1;
            inpt >>= 1;
        }
        return r;
    #endif
}

/*Find all terms for the given n and d, reusing the given array. Loop
 *through all possible monomials of order d starting with variable x_n.
 *Record the values inside of term, copy term into the allTerms array,
 *increment allTerms' index, and recur.*/
void getDegTerms(unsigned int *term, unsigned int index, unsigned int rest_var, unsigned int first_var){
    unsigned int i;
    for(i = 0; first_var - i >= rest_var; i++){
        term[index] = first_var - i;
        if(rest_var > 1)
            getDegTerms(term, index + 1, rest_var - 1, first_var - i - 1);
        else{
            memcpy(allTerms[C], term, (sizeof(int) * index + 1));
            C++;
        }
    }
}

/*Fills in the allTerms array with all terms for polynomials
 *of the given n and d. Will produce these arrays in order
 *of largest degree to the smallest degree, which will be 1.
 *For example: n = d = 3 produces: ((321) (32) (31) (3) (2) (1))*/
void getAllTerms(unsigned int n, unsigned int d){
    //Temporary array to store each term.
    unsigned int term[d];
    unsigned int i;
    for(i = d; i > 0; i--){
        //get all terms for the given degree
        getDegTerms(term, 0, i, n);
    }
}

/*Use the allTerms array to compute an extended input
 *for polynomials of n variables and degree d. This will
 *be done in order of smallest degree to largest.*/
unsigned long long extendInput(unsigned int n, unsigned long long inpt, unsigned int d){
    unsigned long long ret = 0;
    //count the total number of terms excluding the constant.
    unsigned int numTerms = countTerms(n, d) - 2;
    unsigned int i, j;

    //Start with linear monomials.
    unsigned int deg = 1;
    unsigned int count = 0;
    while (deg <= d){
        //There will be (n chooes deg) monomials of this degree.
        for (i = 0; i < binom(n, deg); i++){
            ret <<= 1;
            unsigned int c = 0;
            //each x_j represents the jth bit on the extended input. Thus,
            //raise 2^j-1 for each j in this term.
            for (j = 0; j < deg; j++){
                c += (1 << allTerms[numTerms-count][j]-1);
            }
            ret += evalTerm(inpt, c);
            count++;
        }
        deg++;
    }
   return ret;
}

/* Fill in the array of extended versions of
 * every possible set of inputs */
void allExtendedInputs(unsigned int n, unsigned int d, unsigned int termCount, unsigned long long *extdInpts){
    unsigned long long numInputs = (1 << n);
    unsigned long long inpt;
    for(inpt=0; inpt<numInputs; inpt++){
        extdInpts[inpt] = extendInput(n, inpt, d);
    }
}

/*print the polynomial in a readable format.
 *ie: x_2x_1 + x_2 + x_1
 *Uses the allTerms array for reference on which
 *variables to print for each term.*/
void printPoly(unsigned long long poly, unsigned int n, unsigned int d){
    unsigned int terms = countTerms(n, d);
    unsigned int i, j;
    unsigned int deg = d;
    //to avoid printing unnecessary '+' signs, the first term of a
    //polynomial will be tracked. When the first term is printed,
    //first is set to 0, and further terms will be paired with a '+'
    unsigned int first = 1;
    unsigned int count = 0; //count when all terms have bene pritned.

    while(count < terms){
        //print all terms of degree deg and arity n.
        for (i = 0; i < binom(n, deg); i++){
            unsigned int c = termCoeff(poly, count);
            if (c != 0){
                //check if first term was printed yet.
                if (first){
                    first = 0;
                }
                else{
                    fprintf(fp,  " + ");
                }
                //check if this term has a coeff. other than
                //1. If so, print it.
                if (c != 1)
                    fprintf(fp, "%d",c);

                //print the term first term.
                for (j = 0; j < deg; j++){
                    fprintf(fp,  "x_%d", allTerms[count][j]);
                }
            }
                count++;
        }
        deg--;
    }
    fprintf(fp, "\n");
}


/* Create the table that will contain values for the evaluation
 * of all polynomials [0, TABLE_X - 1] under inputs [0, TABLE_Y - 1].
 * The result is a TABLE_X by TABLE_Y lookup table, that can be used
 * for evaluating polynomials in chunks of T. These constants can be
 * set at the beginning of the code.*/
void fillEvalTable(unsigned int table[TABLE_X][TABLE_Y]){
    unsigned int h, v;

    //fill in first values for polynomials [0, P - 1]
    for (h = 0; h < P; h++){
        for (v = 0; v < TABLE_Y; v++)
            table[h][v] = (h % P) * (v & 1);
    }

    //fill in the rest of the polynomials, using summation with
    //already filled in values.
    for (h = P; h < TABLE_X; h++){
        for (v = 0; v < TABLE_Y; v++){
            table[h][v] = (h % P) * (v & 1) + table[h / P][v >> 1];
        }
    }
}

/* Given the results of the partials table, computer
 * the V-Value of a polynomial. This value can be
 * directly converted to the S-value, without
 * compromising precision and roundoff errors.
 * The V-Value will always be an integer.*/
__inline unsigned long long computeV(int a, int b, int c){
    int t = a;

    if (abs(b) < t)
        t = b;
    if (abs(c) < t)
        t = c;

    a = a - t;
    b = b - t;
    c = c - t;

    if (a == 0)
        return (b + c) * (b + c) - (3 * (b * c));
    else if (b == 0)
        return (a + c) * (a + c) - (3 * (a * c));
    else
        return (a + b) * (a + b) - (3 * (a * b));
}

/* Compute the S-Value (V-Value) of all polynomials and returns
 * the maximum value found. Prints the polynomials that return
 * this maximum in a file of the form  "ndpq_DD.MM.YY.txt" */
unsigned long long evalAllPolys(unsigned int n, unsigned int d){
    unsigned int i, j, k;

    //number of terms for polynomials with arity n, degree d.
    unsigned int termCount = countTerms(n, d);
    //total number of polynomials with arity n, degree d.
    //divide by P because we do not need to look at constants, which
    //will be the final monomial.
    unsigned long long polyCount = countPolys(termCount) / P;

    //count number of inputs we can have
    unsigned long long numInpts = (1 << n);

    //fill in allTerms array with all terms.
    getAllTerms(n, d);
    ///////////////////////////////////////////////////////////////
    //precompute number of 1 bits in each input to use later.
    //record number of 1 bits for each input in inptBits
    unsigned int inptBits[numInpts];

    unsigned long long l;

    for (l = 0; l < numInpts; l++){
        inptBits[l] = countSetBits(l);
    }

    unsigned long long inpt;
    ///////////////////////////////////////////////////////////////
    //Pre-compute the extended versions of all possible inputs
    unsigned long long extdInpts[numInpts];
    allExtendedInputs(n, d, termCount, extdInpts);

    //We will store results for each polynomial in the partials table.
    static unsigned int partials[Q][P];

    //create the lookup table of evaluated polynomials % P^T
    static unsigned int evalTable[TABLE_X][TABLE_Y];
    fillEvalTable(evalTable);

    ////////////////////////////////////////////////////////////////
    //Compute the Base polynomials: Polynomials of arity n and
    //degree d that do not contain any linear terms (1 > degree >= d).
    unsigned long long base = 0;

    for (i = 2; i < d+1; i++){
        base += binom(n, i);
    }

    base = ipow(P, base);
    ////////////////////////////////////////////////////////////////
    //Find the number of sums of linear monomials we can make to
    //base polynomials (see optimization). Record in monoms.
    unsigned int monoms;

    #if P == 2
    monoms = n+1;
    #else
    //P == 3
    monoms = (n + 1) * (n + 2) / 2;
    #endif
    ////////////////////////////////////////////////////////////////
    //Create an array to record all sums to add to base polynomials
    //necessary to check all polynomials disregarding permutations.
    unsigned long long linsums[monoms];
    //index for linsums.
    unsigned int count = 0;

    //Precalculate the sums for adding linear terms to a base poly.
    #if P == 2
    unsigned long long linearSum = 0;

    for (i = 0; i < n+1; i++){
        if (i != 0){
            //Add a linear term
            linearSum = linearSum + ipow(P, i-1) * base;
        }
        linsums[count] = linearSum;
        count++;
     }
    #else
    //P == 3

    //Addition to a base polynomial to sum a set of linear terms.
    unsigned long long linearSum = 0;
    //Addition to a base polynomial to sum a set of linear terms,
    //including possible 2 coefficients on these linear terms.
    unsigned long long linearSumCoeff = 0;

    //The power of P to get to the first linear term.
    unsigned int startPower = termCount - n - 2;

    //Find all sums of linear terms without 2 coefficeints.
    for (i = 0; i < n+1; i++){
            if (i != 0)
                linearSum += ipow(P, i-1) * base;
                linearSumCoeff = linearSum;
            //add in a 2 coefficients as appropriate (see optimization).
            for (j = 0; j < i + 1; j++){
                if (j > 0)
                    linearSumCoeff += ipow(P, startPower + j);

                linsums[count] = linearSumCoeff;
                count++;
            }
     }
    #endif

    unsigned long long maxV = 0; //Maximum V-Value seen so far.
    unsigned int counter = 0; //number of polynomials that have the Maximum

    ////////////////////////////////////////////////////////////////
    //Main Loop. Go through all base polynomials, go through all
    //sums of linear monomials, record results.
    unsigned long long poly = 0;
    unsigned long long polyWithSums = 0;

    for(poly = 0; poly < base; poly++){

        //EV progress check. Print when last x bits all 1
        if ( !(poly & ((1 << 25) - 1) ) )
          printf("poly=%llu, base=%llu;", poly, base);

        polyWithSums = poly;

        //add the ith sum in the linsums array to add
        //a set of linear terms.
        for (i = 0; i < monoms; i++){
            polyWithSums = poly + linsums[i];

            //reset partials array.
            for(j = 0; j < Q; j++){
                for(k = 0; k < P; k++){
                    partials[j][k] = 0;
                }
            }
            //Evaluate this polynomial for all inputs. Fill in
            //partials table.
            for(inpt = 0; inpt < numInpts; inpt++){
                //retrieve the resulting extended input.
                unsigned long long extInpt = extdInpts[inpt];
                unsigned long long polyTemp = polyWithSums;
                //evaluated value of the polynomial for this input.
                unsigned int polyVal = 0;
                do {
                    //look up the value of the chunk in the evalTable
                    polyVal += evalTable[polyTemp % TABLE_X][extInpt & (TABLE_Y-1)];
                    //divide by chunk size to find next polynomial.
                    #if P == 2
                    polyTemp >>= T;
                    #else
                    polyTemp /= TABLE_X;
                    #endif
                    //shift input by T
                    extInpt >>= T;
                }while (polyTemp > 0);

                //Fill in partials table
                partials[inptBits[inpt] % Q][polyVal % P] += 1;
            }
            ////////////////////////////////////////////////////
            //Find V(a,b,c)
            int a, b, c;
            unsigned long long vVal;
            #if P == 3
                a = partials[0][0] - partials[1][0];
                b = partials[0][1] - partials[1][1];
                c = partials[0][2] - partials[1][2];
                vVal = computeV(a,b,c);
            #else
                a = partials[0][0] - partials[0][1];
                b = partials[1][0] - partials[1][1];
                c = partials[2][0] - partials[2][1];
                vVal = computeV(a,b,c);
            #endif

            ////////////////////////////////////////////////////
            if (vVal > maxV){ //new Maximum
                //rewrite the file, so that only the latest Max is printed.
                fclose(fp);
                fp = fopen(name,"w");
                //New Max vVal
                maxV = vVal;
                fprintf(fp, "\n----------------MAX V-VAL: %llu----------------\n\n", maxV);
                //print the polynomial.
                printPoly(polyWithSums, n, d);
                //reset counter to one.
                counter = 1;
            }
             //add to list of polynomials under current max.
            else if (vVal == maxV){
                counter++;
                fprintf(fp, "\n");
                printPoly(polyWithSums, n, d);
            }
        }
    }
    fprintf(fp,  "\nThere were %d instances of the max V-Val\n\n", counter);
    return maxV;
}

/* Set the name of this file*/
void setFileName(int n, int d){
    time_t now;
    struct tm *today;
    char date[10];
    //get current date
    time(&now);
    today = localtime(&now);
    //print it in DD.MM.YY format.
    strftime(date, 15, "%d.%m.%Y", today);

    sprintf(name, "%d%d%d%d_%s.txt",n,d,P,Q,date);
}

int main(int argc, char *argv[]){
    ////////////////////////////////////////////////
    //Set up n and d
    unsigned int n, d;

    if (argc == 3){
        n = atoi(argv[1]);
        d = atoi(argv[2]);
    }
    else{
        printf("Please give values for n and d!\n");
        return 1;
    }
    ////////////////////////////////////////////////

    //set name for the output file
    setFileName(n, d);
    fp = fopen(name,"w");

    printf("I am doing:\nn = %d, d = %d, mod = %d, fmod = %d\n", n, d, P, Q);

    //set up timer to keep track of how long computation takes.
    clock_t start, end;
    double cpu_time_used;
    start = clock();

    //do the computation!
    unsigned long long maxV = evalAllPolys(n, d);

    //stop timer.
    end = clock();
    cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

    //print out maximum V-Value, and time spent computing.
    fprintf(fp, "maximum computed V-value: %llu\n", maxV);
    fprintf(fp, "time computing: %f\n", cpu_time_used);

    //get and print the date and time computation finished.
    time_t now;
    struct tm *today;
    char date[30];
    time(&now);
    today = localtime(&now);
    strftime(date, 30, "%c", today);
    fprintf(fp, "\nFinished on: %s\n", date);

    //print the n,d,P,Q values once more.
    fprintf(fp, "\nn = %d, d = %d, P = %d, Q = %d\n\n", n, d, P, Q);

    return 0;
}