/* Pascal Molin
 * November 2008
 * Build number fields extensions using linear programming.
 * Needs tables of small polynomials (see Pari site)
 *
 * Copyright Pascal Molin
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 *    1. Redistributions of source code must retain the above
 *    copyright notice, this list of conditions and the following disclaimer.
 *
 *    2. Redistributions in binary form must reproduce the above
 *    copyright notice, this list of conditions and the following
 *    disclaimer in the documentation and/or other materials
 *    provided with the distribution.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <glpk.h>
#include "pari/pari.h"
#include "pari/paripriv.h"

ulong VERBOSE = 0;
int RANDOMSEED = 0;
int ERROR = 10; /* percentage of error allowed in mip objective */
int NOCYCLOTOMIC = 1; /* forbid cyclotomic field */
int DISCOBJ = 1; /* use discriminant in objective */
char* obj[]={"degree","degree+1/disc","degree+tanh(disc/100)","log(disc)"};

typedef struct {
  int m; /* number of roots of unity */
  GEN a; /* a table of decomposition of order in residues class */
  int *nbq; /* number of rows taken by residue class */
  int nrows;
  int *poldegrees; /* degrees of polynomials to be considered */
  int primebound;
} badfield;

void
init_struct(badfield *F) {
  int i,j,powi,nrows,o,p,a;
  byteptr pt = diffptr;
  /* init tables */
  int *t = F->poldegrees;
  for(i=0;i<F->m;i++) { *t++ = 0; }
  /* compute residue classes */
  for(i=2;i<=F->m-1;i++) {
    if(cgcd(i,F->m) != 1) continue;
    /* find order */
    powi=i;
    for(o=1;powi != 1;o++) {
      powi = powi*i % F->m;
    }
    gel(F->a,i) = factoru(o); /* factorisation of order */
    F->nbq[i] = lg(gel(gel(F->a,i),1))-1; /* number of rows occupied */
    /* add necessary degrees to poldegrees */
    for(j=1; j <= F->nbq[i]; j++) {
      int p,e;
      p  = (int)gel(gel(gel(F->a,i),1),j);
      e = (int)gel(gel(gel(F->a,i),2),j);
      if (e > F->poldegrees[p]) F->poldegrees[p] = e; 
    }
  }
  /* compute the number of rows */
  nrows = 0;
  for(p=0; p < F->primebound;) {
    NEXT_PRIME_VIADIFF_CHECK(p,pt);
    if(!(F->m % p) || p % F->m == 1) continue; /* if ramified or of order 1 */
    a = p % F->m; /* residue class */
    nrows += F->nbq[a];
  }
  F->nrows = nrows;
}

/* start mip structure */
void
fill_mip(glp_prob *lp, badfield *F) {
  int i,p,a,k;
  byteptr pt = diffptr;
  glp_set_prob_name(lp, "bad field");
  glp_set_obj_name(lp, "degre");
  glp_set_obj_dir(lp, GLP_MIN);
  glp_add_rows(lp,F->nrows);
  /* put row names */
  i = 0;
  for(p=0; p < F->primebound;) {
    NEXT_PRIME_VIADIFF_CHECK(p,pt);
    if(!(F->m % p) || p % F->m == 1) continue; /* if ramified or of order 1 */
    a = p % F->m;
    for(k=1;k<=F->nbq[a];k++) {
      i++;
      char *s = pari_sprintf("p=%i;a=%i;q=%i",p,a,gel(gel(gel(F->a,a),1),k));
      /* printf("ligne %i = %s\n",i,s); */
      glp_set_row_name(lp, i, s);
      glp_set_row_bnds(lp, i, GLP_LO, 1, 1);
    }
  }
  /* to avoid a cyclotomic field, one can add a false condition */
  if(NOCYCLOTOMIC) {
    /* replace last condition */
    glp_set_row_name(lp, i, "bad inertia condition");
    glp_set_row_bnds(lp, i, GLP_UP, 0, 0);
  }
}

/* add a column to the structure corresponding to polynomial T.
 * compute decomposition aver all primes */
void
add_pol_mip(glp_prob *lp, badfield *F, GEN T) {
  int j; /* column in the matrix */
  int i; /* row */
  char *s;
  int ind[F->nrows]; /* indexes of values */
  double val[F->nrows]; /* 1 values */
  int len=0; /* number of non-zero values */
  int a,k,p,deg;
  byteptr pt = diffptr;
  double termobj, disc = (double)itou(absi(RgX_disc(T)));
  /* pari_printf("adding pol Q=%Ps\n",T); */
  j = glp_add_cols(lp,1); /* index of the column */
  s = pari_sprintf("%Ps",T);
  glp_set_col_name(lp,j,s);
  free(s);
  /* type binary */
  glp_set_col_kind(lp, j, GLP_BV);
  /* coefficient for objective :
   * minimise degree and disc */
  deg = degpol(T);
  switch(DISCOBJ) {
    case 1:
      termobj = deg-1/disc;
      break;
    case 2 :
      termobj = deg+2*tanh(disc/100);
      break;
    case 3 :
      termobj = log(disc);
      break;
    default:
      termobj = deg;
  }
  glp_set_obj_coef(lp,j,termobj);
  /* fill the column */
  p = 2;
  i = 0;
  for(p=0; p < F->primebound;) {
    GEN dec, q, e;
    int gcdf;
    NEXT_PRIME_VIADIFF_CHECK(p,pt);
    if( F->m % p == 0 || p % F->m == 1) continue; /* if ramified or of order 1 */
    a = p % F->m; /* residue class */
    dec = gel(FpX_degfact(T,utoi(p)),1);
    /* find gcd */
    gcdf = 0;
    for(k=1;k<lg(dec);k++)
    {
      gcdf = gcdf ? dec[k] : cgcd(gcdf,dec[k]);
      if(gcdf == 1) break;
    }
    q = gel(gel(F->a,a),1); /* prime decomposition of order */
    e = gel(gel(F->a,a),2); /* exponents */
    for(k=1; k <= F->nbq[a]; k++) {
      i++; /* index of the current row */
      /* put coeff if allows order on q */
      if(u_lval(gcdf,(ulong)gel(q,k)) >= (long)gel(e,k)) {
        /* this polynomial gives enough inertia on the q-part for prime p */
        ind[++len] = i;
        val[len] = 1;
      }
    }
  }
  glp_set_mat_col(lp,j,len,ind,val);
}

/* print stats about lp */
void
print_stats(glp_prob *lp) {
  int nr, nc, i;
  nr=glp_get_num_rows(lp);
  nc=glp_get_num_cols(lp);
  printf("number of lines x columns : %i x %i\n",nr,nc);
  for(i=1;i<=nr;i++) {
    printf("line %i : inf = %f ; sup = %f\n",i,
        glp_get_row_lb(lp,i),
        glp_get_row_ub(lp,i));
  }
}

void
output_col(glp_prob *lp, int j) {
  int ind[glp_get_num_rows(lp)];
  double val[glp_get_num_rows(lp)];
  int i,len;
  /* printf("  objective coefficient : %2.0f\n",glp_get_obj_coef(lp,j)); */
  len = glp_get_mat_col(lp,j,ind,val);
  for(i=1;i<=len;i++) {
    printf(" %s  : >= %2.0e, val = %2.0f\n",glp_get_row_name(lp,ind[i]),
        glp_get_row_lb(lp,ind[i]),val[i]);
  }
}

/* retrieve a solution */
void
retrieve_sol(glp_prob *lp) {
  double z,x;
  int j,nofirst=0;
  z = glp_mip_obj_val(lp);
  printf("=== better integer solution found : obj = %.2f\n",z);
  printf("[");
  for(j=1; j<= glp_get_num_cols(lp); j++) {
    x = glp_mip_col_val(lp, j);
    if (x) {
      printf(nofirst++ ? ",%s" : "%s",glp_get_col_name(lp,j));
      if (VERBOSE) output_col(lp,j);
    };
  }
  printf("]\n");
}

/* trap mip execution when better integer solution found */
void
trap_intsols(glp_tree *tree, void *info) {
  switch (glp_ios_reason(tree)) {
    case GLP_IBINGO:
      retrieve_sol(glp_ios_get_prob(tree));
      break;
    default: /* ignore call for other reasons */
      break;
  }
}

void
mipsolve(int m, int primebound, int npols, int tcoeffs, char* path) {
  badfield F;
  int nbq[m];
  int poldegrees[m];
  int j, ncols, e, p;
  long v;
  GEN randmax;
  pari_sp av;
  glp_prob *lp;
  glp_iocp parm;
  pari_init(100000000,2);
  lp = glp_create_prob();
  /* initialise with default control parameters */
  glp_init_iocp(&parm);
  F.nbq = nbq; F.poldegrees = poldegrees;
  F.m = m; F.primebound = primebound;
  F.a=cgetg(m,t_VEC);
  init_struct(&F);
  fill_mip(lp,&F);
  /* construct random polynomials */
  if (RANDOMSEED) setrand(utoi(RANDOMSEED));
  randmax = utoi(tcoeffs);
  pari_var_init();
  v = fetch_user_var("x");
  av=avma;
  /* select primes needed as degree of polynomial */
  for(p=0; p < m; p++) {
    if(!poldegrees[p]) continue;
    /* and exponents */
    for(e=1; e <= poldegrees[p]; e++) {
      /* random polynomial of degree p^e */
      int d = upowuu(p,e);
      if(d<8) { /* use pols in precomputed tables */
        int sig = d%2; /* minimum number of real places */
        for(;sig <=d;sig+=2) {
          avma = av;
          char *s=pari_sprintf("%s/T%i%i.gp",path,d,sig);
          GEN vect = gp_readvec_file(s);
          long nbj = glength(vect) < npols/d ? glength(vect) : npols/d;
          pari_sp av2=avma;
          for(j=1; j < nbj; j++) {
            GEN T;
            avma=av2;
            T = gtopoly(gel(gel(vect,j),2),v);
            add_pol_mip(lp,&F,T);
            ncols++;
          }
        }
      }
      else { /* no precomputed fields for this degree */
        /* number of such pols desired */
        for(j=0; j < npols; j++) {
          avma=av;
          d=d+2; /* size of cgetg for degree d polynomial */
          GEN T = cgetg(d, t_VEC);
          for(;;) {
            int coeff;
            for(coeff=1;coeff<d;coeff++) {
              gel(T,coeff)= subii(randomi(randmax),gmul2n(randmax,-1));
            };
            gel(T,d-1) = gen_1;
            T = gtopolyrev(T,v);
            if(gisirreducible(T)) break;
          }
          add_pol_mip(lp,&F,T);
          ncols ++;
        }
      }
    }
  }
  glp_simplex(lp, NULL);
  /* parm.presolve = GLP_ON; use mip presolver */
  parm.out_frq = 60000; /* print info only every minute */
  parm.mip_gap = (double)ERROR / 100.0 ;
  parm.cb_func = trap_intsols;
  glp_intopt(lp, &parm);
  retrieve_sol(lp);
  glp_delete_prob(lp);
}

  int
main(int argc, char **argv)
{
  /* default values */
  int m = 18, primebound = 100, npols = 1000, tcoeffs = 10;
  char path[100] = "nftables"; /* path for nftables */
  int c;
  extern char *optarg;
  extern int optind;
  extern int opterr;
  extern int optopt;
  opterr=0;
  while ((c = getopt(argc, argv, "vcd:m:p:n:t:r:e:I:")) != EOF) {
    switch (c) {
      case 'v': /* verbosity */
        VERBOSE = 1;
        break;
      case 'm': /* mimics m-th roots of unity */
        sscanf(optarg,"%i",&m);
        break;
      case 'p': /* up to p=primebound */
        sscanf(optarg,"%i",&primebound);
        break;
      case 'n': /* number of polynomials used */
        sscanf(optarg,"%i",&npols);
        break;
      case 't': /* size of random coefficients */
        sscanf(optarg,"%i",&tcoeffs);
        break;
      case 'r': /* random seed */
        sscanf(optarg,"%i",&RANDOMSEED);
        break;
      case 'c': /* allow cyclotomic fields */
        NOCYCLOTOMIC = 0;
        break;
      case 'e': /* error allowed in % */
        sscanf(optarg,"%i",&ERROR); 
        break;
      case 'd': /* choose objective function */
        /* DISCOBJ = 0; */
        sscanf(optarg,"%i",&DISCOBJ); 
        break;
      case 'I' : /* path where the tables of number fields are to be found */
        sscanf(optarg,"%s",path);
        break;
      case '?':
        { (void)fprintf(stderr,
            "usage: findpol_mip [-cvdm][-m roots][-p primebound][-n number of pols]\n\
            [-t size of random coefficients][-r random seed][-e error]\n\
            with -c : can output cyclotomic field\n\
            -v : verbose\n\
            -e n : allow n%% error from optimality\n\
            -n n : use n pols in list (default 1000)\n\
            -I path : path of PARI nf tables\n\
            -d : choose ojective\n\
            0 = only degree\n\
            1 = first degree, then disc (deg-1/disc)\n\
            2 = first degree, then disc (deg+tanh(disc/100))\n\
            3 = disc only (log(disc))\n\
            ");
        exit (2);
        } 
    }
  }
  printf("Searching m=%i, primebound=%i, npols=%i, tcoeffs=%i\n",m,primebound,npols,tcoeffs);
  printf("Objective %i : %s\n",DISCOBJ,obj[DISCOBJ]);
  mipsolve(m,primebound,npols,tcoeffs,path);
  return 0;
}