//use"realroot";
include"mubase.mmx";
include"pencilregular.mmx";
include"curvesurface.mmx";
include "EEA.mmx";
include "solvepolynomial.mmx";
//R:=QQ['t,'x,'y,'z];
//f1:=R<<"t^3+t+1";
//g1:=R<<"t^3-t+1";
//h1:=R<<"t^4-2*t+1";
//f2:=R<<"t^4+2*t+1";
//g2:=R<<"t^6-2*t+2";
//h2:=R<<"t^4-t+1";
//f3:=R<<"t^4+t^3+t+1";
//g3:=R<<"t^5-2*t+2";
//h3:=R<<"t^3-t+1";
//f4:=R<< "t^6+t^3+t^2";
//g4:=R<<"t^6-t^4-t^2";
//h4:=R<<"t^6+t^5+t^4-t-1";
//var0:=R<<"t";
//var1:=R<<"x";
//var2:=R<<"y";
//var3:=R<<"z";
//f5:=R<<"t^4-40*t^3+40*t+1";
//g5:=R<<"t^4+480*t^2+1";
//h5:=R<<"t^4+40*t^3+480*t^2+40*t+1";
//cur1:=R<<"3*t^8-4*t^7+5*t^6-4*t^5+4*t^4-4*t^3-3*t^2+5*t-1";
//cur2:=R<<"-1*t^8-4*t^7+4*t^6+3*t^5-5*t^4+9*t^2-10*t+3";
//cur3:=R<<"2*t^8-t^5-t^2+1";
//f:=[f1,f2,f3];
//f6:=R<<"t^2+t+1";
//g6:=R<<"t^3-2*t+3";
//h6:=R<<"1";
//f9:=R<<"t^5-2*t^3+2";
//g9:=R<<"t^4+3";
//h9:=R<<"1";
//Example 2.1 of paper Laurent and Carlos
//f:=4 * t^10 + 28 * t^9 + 85 * t^8 + 146 * t^7 + 155 * t^6 + 104 * t^5 + 43 * t^4 + 10 * t^3 + t^2=t^2*(2*t+1)^2*(t+1)^6;
//g:=12 * t^7 + 20 * t^6 + 15 * t^5 + 6 * t^4 + t^3=t^3*(2*t+1)^5*(3*t^2+2*t+1);
//h:=-t^10 - 10 * t^9 - 45 * t^8 - 120 * t^7 - 210 * t^6 - 252 * t^5 - 210 * t^4 - 120 * t^3 - 45 * t^2 - 10 * t - 1=-(t+1)^10;
// Singulartity of multiplicity 6: [0,-2,0] cors t=-1.
// Singularity of multiplicity 4: [0,0,-1],[0,0,-9.76*10^{-4}]
//f7:=R<<"4 * t^10 + 28 * t^9 + 85 * t^8 + 146 * t^7 + 155 * t^6 + 104 * t^5 + 43 * t^4 + 10 * t^3 + t^2";
//g7:=R<<"12 * t^7 + 20 * t^6 + 15 * t^5 + 6 * t^4 + t^3";
//h7:=R<<"-t^10 - 10 * t^9 - 45 * t^8 - 120 * t^7 - 210 * t^6 - 252 * t^5 - 210 * t^4 - 120 * t^3 - 45 * t^2 - 10 * t - 1";
// Example of Chen and all
//f8:=R<<"t^7 - 5 * t^6 + 10 * t^5 - 10 * t^4 + 5 * t^3 - t^2";
//g8:=R<<"t^7 - 7 * t^6 + 21 * t^5 - 35 * t^4 + 35 * t^3 - 21 * t^2 + 7 * t - 1";
//h8:=R<<"t^7";
//var:=[var0,var1,var2,var3];
singular_point_plane(L: Vector Generic,var: Vector Generic, k: Int
):Vector Generic=={
  list:Vector Generic:=mubase(L);
  m: Int:= max(degree list[0], degree list[1]);
  M1 : Matrix Generic:=matrix(as_polynomial_POMR(list[0][0,0],var[0]),as_polynomial_POMR(list[0][0,1],var[0]),
as_polynomial_POMR(list[0][0,2],var[0]));
 M2 : Matrix Generic:=matrix(as_polynomial_POMR(list[1][0,0],var[0]),as_polynomial_POMR(list[1][0,1],var[0]), as_polynomial_POMR(list[1][0,2],var[0]));
 f:POMR:=(M1*[var[1],var[2],var[3]])[0];
// mmout<<f <<"\n";
 g:POMR:=(M2*[var[1],var[2],var[3]])[0];
// mmout << g<<"\n";
// if degree(f,0)> degree(g,0) then
   sres_polynomial:Vector Generic:= reverse(sturm_seq(f,g,0));
// else  sres_polynomial:Vector Generic:= reverse(sturm_seq(g,f,0));
// mmout<<sres_polynomial<<"\n";
   n:Int:= #sres_polynomial;
   //mmout<< n<<"\n";
 sres:Vector Generic:= [as_polynomial_POMR(0,var[0])| i in 0..(n-1)];
 //mmout<< type(sres[0]:>Generic) <<"\n";
 for i in 0..(n-1) do 
   sres[i]:=eval(sres_polynomial[i][i],[0,L[0],L[1],L[2]]);
   //mmout<< sres<<"\n";
   if k> min(degree(f,0),degree(g,0)) then 
   {if degree(f,0)=degree(g,0) then return([]);
     if degree(f,0)> degree(g,0)  then
       {if k!=degree(f,0) then return([]);
       if k=degree(f,0) then
     {sres_multiply1: Vector Generic:=[as_polynomial_POMR(0,var[0])| i in 0..(degree(g,0)+1)];
         for i in 0..degree(g,0)+1 do sres_multiply1[i]:=eval(sres_polynomial[n-2][i],[0,L[0],L[1],L[2]]);
         //mmout <<sres_multiply1<<"\n";
sres_multiply2:Vector Generic:= [as_polynomial_POMR(0,var[0])| i in 0..(degree(g,0)+2)]; for i in 0..degree(g,0)+1 do sres_multiply2[i]:=eval(sres_polynomial[n-2][i],[0,L[0],L[1],L[2]]);
sres_multiply2[degree(g,0)+1]:=eval(sres_polynomial[n-1][degree(f,0)],[0,L[0],L[1],L[2]]);
 //mmout <<sres_multiply2<<"\n";
 delta:POMR:=gcd(sres_multiply1);
// if degree(delta,0)>=1 then {
 delta:=delta/gcd(delta,derive(delta,var[0]));
   delta:=delta/gcd(delta,sres_multiply2[degree(g,0)+1]);//}
 //mmout<<delta<<"\n";
       }
     }
      if degree(f,0)< degree(g,0)  then
       {if k!=degree(g,0) then return([]);
       if k=degree(g,0) then
     {sres_multiply1: Vector Generic:=[as_polynomial_POMR(0,var[0])| i in 0..(degree(f,0)+1)];
         for i in 0..degree(f,0)+1 do sres_multiply1[i]:=eval(sres_polynomial[n-1][i],[0,L[0],L[1],L[2]]);
         //mmout <<sres_multiply1<<"\n";
sres_multiply2:Vector Generic:= [as_polynomial_POMR(0,var[0])| i in 0..(degree(f,0)+2)]; for i in 0..degree(f,0)+1 do sres_multiply2[i]:=eval(sres_polynomial[n-1][i],[0,L[0],L[1],L[2]]);
sres_multiply2[degree(f,0)+1]:=eval(sres_polynomial[n-2][degree(g,0)],[0,L[0],L[1],L[2]]);
 //mmout <<sres_multiply2<<"\n";
delta:POMR:=gcd(sres_multiply1);
 //if degree(delta,0)>=1 then {
   delta:=delta/gcd(delta,derive(delta,var[0]));
   delta:=delta/gcd(delta,sres_multiply2[degree(f,0)+1]);//}
 //mmout<<delta<<"\n";
       }
     }
    }
   else 
 //sres[n-1]:=eval(f,[0,L[0],L[1],L[2]]);
 //sres[n]:=eval(g,[0,L[0],L[1],L[2]]);
 {sres_multiply1:Vector Generic:= [as_polynomial_POMR(0,var[0])| i in 0..(k+1)];
 for i in 0.. #sres_multiply1 do
   sres_multiply1[i]:=sres[i];
 sres_multiply2:Vector Generic:= [as_polynomial_POMR(0,var[0])| i in 0..k];
 for i in 0.. #sres_multiply2 do
   sres_multiply2[i]:=sres[i];
// mmout<< gcd(sres_multiply2)<<"\n";
 // mmout<< gcd(sres_multiply1)<<"\n";
 if k< min(degree(f,0),degree(g,0)) then
 {delta:=gcd(sres_multiply2);
    //if degree(delta,0)>=1 then {
   delta:POMR:=delta/gcd(delta,derive(delta,var[0]));
   delta:=delta/gcd(delta,sres[k]);//}
   }
 if k=min(degree(f,0),degree(g,0))then
  if degree(f,0)>= degree(g,0) then {
     delta:POMR:=gcd(sres_multiply2);
      //if degree(delta,0)>=1 then {
   delta:=delta/gcd(delta,derive(delta,var[0]));
  // delta:=delta/gcd([delta, eval(sres_polynomial[n-2][degree(g,0)],[0,L[0],L[1],L[2]]), eval(sres_polynomial[n-1][degree(f,0)],[0,L[0],L[1],L[2]])]);}
   list_coeff:Vector Generic:=[as_polynomial_POMR(0,var[0])| i in 0..k+1];
   for i in 0..(k+1) do
     list_coeff[i]:= eval(sres_polynomial[n-2][i],[0,L[0],L[1],L[2]]);
   delta:=delta/gcd(delta,gcd list_coeff);//}
   }
 if degree(f,0)< degree(g,0) then {
   delta:POMR:=gcd(sres_multiply2);
    //if degree(delta,0)>=1 then {
   delta:=delta/gcd(delta,derive(delta,var[0]));
   delta:=delta/gcd(delta, eval(sres_polynomial[n-1][degree(f,0)],[0,L[0],L[1],L[2]]));//}
   }
    }
 // mmout<< delta<<"\n";
  if eval(delta,[0,0,0,0])=0 then {check:=1;
     point0:Vector Generic:=[eval(L[0],[0,0,0,0]),eval(L[1],[0,0,0,0]),eval(L[2],[0,0,0,0])];}
     //mmout<<point0<<"\n";}
  else check:=0;
 // mmout<<check<<"\n";
  while eval(delta,[0,0,0,0])=0 do delta:=delta/var[0];
  //mmout<<delta<<"\n";
  list_coeff:= coefficients delta; 
  //mmout<<list_coeff<<"\n";
  list_denominator:Vector Generic:=[0|i in 0..#list_coeff];
 // mmout<<list_denominator<<"\n";
 // mmout<<type(list_coeff[0]:>Generic)<<"\n";
  for i in 0..#list_coeff do list_denominator[i]:=denominator(list_coeff[i]);
  //mmout<< list_denominator<<"\n";
  if #list_denominator=1 then delta:=list_denominator[0]*delta;
  else   delta:=lcm(list_denominator)*delta;
//  mmout<<delta<<"\n";
 // list_root_interval:Vector Generic:=solve tensor delta;
  //mmout<<list_root_interval<<"\n";
  //mmout<<list_root_interval<<"\n";
// list_root: Vector Generic:=[0| i in 0..#list_root_interval];
 // for i in 0..#list_root_interval do
  // list_root[i]:=as_floating((lower list_root_interval[i]+ upper list_root_interval[i])/2);
  list_root:Vector Generic:=solve delta;
  if check=1 then list_point:Vector Generic:=[point0];
  else list_point: Vector Generic:=[];
  for i in 0..#list_root do
    list_point:=cons([eval(L[0],[list_root[i],0,0,0]),eval(L[1],[list_root[i],0,0,0]),eval(L[2],[list_root[i],0,0,0])],list_point);
return(delete_dif list_point);
};

lcm(L:Vector Generic):Int=={
  lcm1:Int:=lcm(L[0],L[1]);
  for i in 2..#L do
    lcm1:=lcm(lcm1,L[i]);
  return(lcm1);
};
derive(f:POMR,var:POMR):POMR=={
  list:Vector Generic:=coefficients(f,0);
 // mmout<<list<<"\n";
  list_derive:Vector Generic:=[0| i in 0..#list -1];
 // mmout<<list_derive<<"\n";
  if degree(f,0)=0 then return(0*var^0);
  else {
  for i in 0..#list_derive do
    list_derive[i]:=list[i+1]*(i+1);
  // mmout<<list_derive<<"\n";
  p:POMR:=list_derive[0]*var^0;
  for i in 1..#list_derive do
    p:=p+list_derive[i]*var^i;}
  return(p);
};