/*
    
               CALCUL DE BASES DE HIBERT DE RESEAUX    
               COMPUTATION OF HILBERT BASES OF LATTICES    


    Loi(tre(accent aigu)ma)c Pottier                                              
    
    Projet SAFIR (Systemes Algebriques Formels pour l'Industrie et la Recherche)
    INRIA Sophia Antipolis 
    2004 route des Lucioles
    06565 Valbonne Cedex FRANCE
    Tel: 93 65 78 19
    Fax: 93 65 77 66
    pottier@sophia.inria.fr
    
    Oct. 1996                       
 -------------------------------------------------------------------------------
This code may be compiled by command:

gcc  -I. -O2 -w -o bastat2web bastat2web.c


*/
#include <stdio.h>
#include <math.h>
/***********************************************************************/
/* 
Exposants=monomes=binomes=vecteurs 
 */
/***********************************************************************/
typedef int *exponent;
/***********************************************************************/
/* liste d'exposants  */
/***********************************************************************/
typedef struct lexp
{
    int first;
    struct lexp *rem;
}
lexp;
/* \section{Variables globales} */

int dim; /* dimension des vecteurs */
lexp *base; /* la base en construction */
lexp *baseend;
lexp *lcpair; /* les paires critiques */
lexp *lcpairend;
lexp *lcpair1; /* les paires critiques */

int ncpair_treated;
int npair_elim;
int npair_elim3;
int npair_elim4;
int change;
void baniere()
{
    printf("\n");
    printf("           CALCUL DE BASES DE HIBERT DE RESEAUX    \n");
    printf("           COMPUTATION OF HILBERT BASES OF LATTICES    \n");
    printf("-------------------------------------------------------------------------------\n");
    printf("Loi(tre(accent aigu)ma)c Pottier                                              \n");
    printf("\n");
    printf("Projet SAFIR (Systemes Algebriques Formels pour l'Industrie et la Recherche)\n");
    printf("INRIA Sophia Antipolis \n");
    printf("2004 route des Lucioles\n");
    printf("06565 Valbonne Cedex FRANCE\n");
    printf("Tel: 93 65 78 19\n");
    printf("Fax: 93 65 77 66\n");
    printf("pottier@sophia.inria.fr\n");
    printf("-------------------------------------------------------------------------------\n");
    printf("\n");
    printf("                         Oct. 1996                       \n");
    printf("\n");
    
}

/*  Min de deux entiers */
int min_int(x,y)
     int x,y;
{
    if (x>y)
    return(y);
    else return(x);
}

/*  Max de deux entiers */
int max_int(x,y)
     int x,y;
{
    if (x>y)
    return(x);
    else return(y);
}

/* Copie */
exponent copy_exp(e)
     exponent e;
{
    int i;
    exponent e1;
    e1= (exponent)malloc(sizeof(int)*dim);
    for  (i=0 ; i<dim ; i++)
    {
	e1[i]=e[i];
    }
    return(e1);
}
/*  Impression  */
print_exp(e)
     exponent e;
{
    int i;
    printf("%s","[");
    for  (i=0 ; i<(dim -1); i++)
    {
	printf("%d ",e[i]);
    }
    printf("%d]\n",e[i]);
}

/* Op\'erations sur les exposants.  */
exponent opp_exp(e1)
     exponent e1;
{
    int i;
    for  (i=0 ; i<dim ; i++)
    {
	e1[i]=-e1[i];
    }
    return(e1);
}

equal_exp(e1,e2)
     exponent e1,e2;
{
    int i;
    for  (i=0 ; i<dim ; i++)
    {
	if (!(e1[i]==e2[i]))
	return(0);
    }
    return(1);
}

exponent add_exp(e1,e2)
     exponent e1,e2;
{
    int i;
    for  (i=0 ; i<dim ; i++)
    {
	e1[i]=e1[i]+e2[i];
    }
    return(e1);
}

exponent sub_exp(e1,e2)
     exponent e1,e2;
{
    int i;
    for  (i=0 ; i<dim ; i++)
    {
	e1[i]=e1[i]-e2[i];
    }
    return(e1);
}
/*  Soustrait q fois e2 de e1.  */
exponent subq_exp(e1,e2,q)
     exponent e1,e2;
     int q;
{
    int i;
    if (!(q==0))
    for  (i=0 ; i<dim ; i++)
    {
	e1[i]=e1[i]-q*e2[i];
    }
    return(e1);
}

trivial_exp(e)
     exponent e;
{
    int i;
    for  (i=0 ; i<dim ; i++)
    {
	if (!(e[i]==0))
	return(0);
    }
    return(1);
}

normalize_exp(exponent e)
{
    int i;
    i=0;
    while (i<dim && e[i]==0)
	i++;
    if (e[i]<0)
	opp_exp(e);
    return e;
}
/* \section{Fonctions sur les listes d'exposants} */

/* Cr\'eation */

lexp *make_lexp(exponent b,lexp *rem)
{
    lexp *l;
    l=(lexp *)malloc(sizeof(lexp));
    l->first = b;
    l->rem = rem;
    return(l);
}

/* Impression des listes de bin\^omes */

print_lexp(s)
     lexp *s;
{
    printf("(\n");
    print_lexp_rec(s);
    printf(")\n");
}
  
print_lexp_rec(s)
     lexp *s;
{
    if (s)
    {
	print_exp(s->first);
	print_lexp_rec(s->rem);
    }
}

length_lexp(lcp)
     lexp *lcp;
{
    int res;
    res=0;
    while (lcp)
    {
	res++;
	lcp=lcp->rem;
    }
    return(res);
}

/* Divisions. */


/*  Quotient de e1 par e2.
 rend q max tel que e1+ - q*e2+ >= 0 et e1- + -q*e2- >= 0 
*/

int quo_exp_pos (e1,e2)
     exponent e1,e2;
{
    int i,q,qmin;
    qmin=-1;        
    for  (i=0 ; i<dim ; i++)
    {
	if (e2[i]>0 && e1[i]<e2[i]) return(0);
	if (e2[i]<0 && e1[i]>e2[i]) return(0);
	if (e2[i]!= 0)
	    {
		q=e1[i]/e2[i];
		if (qmin==-1) qmin=q; else qmin=min_int(qmin,q);
	    }
    }
    return max_int(qmin,0);
}
int quo_exp_neg (e1,e2)
     exponent e1,e2;
{
    int i,q,qmin;
    qmin=-1;        
    for  (i=0 ; i<dim ; i++)
    {
	if (e2[i]<0 && e1[i]<-e2[i]) return(0);
	if (e2[i]>0 && e1[i]>-e2[i]) return(0);
	if (e2[i]!= 0)
	    {
		q=e1[i]/(-e2[i]);
		if (qmin==-1) qmin=q; else qmin=min_int(qmin,q);
	    }
    }
    return -max_int(qmin,0);
}

/* Physique */
rem_exp_lexp(exponent v,lexp *lv)
{
    exponent v1;
    int q;
    change=0;
    while (lv)
	{
	    v1=lv->first;
	    lv=lv->rem;
	    q=quo_exp_pos(v,v1);
	    if (q==0) q=quo_exp_neg(v,v1);
	    if (q!=0)
		{
		    change=1;
		    subq_exp(v,v1,q);
		}
	}
}
/* Physique */
rem_exp_exp(exponent v,exponent v1)
{
    int q;
    q=quo_exp_pos(v,v1);
    if (q==0) q=quo_exp_neg(v,v1);
    if (q!=0) subq_exp(v,v1,q);
}

/* Paires critiques. */

same_orthant(exponent e1,exponent e2)
{
    int i;
    for(i=0;i<dim;i++)
	{
	    if (e1[i]>0 && e2[i]<0) return 0;
	    if (e1[i]<0 && e2[i]>0) return 0;
	}
    return 1;
}
opp_orthant(exponent e1,exponent e2)
{
    int i;
    for(i=0;i<dim;i++)
	{
	    if (e1[i]<0 && e2[i]<0) return 0;
	    if (e1[i]>0 && e2[i]>0) return 0;
	}
    return 1;
}
member_lexp(exponent b,lexp *l)
{
    while (l)
    {
	if (equal_exp(b,l->first)) return(1);
	l=l->rem;
    }
    return(0);
}

void add_cpair(exponent cp)
{
    /* ESSAI 1 qui a l'air concluant (gain de 5%):
       on reduit la paire par la base avant de l'ajouter */
    rem_exp_lexp(cp,base);
    /* ESSAI 2 non concluant: on la reduit aussi par les paires existantes */
    /* if (!trivial_exp(cp))  rem_exp_lexp(cp,base);*/
    if (trivial_exp(cp))
	npair_elim3++;
    else
	{
	    if (lcpair==0)
		{
		    lcpair=make_lexp(cp,0);
		    lcpairend=lcpair;
		}
	    else
		{
		    lcpairend->rem=make_lexp(cp,0);
		    lcpairend=lcpairend->rem;
		}
	}
}

void add_cpairs(exponent v1, exponent v2)
{
    lexp *lv;
    exponent v3;
    int elimp,elimm;
    lv=base;
    elimp=same_orthant(v1,v2);
    elimm=opp_orthant(v1,v2);
    while (lv && (elimp==0 || elimm==0))
	{
	    v3=lv->first;
	    lv=lv->rem;
	    if (v3!=v1 && v3!=v2)
		{
		    if (elimm==0 &&
			(betweenm1(v1,v2,v3) || betweenm2(v1,v2,v3)))
			{
			    npair_elim++;
			    elimm=1;
			}
		    if (elimp==0 &&
			(betweenp1(v1,v2,v3) || betweenp2(v1,v2,v3)))
			{
			    npair_elim++;
			    elimp=1;
			}
		}
	}
    if (elimp==0)
	add_cpair(add_exp(copy_exp(v1),v2));
    if (elimm==0)
	add_cpair(sub_exp(copy_exp(v1),v2));
}
  
betweenm1(exponent v1, exponent v2, exponent v3)
{
    int i;
    for (i=0;i<dim;i++)
	    if (v3[i]>v1[i])
		{
		    if (v3[i]>v2[i]) return 0;
		}
	    else
		if (v3[i]!=v1[i] && v3[i]<v2[i]) return 0;
    return 1;
}
betweenm2(exponent v1, exponent v2, exponent v3)
{
    int i;
    for (i=0;i<dim;i++)
	if (v3[i]>-v1[i])
	   {
	       if (v3[i]>-v2[i]) return 0;
	   }
        else
	   if (v3[i]!=-v1[i] && v3[i]<-v2[i]) return 0;
    return 1;
}
betweenp1(exponent v1, exponent v2, exponent v3)
{
    int i;
    for (i=0;i<dim;i++)
	if (v3[i]>v1[i])
	   {
	       if (v3[i]>-v2[i]) return 0;
	   }
        else
	   if (v3[i]!=v1[i] && v3[i]<-v2[i]) return 0;
    return 1;
}
betweenp2(exponent v1, exponent v2, exponent v3)
{
    int i;
    for (i=0;i<dim;i++)
	if (v3[i]>-v1[i])
	   {
	       if (v3[i]>v2[i]) return 0;
	   }
        else
	   if (v3[i]!=-v1[i] && v3[i]<v2[i]) return 0;
    return 1;
}

void lcpair_exp_lexp(exponent v,lexp *l)
{
    while (l)
	{
	    add_cpairs(v,l->first);
	    l=l->rem;
	}
}
void lcpair_lexp(lexp *l)
{
    while (l)
	{
	    lcpair_exp_lexp(l->first,l->rem);
	    l=l->rem;
	}
}
/************ ESSAI3 non concluant:
ajout d'un vecteur a la base:
calcul de toutes les paires (seulement le critere 1), interreduction entre elles, et ajout aux anciennes paires.
*/
void lcpair_exp_lexp1(exponent v,lexp *l)
{
    lcpair1=0;
    while (l)
	{
	    add_cpairs1(v,l->first);
	    l=l->rem;
	}
    lcpair1=interred2(0,lcpair1);
    lcpair=append_lexp(lcpair,lcpair1);
    lcpairend=lcpair;
    while (lcpairend->rem) lcpairend=lcpairend->rem;
}
void add_cpairs1(exponent v1, exponent v2)
{
    if (same_orthant(v1,v2)==0)
	lcpair1=make_lexp(add_exp(copy_exp(v1),v2),lcpair1);
    if (opp_orthant(v1,v2)==0)
	lcpair1=make_lexp(sub_exp(copy_exp(v1),v2),lcpair1);
}  
/*******************************/

/* Tri lex decroissant d'une liste d'exposants*/

lexp *qsort_lexp(l)
     lexp *l;
{
    int x,y;
    lexp *l1,*l2;
    l1=0;
    l2=0;
    if (length_lexp(l)<2)
    return(l);
    else
    {
	x=l->first;
	l=l->rem;
	while (l)
	{
	    y=l->first;
	    l=l->rem;
	    if (lex_exp(x,y)==-1)
	    l1=make_lexp(y,l1);
	    else
	    l2=make_lexp(y,l2);
	}
	return(append_lexp(qsort_lexp(l1),
			    make_lexp(x,qsort_lexp(l2))));
    }
}
lexp *append_lexp(l1,l2)
     lexp *l1,*l2;
{
    lexp *laux;
    laux=l1;
    if (l1==0)
    return(l2);
    else
    {
	while (laux->rem != 0)
	{
	    laux=laux->rem;
	}
	laux->rem=l2;
	return(l1);
    }
}

lexp *append1_lexp(lexp *l1, exponent v)
{
    lexp *laux;
    laux=l1;
    if (l1==0)
    return(make_lexp(v,0));
    else
    {
	while (laux->rem != 0)
	{
	    laux=laux->rem;
	}
	laux->rem=make_lexp(v,0);
	return(l1);
    }
}
/* Ordre lexicographique sur les parties positives: rend -1 0 ou 1  */

int lex_exp(e1,e2)
     exponent e1,e2;
{
    int k;
    k=0;
    while (k<dim)
    {
	if (e1[k]>e2[k])
	return(1);
	else
	{
	    if (e1[k]<e2[k])
	    return(-1);
	    else k++;
	}
    }
    return(0);
}
/* interreduction: lv1 est deja interreduite*/
interred2(lexp *l1, lexp *l2)
{
    exponent v1,v2;
    lexp *l11,*l12,*lv2;
    while (l2)
	{
	    v2=l2->first;
	    l2=l2->rem;
	    rem_exp_lexp(v2,l1);
	    if (!trivial_exp(v2))
		{
		    normalize_exp(v2);
		    l11=0;
		    l12=0;
		    lv2=make_lexp(v2,0);
		    while (l1)
			{
			    v1=l1->first;
			    l1=l1->rem;
			    change=0;
			    rem_exp_lexp(v1,lv2);
			    if (change==1)
				l12=make_lexp(v1,l12);
			    else
				l11=make_lexp(v1,l11);
			}
		    l1=append_lexp(l11,lv2);
		    l2=append_lexp(l2,l12);
		}
	}
    return l1;
}    
    
add_base(exponent v)
{
    baseend->rem=make_lexp(v,0);
    baseend=baseend->rem;
}
/* \subsection{Algorithme de completion.} */
/* interreduction tous les 100 nouveaux vecteurs */
lexp *basta()
{
    exponent cp;
    lexp *laux,*laux1;
    int i,new_vect;
    base=interred2(0,base);
    base=qsort_lexp(base);    
    baseend=base; while (baseend->rem !=0) baseend=baseend->rem;
    lcpair=0;
    lcpairend=0;
    lcpair_lexp(base);
    ncpair_treated=0;
    npair_elim=0;
    npair_elim3=0;
    new_vect=0;
    while (lcpair)
    {
/*	print_lexp(base);
	print_lexp(lcpair);*/
	cp=lcpair->first;
	lcpair=lcpair->rem;
	ncpair_treated++;
	/* ESSAI 1.1: concluant (10%), les paires de lcpair sont reduites 
	   pour la base, et non nulles
	   rem_exp_lexp(cp,base);
	   if (!trivial_exp(cp))*/
	{
		printf("[%d,%d]",length_lexp(base),length_lexp(lcpair));
		new_vect++;
		normalize_exp(cp);
		print_exp(cp);
		/* ESSAI 4:(???) interreduction tous les 100 vecteurs ajoutes */
		if /*(new_vect%100==0)*/ (0)
		    {
			printf("  INTERREDUCTION  ");
			add_base(cp);
			base=interred2(0,base);
			base=qsort_lexp(base);
			baseend=base;
			while (baseend->rem !=0) baseend=baseend->rem;
			lcpair=0;
			lcpairend=0;
			lcpair_lexp(base);
		    }
		else
		    {
			lcpair_exp_lexp(cp,base);
			add_base(cp);
	      	/* ESSAI 1.1: les paires de lcpair sont reduites pour la base,
		 il suffit de les reduire par le nouveau vecteur,
		 et de virer les nulles */
			laux=make_lexp(0,lcpair);
			laux1=laux;
			while (laux->rem)
			    {
				rem_exp_exp(laux->rem->first,cp);
				if (trivial_exp(laux->rem->first))
				    {
					laux->rem=laux->rem->rem;
					npair_elim3++;
				    }
				else
				    laux=laux->rem;
			    }
			lcpair=laux1->rem;
			if (lcpair)
			    {
				lcpairend=lcpair;
				while (lcpairend->rem) lcpairend=lcpairend->rem;
			    }
		    }
	}
    }
    base=interred2(0,base);
    base=qsort_lexp(base);
}
/* Lancement */
main()
{
    int i,rank,r;
    exponent b;

    baniere();
    printf("Dimension des vecteurs (Dimension of vectors):\n");
    scanf("%d",&dim);
    printf("Nombre de generateurs(Number of generators):\n");
    scanf("%d",&rank);
    printf("Donnez les generateurs, une ligne par generateur (Give generators, one row for each generator):\n");
    base=0;
    for (r=1; r<=rank; r++)
    {
	b=(exponent)malloc(sizeof(int)*dim);
	for (i=0;i<dim; i++)
	{
	    scanf("%d",&b[i]);
	}
	base=make_lexp(normalize_exp(b),base);
    }
    basta();
    printf("\nBase de Hilbert(Hilbert basis):\n");
    print_lexp(base);
    printf("Nombre de vecteurs dans la base (Number of vectors in the basis): %d\n", length_lexp(base));
}