function MinimalPresentations3S(n1,n2,n3,P) 
%==========================================================================
% - Función que da la dimensión de inmersión del semigrupo numérico generado
% por n1,n2,n3.
% - En caso de ser 3, calcula los generadores de la presentación minimal del
% semigrupo mediante los cálculos de Herzog.
% - Además da los generadores minimales del ideal asociado a la curva
% monomial de parámetros n1,n2,n3. 
% - Si P=1, entonces hace una gráfica en 3D donde se muestra la curva como
% intersección de las hipersuperfícies encontradas. 
%==========================================================================
tic
[c,r,CD]=Coeficientesci(n1,n2,n3);
if CD==1
    if r(2,1)==0
        Caso=2;
      % Generadores de la presentación minimal
        fprintf('Generadores de la presentación minimal:\n');
        fprintf('\t[(0,0,%d),(0,%d,0)],[(%d,0,0),(0,%d,%d)]\n',c(3),c(2),c(1),r(1,2),r(1,3));
      % Generadores del ideal asociado  
        p1=r(1,2)+floor(r(1,3)/c(3))*c(3);
        p2=r(1,3)-floor(r(1,3)/c(3))*c(3);
        fprintf('Generadores de I_S:\n');
        fprintf('\tI_S=(F_1=y^%dz^%d-x^%d,F_2=z^%d-y^%d)\n',p1,p2,c(1),c(3),c(2));
        if P==1
          % Dibujamos la curva monomial como intersección de las
          % hipersuperfícies obtenidas
            f = @(x,y,z) y.^p1.*z.^p2-x.^c(1);
            g = @(x,y,z) z.^c(3)-y.^c(2);
            h = 0;
            PlotMonomialCurves(Caso,n1,n2,n3,f,g,h)
        end
    elseif r(3,2)==0
        Caso=2;
      % Generadores de la presentación minimal
        fprintf('Generadores de la presentación minimal:\n');
        fprintf('\t[(%d,0,0),(0,0,%d)],[(0,%d,0),(%d,0,%d)]\n',c(1),c(3),r(2,1),c(2),r(2,3));
      % Generadores del ideal asociado  
        p1=r(2,1)+floor(r(2,1)/c(1))*c(1);
        p2=r(1,3)-floor(r(2,1)/c(1))*c(3);
        fprintf('Generadores de I_S:\n');
        fprintf('\tI_S=(F_1=x^%dz^%d-y^%d,F_2=x^%d-z^%d)\n',p1,p2,c(2),c(1),c(3));
      if P==1
          % Dibujamos la curva monomial como intersección de las
          % hipersuperfícies obtenidas
            f = @(x,y,z) x.^p1.*z.^p2-y.^c(2);
            g = @(x,y,z) x.^c(1)-z.^c(3);
            h = 0;
            PlotMonomialCurves(Caso,n1,n2,n3,f,g,h)
      end
    elseif r(1,3)==0
        Caso=2;
      % Generadores de la presentación minimal
        fprintf('Generadores de la presentación minimal:\n');
        fprintf('\t[(0,%d,0),(%d,0,0)],[(0,0,%d),(%d,%d,0)]\n',c(2),c(1),c(3),r(3,1),r(3,2));
      % Generadores del ideal asociado  
        p1=r(3,1)+floor(r(3,2)/c(2))*c(1);
        p2=r(3,2)-floor(r(3,2)/c(2))*c(2);
        fprintf('Generadores de I_S:\n');
        fprintf('\tI_S=(F_1=x^%dy^%d-z^%d,F_2=y^%d-x^%d)\n',p1,p2,c(3),c(2),c(1));
      if P==1
          % Dibujamos la curva monomial como intersección de las
          % hipersuperfícies obtenidas
            f = @(x,y,z) x.^p1.*y.^p2-z.^c(3);
            g = @(x,y,z) y.^c(2)-x.^c(1);
            h = 0;
            PlotMonomialCurves(Caso,n1,n2,n3,f,g,h)
      end
    else
      % Caso NO Intersección completa
        Caso=1;
        fprintf('Generadores de la presentación minimal:\n');
        fprintf('\t[(0,%d,%d),(%d,0,0)],[(%d,0,%d),(0,%d,0)],[(0,0,%d),(%d,%d,0)]\n',r(1,2),r(1,3),c(1),r(2,1),r(2,3),c(2),c(3),r(3,1),r(3,2));
      % Generadores del ideal asociado  
        fprintf('Generadores de I_S:\n');
        fprintf('\tI_S=(F_1=y^%dz^%d-x^%d,F_2=x^%dz^%d-y^%d,F_3=x^%dy^%d-z^%d)\n',r(1,2),r(1,3),c(1),r(2,1),r(2,3),c(2),r(3,1),r(3,2),c(3));
      if P==1
          % Dibujamos la curva monomial como intersección de las
          % hipersuperfícies obtenidas
            f = @(x,y,z) y.^r(1,2).*z.^r(1,3)-x.^c(1);
            g = @(x,y,z) x.^r(2,1).*z.^r(2,3)-y.^c(2);
            h = @(x,y,z) x.^r(3,1).*y.^r(3,2)-z.^c(3);
            PlotMonomialCurves(Caso,n1,n2,n3,f,g,h)
      end
    end
end
toc
end

function [c,r,CD]=Coeficientesci(n1,n2,n3)
% Declaración de variables
  c=zeros(3,3);
  mu=zeros(1,3);
  r=zeros(3,3);
% Ordenamos n1<n2<n3
  n=sort([n1 n2 n3]);
% Comprobamos mcd(n1,n2,n3)=1
  [d12,s,t]=gcd(n1,n2);
  if 1~=gcd(d12,n3)
     fprintf('¡ERROR: No se puede generar el semigrupo numérico!')
     CD=0;
     return;
  elseif n1==d12
      e=2;
      fprintf('¡Error: e(<%d,%d,%d>)=%d!\n',n(1),n(2),n(3),e);
      CD=0;
      return;
  end
% Comprobamos que e(S)=3, mirando si n3 está en <n1,n2>
  if d12>1
      e=3;
      fprintf('e(<%d,%d,%d>)=%d\n',n(1),n(2),n(3),e);
      CD=1;
  elseif n(3)>n(1)*n(2)-n(1)-n(2)
      e=2;
      fprintf('¡Error: e(<%d,%d,%d>)=%d!\n',n(1),n(2),n(3),e);
      CD=0;
      return;
  else
      F1=-n(3)*(s)/(n(2));
      F2=n(3)*(t)/(n(1));
      if F1~=floor(F1) && floor(F1)==floor(F2)
          e=3;
          fprintf('e(<%d,%d,%d>)=%d\n',n(1),n(2),n(3),e);
          CD=1;
      else
          e=2;
          fprintf('¡Error: e(<%d,%d,%d>)=%d!\n',n(1),n(2),n(3),e);
          CD=0;
          return;
      end
  end
  if e==3
    % Control de dimensión de inmersión=3  
      CD=1;
    % Si e(S)=3, calculamos los coeficientes ci y rij
      for i=1:3
          j=mod(i+1,3);
          if j==0
              j=3;
          end
          k=mod(i+2,3);
          if k==0
              k=3;
          end
          dij=gcd(n(i),n(j));
          dik=gcd(n(i),n(k));
          M=min(n(i)*n(j)/dij,n(i)*n(k)/dik);
          [djk,s,t]=gcd(n(j),n(k));
          for lambda=2:M
              if mod(lambda*n(i),djk)==0
                  F1=-lambda*n(i)*(s)/(n(k));
                  F2=lambda*n(i)*(t)/(n(j));
                    if F1==floor(F1)
                      c(i)=lambda;
                      mu(i)=floor(F1);
                      break;
                  elseif floor(F1)<floor(F2)
                      c(i)=lambda;
                      mu(i)=floor(F1)+1;
                      break;
                    end
              end
          end
          x0=c(i)*n(i)/djk*s;
          y0=c(i)*n(i)/djk*t;
          r(i,j)=x0+mu(i)*n(k)/djk;
          r(i,k)=y0-mu(i)*n(j)/djk;
      end
  end
end

function PlotMonomialCurves(Caso,n1,n2,n3,f,g,h)
% Declaración de variables
  t0=-1;
  tf=1;
  M=10000;
  x=zeros(1,M);
  y=zeros(1,M);
  z=zeros(1,M);
  t = t0:(tf-t0)/M:tf;
% Calculamos la curva monomial
  for i=1:M
      x(i)=t(i)^n1;
      y(i)=t(i)^n2;
      z(i)=t(i)^n3;
  end
  if Caso==2
      figure(1)
      hold on
      fimplicit3(f,[-1 1],'red','FaceAlpha',0.7,'EdgeColor','white','MeshDensity',75)
      fimplicit3(g,[-1 1],'FaceAlpha',0.7,'EdgeColor','none','MeshDensity',75)
      plot3(x,y,z,'-o','Color','black','MarkerSize',10,'LineWidth',10.5)
      legend('\bf F_1','\bf F_2','\bf Curva monomial','Position',[0.2 0.6 0.2 0.3])
      hold off
  end
  if Caso==1
      figure(1)
      hold on
      fimplicit3(f,[-1 1],'FaceAlpha',0.7,'EdgeColor','cyan','MeshDensity',50)
      fimplicit3(g,[-1 1],'red','FaceAlpha',0.7,'EdgeColor','white','MeshDensity',50)
      fimplicit3(h,[-1 1],'MeshDensity',50)
      plot3(x,y,z,'-o','Color','black','LineWidth',10.5)
      legend('\bf F_1','\bf F_2','\bf F_3','\bf Curva monomial','Position',[0.2 0.6 0.2 0.3])
      hold off
  end
end