(* ::Package:: *)

BeginPackage["OmegaPrimality`"];


Authors="J. I. Garc\[IAcute]a Garc\[IAcute]a and A. Vigneron Tenorio,\nDpto. Matem\[AAcute]ticas, 
Universidad de C\[AAcute]diz
ignacio.garcia@uca.es, alberto.vigneron@uca.es\n";

Commands="The package defines the commands: 
OmegaPrimalityOfElemAffSG, OmegaPrimalityOfAffSG, OmegaPrimalityOfElemSG, OmegaPrimalityOfSG.
Try command::usage for more help";

OmegaPrimalityOfElemAffSG::usage="OmegaPrimalityOfElemAffSG[a,sg], 
if a belongs to the affine semigroup generated by sg, 
it computes the \[Omega]-primality of the element.
Example: 
	OmegaPrimalityOfElemAffSG[{154, 118}, {{5, 3}, {5, 11}, {2, 7}, {11, 4}}] 
computes the \[Omega]-primality of the element {154,118} in the affine semigroup generated by {5, 3}, {5, 11}, {2, 7}, {11, 4}.
Example: 
	OmegaPrimalityOfElemAffSG[{10000},{{115},{212},{333},{571}}]
computes the \[Omega]-primality of 10000 in the numerical semigroup generated by 115, 212, 333, 571.
Example: 
	OmegaPrimalityOfElemAffSG[{41},{{5},{12},{19}}]
computes the \[Omega]-primality of 41 in the numerical semigroup generated by 5, 12, 19.";

OmegaPrimalityOfAffSG::usage="OmegaPrimalityOfAffSG[sg], 
it returns the \[Omega]-primality of the affine semigroup generated by sg.
Example: 
	OmegaPrimalityOfAffSG[{{5, 3}, {5, 11}, {2, 7}, {11, 4}}] 
computes the \[Omega]-primality of the affine semigroup generated by {5, 3}, {5, 11}, {2, 7}, {11, 4}.
Example: 
	OmegaPrimalityOfAffSG[{{115},{212},{333},{571}}] 
computes de \[Omega]-primality of the semigroup numerical semigroup generated by 115, 212, 333, 571.
Example: 
	OmegaPrimalityOfAffSG[{{5},{12},{19}}]
";

relatorsAffSG::usage="relatorsAffSG[sg] where sg is a system of generators of the affine semigroup 
computes a set of relators the semigroup generated by sg.
Example: 
	relatorsAffSG[{{5, 3}, {5, 11}, {2, 7}, {11, 4}}]
Example: 
	relatorsAffSG[{{115},{212},{333},{571}}]
Example: 
	relatorsAffSG[{{5},{12},{19}}]";

OmegaPrimalityOfElemSG::usage="OmegaPrimalityOfElemSG[\[Alpha],\[Rho]] where \[Alpha] is a element of the semigroup 
and \[Rho] a presentation. Computes the \[Omega]-primality of \[Alpha] in the semigroup defined by \[Rho].
Example: 
	r1=relatorsAffSG[{{5, 3}, {5, 11}, {2, 7}, {11, 4}}];
	OmegaPrimalityOfElemSG[{3,5,2,10},r1]
Example:
	r2=relatorsAffSG[{{115},{212},{333},{571}}];
	OmegaPrimalityOfElemSG[{73,6,1,0},r2]
Example:
	r3=relatorsAffSG[{{5},{12},{19}}];
	OmegaPrimalityOfElemSG[{2,1,1},r3]";

OmegaPrimalityOfSG::usage="OmegaPrimalityOfSG[\[Rho]], where \[Rho] is a presentation of the semigroup.
Example:
	r1=relatorsAffSG[{{5, 3}, {5, 11}, {2, 7}, {11, 4}}];
	OmegaPrimalityOfSG[r1]
Example:
	r2=relatorsAffSG[{{115},{212},{333},{571}}];
	OmegaPrimalityOfSG[r2]
Example:
	r3=relatorsAffSG[{{5},{12},{19}}];
	OmegaPrimalityOfSG[r3]";

(*x::usage="";*)

Print[Authors<>"\n"<>Commands];


Begin["`Private`"];


parABinomio[l_]:=Module[{i,pi=1,pd=1,dim=Length[l[[1]]]},
For[i=1,i<=dim,i++,
pi=pi Subscript[x, i]^l[[1,i]];
pd=pd Subscript[x, i]^l[[2,i]];
];
Return[pi-pd];
];

monomioAPunto[monomio_,numIndet_]:=Table[Exponent[monomio,Subscript[x, i]],{i,1,numIndet}];

puntoAMonomio[l_]:=Module[{i,p=1},
For[i=1,i<=Length[l],i++,
p=p Subscript[x, i]^l[[i]]
];
Return[p];
];

listaVars[n_]:=Table[Subscript[x, i],{i,n}];

listaVars[n_,t_]:=Table[Subscript[t, i],{i,n}];

binomio2Par[bin_,dim_]:=Module[{aux},
aux=MonomialList[bin,listaVars[dim],"DegreeReverseLexicographic"];
Return[monomioAPunto[#,dim]& /@ aux];
];
ReesElements[m_,nvars_]:=Module[{lbinomios={},i},
For[i=1,i<=nvars,i++,
lbinomios=lbinomios~Join~{m-m Subscript[x, i]}
];
Return[lbinomios];
];

esMinimal[m_,bg_,nv_,mu_]:=Module[{maux},
For[i=1,i<=nv,i++,
If[Exponent[m,Subscript[x, i]]>=1,
maux=PolynomialReduce[m/Subscript[x, i],bg,listaVars[nv],MonomialOrder->DegreeReverseLexicographic][[2]];
If[maux==mu,
Return[False];
];
];
];
Return[True];
];


Algoritmo16[bgDel4ti2_,monomio_,nv_]:=Module[{mu,A={},Aaux={},B={},i,alpha,beta,C={},bgRees,rees,pares,sistgen},
rees=ReesElements[monomio,nv];
sistgen=rees~Join~bgDel4ti2;
bgRees= GroebnerBasis[sistgen,listaVars[nv ],MonomialOrder->DegreeReverseLexicographic] ;
pares=MonomialList[#,listaVars[nv],DegreeReverseLexicographic] & [bgRees];
mu=PolynomialReduce[monomio,bgRees,listaVars[nv],MonomialOrder->DegreeReverseLexicographic][[2]];
Aaux={mu};
While[A!=Aaux,
(*Print["."];*)
A=Aaux;
B={};
For[i=1,i<=Length[A],i++,
For[j=1,j<=Length[bgRees],j++,
{alpha,beta}=pares[[j]];
B=B\[Union]{PolynomialLCM[A[[i]],beta]*alpha/beta};
];
];
B=GroebnerBasis[B,listaVars[nv],MonomialOrder->DegreeReverseLexicographic];
C=Select[B,(esMinimal[#,bgRees,nv,mu])&];
Aaux=A \[Union] C;
];
Print["Length of output of Alg16=",Length[A]];
Return[A];
];


baseDeSGAfin[l_]:=Module[{dim,bg={},i,j,aux,lu={},base},
dim=Length[l[[1]]];
For[i=1,i<=Length[l],i++,
lu=lu\[Union]{Subscript[u, i]};
aux=1;
For[j=1,j<=dim,j++,aux=aux*Subscript[u, j]^l[[i,j]]];
bg=bg\[Union]{aux-Subscript[x, i]};
];
base=GroebnerBasis[ bg, listaVars[dim]\[Union]lu, lu];
Return[base];
];


relatorsAffSG[l_]:=binomio2Par[#,Length[l]]& /@ baseDeSGAfin[l];


OmegaPrimalityOfElemSG[x_,rel_]:=Module[{bg,ng=Length[x],aux},
bg=parABinomio /@ rel;
monomio=puntoAMonomio[x];
aux=monomioAPunto[#,ng]&  /@ Algoritmo16[bg,monomio,ng];
omega=Max[(Plus @@ # )& /@ aux];
Return[omega];
];

OmegaPrimalityOfSG[rel_]:=Module[{ng=Length[rel[[1,1]]],aux},
bg=parABinomio /@ rel;
aux=Max[(Plus @@ # )& /@ (monomioAPunto[#,ng]&  /@ Algoritmo16[bg,Subscript[x, #],ng])] & /@ Table[i,{i,ng}];
Print["\[Omega]-primalities of the generators: ",aux];
Return[Max[aux]];
];


OmegaPrimalityOfElemAffSG[x_,sistGen_]:=Module[{ng=Length[sistGen],aux,bg,omega,sol},
aux=FindInstance[
(x==listaVars[ng,t].sistGen)&& (And @@ ((#>=0)&/@listaVars[4,t])),listaVars[4,t],Integers
];
If[aux!={},
sol=listaVars[ng,t]/.aux[[1]];
Print["The expresion of the element is ",sol];
bg=baseDeSGAfin[sistGen];
omega=Max[(Plus @@ # )& /@ (monomioAPunto[#,ng]&  /@ Algoritmo16[bg,puntoAMonomio[sol],ng])];
Return[omega],
Print["No solution found, it seems that " ,x ," is not an element of the semigroup"];
Return[];
]
];

OmegaPrimalityOfAffSG[sistGen_]:=Module[{bg,ng=Length[sistGen],aux},
bg=baseDeSGAfin[sistGen];
aux=Max[(Plus @@ # )& /@ (monomioAPunto[#,ng]&  /@ Algoritmo16[bg,Subscript[x, #],ng])] & /@ Table[i,{i,ng}];
Print["\[Omega]-primalities of the generators: ",aux];
Return[Max[aux]];
];


End[];
EndPackage[];
