/*****************************************************************/
/* A SAT encoding of Erdos's discrepancy problem                 */
/*   based on Sinz's sequential counter                          */
/* Copyright (C) 2014 by Boris Konev & Liverpool University      */
/*****************************************************************/

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


#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))

long int discrepancy, length;
long int verbose = 0;
long int useCompleteMultiplicativity = 0;
long int useMultiplicativity = 0;
long int useOne = 0;



// return the largest divisor of n, which does not exceed sqrt(n)
// this is convenient as for n=m^2 split(n) is m
long int split(long int n) {
    long int i = sqrt(n);
    for (; ; i--) {
        if (n%i == 0) {
            return i;
        }
    }
}




void printMultiplicativityRulesAux(long int m, long int n, long int limit) {
    if (m*n > limit) return;
    // else
    if(n>1) { 
        printf("X%ld X%ld X%ld \n",   m, n, m*n);
        printf("~X%ld ~X%ld X%ld \n", m, n, m*n);
        printf("~X%ld X%ld ~X%ld \n", m, n, m*n);
        printf("X%ld ~X%ld ~X%ld \n", m, n, m*n);
    }
    printMultiplicativityRulesAux(2*m-n,m, limit);
    printMultiplicativityRulesAux(2*m+n,m, limit);
    printMultiplicativityRulesAux(m+2*n, n, limit);
}

void printMultiplicativityRules(long int n) {
    printf("~ multiplicativity\n");
    printMultiplicativityRulesAux(2,1, n);
    printMultiplicativityRulesAux(3,1, n);
}

void printCompleteMultiplicativityRules(long int n) {
    printf("~ complete multiplicativity\n");
    printf("X1 \n");
    long int i = 4;
    for( ; i<=n; i++) {
        long int j = split(i);
        //printf("i = %ld, j = %ld\n", i, j);
        if (j==1) continue;
        long int k = i/j;
        if (j!=k) {
            printf("X%ld X%ld X%ld \n", j, k, i);
            printf("~X%ld ~X%ld X%ld \n", j, k, i);
            printf("~X%ld X%ld ~X%ld \n", j, k, i);
            printf("X%ld ~X%ld ~X%ld \n", j, k, i);
        }
        else {
            printf("X%ld \n",  i);
        }
    }
}
/* An alternative implementation of */
/* multiplicativity. Can run faster but doesn't reduce*/
/* the unsatisfiability certificate in size */
/*
int deSquare(int n) {
    int i=2;
    int isq, nsqrt;

    nsqrt = sqrt(n);
    while (i <= nsqrt)  {
        isq = i*i;
        if (n%isq == 0) {
            n = n/isq;
            nsqrt = sqrt(n);
            // not moving i in case if n is divisable by i^4.
        }
        else { 
            i++;
        }
    }
    return n;
}

void printCompleteMultiplicativityRules(int n) {
    printf("~ multiplicativity condition\n");
    printf("X1 \n");
    int i = 4;
    for(i; i<=n; i++) {
        int iprim = deSquare(i);
        //printf("i = %ld, iprim = %ld\n", i, iprim);
        if(iprim == 1) { // i is a square number
            printf("X%ld \n",  i); // must be 1
        }
        else if(iprim < i) {
            printf("~X%ld X%ld \n", i, iprim);
            printf("X%ld ~X%ld \n", i, iprim);
        }
        else {
            int j = split(i);
            //printf("i = %ld, j = %ld\n", i, j);
            if (j==1) {
                //printf("%ld is prime\n", i);
                continue;
            }
            int k = i/j;
            if (j!=k) {
                printf("X%ld X%ld X%ld \n", j, k, i);
                printf("~X%ld ~X%ld X%ld \n", j, k, i);
                printf("~X%ld X%ld ~X%ld \n", j, k, i);
                printf("X%ld ~X%ld ~X%ld \n", j, k, i);
            }
        }
    }
}
*/

// returns 1 if s^k_j is always true and 0 otherwise
int isTrue (long int k, long int j) {
    //printf(" (isTrue %ld, %ld..", k,j);
    if((j>2*k-2+discrepancy) ||(k==0)) {
        //printf("returning 1) ");
        return 1;
    }
    else {
        //printf("returning 0) ");
        return 0;
    }
}

// returns 1 if s^k_j is always false and 0 otherwise
int isFalse (long int k, long int j) {
    //printf(" (isFalse %ld, %ld..", k,j);
    if((j<2*k-discrepancy) || ((k>0)&&(j<k))) {
        //printf("returning 1) ");
        return 1;
    }
    else {
        //printf("returning 0) ");
        return 0;
    }
}

// prints a CNF formula enforcing that the 
// sequence x_d, x_{2d},..., x_{l/d} is C-bounded
void cBounded (long int d, long int l) {
    long int n = l / d;

    // a little optimisation:
    // notice that the sum of n elements is even when n is even
    // and odd when n is odd
    //
    // Thus, we can shorten the sequence by 1 if n is odd and C is even
    // or vice versa
    //
    if(n%2 == 0) {
        // n is even
        // check if  discrepancy is odd
        if(discrepancy % 2 == 0) 
            n--;
    }
    else {
        // n is odd
        // check if discrepancy is even
        if(discrepancy % 2 != 0) 
            n--;
    }

    long int k,j;
    // now, print clauses
    for (k = 1; k <= n; k++) {
        //for(j = 0; j <= n; 
        for(j = MAX(1, 2*k-discrepancy-1);  // a conservative estimate
                j <= MIN(n, 2*k-2+discrepancy+1); 
                    j++) {
                //        printf("k=%ld, j=%ld\n",k,j);

            // s^{k}_{j}  <-->  s^{k}_{j-1} \/ (s^{k-1}_{j-1} /\ x_{j})
            //
            // RHS -> LHS
            // converts to
            //
            // s^{k}_{j-1}  --> s^{k}_{j}    
            // s^{k-1}_{j-1} /\ x_{j} --> s^{k}_{j}    
            //
            //   in CNF:
            //   -s^{k}_{j-1}  \/ s^{k}_{j}    
            if(!isFalse(k,j-1) && !isTrue(k,j)) {
                if(verbose) printf("1  k=%ld, j=%ld ", k, j);
                if(!isTrue(k,j-1)) {
                    printf("~%ldS%ld.%ld ", d, k, j-1);
                }
                if(!isFalse(k,j)) {
                    printf("%ldS%ld.%ld", d, k, j);
                }
                printf("\n");
            }
            //   -s^{k-1}_{j-1} \/ -x_{j} \/ s^{k}_{j}    
            if(!isFalse(k-1,j-1) && !isTrue(k,j)) {
                if(verbose) printf("2  k=%ld, j=%ld ", k, j);
                if(!isTrue(k-1,j-1)) {
                    printf("~%ldS%ld.%ld ", d, k-1, j-1);
                }
                //
                printf("~X%ld ", d*(j));
                //
                if(!isFalse(k,j)) {
                    printf("%ldS%ld.%ld", d, k, j);
                }
                printf("\n");
            }
            // the other direction (LHS->RHS)
            // s^{k}_{j}  -->  s^{k}_{j-1} \/ s^{k-1}_{j-1}
            // s^{k}_{j}  -->  s^{k}_{j-1} \/ x_{j}
            //
            //   in CNF:
            //   -s^{k}_{j}  \/  s^{k}_{j-1} \/ s^{k-1}_{j-1}
            if(!isFalse(k,j) && !isTrue(k,j-1) && !isTrue(k-1,j-1)) {
                if(verbose) printf("3  k=%ld, j=%ld ", k, j);
                if(!isTrue(k,j)) {
                    printf("~%ldS%ld.%ld ", d, k, j);
                }
                if(!isFalse(k,j-1)) {
                    printf("%ldS%ld.%ld ", d, k, j-1);
                }
                if(!isFalse(k-1,j-1)) {
                    printf("%ldS%ld.%ld ", d, k-1, j-1);
                }
                printf("\n");
            }
            //   -s^{k}_{j}  \/  s^{k}_{j-1} \/ x_{j}
            if(!isFalse(k,j) && !isTrue(k,j-1)) {
                if(verbose) printf("4  k=%ld, j=%ld ", k, j);
                if(!isTrue(k,j)) {
                    printf("~%ldS%ld.%ld ", d, k, j);
                }
                if(!isFalse(k,j-1)) {
                    printf("%ldS%ld.%ld ", d, k, j-1);
                }
                printf("X%ld\n", d*(j));
            }
            //printf("end\n");
        }
    }
}

// Print the usage message
void usage(char* s) {
        fprintf(stderr,"Usage: %s [-m] [-c] [-X num] [-noX num] discrepancy length\n",s);
        printf("\twhere\t-m adds multiplicativity restrictions\n");
        printf("\t\t-c adds complete multiplicativity restrictions\n");
        printf("\t\t-X n sets n\'s x to be +1\n");
        printf("\t\t-noX n sets n\'s x to be -1\n");
        printf("\t\t-o just the sequence with difference one is required to be cbounded\n");
}

// The entry point
int main(int argc, char* argv[])
{
    long int a = 1; // position in argv[]
    if(argc<3) {
        usage(argv[0]); // at least 2 args should be provided: C and n
        exit(-1);
     }

    // processing flags
    while(strncmp(argv[a],"-",1)==0) {
        if(strncmp(argv[a],"-v",2)==0) { 
            verbose = 1;
            a++;
        }
        if(strncmp(argv[a],"-c",2)==0) { 
            useCompleteMultiplicativity = 1;
            a++;
        }
        if(strncmp(argv[a],"-m",2)==0) { 
            useMultiplicativity = 1;
            a++;
        }
        if(strncmp(argv[a],"-o",2)==0) { 
            useOne = 1;
            a++;
        }
        // sets the value of a particular X to true
        if(strncmp(argv[a],"-X",2)==0) { 
            a++;
            long int x;
            if(sscanf(argv[a],"%ld",&x)!=1) {
                usage(argv[0]);
                exit(-1);
            }
            printf("X%ld\n", x);
            a++;
        }
        // sets the value of a particular X to false
        if(strncmp(argv[a],"-noX",2)==0) { 
            a++;
            long int x;
            if(sscanf(argv[a],"%ld",&x)!=1) {
                usage(argv[0]);
                exit(-1);
            }
            printf("~X%ld\n", x);
            a++;
        }
    }
    // scanning for the discrepancy and length
    //
    if(sscanf(argv[a], "%ld", &discrepancy)!=1||
            sscanf(argv[a+1], "%ld", &length)!=1) {
        usage(argv[0]);
        exit(-1);
     }

    // any string starting with "~ " (followed by space) is a comment
    printf("~ EDP(%ld, %ld)\n",   discrepancy, length);

    if(useCompleteMultiplicativity) {
        //printCompleteMultiplicativityRules(127600);
        printCompleteMultiplicativityRules(length);
    }

    if(useMultiplicativity) {
        //printMultiplicativityRules(127600);
        printMultiplicativityRules(length);
    }

    printf("~ discrepancy\n");
    long int k;
    if(useOne) {
        // to be used only in conjunction with -c
        // as for completely multiplicative sequences cboundednes of the
        // "main" sequence implies cboundedness of all subsequences
        // with indices forming a homogeneous progression
        // This does not lead to any noticible reductions in the
        // unsatisfiability certificate size, though..
        cBounded(1, length);
    } else {
        for (k = 1; (discrepancy+1)*k <= length; k++) {
            cBounded (k, length);
        }
    }
}