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

#define AJM(x,y) ajm[x][y] = ajm[y][x]

typedef enum { kNone, kStar, kMax, kBranch } Mode;

static void xit(int rc, const char *msg, ...) {
    va_list v;
    va_start(v, msg);
    vfprintf(rc ? stderr : stdout, msg, v);
    exit(rc);
    va_end(v);
}

static int ij2var(int i, int j, int vert) {
    int ma = i > j ? i : j;
    int mi = i < j ? i : j;
    int no = (2*vert - 1 - mi)*mi/2 + ma - mi;
    return no;
}

static int product(int **ajm, int *sum, int i, int j, int m, int n) {
    if (i == j || m == n)
        return 0;
    if ((i == m && j == n) || (i == n && j == m)) {
        if (ajm[i][j] < 0) {
            printf("+1 x%d ",-ajm[i][j]);
            return 1;
        } else {
            *sum += ajm[i][j];
            return 0;
        }
    }
    if (ajm[i][j] < 0) {
        if (ajm[m][n] < 0) {
            printf("+1 x%d x%d ", -ajm[i][j], -ajm[m][n]);
            return 1;
        } else {
            if (ajm[m][n]) {
                printf("+1 x%d ", -ajm[i][j]);
                return 1;
            }
        }
    } else {
        if (ajm[m][n] < 0) {
            if (ajm[i][j]) {
                printf("+1 x%d ", -ajm[m][n]);
                return 1;
            }
        } else {
            *sum += ajm[i][j]*ajm[m][n];
        }
    }
    return 0;
}

static int fill_matrix(int **ajm, int verts, int valence, int adjac, int ndjac, Mode mode, int branch) {
    int i, j, cnt;
    for(i=0; i<verts; i++) {
        ajm[i][i] = 0;
    }
    if (mode == kNone)
        goto count;

    for(i=1; i<verts; i++) {
        ajm[0][i] = ajm[i][0] = i <= valence;
    }
    if (mode == kStar)
        goto count;

    if (verts == 99 && valence == 14 && adjac == 1 && ndjac == 2) {
        int k;
 
        for (k=15, i=1; i <= 13; i += 2) {
           AJM(i, i+1) = 1;
           for (j=i+2; j <= 14; j++) {
               if (k >= 98)
                   xit(4, "Miscalculated_1 k = %d\n", k);

               AJM(i, k) = 1;
               AJM(k, j) = 1;
               k++;
               AJM(i+1, k) = 1;
               AJM(k, j) = 1;
               k++;
           }
        }
        if (k != 99)
            xit(4, "Miscalculated_2 k = %d\n", k);
        for (i=0; i<=14; i++) {
            for (j=0; j<99; j++) {
                if (ajm[i][j] < 0)
                    AJM(i,j) = 0;
            }
        }
    } 
    if (verts == 65 && valence == 32 && adjac == 15 && ndjac == 16) {
        AJM(0,1) = 1;
        AJM(0,2) = 0; // partly undo the star
        for(i=3; i<=17; i++) {
            AJM(0,i) = 1;
            AJM(1,i) = 1;
        }
        for(i=0; i<16; i++) {
            AJM(0, 18+i) = 1;
            AJM(1, 34+i) = 1;
        }
        AJM(2, 18) = 1;
        AJM(2, 34) = 1;
        for(i=0; i<verts; i++) {
            if (ajm[0][i] < 0)
                AJM(0,i) = 0;
            if (ajm[1][i] < 0)
                AJM(1,i) = 0;
        }
    }
 
    if (mode == kMax) // special case
        goto count;
    if (verts == 99 && valence == 14 && adjac == 1 && ndjac == 2) {
        #define AJMK(x,y) AJM((2*x+13), (2*y+13))    
        switch(branch) {
            case 0: AJMK(1,2) = 1; AJMK(3,4) = 1; AJMK(5,6) = 1; AJMK(7,8) = 1; AJMK(9,10) = 1; AJMK(11,12) = 1; break;
            case 1: AJMK(1,3) = 1; AJMK(2,4) = 1; AJMK(5,6) = 1; AJMK(7,8) = 1; AJMK(9,10) = 1; AJMK(11,12) = 1; break;
            case 2: AJMK(1,3) = 1; AJMK(2,4) = 1; AJMK(5,7) = 1; AJMK(6,8) = 1; AJMK(9,10) = 1; AJMK(11,12) = 1; break;
            case 3: AJMK(1,3) = 1; AJMK(2,4) = 1; AJMK(5,7) = 1; AJMK(6,8) = 1; AJMK(9,11) = 1; AJMK(10,12) = 1; break;
           
            case 4: AJMK(1,6) = 1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(7,8) = 1; AJMK(9,10) = 1; AJMK(11,12) = 1; break;
            case 5: AJMK(1,6) = 1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(7,9) = 1; AJMK(8,10) = 1; AJMK(11,12) = 1; break;
            case 6: AJMK(1,6) = 1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(7,12) =1; AJMK(8, 9) = 1; AJMK(10,11) = 1; break;
            
            case 7: AJMK(1,8) = 1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(6,7) = 1; AJMK(9,10) = 1; AJMK(11,12) = 1; break;
            case 8: AJMK(1,8) = 1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(6,7) = 1; AJMK(9,11) = 1; AJMK(10,12) = 1; break;
            
            case 9: AJMK(1,10) =1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(6,7) = 1; AJMK(8, 9) = 1; AJMK(11,12) = 1; break;
            case 10:AJMK(1,12) =1; AJMK(2,3) = 1; AJMK(4,5) = 1; AJMK(6,7) = 1; AJMK(8, 9) = 1; AJMK(10,11) = 1; break;
            default: xit(1, "--branch must be 0..10 for this graph.");
        }
        for (i=1; i <= 12; i++) {
            for (j=i+1; j <=12; j++)
                if (ajm[2*i+13][2*j+13] < 0)
                    AJMK(i,j) = 0;
        }    
    
    } else if (verts == 65 && valence == 32 && adjac == 15 && ndjac == 16) {
        if (branch < 0 || branch > 15)
            xit(1, "--branch must be 0..15 for this graph.");
        for (i=0; i < branch; i++) {
            AJM(2,i+3) = 1;
            AJM(2,i+50) = 1; 
        }
        for (i=0; i < 16-branch; i++) {
            AJM(2,i+18) = 1;
            AJM(2,i+34) = 1; 
        }
        for(i=0; i<verts; i++) {
            if (ajm[2][i] < 0)
                AJM(2,i) = 0;
        }
    }

  count:;  
    cnt = 0;
    for(i=0; i < verts; i++) {
        for (j=i+1; j < verts; j++) {
            if (ajm[i][j] == -1) {
                ajm[i][j] = ajm[j][i] = --cnt;
            }
        }
    }
    return -cnt;
}

int main( int argc, char *argv[] ) {
    int i, j, k, verts=0, valence=0, adjac=0, ndjac=0;
    int vars=0, constrs=0, eqns=0, prods=0;
    int **ajm;
    Mode mode = kMax;
    int branch = 0;
    int no_subst = 0;
    int num_count = 0;
    
    for (i = 1; i < argc; i++) {
        if (!strcmp(argv[i],"--help")) {
            xit(0, "Usage: srg2opb --srg=VERTICES,VALENCE,ADJACENT,NONADJACENT [--level=none|star|max|NUMBER]\n"
                   "               [--no-subst]\n"
                   "Generate *.opb description file for a given SRG\n");
        } else if (!strcmp(argv[i],"--version")) {
            xit(0, "1.0\n");
        } else if (!strcmp(argv[i],"--no-subst")) {
            no_subst = 1;
        } else if (!strncmp(argv[i],"--level=", 8)) { 
            char *p=argv[i]+8;
            if (!strcmp(p, "none"))
                mode = kNone;
            else if (!strcmp(p, "star"))
                mode = kStar;
            else if (!strcmp(p, "max"))
                mode = kMax;
            else {
                for(j=0; p[j] >= '0' && p[j] <= '9'; j++)
                    ;
                if (j && p[j]==0) {
                    mode = kBranch;
                    branch = atoi(p);
                } else {
                    xit(1, "Wrong argument: %s\n", p);
                }
            }
        } else if (!strncmp(argv[i],"--srg=", 6)) {
            if(sscanf(argv[i]+6, "%d,%d,%d,%d", &verts, &valence, &adjac, &ndjac) != 4)
                xit(1, "Wrong --srg= format.\n");
        } else {
            xit(1, "Unknown argument: %s\n", argv[i]);
        }
    }

    if (verts==0 || valence == 0 || (verts-valence-1)*ndjac != (valence-adjac-1)*valence)
        xit(1, "Contradictary SRG params: (%d, %d, %d, %d)\n", verts, valence, adjac, ndjac);

    ajm = malloc(sizeof(int *)*verts);
    for (i=0; i<verts; i++) {
        ajm[i] = malloc(sizeof(int)*verts);
        memset(ajm[i], 0xff, sizeof(int)*verts);
    }
    vars = fill_matrix(ajm, verts, valence, adjac, ndjac, mode, branch);
    for(i=0; i<verts; i++) {
        for(j=i+1; j<verts; j++) { 
            int have_constr = 0;
            for(k=0; k<verts; k++) {
                if (ajm[i][k] < 0) {
                    if (ajm[k][j] < 0) { // product
                        prods++;
                    } else {
                        have_constr |= 1;
                    }
                } else {
                    if (ajm[k][j] < 0) {
                         have_constr |= 1;
                    }
                }
            }
            constrs += have_constr;
        }
        { int cnt = 0;
          int sum = 0;
          for (j=0; j<verts; j++) {
              if (ajm[i][j] < 0)
                  cnt++;
              else
                  sum += ajm[i][j];
          }
          if (cnt)
              eqns++;
          else if (sum != valence)
              xit(4, "Miscalculated line=%d\n",i);
        }
    }

    for(i=0; i<verts;i++)
        for(j=i+1; j<verts; j++)
            num_count += (ajm[i][j] >= 0);

    if (no_subst)
         printf("* #variable= %d #constraint= %d #product= %d sizeproduct= %d\n"
            ,verts*(verts-1)/2, verts*(verts+1)/2 + num_count, verts*(verts-1)*(verts-2)/2, verts*(verts-1)*(verts-2));
    else
        printf("* #variable= %d #constraint= %d #product= %d sizeproduct= %d\n"
            ,vars, constrs+eqns, prods, 2*prods);
        
    printf("* ");
    for (i=0; i<argc; i++)
        printf("%s ", argv[i]);
    printf("\n");
    
    for(i=0; i<verts; i++) {
        printf("* ");
        for(j=0; j<verts; j++) { 
            if(ajm[i][j] >= 0)
                printf("%d", ajm[i][j]);
            else
                printf("?");
        }
        printf("\n");
    }
    
    if (no_subst) {
        for(i=0; i < verts; i++)  // constants
            for(j=i+1; j<verts; j++)
                if (ajm[i][j] >= 0)
                    printf("+1 x%d = %+d;\n", ij2var(i, j, verts), ajm[i][j]);

        for (i=0; i<verts; i++) { // linear constrains
            for(j=0; j < verts; j++) {
                if (i != j)
                    printf("+1 x%d ",ij2var(i, j, verts));
            }
            printf("= %+d;\n", valence);
        }
        
        for (i=0; i<verts; i++) { // non-linear constrains
            for(j=i+1; j<verts; j++) {
                for(k=0; k<verts; k++) {
                    if (k != i && k != j)
                        printf("+1 x%d x%d ", ij2var(i, k, verts), ij2var(j, k, verts));                 
                }
                printf("%+d x%d = %+d;\n", ndjac-adjac, ij2var(i, j, verts), ndjac);
            }
        }
    } else {
        for (i=0; i<verts; i++) { // linear
            for(j=0; j<verts; j++) {
                if(ajm[i][j] < 0)
                    break;
            }
            if (j < verts) {
                int sum = 0;
                for(j=0; j<verts; j++) {
                    if (ajm[i][j] >= 0)
                        sum += ajm[i][j];
                    else
                        printf("+1 x%d ",-ajm[i][j]);
                }
                printf("= %+d;\n", valence-sum);
            }
        }
        for (i=0; i<verts; i++) { // non-linear
            for(j=i+1; j<verts; j++) {
                // printf("* %d %d\n", i, j);
      
                int var=0, sum = 0;   
                for(k=0; k<verts; k++) {
                    var |= product(ajm, &sum, i,k,k,j);
                }
                if (var) {
                    if (ajm[i][j] < 0)
                        printf("%+d x%d = %+d;\n", ndjac-adjac, -ajm[i][j], ndjac-sum);
                    else
                        printf(" = %+d;\n", ndjac-sum-(ndjac-adjac)*ajm[i][j]);
                }
            }
        }
    }
    for (i=0; i<verts; i++)
        free(ajm[i]);
    free(ajm);
    return 0;
}