(* ::Package:: *)

BeginPackage["CircleSG`"]
(* Authors: 
I. Garc\'{i}a-Garc\'{i}a <ignacio.garcia@uca.es>, 
M.A. Moreno-Fr\'{i}as <mariangeles.moreno@uca.es>, 
A. S\'{a}nchez-R.-Navarro <alfredo.sanchez@uca.es>, 
A. Vigneron-Tenorio <alberto.vigneron@uca.es>

We would like to thank P. A. Garc\'{i}a-S\'{a}nchez for his comments and suggestions.
*)
CircleSG::usage="CircleSG[a,b,r,view], if exists, computes the minimal system of generators of the semigroup associated to the circle with center (a, b) and radius r (a, b and r must belong to the rational number set). If view is True, a plot will be rendered. This program does not yet computes the minimal system of generators if the circle intersects both axis in one or more points. If you detect any error, please notify us."
Begin["Private`"]
CurveIntersections[a_, b_, r_] := Module[{},
If[a!= 0&&b!= 0,
Return[{{(a^3+a b^2-a r^2-Sqrt[a^2 b^2 r^2+b^4 r^2-b^2 r^4])/(a^2+b^2), -(1/b)(-a^2-b^2+a^4/(a^2+b^2)+(a^2 b^2)/(a^2+b^2)+r^2-(a^2 r^2)/(a^2+b^2)-(a Sqrt[a^2 b^2 r^2+b^4 r^2-b^2 r^4])/(a^2+b^2))},{(a^3+a b^2-a r^2+Sqrt[a^2 b^2 r^2+b^4 r^2-b^2 r^4])/(a^2+b^2),-(1/b)(-a^2-b^2+a^4/(a^2+b^2)+(a^2 b^2)/(a^2+b^2)+r^2-(a^2 r^2)/(a^2+b^2)+(a Sqrt[a^2 b^2 r^2+b^4 r^2-b^2 r^4])/(a^2+b^2))}}];,
If[a!= 0,
Return[{{(a^2-r^2)/a,-((r Sqrt[a^2-r^2])/a)},{(a^2-r^2)/a,(r Sqrt[a^2-r^2])/a}}];,
Return[{{-((r Sqrt[b^2-r^2])/b),(b^2-r^2)/b},{(r Sqrt[b^2-r^2])/b,(b^2-r^2)/b}}];
];
];
];
FactorMul[s_, g_]:=Module[
{resto, cocie, restoy, cociey},
If[g[[1]]!= 0,
(*then*) 
resto=Mod[s[[1]],g[[1]]];
cocie=s[[1]]/g[[1]];,
(*else*)
resto=-1;
];
If[g[[2]]!= 0,
(*then*)
restoy=Mod[s[[2]],g[[2]]];
cociey=s[[2]]/g[[2]];
If[cocie!= cociey, 
(*then*)
resto=restoy=-1;
];,
restoy=-1;
];
If[resto<0, resto = restoy; cocie=cociey];
Return[{resto, cocie}];
];
GetB1B2[s_, g_] := Module[
{lmbdt,done,kt,p,resto,cocie,i,k,t, gp,gg,B1,B11},
lmbdt=0;
done=False;
kt={};
For[p=1, !done && p<=Length[g], p++,
{resto, cocie}=FactorMul[s[[1]], g[[p]]];
If[resto== 0,
(*then*) AppendTo[kt, cocie];
i=2;
For[t=2, t<= Length[s],t++,
{resto, cocie}=FactorMul[s[[t]], g[[p]]];
If[resto== 0,(*then*) AppendTo[kt, cocie]; i=i+1;];
];
lmbdt=Max[kt];
gp=g[[p]];
gg=Delete[g, {p}]; done = True;
];
];
B1 = s;
B1=Union[B1, gg];
B11={};
Do[
Do[
If[!MemberQ[ kt, k],(*then*) AppendTo[B11, gg[[i]]+k*gp];];,
{k, 1, lmbdt-1}
];,
{i,1,Length[gg]}
];
Return[{B1,B11}];
];
GetB1B2var[s_, g_] := Module[
{lmbdt,done,kt,p,resto,cocie,i,k,t, gp,gg,B1,B11},
lmbdt={};
done=False;
kt={};
lmbdt=Max[s[[All,1]]];
If[MemberQ[g,{1,0}],
(*then*)
gg=Complement[g,{{1,0}}];
gp={1,0};,
(*else*)
Print["GetB1B2var>> Error, no esta el {1,0}."];
Abort[];
];
B1 = s;
B1=Union[B1, gg];
B11={};
Do[
Do[
If[!MemberQ[ kt, k],(*then*) AppendTo[B11, gg[[i]]+k*gp];];,
{k, 1, lmbdt-1}
];,
{i,1,Length[gg]}
];
Return[{B1,B11}];
];
EveMinimal[xc_, lc_]:=Module[{i,il},
For[il=1,il<=Length[lc],il++,
For[i=1,i<=Length[lc[[il]]],i++,
If[lc[[il,i]]>xc[[i]], Break[];];
];
If[i== Length[lc[[il]]]+1,
Return[False];
];
];
Return[True];
];
Evelyne[A_, b_]:=Module[{m,n,w,nohom,mZero,mZerom,mbA, mBlocked,vBlocked,x,myStack,kstp,Atx,i,esc,S},
m=Length[A];
n=w=Length[A[[1]]];
mZero=ConstantArray[0,{w,1}];
mZerom=ConstantArray[0,{m,1}];
If[b== mZerom,
nohom=False;
mbA = A;,
w=w+1; 
nohom=True;
mbA=ArrayFlatten[{{A,-b}}];
];
mZerom = Flatten[mZerom];
myStack={{0},ConstantArray[0,{w,w}]};
mBlocked = ConstantArray[True,{w,w}];
mBlocked[[1]]=ConstantArray[False,{1,w}][[1]];
x=Table[0,{i,w}];
x[[1]] = 1;
myStack[[1]]=1;
myStack[[2]][[myStack[[1]]]] = x;
For[i=2,i<= w,i++,
x[[i-1]]--; x[[i]]++;
myStack[[2]][[++myStack[[1]]]] = x;
mBlocked[[i]][[i;;w]]=Table[False,{j,w-i+1}];
];
If[nohom,mBlocked[[w,w]] = True];
kstp=0;
S={};
While[myStack[[1]]>0,
kstp++;
x=myStack[[2]][[myStack[[1]]--]];
Atx=mbA.x;
If[Atx== mZerom,
(* then *)
AppendTo[S,x[[1;;n]]];
,
(* else *)
vBlocked=mBlocked[[myStack[[1]]+1]];
For[i=1,i<= w,i++,
If[!vBlocked[[i]],
esc =Atx.mbA[[All,i]];
If[esc<0,
x[[i]]++;
If[EveMinimal[x[[1;;n]],S],
myStack[[2]][[++myStack[[1]]]]=x;
mBlocked[[myStack[[1]]]]=vBlocked;
If[nohom  && i==w, 
mBlocked[[myStack[[1]],w]]=True;
];
];
vBlocked[[i]]=True;
x[[i]]--;
];
];
];
];
];
Return[S];
];
CalcCota[A_,b_]:=Module[{lA=A, lb=b,vCota,m,n,i,j,vauxmax},
m=Length[lA];
If[m!= Length[lb],
Print["CalcCota>> El tama\[CapitalATilde]\[PlusMinus]o de las matrices es distinto."]; Abort[];
];
n=Length[lA[[1]]];
vCota={};
i=1;
vauxmax= Max[lb];
While[i<=n ,
vCota=Append[vCota, vauxmax];
For[j=1,j<m,j++,
If[lA[[j,i]]!= 0,
vCota[[i]]=Min[vCota[[i]],Floor[lb[[j,1]]/lA[[j,i]]]];
];
];
i++;
];
Return[vCota];
];
Evelyne1sc[A_, b_]:=Module[{m,n,w,nohom,mZero,mZerom,mbA,vCota,myStack,mBlocked,vBlocked,x,i,j,kstp,Atx,esc},
m=Length[A];
n=w=Length[A[[1]]];
mZero=ConstantArray[0,{w,1}];
mZerom=Flatten[ConstantArray[0,{m,1}]];
If[b== mZero,
Print["Evelyne1sc>> Esto no funciona con t.i. a 0."];
Abort[];
];
nohom=True;
w=w+1;
mbA=ArrayFlatten[{{A,-b}}];
vCota=CalcCota[A,b];
vCota=Append[vCota,1];
myStack={{0},ConstantArray[0,{w,w}]};
mBlocked = ConstantArray[True,{w,w}];
mBlocked[[1]]=ConstantArray[False,{1,w}][[1]];
x=Table[0,{i,w}];
x[[1]] = 1;
myStack[[1]]=1;
myStack[[2]][[myStack[[1]]]] = x;
For[i=2,i<= w,i++,
x[[i-1]]--; x[[i]]++;
myStack[[2]][[++myStack[[1]]]] = x;
mBlocked[[i]][[i;;w]]=Table[False,{j,w-i+1}];
];
If[nohom,mBlocked[[w,w]] = True];
kstp=0;
While[myStack[[1]]>0,
kstp++;
x=myStack[[2]][[myStack[[1]]--]];
Atx=mbA.x;
If[Atx== mZerom,
Return[x[[1;;n]]];,
vBlocked=mBlocked[[myStack[[1]]+1]];
For[i=1,i<= w,i++,
If[!vBlocked[[i]],
esc =Atx.mbA[[All,i]];
If[esc<0 && x[[i]]<=vCota[[i]],
x[[i]]++;
myStack[[2]][[++myStack[[1]]]]=x;
mBlocked[[myStack[[1]]]]=vBlocked;
If[nohom  && i==w, 
mBlocked[[myStack[[1]],w]]=True;
];
vBlocked[[i]]=True;
x[[i]]--;
];
];
];
];
];
];
GetMinGenSysCon[Lops_]:=Module[
{lLops=Lops , nuelin, result, ii, TI, facLops, Ac, bti, ToS},
nuelin = Length[lLops];
result={};
ii=1;
While [ii<=Length[lLops],
TI=lLops[[ii]];
lLops=Delete[lLops,ii];
facLops=Select[lLops,#[[1]]<=TI[[1]] && #[[2]]<=TI[[2]]&];
If[Length[facLops]>0,
Ac=Transpose[facLops];
bti=Transpose[{TI}];
ToS=Evelyne1sc[Ac,bti];
If[ToS==Null,
(*then*)result=Append[result,TI];lLops=Insert[lLops, TI, ii];,
(*else*) ii=ii-1; 
];,
(*else*)
result=Append[result,TI];
lLops=Insert[lLops, TI, ii];
];
ii=ii+1;
];
Return[result];
];
GetMinGenSysF1[Lops_, ptau1_,ptau2_, s1_, s2_]:=Module[
{lLops=Lops , nuelin, result, ii, piv, dif, des, facLops, pend, done},
nuelin = Length[lLops];
ii=1;
While [ii<=nuelin,
des=lLops[[ii]];
facLops=Select[lLops,#[[1]]<=  des[[1]] && #[[2]]<=  des[[2]] && (#[[1]]!=   des[[1]] || #[[2]]!=   des[[2]])&];
done=False;
piv=1;
While[!done && piv<=Length[facLops],
dif =des-facLops[[piv]];
If[dif[[1]]!= 0,
(* then *)
pend=Simplify[dif[[2]]/dif[[1]]];
If[pend<  ptau1 && pend> ptau2 ,(* then *)done=True;];
If[(pend== ptau1 && (dif[[1]]*s1[[2]]- dif[[2]]*s1[[1]]== 0 && Element[dif[[1]]/s1[[1]],Integers] )) || (pend== ptau2 && (dif[[1]]*s2[[2]]- dif[[2]]*s2[[1]]== 0 && Element[dif[[1]]/s2[[1]],Integers] )), done=True;];
];
piv++;
];
If[!done,(* then *)ii++;,(* else *)lLops=Delete[lLops,ii];nuelin--;];
];
Return[lLops];
];
GetMinGenSysF2[Lops_, ptau1_, s1_,a_, b_, c_]:=Module[
{lLops=Lops , nuelin, result, ii, piv, dif,des, facLops, pend, done},
nuelin = Length[lLops];
ii=1;
While [ii<=nuelin,
des=lLops[[ii]];
facLops=Select[lLops,#[[1]]<=  des[[1]] && #[[2]]<=  des[[2]] && (#[[1]]!=   des[[1]] || #[[2]]!=   des[[2]])&];
done=False;
piv=1;
While[!done && piv<=Length[facLops],
dif=des- facLops[[piv]];
If[dif[[1]]!= 0,
(* then *)
pend=Simplify[(dif[[2]]/dif[[1]])];
If[(pend== 0 && (Mod[a *dif[[1]],b]<= c*dif[[1]]==1)) || (0<pend<  ptau1)  || (pend== ptau1 && (dif[[1]]*s1[[2]]- dif[[2]]*s1[[1]]== 0 && Element[dif[[1]]/s1[[1]],Integers] )),
(*then*) done=True;
];
];
piv++;
];
If[!done,(* then *)ii++;,(* else *)lLops=Delete[lLops,ii];nuelin--;];
];
Return[lLops];
];
GetMinGenSysFi1[Lops_, jd_, ptau1_, ptau2_,s1_,s2_]:=Module[
{lLops=Lops , nuelin, result, ii, piv, dif, des, facLops, pend, done},
nuelin = Length[lLops];
ii=1;
While [ii<=nuelin,
des=lLops[[ii]];
facLops=Select[lLops,#[[1]]<=  des[[1]] && #[[2]]<=  des[[2]] && (#[[1]]!=   des[[1]] || #[[2]]!=   des[[2]])&];
done=False;
piv=1;
While[!done && piv<=Length[facLops],
dif =des-facLops[[piv]];
If[!MemberQ[jd, dif],
If[dif[[1]]!= 0,
(* then *)
pend=Simplify[dif[[2]]/dif[[1]]];
If[pend< ptau1 && pend> ptau2 ,(* then *)done=True;];
If[(pend== ptau1 && (dif[[1]]*s1[[2]]- dif[[2]]*s1[[1]]== 0 && Element[dif[[1]]/s1[[1]],Integers] )) || (pend== ptau2 && (dif[[1]]*s2[[2]]- dif[[2]]*s2[[1]]== 0 && Element[dif[[1]]/s2[[1]],Integers] )), done=True;];
];
];
piv++;
];
If[!done,(* then *)ii++;,(* else *)lLops=Delete[lLops,ii];nuelin--;];
];
Return[lLops];
];
GetMinGenSysFi2[Lops_, jd_, ptau1_, s1_,a_, b_, c_]:=Module[
{lLops=Lops , nuelin, result, ii, piv, dif, des, facLops, pend, done},
nuelin = Length[lLops];
ii=1;
While [ii<=nuelin,
des=lLops[[ii]];
facLops=Select[lLops,#[[1]]<=  des[[1]] && #[[2]]<=  des[[2]] && (#[[1]]!=   des[[1]] || #[[2]]!=   des[[2]])&];
done=False;
piv=1;
While[!done && piv<=Length[facLops],
dif =des-facLops[[piv]];
If[!MemberQ[jd, dif],
If[dif[[1]]!= 0,
(* then *)
pend=Simplify[dif[[2]]/dif[[1]]];
If[pend< ptau1 && pend> 0,(* then *)done=True;];
If[(dif[[2]]== 0 &&  (Mod[a *dif[[1]],b]<= c*dif[[1]])== 1), done=True;];
];
];
piv++;
];
If[!done,(* then *)ii++;,(* else *)lLops=Delete[lLops,ii];nuelin--;];
];
Return[lLops];
];
GetMinGenSysNDC[Lops_,nc_, a_,b_,r_]:=Module[
{lLops=Lops,rc,ambc,rc1,dx,dx1,dy,dy1,dif,nuelin,ii,des,facLops, done, piv,k},
nuelin = Length[lLops];
ii=1;
ambc=a*a+b*b;
While [ii<=nuelin,
des=lLops[[ii]];
facLops=Select[lLops,#[[1]]<=  des[[1]] && #[[2]]<=  des[[2]] && (#[[1]]!=   des[[1]] || #[[2]]!=   des[[2]])&];
done=False;
piv=1;
While[!done && piv<=Length[facLops],
dif=des-facLops[[piv]];
k=Floor[Sqrt[(dif[[1]]*dif[[1]]+dif[[2]]*dif[[2]])/ambc]];
dx=dif[[1]]-k*a; dy=dif[[2]]-k*b;
dx1=dif[[1]]-(k+1)*a; dy1=dif[[2]]-(k+1)*b;
rc=(k*r)*(k*r); rc1=((k+1)*r)*((k+1)*r);
If[(dx*dx+dy*dy)<= rc || (dx1*dx1+dy1*dy1)<= rc1,
(*then*) done=True;,
piv++;
];
];
If[!done,(* then *)ii++;,(* else *)lLops=Delete[lLops,ii];nuelin--;];
];
Return[lLops];
];
Int2InCone[m1_,m2_,p1_,p2_]:=Module[
{malla, m, limx,limy,x,y,c},
malla={};
If[m2!= 0,
limy=p2[[2]];
For[y=1,y<=limy, y++,
limx=Floor[y/m2];
x=Ceiling[y/m1];
If[x*m1== y,x++];
For[x=x, x<limx, x++,
malla=Append[malla, {x,y}];
];
If[x*m2<y,malla=Append[malla,{x,y}];]; 
];
];
m=p2[[2]]-p1[[2]];
If[m!= 0,
m=(p2[[1]]-p1[[1]])/m;
];
limy=p1[[2]];
c=-m*p1[[2]]+p1[[1]];
For[y=Max[1,Ceiling[p2[[2]]]],y<=limy, y++,
limx=Floor[m*y+c];
x=Ceiling[y/m1];
If[x*m1== y,x++];
For[x=x, x<limx, x++,
malla=Append[malla, {x,y}];
];
If[(m*y-x+c)*Sign[c]> 0,malla=Append[malla,{x,y}];];
];
Return[malla];
];
Int2InConeOnAbcis[m1_,p1_,p2_]:=Module[
{malla, m, limx,limy,x,y,c},
malla={};
m=p2[[2]]-p1[[2]];
If[m!= 0,
(*then*)m=(p2[[1]]-p1[[1]])/m;
];
limy=p1[[2]];
c=-m*p1[[2]]+p1[[1]];
For[y=1,y<=limy, y++,
limx=Floor[m*y+c];
x=Ceiling[y/m1];
If[x*m1== y,x++];
For[x=x, x<limx, x++,
malla=Append[malla, {x,y}];
];
If[(m*y-x+c)*Sign[c]> 0,malla=Append[malla,{x,y}];];
];
Return[malla];
];
Int2InCircle[x0_,y0_,r0_]:=Module[
{xa=Floor[x0],ya=Floor[y0],ra=Floor[r0]+2,ta},
Return[Select[Flatten[Table[{xx,yy},{xx,xa-ra,xa+ra},{yy,ya-ra,ya+ra}],1],(#[[1]]-x0)^2+(#[[2]]-y0)^2<=r0^2&]];
];
NearRational[x_]:=Module[{res},
If[Element[x,Rationals],
Return[x];,
res=Rationalize[N[x],10^(Round[Log[10,FractionalPart[N[x]]]]-1)];
If[res==0, Print["NearRational>> El racional cercano no se calculo como era esperable."];,
Return[res];
];
];
];
MinGenPGS[l_,a_,b_,c_]:= (* Generadores minimales de semigrupo proporcionalmente generado, de Pedro A. Garc\[CapitalATilde]\[DiscretionaryHyphen]a S\[CapitalATilde]\[DownExclamation]nchez *)Select[l, Function[y,Length[Select[Select[l,(#<y)&],Function[x,Mod[a (y-x),b]<= c (y-x)]]]==0]];
TriOri[p1_,p2_,p3_]:=
Sign[(p1[[1]]-p3[[1]])*(p2[[2]]-p3[[2]])-(p1[[2]]-p3[[2]])*(p2[[1]]-p3[[1]])];
SqTri[p1_,p2_,p3_,c1_,r1_,c2_,r2_,m1_,m2_]:=Module[
{cdl,cur,sqpts,fndis},
cdl={Floor[Min[p1[[1]],p2[[1]],p3[[1]]]], Floor[Min[p1[[2]],p2[[2]],p3[[2]]]]};
cur={Ceiling[Max[p1[[1]],p2[[1]],p3[[1]]]], Ceiling[Max[p1[[2]],p2[[2]],p3[[2]]]]};
sqpts=Flatten[Table[{i,j},{i,cdl[[1]], cur[[1]],1}, {j,cdl[[2]],cur[[2]],1}],1];
sqpts=Select[sqpts,#[[2]]>= 0&];
fndis=Function[t,If[TriOri[p1,p2,p3]>= 0,If[TriOri[p1,p2,t]>=0 && TriOri[p2,p3,t]>= 0&&TriOri[p3,p1,t]>= 0,True,False],If[TriOri[p1,p2,t]<0 && TriOri[p2,p3,t]< 0&&TriOri[p3,p1,t]< 0,True,False]]&& EuclideanDistance[c1,t]> r1 && EuclideanDistance[c2,t]> r2  &&t[[2]]/t[[1]]!= m1 && t[[2]]/t[[1]]!= m2];
Return[Select[sqpts,fndis]];
];

CircleSG[a0_, b0_, r0_, view_:False] := Module[
{a,b,r,symbis,H,m1,m2,negslope, alfa, beta, t,t0,alfa0,beta0,dab,dnr,alfa1,beta1, Ac, bti,CpG, CpG1, cp, ap, bp, CpG1min,CpG2, C4,C2,g1,g2,F,FDT,s1,s2,
B1, B1p,B1B2, tA1,tA2,tB1,tB2,d1,d2,d,x,y,i0,x0,y0,x1,y1,i1,i2,i3,sig,lcm, plasta,plastap, incir},
a=Rationalize[a0];b=Rationalize[b0];r=Rationalize[r0];
If[r<=0, Print["Radius must be bigger than 0"]; Return[];];
If[a^2+b^2 < r^2, Return[{{1,0},{0,1}}];];
If[a+r<=0 || b+r<=0, Return[{{0, 0}}];];
If[((a - r < 0 && b - r < 0) || (a - r < 0 && 
       b - r == 0) || (a - r == 0 && b - r == 0) || (a - r == 0 && 
       b - r < 0) ) && (b + r > 0 && a + r > 0), 
Print["This case is not implemented yet, see ?CircleSG."];
Return[];
];
If[a - r <= 0  && b - r > 0,
symbis=True;
x=a;a=b;b=x; Clear[x];,
symbis=False;
];
If[!Element[a,Rationals]||!Element[b,Rationals]||!Element[r,Rationals],
(*then*)  MessageDialog ["Los par\[CapitalATilde]\[DownExclamation]metros no son racionales."];Return[];];
H=CurveIntersections[a,b,r];
If[H[[1,1]]!=0,
(*then*)
m1=Simplify[H[[1,2]]/H[[1,1]]];
If[! Element[m1, Rationals],
(*then*)  Print["The semigroup associated to the circle with center (", a,", ",b, ") and radius (", r, ") is not finitely generated"]; Return[];];
];
If[H[[2,1]]!=0,
(*then*) 
m2=Simplify[H[[2,2]]/H[[2,1]]];
If[! Element[m2, Rationals],
(*then*)  Print["The semigroup associated to the circle with center (", a,", ",b, ") and radius (", r, ") is not finitely generated"]; Return[];];
];
negslope = False;
If[m2<0, negslope = True;];
If[negslope, 
(*then*)
alfa=a-Sqrt[-b^2+r^2];beta=a+Sqrt[-b^2+r^2];
t=Minimize[{t0,t0*beta>=(t0+1)*alfa&& t0>=0&& Element[t0, Integers]},{t0}];
alfa0=Max[Table[Floor[i*alfa]/i,{i,t[[1]]+1}]];
beta0=Min[Table[(Floor[i*beta]+1)/i,{i,t[[1]]}]];
dab=Min[Abs[alfa-alfa0],Abs[beta-beta0]];
dnr=NearRational[dab/2];
alfa1=alfa0+dnr; beta1=beta0-dnr;
If[alfa1>alfa || beta1<beta, Print["CircleSG>> Ooops, esto no deberia haber ocurrido..."]; Return[];];
Ac={{Denominator[alfa1], -Numerator[alfa1], -1,  0},{Denominator[beta1], -Numerator[beta1], 0,  1}};
bti={{0},{0}};
CpG=Evelyne[Ac,bti];
If[CpG== {},
Print["CircleSG>> El cono para obtener los de y=0 no tiene generadores...."];
Return[];
];
CpG1= Sort[Union[CpG[[All,1]]]];
cp=LCM[Denominator[beta1/(beta1-alfa1)],Denominator[alfa1*beta1/(beta1-alfa1)]];
ap=cp*beta1/(beta1-alfa1); bp=cp*alfa1*beta1/(beta1-alfa1);
CpG1min = MinGenPGS[CpG1, ap,bp,cp];
CpG2=Append[{#},0]&/@CpG1min;
Ac={{Numerator[m1], -Denominator[m1], -1,  0},{0, 1, 0, -1 }};,
(* else *)
Ac={{Numerator[m1], -Denominator[m1], -1,  0},{Numerator[m2], -Denominator[m2], 0,  1}};
];
bti={{0},{0}};
C4=Evelyne[Ac,bti];
C2=C4[[All,1;;2]];
If[negslope,
(*then*)
C2=Complement[C2, {g1=H[[1]]/GCD[H[[1,1]],H[[1,2]]], g2={1,0}}];,
(*else*)
C2=Complement[C2, {g1=H[[1]]/GCD[H[[1,1]],H[[1,2]]], g2=H[[2]]/GCD[H[[2,1]],H[[2,2]]]}];
];
C2=Prepend[C2, g2]; C2=Prepend[C2, g1];


F=GetMinGenSysCon[C2];
FDT=GetMinGenSysF1[C2, m1, m2];

If[negslope,
(*then*)
s1 = Apply[LCM, Denominator[ H[[1]]]]*H[[1]];
B1p=GetB1B2[{s1}, F];
B1p=GetMinGenSysF2[Union[B1p[[1]],B1p[[2]]], m1, s1, ap, bp, cp];
B1p=GetB1B2var[CpG2,B1p];
B1p=GetMinGenSysF2[Union[B1p[[1]],B1p[[2]]],m1, s1, ap, bp, cp];,
(* else *)
s1 = Apply[LCM, Denominator[ H[[1]]]]*H[[1]]; s2 =Apply[LCM, Denominator[ H[[2]]]]*H[[2]];
B1p=GetB1B2[{s1}, F];
B1p=GetMinGenSysF1[Union[B1p[[1]],B1p[[2]]],m1, m2, s1, s2];
B1p=GetB1B2[{s2},B1p];
B1p=GetMinGenSysF1[Union[B1p[[1]],B1p[[2]]],m1, m2, s1, s2];
];
tA1=Numerator[m1]; tB1=Denominator[m1];
d1=Minimize[{(tA1*x-tB1*y)/Sqrt[tA1^2+tB1^2],tA1*x-tB1*y>0 && tB1*x+tA1*y-(tB1*s1[[1]]+tA1*s1[[2]])>=0 && tB1*x+tA1*y-(tB1*2*s1[[1]]+tA1*2*s1[[2]])<=0 && tA1*x-tB1*y-Sqrt[tA1^2+tB1^2]*Sqrt[2] <=0 && x>0 && y>0 && Element[x, Integers] && Element[y, Integers]},{x,y}];
d=d1[[1]];

tA2=Numerator[m2]; tB2=Denominator[m2];
If[!negslope,
(*then*)
d2=Minimize[{-(tA2*x-tB2*y)/Sqrt[tA2^2+tB2^2],tA2*x-tB2*y<0 && tB2*x+tA2*y-(tB2*s2[[1]]+tA2*s2[[2]])>=0 && tB2*x+tA2*y-(tB2*2*s2[[1]]+tA2*2*s2[[2]])<=0 && tA2*x-tB2*y+Sqrt[tA2^2+tB2^2]*Sqrt[2] >=0 && x>0 && y>0 && Element[x, Integers] && Element[y, Integers]},{x,y}];
d=Min[d,d2[[1]]];
];
x0=(a^3+a b^2+2 a^3 i0+2 a b^2 i0-a r^2-2 a i0 r^2+Sqrt[-b^2 (a^2+b^2-r^2) (a^2+b^2-(r+2 i0 r)^2)])/(2 (a^2+b^2));
y0=(a^2 b^2 (1+2 i0)+b^2 (1+2 i0) (b^2-r^2)-a Sqrt[-b^2 (a^2+b^2-r^2) (a^2+b^2-(r+2 i0 r)^2)])/(2 b (a^2+b^2));
x1=(a^3+a b^2+2 a^3 i0+2 a b^2 i0-a r^2-2 a i0 r^2-Sqrt[-b^2 (a^2+b^2-r^2) (a^2+b^2-(r+2 i0 r)^2)])/(2 (a^2+b^2));
y1=(a^2 b^2 (1+2 i0)+b^2 (1+2 i0) (b^2-r^2)+a Sqrt[-b^2 (a^2+b^2-r^2) (a^2+b^2-(r+2 i0 r)^2)])/(2 b (a^2+b^2));
i1=Minimize[{i0,(-(tA2*x0-tB2*y0)/Sqrt[tA2^2+tB2^2])<= d && x0>= 0 &&  i0>= 0 && Element[i0, Integers]},{i0}];
i3=Minimize[{i0,x0>= 0 &&  i0>= 0 && Element[i0, Integers]},{i0}];
i3=i3[[1]];
plasta={};
Do[
i0=i;
AppendTo[plasta,SqTri[i*H[[1]],(i+1)*H[[1]],{x1,y1},i*{a,b},i*r,(i+1)*{a,b},(i+1)*r ,m1,m2]];
AppendTo[plasta, SqTri[i*H[[2]],(i+1)*H[[2]],{x0,y0},i*{a,b},i*r,(i+1)*{a,b},(i+1)*r , m1,m2]];,
{i,i3,i1[[1]]-1}
];
plasta=Flatten[plasta,1];

If[negslope,
(*then*)
plastap=Int2InConeOnAbcis[m1,i3*H[[1]], {i3*alfa1,0}];,
(*else*)
plastap=Int2InCone[m1, m2, i3*H[[1]], i3*H[[2]]];
];
For[i=1,i<=i3,i++,
incir=Int2InCircle[a*i,b*i,r*i];
plastap=Complement[plastap, incir];
];
plasta=Union[plasta,plastap];

B2p=B1p;
CoEl = plasta;
CoEld={};
i=1;
While[Length[CoEl]> 0 && i<= Length[CoEl],
If[MemberQ[B2p, CoEl[[i]]],
(*then*)
B2p=Complement[B2p,{CoEl[[i]]}];
CoEld=Union[CoEld,{CoEl[[i]]}];
B2pB=Union[(CoEl[[i]]+#)&/@B2p, {2*CoEl[[i]]}, {3*CoEl[[i]]}];
If[negslope,
B2p=GetMinGenSysFi2[Union[B2p, B2pB], CoEld, m1, s1, ap, bp, cp];,
B2p=GetMinGenSysFi1[Union[B2p, B2pB], CoEld, m1, m2, s1, s2];
];
CoEl=Delete[CoEl,i];
i=1;,
(*else*) 
i++
];
];

If[symbis,
x=a;a=b;b=x;
plasta={#[[2]],#[[1]]}&/@plasta;
B2p={#[[2]],#[[1]]}&/@B2p;
H={#[[2]],#[[1]]}&/@H;
H={H[[2]],H[[1]]};
If[H[[1,1]]!=0,
m1=Simplify[H[[1,2]]/H[[1,1]]];,
m1=0;
];
If[H[[2,1]]!=0,
m2=Simplify[H[[2,2]]/H[[2,1]]];,
m2=0;
];
];

If[view,
If[plasta=={},
Print[Show[
Graphics[Line[{{0,0}, i1[[1]]*H[[1]]}]], Graphics[Line[{{0,0}, i1[[1]]*H[[2]]}]], Graphics[Table[Circle[{a*i,  b*i}, r*i], {i,1,i1[[1]]+2}]],
ListPlot[B2p,PlotStyle->Directive[Red,Thick,PointSize[0.005]]],
Axes-> True,AxesOrigin -> {0,0},PlotRange->{{0,All},{0,All}},PlotRangeClipping->True
]],
Print[Show[
Graphics[Line[{{0,0}, i1[[1]]*H[[1]]}]], Graphics[Line[{{0,0}, i1[[1]]*H[[2]]}]], 
ListPlot[plasta,PlotStyle->Directive[Blue,Thick,PointSize[0.005]]],
Graphics[Table[Circle[{a*i,  b*i}, r*i], {i,1,i1[[1]]+2}]],
ListPlot[B2p,PlotStyle->Directive[Red,Thick,PointSize[0.005]]],
Axes-> True,AxesOrigin -> {0,0},PlotRange->{{0,All},{0,All}},PlotRangeClipping->True
]]
]
];

Return[B2p];
];
End[]
EndPackage[]



