#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <assert.h>
#include <stdarg.h>
#include <stdbool.h>
#include <sys/time.h>

#include <netdb.h>
#include <netinet/in.h>


/* Search for a graph one edge at a time */


typedef unsigned long long ull;
typedef const unsigned long long cull;

typedef struct tag_Decode {
    unsigned short i, j;
} Decode;

typedef struct tag_Match {
    int i, j, n, to;
} Match;

static inline void set2_set(ull *s, int i) {
    if (i < 64)
        s[0] |= 1ull << i;
    else
        s[1] |= 1ull << (i-64);
}

static inline ull set2_test(ull *s, int i) {
    if (i < 64)
        return s[0] & (1ull << i);
    else
        return s[1] & (1ull << (i-64));
}
                 
static inline int bitcount_1(cull *p, int i) {
    return __builtin_popcountll(p[i]);
}

static inline int bitcount_2(cull *p, int i) {
    return __builtin_popcountll(p[2*i]) +  __builtin_popcountll(p[2*i + 1]);
}

static inline int bitcount_n(cull *p, int n, int i) {
    int j, e = 0;
    for (j=0; j<n; j++)
        e += __builtin_popcountll(p[i*n + j]);
    return e;
}


static inline int bitdiff_1(cull *p, int i, int j, int b) {
    return abs(__builtin_popcountll(p[i] & p[j]) - b);
}

static inline int bitdiff_2(cull *p, int i, int j, int b) {
    return abs(__builtin_popcountll(p[2*i] & p[2*j]) + __builtin_popcountll(p[2*i + 1] & p[2*j + 1]) - b);
}

static inline int bitdiff_n(cull *p, int n, int i, int j, int b) {
    int k, e = 0;
    for(k=0; k<n; k++)
        e += __builtin_popcountll(p[i*n+k] & p[j*n+k]); 
    return abs(e - b);
}


static inline void toggle_edge_1(ull *p, int i, int j) {
    p[i] ^= 1ull<<j;
    p[j] ^= 1ull<<i;
}

static inline void toggle_edge_2(ull *p, int i, int j) {
    assert(j < i);
    if(i < 64) {
            p[2*i] ^= 1ull<<j;
            p[2*j] ^= 1ull<<i;
    } else {
        if (j < 64) {
            p[2*i] ^= 1ull<<j;
            p[2*j+1] ^= 1ull<<(i-64);
        } else {
            p[2*i+1] ^= 1ull<<(j-64);
            p[2*j+1] ^= 1ull<<(i-64);
        }
    }
}

static inline void toggle_edge_2u(ull *p, int i, int j) {
    if (j < i)
        return toggle_edge_2(p, i, j);
    else
        return toggle_edge_2(p, j, i);
}

static inline void toggle_edge_n(ull *p, int n, int i, int j) {
    p[i*n + j/64] ^= 1ull<<(j%64);
    p[j*n + i/64] ^= 1ull<<(i%64);
}

static void make_edge_n(ull *p, int n, int i, int j) {
    p[i*n + j/64] |= 1ull<<(j%64);
    p[j*n + i/64] |= 1ull<<(i%64);
}

static void make_edge_n2(ull *p, ull *q, int n, int i, int j) {
    make_edge_n(p, n, i, j);
    make_edge_n(q, n, i, j);
}

static inline ull has_edge_1(cull *p, int i, int j) {
    return p[i] & 1ull<<j;
}

static inline ull has_edge_2(cull *p, int i, int j) {
    if (j < 64) {
        return p[2*i] & 1ull<<j;
    } else {
        return p[2*i+1] & 1ull<<(j-64);
    }
}

static inline ull has_edge_n(const ull *p, int n, int i, int j) {
    return p[n*i + j/64] & 1ull<<(j%64);
}



#define CASE_OP(x)  case x:  e += bitdiff_2(p, i, x, 1 + !has_edge_2(p, i, x));
#define MULTICASE_OP(x) CASE_OP(x-1) CASE_OP(x-2) CASE_OP(x-3) CASE_OP(x-4) CASE_OP(x-5) \
                        CASE_OP(x-6) CASE_OP(x-7) CASE_OP(x-8) CASE_OP(x-9) CASE_OP(x-10)

static int error_count_99(cull *p) {
    int i;
    int e = 0;
    for (i=0; i < 99; i++) {
        switch(i-1) {
            MULTICASE_OP(100)
            MULTICASE_OP( 90)
            MULTICASE_OP( 80)
            MULTICASE_OP( 70)
            MULTICASE_OP( 60)
            MULTICASE_OP( 50)
            MULTICASE_OP( 40)
            MULTICASE_OP( 30)
            MULTICASE_OP( 20)
            MULTICASE_OP( 10)
        }
    }
    return e; 
}
#undef CASE_OP

#define CASE_OP(x)  case x: e += bitdiff_2(p, i, x, has_edge_2(p, i, x) ? adjac : ndjac);
static int error_count_2(cull *p, int verts, int adjac, int ndjac) {
    int i;
    int e = 0;
    for (i=0; i < verts; i++) {
        switch(i-1) {
            MULTICASE_OP(130)
            MULTICASE_OP(120)
            MULTICASE_OP(110)
            MULTICASE_OP(100)
            MULTICASE_OP( 90)
            MULTICASE_OP( 80)
            MULTICASE_OP( 70)
            MULTICASE_OP( 60)
            MULTICASE_OP( 50)
            MULTICASE_OP( 40)
            MULTICASE_OP( 30)
            MULTICASE_OP( 20)
            MULTICASE_OP( 10)
        }
    }
    return e; 
}
#undef CASE_OP

#define CASE_OP(x)  case x: e += bitdiff_1(p, i, x, has_edge_1(p, i, x) ? adjac : ndjac);
static int error_count_1(cull *p, int verts, int adjac, int ndjac) {
    int i;
    int e = 0;
    for (i=0; i < verts; i++) {
        switch(i-1) {
            MULTICASE_OP( 70)
            MULTICASE_OP( 60)
            MULTICASE_OP( 50)
            MULTICASE_OP( 40)
            MULTICASE_OP( 30)
            MULTICASE_OP( 20)
            MULTICASE_OP( 10)
        }
    }
    return e; 
}
#undef CASE_OP

static int error_count_n(cull *p, int n, int verts, int adjac, int ndjac) {
    if (n == 1)
        return error_count_1(p, verts, adjac, ndjac);
    else if (n == 2) {
        if (verts == 99 || adjac == 1 || ndjac == 2)
            return error_count_99(p);
        else 
            return error_count_2(p, verts, adjac, ndjac);
    } else {
        int i, j;
        int e = 0;
        for (i=0; i < verts; i++) {
            for (j = 0; j < i; j++)
                e += bitdiff_n(p, n, i, j, has_edge_n(p, n, i, j) ? adjac : ndjac); 
        }
        return e;
    } 
}


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 check_matrix_eqn(unsigned long long *r, int verts, int valence, int adjac, int ndjac) {
    int i, j, k;
    int wcount = (verts + 63)/64;
    
    for (i=0; i<verts; i++) {
        for (j=0; j<verts; j++) {
            int a1 = ndjac, a2 = 0;
            for (k=0; k < wcount; k++)
                a2 += __builtin_popcountll(r[i*wcount + k] & r[j*wcount + k]);
            if (i == j) 
                a1 += valence - ndjac;
            if (has_edge_n(r, wcount, i, j))
                a1 += adjac - ndjac;           
            if (a1 != a2) {
                printf("# Bad matrix\n");
                return 1;
            }
        }
    }
    return 0;
}

static char * strnchr(char *buf, int len, char c) {
    while(len-- > 0 && *buf != c && *buf != 0)
        buf++;
    return *buf == c ? buf : 0;
}


static int socket_connect(char *hostname, int port, char *msg, int len) {
    int fd;
    struct sockaddr_in serv_addr;
    struct hostent *server;
    
    /* Create a socket point */
    fd = socket(AF_INET, SOCK_STREAM, 0);
    if (fd < 0)
        xit(1, "ERROR opening socket\n");

    server = gethostbyname(hostname);
    if (server == NULL)
        xit(1, "No such host: %s\n", hostname);
  
    memset(&serv_addr, 0, sizeof(serv_addr));
    serv_addr.sin_family = AF_INET;
    memcpy(&serv_addr.sin_addr.s_addr, server->h_addr, server->h_length);
    serv_addr.sin_port = htons(port);

    /* Now connect to the server */
    if (connect(fd, (struct sockaddr*)&serv_addr, sizeof(serv_addr)) < 0)
        xit(1, "ERROR connecting\n");
    
    if(write(fd, msg, len) != len)
        xit(1, "ERROR writing to socket");
    
    return fd;
}

/*
static int in_array(int *array, int a, int cnt) {
    int i;
    for(i=0; i<cnt; i++) {
        if (a == array[i])
            return 1;
    }
    return 0;
}
*/


static int neighbours_2(int *buf, cull *p, int i, int j) {
    int k, *q=buf;
    ull a = p[2*i] & p[2*j];
    for (k=0; k<64; k++) {
        if (a & (1ull<<k))
            *q++ = k;
    }
    a = p[2*i + 1] & p[2*j + 1];
    for (k=0; k<64; k++) {
        if (a & (1ull<<k))
            *q++ = k+64;
    }
    return q-buf;
}


static int rank_proc(const void *p, const void *q) {
    return ((const int *)p)[1] - ((const int *)q)[1];
}

static int match_proc(const void *p, const void *q) {
    return ((const Match *)p)->n - ((const Match *)q)->n;
}

int main(int argc, char **argv) {
    int i, j;
    unsigned long long iter, iters=0;
    unsigned long long acc_error=0, acc_error2=0;
    unsigned case_no=0;
    unsigned tries=0;
    unsigned debug=0;
    int var_edges = 0;
    double foo_factor = 1, foo_exp_beg = 1, foo_exp_end = 1;
    unsigned int seed = ~0;
    char *fname = 0;
    char *portstr = 0;
    char *hostname = 0;
    int portno = 0;
    int no_graph = 0;
    int readjust = 0;
    int temperature, best_temperature, initial_temperature, max_temperature;
    int verts=0, valence=0, adjac=0, ndjac=0, words;
    ull *graph, *ngraph, *egraph, *best_graph;
    Decode *dec;
    bool more_edges=false;
    enum {kAnnealing, kSearch} mode = kAnnealing;
    int n_effect; // Effect of %n on the number of scanned fields.

    { int dummy1, dummy2;
      n_effect = sscanf("1", "%d%n", &dummy1, &dummy2) - 1;
    }
    
    for (i = 1; i < argc; i++) {
        if (!strcmp(argv[i],"--help")) {
            xit(0, "srg [--foo-factor=FLOAT] [--foo-exp=FLOAT] [--seed=INT] [--iters=INT] [--in=NAME]\n"
            "    [--more-edges] [--srg=INT,INT,INT,INT]\n");
        } else if (!strcmp(argv[i],"--version")) {
            xit(0, "1.0\n");
        } else if (!strcmp(argv[i],"--no-graph")) {
            no_graph = 1;
        } else if (!strncmp(argv[i],"--foo-factor=", 13)) { 
            sscanf(argv[i]+13, "%lf", &foo_factor);
        } else if (!strncmp(argv[i],"--foo-exp-beg=", 14)) {
            sscanf(argv[i]+14, "%lf", &foo_exp_beg);
        } else if (!strncmp(argv[i],"--foo-exp-end=", 14)) {
            sscanf(argv[i]+14, "%lf", &foo_exp_end);
        } else if (!strncmp(argv[i],"--seed=", 7)) {
            sscanf(argv[i]+7, "%u", &seed);
        } else if (!strncmp(argv[i],"--case=", 7)) {
            sscanf(argv[i]+7, "%u", &case_no);
        } else if (!strcmp(argv[i], "--debug")) {
            debug = 1;
        } else if (!strncmp(argv[i],"--iters=", 8)) {
            char c;
            if (sscanf(argv[i]+8, "%llu%c", &iters, &c) == 2) {
                switch(c) {
                case 'K': case 'k':
                    iters *= 1000;
                    break;
                case 'M': case 'm':
                    iters *= 1000000;
                    break;
                case 'G': case 'g':
                    iters *= 1000000000;
                    break;
                default:
                    xit(1, "Unknown multiplier: %c\n", c);
                }
            }
        } else if (!strncmp(argv[i],"--in=", 5)) {
            fname = argv[i]+5;
            portstr = strchr(fname,':');
            if (portstr) {
                if (strlen(portstr) < 2)
                    portno = 5000;
                else
                    portno = atoi(portstr + 1);
                hostname = malloc(portstr-fname+1);
                if (!hostname)
                    xit(1, "malloc failed\n");
                memcpy(hostname, fname, portstr-fname);
                hostname[portstr-fname] = 0;
            }
        } 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 if (!strcmp(argv[i], "--more-edges")) {
            more_edges = true;
        } else if (!strncmp(argv[i], "--mode=", 7)) {
            if (!strcmp(argv[i] + 7, "annealing"))
                mode = kAnnealing;
            else if (!strcmp(argv[i] + 7, "search"))
                mode = kSearch;
            else
                xit(1, "Unknown mode: %s\n", argv[i]+7);
        } else if (!strncmp(argv[i], "--readjust=", 11)) {
            sscanf(argv[i]+11, "%d", &readjust);
        } else {
            xit(1, "Unknown argument: %s\n", argv[i]);
        }
    }
    if (seed == ~0) {
        struct timeval tv;
        if (gettimeofday(&tv, 0))
            xit(1, "gettimeofday() failed\n");
        seed = tv.tv_sec ^ tv.tv_usec;
    }
    srandom(seed);
   
    if (fname == 0) {
        if (verts <=0 || valence <= 0)
            xit(2, "Missing --srg= option\n");

        if ((verts-valence-1)*ndjac != (valence-adjac-1)*valence)
            xit(2, "Contradictary SRG params.\n");

        // Allocating and zeroing the graph
        words = (verts+63)/64;
        graph = (ull *)calloc(words * verts, sizeof(ull));  // workspace
        ngraph = (ull *)calloc(words * verts, sizeof(ull)); // known non-edges
        egraph = (ull *)calloc(words * verts, sizeof(ull)); // known edges
        best_graph = (ull *)calloc(words * verts, sizeof(ull)); // best graph found
        dec = (Decode *)calloc(verts*(verts-1)/2, sizeof(Decode));

        if (graph == 0 || ngraph == 0 || egraph == 0 || best_graph == 0 || dec == 0)
            xit(3, "Graph allocation error.\n"); 

        // Draw known edges
        for (i=0; i < verts; i++)   // no loops
            make_edge_n(ngraph, words, i, i);
          
        if (more_edges && verts==65 && valence==32 && adjac==15 && ndjac==16) {

            make_edge_n2(graph, egraph, words, 0, 1);

            for (i=2; i <= 16; i++) {
                make_edge_n2(graph, egraph, words, 0, i);
                make_edge_n2(graph, egraph, words, 1, i);
            }
            for (i=0; i < 16; i++) {
                make_edge_n2(graph, egraph, words, 0, i+17);
                make_edge_n2(graph, egraph, words, 1, i+33);
            }
            for (i=17; i < 65; i++) {
                if (!has_edge_n(graph, words, 0, i))
                    make_edge_n(ngraph, words, 0, i);
                if (!has_edge_n(graph, words, 1, i))
                    make_edge_n(ngraph, words, 1, i);
            }
            make_edge_n2(graph, egraph, words, 17, 49); // always there, may be set aganin later
            make_edge_n2(graph, egraph, words, 33, 49); // ditto
            if (case_no >= 0 && case_no < 16) {
                for (i=0; i < case_no; i++) {
                    make_edge_n2(graph, egraph, words, i+2, 49);
                    make_edge_n2(graph, egraph, words, 49, i+50);
                }
                for (i=0; i < 16-case_no; i++) {
                    make_edge_n2(graph, egraph, words, i+17, 49);
                    make_edge_n2(graph, egraph, words, i+33, 49);
                }
                for (i=0; i < 65; i++) {
                    if (!has_edge_n(graph, words, 49, i))
                        make_edge_n(ngraph, words, 49, i);
                }
            }

        } else if (more_edges && verts==99 && valence==14 && adjac==1 && ndjac==2) {
            int k;
            int ac[14];
            for (i=1; i <= 14; i++)
                make_edge_n2(graph, egraph, words, 0, i);
            for (i=1; i <= 13; i+=2)
                make_edge_n2(graph, egraph, words, i, i+1);
            for (k=15, i=1; i <= 14; i++) {
                for (j=1; j<i; j++) {
                    if (!has_edge_n(graph, words, i, j)) {
                        make_edge_n2(graph, egraph, words, i, k);
                        make_edge_n2(graph, egraph, words, j, k);
                        if (j==1)
                            ac[i-2] = k;
                        k += 1;
                    }
                }
            }
            if (k != 99)
                xit(4, "Miscalculated k = %d\n", k);
            for (i=0; i<=14; i++) {
                for (j=0; j<99; j++) {
                    if (!has_edge_n(graph, words, i, j))
                        make_edge_n(ngraph, words, i, j);
                }
            }
            if (case_no >= 0) {
                #define AJMK(x,y) make_edge_n2(graph, egraph, words, ac[x], ac[y])
                switch(case_no) {
                    case 0: AJMK(1,2); AJMK(3,4); AJMK(5,6); AJMK(7,8); AJMK(9,10); AJMK(11,12); break;
                    case 1: AJMK(1,3); AJMK(2,4); AJMK(5,6); AJMK(7,8); AJMK(9,10); AJMK(11,12); break;
                    case 2: AJMK(1,3); AJMK(2,4); AJMK(5,7); AJMK(6,8); AJMK(9,10); AJMK(11,12); break;
                    case 3: AJMK(1,3); AJMK(2,4); AJMK(5,7); AJMK(6,8); AJMK(9,11); AJMK(10,12); break;
                   
                    case 4: AJMK(1,6); AJMK(2,3); AJMK(4,5); AJMK(7,8); AJMK(9,10); AJMK(11,12); break;
                    case 5: AJMK(1,6); AJMK(2,3); AJMK(4,5); AJMK(7,9); AJMK(8,10); AJMK(11,12); break;
                    case 6: AJMK(1,6); AJMK(2,3); AJMK(4,5); AJMK(7,12); AJMK(8, 9); AJMK(10,11); break;
                    
                    case 7: AJMK(1,8); AJMK(2,3); AJMK(4,5); AJMK(6,7); AJMK(9,10); AJMK(11,12); break;
                    case 8: AJMK(1,8); AJMK(2,3); AJMK(4,5); AJMK(6,7); AJMK(9,11); AJMK(10,12); break;
                    
                    case 9: AJMK(1,10); AJMK(2,3); AJMK(4,5); AJMK(6,7); AJMK(8, 9); AJMK(11,12); break;
                    case 10:AJMK(1,12); AJMK(2,3); AJMK(4,5); AJMK(6,7); AJMK(8, 9); AJMK(10,11); break;
                    default: xit(1, "--case must be 0..10 for this graph.");
                }
                for (i=1; i<=12; i++) {
                    for(j=1; j<i; j++) {
                        if (!has_edge_n(graph, words, ac[i], ac[j]))
                            make_edge_n(ngraph, words, ac[i], ac[j]);
                    }
                }
            }           
        } else {
            for (i=1; i < verts; i++) { 
                if (i <= valence) {
                    make_edge_n2(graph, egraph, words, 0, i);
                } else {
                    make_edge_n(ngraph, words, 0, i);
                }
            }
        }            
       
        // Set up the decoder
        for (i = 0; i < verts; i++) {
           for (j = 0; j < i; j++) {
               if (!has_edge_n(graph, words, i, j) && !has_edge_n(ngraph, words, i, j)) {
                   dec[var_edges].i = i;
                   dec[var_edges].j = j;
                   var_edges += 1;
               }
           }
        }
        
    } else {
        char *inbuf = 0, *pbuf, *pend;
        size_t insize=0;
        const size_t page_sz = 2048;
        
        if (portstr) { // host:port
            int inbuf_size = 0;
            int fd = socket_connect(hostname, portno, "# SEEDME", 8);

            for(;;) {
                size_t csz;
                
                if (inbuf_size - insize < page_sz/2) {
                    inbuf_size += page_sz;
                    inbuf=realloc(inbuf, inbuf_size);
                }
                if (inbuf == 0)
                    xit(1, "Buffer re-allocation error\n");
            
                csz = read(fd, inbuf+insize, inbuf_size-insize);
                if (csz < 0)
                    xit(1, "Error reading from socket\n");
                insize += csz;
                if (csz == 0)
                    break;
            }
            shutdown(fd, SHUT_RDWR);
            close(fd);
        } else {
            FILE *f;
            if (!strcmp(fname, "-"))
                f = stdin;
            else {
                f = fopen(fname, "rb");
                if (!f)
                    xit(1, "Can't open %s\n", fname);
            }
            for(;;) {
                size_t csz;
                
                inbuf=realloc(inbuf, insize+page_sz);
                if (!inbuf)
                    xit(1, "Buffer re-allocation error\n");
            
                csz = fread(inbuf+insize, 1, page_sz, f);
                if (csz != page_sz) {
                    insize += csz;
                    break;
                }
                insize += page_sz;
            }
            if (f != stdin)
                fclose(f);
        }    
  
        pend = inbuf + insize;
        pbuf = inbuf; // Start from the beginning
    
        while (pbuf[0] == '#') { // Process the header
            char *pend = strnchr(pbuf, insize-(pbuf-inbuf), '\n');
            if (pend == 0)
                xit(1, "Missing end of line\n");
            *pend=0;
            if (!memcmp(pbuf, "# SRG(", 6)) {
                if(sscanf(pbuf+6, "%d,%d,%d,%d", &verts, &valence, &adjac, &ndjac) != 4)
                    xit(1, "Wrong SRG comment.\n");
            } else if (!memcmp(pbuf, "# Case", 6)) {
                if(sscanf(pbuf+6, "%u", &case_no) != 1)
                    xit(1, "Wrong Case comment.\n");
            } else if (!memcmp(pbuf, "# Tries", 7)) {
                if(sscanf(pbuf+7, "%u", &tries) != 1)
                    xit(1, "Wrong Case= command.\n");
            }
            pbuf = pend+1;
        }

        if ((verts-valence-1)*ndjac != (valence-adjac-1)*valence)
            xit(2, "Contradictary SRG params.\n");

        // Allocating and zeroing the graph
        words = (verts+63)/64;
        graph = (ull *)calloc(words * verts, sizeof(ull));  // workspace
        ngraph = (ull *)calloc(words * verts, sizeof(ull)); // known non-edges
        egraph = (ull *)calloc(words * verts, sizeof(ull)); // known edges
        best_graph = (ull *)calloc(words * verts, sizeof(ull)); // best graph found
        dec = (Decode *)calloc(verts*(verts-1)/2, sizeof(Decode));

        if (graph == 0 || ngraph == 0 || egraph == 0 || best_graph == 0 || dec == 0)
            xit(3, "Graph allocation error.\n"); 

        for(i=0; i < verts; i++) {
            int v, n;
            if (sscanf(pbuf, "%d%n", &v, &n) != 1 + n_effect || v != i) {
                xit(1, "Incorrect counter %d at the line %d\n", v, i);
            }
            pbuf += n;
            for(j=0; j < verts; j++) {
                while(pbuf < pend && isspace(*pbuf))
                    pbuf++;
                if (pbuf >= pend)
                    xit(1,"Unexpected end of block\n");
                switch (*pbuf) {
                    case '-':
                        if (readjust == 0 || i == j)
                            make_edge_n(ngraph, words, i, j);
                        break;
                    case '+':
                        if (readjust == 0)
                            make_edge_n2(graph, egraph, words, i, j);
                        else
                            make_edge_n(graph, words, i, j);
                        break;
                    case '0':
                        break;
                    case '1':
                        make_edge_n(graph, words, i, j);
                        break;
                    default:
                        xit(1, "Incorrect character at (%d,%d)\n", i, j);
                }
                pbuf += 1;
            }
        }
        free (inbuf);
        inbuf = pbuf = pend = 0; // for safety
        
        // Set up the decoder
        for (i = 0; i < verts; i++) {
            for (j = 0; j < i; j++) {
                if (!has_edge_n(egraph, words, i, j) && !has_edge_n(ngraph, words, i, j)) {
                    dec[var_edges].i = i;
                    dec[var_edges].j = j;
                    var_edges += 1;
                }
            }
        }
    }

    max_temperature = best_temperature = initial_temperature = temperature = error_count_n(graph, words, verts, adjac, ndjac);
    if (mode == kAnnealing) { 
        int jj;
        int t = temperature;
        double diff = foo_exp_end - foo_exp_beg;
        memcpy(best_graph, graph, words*verts*sizeof(ull));
      
        if (readjust) {
            
            int xoox = random() % verts;

            for (jj=0; jj<3; jj++) {

            if (verts==99 && valence==14 && adjac==1 && ndjac==2) {
                int star[98];
                int istar = 0;
                for(i=0; i<99; i++) {
                    if (i != xoox && has_edge_2(graph, xoox, i))
                        star[istar++] = i;
                }
                if (debug) {
                    fprintf(stderr, "Center: %d\n", xoox);
                    fprintf(stderr, "Star edges: ");
                    for (i=0; i<istar; i++)
                        fprintf(stderr, "%d ", star[i]);
                    fprintf(stderr, "\nError: %d\n", error_count_n(graph, words, verts, adjac, ndjac));
                }

                if (debug)
                   fprintf(stderr, "Starred: ");
                while(istar < 14) {
                    int imin = 0;
                    int emin = 100000000;
                    for (i=0; i<99; i++) {
                        if (i != xoox && !has_edge_2(graph, xoox, i)) {
                            int e;
                            toggle_edge_2u(graph, xoox, i);
                            e = error_count_99(graph);
                            if (e < emin) {
                                imin = i;
                                emin = e;
                            }
                            toggle_edge_2u(graph, xoox, i);
                        }
                    }
                    toggle_edge_2u(graph, xoox, imin);
                    star[istar++] = imin;
                    if (debug)
                       fprintf(stderr, "%d ", imin);
                }
                if (debug)
                    fprintf(stderr, "\nError: %d\n", error_count_n(graph, words, verts, adjac, ndjac));

                if (debug)
                   fprintf(stderr, "Un-starred: ");
                while (istar > 14) {
                    int imin = 0;
                    int emin = 100000000;
                    for (i=0; i<istar; i++) {
                        int e;
                        toggle_edge_2u(graph, xoox, star[i]);
                        e = error_count_99(graph);
                        if (e < emin) {
                            imin = i;
                            emin = e;
                        }
                        toggle_edge_2u(graph, xoox, star[i]);
                    }
                    toggle_edge_2u(graph, xoox, star[imin]);
                    memmove(star+imin, star+imin+1, (istar-imin-1)*sizeof(*star));
                    istar--;
                    if (debug)
                       fprintf(stderr, "%d ", star[imin]);
                }
                if (debug)
                    fprintf(stderr, "\nError: %d\n", error_count_n(graph, words, verts, adjac, ndjac));
  
                for (i=0,j=14; i<j; ) {
                    int k, e = 0;
                    if (debug)
                        fprintf(stderr, "\nProcessing edge: %d\n", star[i]);
                    
                    for (k=i+1; k<j; k++) {
                        if (has_edge_2(graph, star[i], star[k])) {
                            if (e) {
                                toggle_edge_2u(graph, star[i], star[k]);
                                if (debug)
                                    fprintf(stderr, "Delete bridge: %d %d\n", star[i], star[k]);
                            } else {
                                e = k;
                                if (debug)
                                    fprintf(stderr, "Keep bridge: %d %d\n", star[i], star[k]);
                            }
                        }
                    }
                    if (e) {
                        int t = star[i+1];
                        star[i+1] = star[e];
                        star[e] = t;
                        for (k=0; k<14; k++) {
                            if (k != i && has_edge_2(graph, star[i+1], star[k])) {
                                toggle_edge_2u(graph, star[i+1], star[k]);
                                if (debug)
                                    fprintf(stderr, "Delete co-bridge: %d %d\n", star[i+1], star[k]);
                            }
                        }
                        i += 2;
                    } else {
                        j -= 1;
                        int t = star[i];
                        star[i] = star[j];
                        star[j] = t;
                    }
                }
                if (debug)
                    fprintf(stderr, "Error: %d\n", error_count_n(graph, words, verts, adjac, ndjac));

                if (j & 1)
                    xit(1, "Odd number of unconnected nodes = %d\n", j);
                for (; j<14; j+=2) {
                    toggle_edge_2u(graph, star[j], star[j+1]);
                    if (debug)
                        fprintf(stderr, "Add bridge: %d %d\n", star[j], star[j+1]);
                }
                if (debug)
                    fprintf(stderr, "\nError: %d\n", error_count_n(graph, words, verts, adjac, ndjac));

            
                // testing
                for (i=0; i<14; i+=2) {
                    if(!has_edge_2(graph, xoox, star[i])) {
                        fprintf(stderr, "Missing edge: %d %d\n", xoox, star[i]);
                        exit(1);
                    }
                    if(!has_edge_2(graph, xoox, star[i+1])) {
                        fprintf(stderr, "Missing edge: %d %d\n", xoox, star[i+1]);
                        exit(1);
                    }
                    if(!has_edge_2(graph, star[i], star[i+1])) {
                        fprintf(stderr, "Missing edge: %d %d\n", star[i], star[i+1]);
                        exit(1);
                    }
                    for(j=i+2; j<14; j++) {
                        if(has_edge_2(graph, star[i], star[j])) {
                            fprintf(stderr, "Extra edge: %d %d\n", star[i], star[j]);
                            exit(1);
                        }
                    }
                }

                // corona
                for (i=0; i<14; i+=2) {
                    if (bitdiff_2(graph, star[i], star[i+1], 1)) {
                        int buf[99];
                        int n = neighbours_2(buf, graph, star[i], star[i+1]);
                        for (j=0; j<n; j++) {
                            if (buf[j] != xoox) {
                                toggle_edge_2u(graph, buf[j], star[i]);
                                toggle_edge_2u(graph, buf[j], star[i+1]);
                                if (debug)
                                    fprintf(stderr, "Nix corona over bridge: %d-%d-%d\n", star[i], buf[j], star[i+1]);
                            }
                        }
                    }
                }
                if (debug)
                    fprintf(stderr, "Error: %d\n", error_count_n(graph, words, verts, adjac, ndjac));
                
                { Match match[84];
                  int k, buf[99];
                  ull set[2] = {0, 0};
                  
                  for (k=i=0; i<14; i++) {
                      for (j=i+1; j<14; j++) {
                          if ((i & 1) == 0 && j == i+1)
                              continue;
                          match[k].i = i;
                          match[k].j = j;
                          match[k].n = bitdiff_2(graph, star[i], star[i+1], 0);
                          match[k].to = -1;
                          k++;
                      }
                  }
                  if (k != 84) {
                      fprintf(stderr, "Miscalculated number of pairs =%d\n", k);
                      exit(1);
                  }
                  qsort(match, 84, sizeof(*match), match_proc);
                  if (debug) {
                      fprintf(stderr, "match = (\n");
                      for (i=0; i< 84; i++)
                          fprintf(stderr, "(%d %d %d)%c", match[i].i, match[i].j, match[i].n, i%12 == 11 ? '\n' : ' ');
                      fprintf(stderr, ")\n");
                  }
                  if (match[0].n <= 0) {
                      fprintf(stderr, "No link to center for: %d %d\n", match[0].i, match[0].j);
                      exit(1);
                  }
                  for (i=0; i<84; i++) {
                      int n = neighbours_2(buf, graph, star[match[i].i], star[match[i].j]);
                      int xoox_cnt = 0;
                      for (j=0; j<n; j++) {
                          if (buf[j] == xoox)
                              xoox_cnt += 1;
                      }
                      if (xoox_cnt == 0) {
                          fprintf(stderr, "No links to the center from %d %d\n", star[match[i].i], star[match[i].j]);
                          exit(1);
                      }
                      if (xoox_cnt > 1) {
                          fprintf(stderr, "Too many links to the center from %d %d\n", star[match[i].i], star[match[i].j]);
                          exit(1);
                      }
                      for (j=0; j<n; j++) {
                          if (buf[j] != xoox) {
                              if (set2_test(set, buf[j])) { // already used
                                  toggle_edge_2u(graph, star[match[i].i], buf[j]);
                                  toggle_edge_2u(graph, star[match[i].j], buf[j]);
                                  if (debug) {
                                      fprintf(stderr, "Nix corona to used point: %d-%d-%d\n", star[match[i].i], buf[j], star[match[i].j]);
                                  }
                              } else {
                                  if (match[i].to < 0) {
                                      match[i].to = buf[j];
                                      set2_set(set, buf[j]);
                                  } else { // disconnect extra points
                                      toggle_edge_2u(graph, star[match[i].i], buf[j]);
                                      toggle_edge_2u(graph, star[match[i].j], buf[j]);
                                      if (debug) {
                                          fprintf(stderr, "Nix double corona: %d-%d-%d\n", star[match[i].i], buf[j], star[match[i].j]);
                                      }
                                  }
                              }
                          }
                      }
                  }
                  if (debug)
                      fprintf(stderr, "Error: %d\n", error_count_n(graph, words, verts, adjac, ndjac));
                  set2_set(set, xoox);
                  for (i=0; i<14; i++)
                      set2_set(set, star[i]);
                  for(i=0, j=-1; i<84; i++) {
                      if (match[i].to < 0) {
                          while (++j < 99 && set2_test(set, j) != 0)
                              ;
                          if (j >= 99) {
                              fprintf(stderr, "Run out of free nodes.\n");
                              exit(1);
                          }
                          if (!has_edge_2(graph, star[match[i].i], j))
                              toggle_edge_2u(graph, star[match[i].i], j);
                          if (!has_edge_2(graph, star[match[i].j], j))
                              toggle_edge_2u(graph, star[match[i].j], j);
                          if (debug)
                              fprintf(stderr, "Making corona: %d-%d-%d\n", star[match[i].i], j, star[match[i].j]);
                      }
                  }
                  if (debug)
                      fprintf(stderr, "Error: %d\n", error_count_n(graph, words, verts, adjac, ndjac));
              }
            }
            // Always accept results of readjustment 
            temperature = error_count_n(graph, words, verts, adjac, ndjac);
            if (temperature > max_temperature)
                max_temperature = temperature;
        }
        
        }
        
        
        for(iter = 0; t && iter < iters; iter++) {  
            int a = random() % var_edges;

            toggle_edge_n(graph, words, dec[a].i, dec[a].j);
            t = error_count_n(graph, words, verts, adjac, ndjac);
            if (iter == iters/2) {
                // search for the best graph in the 2nd half of the process
                if (best_temperature >= initial_temperature) {
                    best_temperature = t;
                    memcpy(best_graph, graph, words*verts*sizeof(ull));
                }
            }
            if (t < temperature) {
                temperature = t;
                if (t < best_temperature) {
                    best_temperature = t;
                    memcpy(best_graph, graph, words*verts*sizeof(ull));
                }
            } else {
                double p, k;
                k = (double)iter/iters;
                p = foo_factor*exp((foo_exp_beg + k*diff)*(temperature - t));
                
                if (p*RAND_MAX > random()) {
                    temperature = t;
                    if (t > max_temperature) // for statistics
                        max_temperature = t;
                } else {
                    toggle_edge_n(graph, words, dec[a].i, dec[a].j);
                }
            }
            acc_error += t - best_temperature;
            acc_error2 += (t - best_temperature)*(t - best_temperature);
        }
    } else if (mode == kSearch) {
        int elems[64];
        int prev_gray = 1<<(iters-1);
        int pow_iters = iters;
        int i;
       
        
        if (iters > 40)
            xit(1,"iters must be 0..40 to perform 2**iters tries\n");
        iter=iters = 1ull<<iters;
 
 #if 1
        { // Pick least fit edges.
          int oe = error_count_n(graph, words, verts, adjac, ndjac);
          int *erank = malloc(2*var_edges*sizeof(int));
          for (i=0; i<var_edges; i++) {
              int ne;
              toggle_edge_n(graph, words, dec[i].i, dec[i].j);
              ne = error_count_n(graph, words, verts, adjac, ndjac);
              toggle_edge_n(graph, words, dec[i].i, dec[i].j);
              erank[2*i] = i;
              erank[2*i+1] = ne - oe;
          }
          qsort(erank, var_edges, 2*sizeof(int), rank_proc);
          for (i=0; i < pow_iters; i++)
              elems[i] = erank[2*i];
          free(erank);
        }
#else       
        // Pick random edges.
        for (i=0; i < pow_iters; i++) {
            int a;
            do {
                a = random() % var_edges;
            } while(in_array(elems, a, i));
            elems[i] = a; 
        }
#endif

        memcpy(best_graph, graph, words*verts*sizeof(ull));
        for (iter=0; iter < iters; iter++) {
            int t;
            ull gray = iter ^ (iter >> 1);
            int j = __builtin_ctz(prev_gray ^ gray);
            toggle_edge_n(graph, words, dec[elems[j]].i, dec[elems[j]].j);
            t = error_count_n(graph, words, verts, adjac, ndjac);
            if (t > max_temperature) // for statistics
                max_temperature = t;
            else if (t < best_temperature) {
                best_temperature = t;
                memcpy(best_graph, graph, words*verts*sizeof(ull));
            }
            acc_error += t - best_temperature;
            acc_error2 += (t - best_temperature)*(t - best_temperature);
            prev_gray=gray;
        }
    } 
   
    if (best_temperature == 0) {
        check_matrix_eqn(best_graph, verts, valence, adjac, ndjac);
    }
    
    { // Output
      int len, header_len;
      char *p, *buf = malloc(verts*(verts + 6) + 1000);
      int end_temp = error_count_n(graph, words, verts, adjac, ndjac);
      if (buf == 0)
          xit(1, "Can't allocate output buffer\n");
      
      p = buf;
      p += sprintf(p, "# SRG(%d, %d, %d, %d)\n", verts, valence, adjac, ndjac);
      p += sprintf(p, "# foo-factor=%lf foo-exp=%lf - %lf\n", foo_factor, foo_exp_beg, foo_exp_end);
      p += sprintf(p, "# Error %d => %d\n", initial_temperature, best_temperature);
      p += sprintf(p, "# Final error= %d\n", end_temp);
      p += sprintf(p, "# Max error= %d\n", max_temperature);
      p += sprintf(p, "# Final iter=%llu, %lf, iters=%llu\n", iter, iters ? (double)iter/iters : 0., iters);
      p += sprintf(p, "# Case %u\n", case_no);
      p += sprintf(p, "# Tries %u\n", tries+1);
      p += sprintf(p, "# Neighborhood %lf +/- %lf\n"
        ,(double)acc_error/iters, sqrt((double)acc_error2/iters 
                                        -  ((double)acc_error/iters)*((double)acc_error/iters)));
      header_len = p - buf;
      
      for(i=0; i < verts; i++) {
         p += sprintf(p, verts <= 100 ? "%02d " : verts <= 1000 ? "%03d " : "%04d ", i);
         for (j=0; j < verts; j++)
         if (has_edge_n(ngraph, words, i, j))
             *p++ = '-';
         else if (has_edge_n(egraph, words, i, j))
             *p++ = '+';
         else
             *p++ = has_edge_n(best_graph, words, i, j) ? '1' : '0';
         *p++ = '\n';
      }

      len = p-buf;
      if (portstr) {
          int fd = socket_connect(hostname, portno, "# RESULT", 8);
          int sent = 0;
          while(sent < len) {
              int n = write(fd, buf + sent, len - sent);
              if (n < 0)
                xit(1, "ERROR writing to socket");
              sent += n;
          }
          shutdown(fd, SHUT_RDWR);
          close(fd);
      }
      fwrite(buf, 1, no_graph ? header_len : len, stdout);
      free(buf);
    }
   
    free(hostname);
    free(graph);
    free(ngraph);
    free(egraph);
    free(best_graph);
    free(dec);
    return 0;
}