#include "cado.h"
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <string.h>
#include <inttypes.h>

#include "utils.h"

#define DEBUG 0
#define MAPLE 0

int
checkVector(int *vec, int ncols)
{
    int ok, i;

    ok = 1;
    for(i = 0; i < ncols; i++)
	if(vec[i] & 1){
	    ok = 0;
	    break;
	}
#if 0
    if(ok)
	printf(" -> y\n");
    else
	printf(" -> n (%d)\n", i);
#endif
    return ok;
}

// str is a row coming from purgedfile
void
str2Vec(int *vec, char *str)
{
    char *t = str;
    int j, k = 0, nc;

    // skip first integer which is the index
    for(; *t != ' '; t++);
#if 1
    sscanf(t+1, "%d", &nc);
    for(j = 0; j < nc; j++){
	// t is on a "space"
	t++;
	// t begins on a number
	sscanf(t, PURGE_INT_FORMAT, &k);
	vec[k]++;
	if(j == (nc-1))
	    // no need to read further
	    break;
	// go to end of number
	for(; *t != ' '; t++);
    }
#else
    // skip second integer which is the number of primes
    for(++t; *t != ' '; t++);
    t++;
    while(1){
	if((*t == '\n') || (*t == ' ')){
	    // new integer read
	    vec[k]++;
	    k = 0;
	    if(*t == '\n')
		break;
	}
	else
	    k = k*10+(*t - '0');
	t++;
    }
#endif
}

void
treatRationalRelation(hashtable_t *H, relation_t rel)
{
    int j, h;
    uint64_t minus2 = (H->need64) ? (uint64_t) (-2) : (uint32_t) (-2);

#if MAPLE >= 1
    fprintf(stderr, "P:=P * (%ld-m*%lu):\n", rel.a, rel.b);
#endif
    for(j = 0; j < rel.nb_rp; j++){
#if MAPLE >= 1
	fprintf(stderr, "R2:=R2*%ld^%d:\n", rel.rp[j].p, rel.rp[j].e);
#endif
	h = getHashAddr(H, rel.rp[j].p, minus2);
	if(H->hashcount[h] == 0){
	    // new empty place
            SET_HASH_P(H,h,rel.rp[j].p);
            SET_HASH_R(H,h,minus2);
	}
	H->hashcount[h] += rel.rp[j].e;
    }
}

void
computefree(relation_t *rel)
{
    int i;

    for(i = 0; i < rel->nb_ap; i++){
	rel->ap[i].r = rel->ap[i].p;
	rel->ap[i].p = rel->a;
	rel->ap[i].e = 1;
    }
}

void
treatAlgebraicRelation(FILE *algfile, hashtable_t *H, relation_t rel)
{
    int j, h;

    if(rel.b > 0)
	computeroots(&rel);
    else
	computefree(&rel);
    fprintf(algfile, "%ld %lu\n", rel.a, rel.b);
#if MAPLE >= 1
    fprintf(stderr, "NORM:=NORM*mynorm(f, %ld, %lu):\n", rel.a, rel.b);
    fprintf(stderr, "AX:=AX*(%ld-%lu*X):\n", rel.a, rel.b);
#endif
    for(j = 0; j < rel.nb_ap; j++){
#if MAPLE >= 2
	fprintf(stderr, "A2:=A2*[%ld, %ld]^%d:\n",
		rel.ap[j].p, rel.ap[j].r, rel.ap[j].e);
#endif
	h = getHashAddr(H, rel.ap[j].p, rel.ap[j].r);
	if(H->hashcount[h] == 0){
	    // new empty place
            SET_HASH_P(H,h,rel.ap[j].p);
            SET_HASH_R(H,h,rel.ap[j].r);
	}
	H->hashcount[h] += rel.ap[j].e;
    }
}

void
finishRationalSqrt(FILE *ratfile, hashtable_t *H, cado_poly pol)
{
    mpz_t prod;
    int j;
    unsigned int i;
    uint64_t minus2 = (H->need64) ? (uint64_t) (-2) : (uint32_t) (-2);

    mpz_init_set_ui(prod, 1);
#if MAPLE >= 1
    fprintf(stderr, "R:=1;\n");
#endif
    for(i = 0; i < H->hashmod; i++){
	if(H->hashcount[i] > 0){
          if (GET_HASH_R(H,i) != minus2)
		continue;
	    if ((H->hashcount[i] & 1)) {
	        fprintf(stderr, "  Odd valuation! At rational prime %"PRIi64"\n",
                        GET_HASH_P(H,i));
		exit(1);
	    }
#if MAPLE >= 1
	    fprintf(stderr, "R:=R*%ld^%d:\n",
		    H->hashtab_p[i], H->hashcount[i]>>1);
#endif
	    // TODO: do better
	    for(j = 0; j < (H->hashcount[i]>>1); j++)
              mpz_mul_ui(prod, prod, GET_HASH_P(H,i));
	    mpz_mod(prod, prod, pol->n);
	}
    }
#if DEBUG >= 1
    gmp_fprintf(stderr, "prod:=%Zd;\n", prod);
    fprintf(stderr, "# We print the squareroot of the rational side...\n");
#endif
    // TODO: humf, I should say...!
    gmp_fprintf(ratfile, "%Zd\n", prod);
    mpz_clear(prod);
}

// returns the sign of m1*a+m2*b
int
treatSign(relation_t rel, cado_poly pol)
{
    mpz_t tmp1, tmp2;
    int s;

    /* first perform a quick check */
    s = (rel.a > 0) ? mpz_sgn (pol->g[1]) : -mpz_sgn (pol->g[1]);
    if (mpz_sgn (pol->g[0]) == s)
      return s;

    mpz_init(tmp1);
    mpz_mul_si(tmp1, pol->g[1], rel.a);
    mpz_init(tmp2);
    mpz_mul_ui(tmp2, pol->g[0], rel.b);
    mpz_add(tmp1, tmp1, tmp2);
    s = (mpz_sgn(tmp1) >= 0 ? 1 : -1);
    mpz_clear(tmp1);
    mpz_clear(tmp2);

    return s;
}

// RETURN VALUE: -1 if file exhausted.
//                0 if product is negative...
//                1 if product is positive...
//
// If check == 0, then we don't need to read purgedfile again and again
// so we store the interesing values in the vec array (ohhhhhh!).
int
treatDep(char *ratname, char *algname, char *relname, char *purgedname, char *indexname, FILE *kerfile, cado_poly pol, int nrows, int ncols, char *small_row_used, int small_nrows, hashtable_t *H, int nlimbs, char *rel_used, int *vec, int rora, int verbose, int check)
{
    FILE *ratfile, *algfile, *relfile, *indexfile, *purgedfile = NULL;
    relation_t rel;
    uint64_t w;
    int ret, i, j, nrel, r, irel, nr, sg, ind;
    char str[1024];

    memset(small_row_used, 0, small_nrows * sizeof(char));
    // now use this dep
    for(i = 0; i < nlimbs; ++i){
	ret = fscanf (kerfile, "%" SCNx64, &w);
	if(ret == -1)
	    return ret;
	ASSERT (ret == 1);
	if (verbose)
	    fprintf(stderr, "w=%" PRIx64 "\n", w);
	for(j = 0; j < 64; ++j){
	    if(w & 1UL){
		ind = (i * 64) + j;
		if(verbose)
		    fprintf(stderr, "+R_%d\n", ind);
		small_row_used[ind] = 1;
	    }
	    w >>= 1;
	}
    }
    // now map to the rels of the purged matrix
    memset(rel_used, 0, nrows * sizeof(char));
    indexfile = gzip_open(indexname, "r");
    fscanf(indexfile, "%d %d", &i, &j); // skip first line
    for(i = 0; i < small_nrows; i++){
	fscanf(indexfile, "%d", &nrel);
	for(j = 0; j < nrel; j++){
	    fscanf(indexfile, PURGE_INT_FORMAT, &r);
	    if(small_row_used[i]){
#if DEBUG >= 1
		fprintf(stderr, "# Small[%d] -> %d\n", i, r);
#endif
#if DEBUG >= 1
		if(rel_used[r])
		    fprintf(stderr, "WARNING: flipping rel_used[%d]\n", r);
#endif
		rel_used[r] ^= 1;
	    }
	}
    }
    gzip_close(indexfile, indexname);
#if MAPLE >= 1
    if((rora == 1) || (rora == 3))
	fprintf(stderr, "R2:=1; P:=1;\n");
    if((rora == 2) || (rora == 3)){
	fprintf(stderr, "A2:=1;\n");
	fprintf(stderr, "AX:=1;\n");
	fprintf(stderr, "NORM:=1;\n");
    }
#endif
    if(check)
	memset(vec, 0, ncols * sizeof(int));
    // FIXME: sg should be removed...
    sg = 1;
    ratfile = fopen(ratname, "w");
    algfile = fopen(algname, "w");
    // now really read the purged matrix in
    if(check){
	purgedfile = gzip_open(purgedname, "r");
	fgets(str, 1024, purgedfile); // skip first line
    }
    // we assume purgedfile is stored in increasing order of the indices
    // of the real relations, so that one pass in the rels file is needed...!
    relfile = gzip_open(relname, "r");
    irel = 0; // we are ready to read relation irel
    for(i = 0; i < nrows; i++){
	if(check)
	    fgets(str, 1024, purgedfile);
	if(rel_used[i]){
	    if(check)
		sscanf(str, "%d", &nr);
	    else
		nr = vec[i];
#if DEBUG >= 1
	    fprintf(stderr, "Reading in rel %d of index %d\n", i, nr);
#endif
	    if(check)
		str2Vec(vec, str);
            skip_relations_in_file(relfile, nr - irel);
            fread_relation(relfile, &rel);
	    irel = nr+1;
	    if((rora == 1) || (rora == 3))
		treatRationalRelation(H, rel);
	    if((rora == 2) || (rora == 3))
		treatAlgebraicRelation(algfile, H, rel);
	    sg *= treatSign(rel, pol);
	    clear_relation(&rel);
	}
    }
    gzip_close(relfile, relname);
    ASSERT(!check || checkVector(vec, ncols));
    if(sg == -1)
	fprintf(stderr, "prod(a-b*m) < 0\n");
    else{
	if((rora == 1) || (rora == 3))
	    finishRationalSqrt(ratfile, H, pol);
    }
    fclose(ratfile);
    fclose(algfile);
    if(check)
	gzip_close(purgedfile, purgedname);
    return (sg == -1 ? 0 : 1);
}

void
SqrtWithIndexAll(char *prefix, char *relname, char *purgedname, char *indexname, FILE *kerfile, cado_poly pol, int rora, int ndepmin, int ndepmax, int verbose, int check)
{
    FILE *indexfile, *purgedfile;
    char ratname[200], algname[200], str[1024];
    hashtable_t H;
    uint64_t w;
    int i, j, ret, nlimbs, nrows, ncols, small_nrows, small_ncols;
    char *small_row_used, *rel_used;
    int *vec; // useful to check dependency relation in the purged matrix
    int need64 = (pol->lpba > 32) || (pol->lpbr > 32);

    purgedfile = gzip_open(purgedname, "r");
    fscanf(purgedfile, "%d %d", &nrows, &ncols);

    indexfile = gzip_open(indexname, "r");
    fscanf(indexfile, "%d %d", &small_nrows, &small_ncols);
    gzip_close(indexfile, indexname);

    nlimbs = ((small_nrows - 1) / 64) + 1;
    // first read used rows in the small matrix
    small_row_used = (char *)malloc(small_nrows * sizeof(char));
    rel_used = (char *)malloc(nrows * sizeof(char));
    if(check)
	vec = (int *)malloc(ncols * sizeof(int));
    else{
	// ohhhhhhhhhhhhh!
	vec = (int *)malloc(nrows * sizeof(int));
	// get rid of end of first line
	fgets(str, 1024, purgedfile);
	for(i = 0; i < nrows; i++){
	    fgets(str, 1024, purgedfile);
	    sscanf(str, "%d", vec+i);
	}
    }
    gzip_close(purgedfile, purgedname);

    // skip first ndepmin-1 relations
    for(j = 0; j < ndepmin; j++)
	for(i = 0; i < nlimbs; ++i)
	    ret = fscanf (kerfile, "%" SCNx64, &w);

    // use a hash table to rebuild P2
    hashInit(&H, nrows, 1, need64);
    while(1){
	if(ndepmin >= ndepmax)
	    break;
	sprintf(ratname, "%s.rat.%03d", prefix, ndepmin);
	sprintf(algname, "%s.alg.%03d", prefix, ndepmin);
	ret = treatDep(ratname, algname, relname, purgedname, indexname, kerfile, pol, nrows, ncols, small_row_used, small_nrows, &H, nlimbs, rel_used, vec, rora, verbose, check);
	if(ret == -1)
	    break;
	fprintf(stderr, "# Treated dependency #%d at %2.2lf\n",
		ndepmin, seconds());
	hashClear(&H);
	ndepmin++;
    }
    hashFree(&H);
    free(small_row_used);
    free(rel_used);
    free(vec);
}

int main(int argc, char *argv[])
{
    char *relname, *purgedname, *indexname, *kername, *polyname;
    FILE *kerfile;
    cado_poly pol;
    int verbose = 1, ndepmin, ndepmax, rora, ret, i, check;

    if(argc != 10){
	fprintf(stderr, "Usage: %s relname purgedname indexname", argv[0]);
	fprintf(stderr, " kername polyname ndepmin ndepmax r|a prefix\n");
	fprintf(stderr, "Dependency relation i will be put in files prefix.i\n");
	return 0;
    }

    /* print the command line */
    fprintf (stderr, "%s.r%s", argv[0], CADO_REV);
    for (i = 1; i < argc; i++)
      fprintf (stderr, " %s", argv[i]);
    fprintf (stderr, "\n");

    relname = argv[1];
    purgedname = argv[2];
    indexname = argv[3];
    kername = argv[4];
    polyname = argv[5];
    ndepmin = atoi(argv[6]);
    ndepmax = atoi(argv[7]);

    kerfile = fopen(kername, "r");

    cado_poly_init(pol);
    ret = cado_poly_read(pol, polyname);
    ASSERT (ret);

    if(!strcmp(argv[8], "r"))
	rora = 1;
    else if(!strcmp(argv[8], "a"))
	rora = 2;
    else if(!strcmp(argv[8], "ar") || !strcmp(argv[8], "ra"))
	rora = 3;
    else
      {
        fprintf (stderr, "Error, 8th argument must be r, a, ar or ra\n");
        exit (1);
      }
    verbose = 0;
    check = 0;
    SqrtWithIndexAll(argv[9], relname, purgedname, indexname, kerfile, pol, rora, ndepmin, ndepmax, verbose, check);

    fclose(kerfile);

    return 0;
}