# The Nearest Singular Algorithm # 
# Input: Monic polynomial f  in C[x] of degree m, multiplicity kk #
# Output: Monic polynomial h in C[x] and a root c in C such that c is #
#          a multiple root of  h with given multiplicity  


nspa:=proc(f,kk)
local
q1,q2,p2,i,m,dq1,dp2,p22,Rp,Ip,delta,GL,RL1,j,k,cc,h,Nm,a,b,RL2,cc1,p3,p33,h2,q22,q3,p4,p44,h3,GL2,
NM,CC1,CC,Nmm,d,d1,q4,p5,t;
 m:=degree(f);
NM:=1000;
CC1:=0;
CC:=0;
 q1:=0;
 for i  to m do 
  q1:=q1+(c*c1)^(i-1);
 od;
 dq1:=diff(q1,c);
 q2:=q1*diff(dq1,c1)-dq1*subs({c=c1,c1=c},dq1);
 p2:=subs(x=c,q1*diff(f,x)-dq1*f);
 dp2:=diff(p2,c);
 if kk=2 then 
    p22:=expand(subs({c=a+I*b,c1=a-I*b},p2));
    p3:=simplify((q2*diff(p2,c)-p2*diff(q2,c))/q1);
    d:=q1^4*p3*subs({c=c1,c1=c},evalc(conjugate(p3)))-q2^4*subs(x=c,f)*subs(x=c1,evalc(conjugate(f)));
    Rp:=coeff(p22,I,0);
    Ip:=coeff(p22,I,1);
    Nm:=subs(x=c,f)*subs(x=c1,evalc(conjugate(f)))/q1;
 elif kk=3 then
     p3:=simplify((q2*diff(p2,c)-p2*diff(q2,c))/q1);
     q3:=simplify((q2*diff(diff(q2,c),c1)-diff(q2,c)*diff(q2,c1))/q1);
     p4:=simplify((q3*diff(p3,c)-p3*diff(q3,c))/q2);
     q22:=q1*diff(subs(c=x,q1),c1)-diff(q1,c1)*subs(c=x,q1);
    p33:=expand(subs({c=a+I*b,c1=a-I*b},p3));
    d:=q2^4*p4*subs({c=c1,c1=c},evalc(conjugate(p4)))-q3^4*p2*subs(x=c1,evalc(conjugate(p2)));
    Rp:=coeff(p33,I,0);
    Ip:=coeff(p33,I,1);
    Nm:=subs(x=c,f)*subs(x=c1,evalc(conjugate(f)))/q1+p2*subs({c=c1,c1=c},p2)/(q1*q2);
 
 elif kk=4 then 
 
    p3:=simplify((q2*diff(p2,c)-p2*diff(q2,c))/q1);
    q3:=simplify((q2*diff(diff(q2,c),c1)-diff(q2,c)*diff(q2,c1))/q1);
    p4:=simplify((q3*diff(p3,c)-p3*diff(q3,c))/q2);
    q22:=q1*diff(subs(c=x,q1),c1)-diff(q1,c1)*subs(c=x,q1);
    q4:=simplify((q3*diff(diff(q3,c),c1)-diff(q3,c)*diff(q3,c1))/q2);
    p5:=simplify((q4*diff(p4,c)-p4*diff(q4,c))/q3);
    d:=q3^4*p5*subs({c=c1,c1=c},evalc(conjugate(p5)))-q4^4*p3*subs(x=c1,evalc(conjugate(p3)));
    p44:=expand(subs({c=a+I*b,c1=a-I*b},p4));
    Rp:=coeff(p44,I,0);
    Ip:=coeff(p44,I,1);
    Nm:=subs(x=c,f)*subs(x=c1,evalc(conjugate(f)))/q1+p2*subs({c=c1,c1=c},p2)/(q1*q2)+ p3*subs({c=c1,c1=c},p3)/(q2*q3);
else
   ERROR(`failture for k>4`);
fi;
   # delta:= e+length(4*m^8*(1+2*norm(f,infinity))^(4*m))+1;
    GL:=[[factor(Rp),factor(resultant(Rp,Ip,a))]];
    #lprint(`nop of GL`, GL);
    #readlib(realroot);
    for i from 1 to nops(GL) do
      #lprint(op(2,op(i,GL)));
     # RL1:=realroot(op(2,op(i,GL)),10^(-delta*degree(op(1,op(i,GL)),a)));
      GL2:=op(2,op(i,GL));
      RL1:=[];
      if type(GL2,'`+`' )=true then 
         RL1:=[fsolve(evalf(GL2),b,real)];
      else
         for t from 1 to nops(GL2) do 
           if degree(op(t,GL2))>0 then 
               RL1:={op(RL1), fsolve(evalf(op(t,GL2)),b,real)};
            fi;
          od;
          RL1:=[op(RL1)];
       fi;
       
       
      #RL1:=[fsolve(op(2,op(i,GL)),b,real)];
     # RL1:=[fsolve(evalf(GL2),b,real)];
      #lprint(`RL1=`,evalf(RL1));
      for j  from 1 to nops(RL1) do
        #lprint(`check`, op(1,op(i,GL)),evalf(subs(b=1/2*(op(1,op(j,RL1))+op(2,op(j,RL1))),op(1,op(i,GL)))));
       # RL2:=realroot(subs(b=1/2*(op(1,op(j,RL1))+op(2,op(j,RL1))),op(1,op(i,GL))),10^(-delta));
        RL2:=[fsolve(subs(b=op(j,RL1),op(1,op(i,GL))),a,real)];       
        #lprint(`RL2=`,evalf(RL2));
        for k from 1 to nops(RL2) do 
          # d1:=evalf(subs({b=1/2*(op(1,op(j,RL1))+op(2,op(j,RL1))),a=1/2*(op(1,op(k,RL2))+op(2,op(k,RL2)))},subs({c=a+I*b,c1=a-I*b},d )));
          d1:=evalf(subs({b=op(j,RL1),a=op(k,RL2)},subs({c=a+I*b,c1=a-I*b},d)));
           #lprint(`d1=`,d1);
          #lprint(`d21=`,evalf(d21));
          if coeff(d1,I,0) >= 0 then
             cc:=op(k,RL2)+op(j,RL1)*I;
            # cc:=1/2*(op(1,op(k,RL2))+op(2,op(k,RL2)))+1/2*(op(1,op(j,RL1))+op(2,op(j,RL1)))*I;
             #lprint(`cc=`, evalf(cc));
             cc1:=evalc(conjugate(cc));
             #lprint(`cc1=`,evalf(cc1));
             #Nm:=subs(x=cc,f)*subs(x=cc1,evalc(conjugate(f)))/subs({c=cc,c1=cc1},q1);
             #lprint(`Nm=`, Nm);
             Nmm:=subs(c=cc, c1=cc1, Nm);
             #lprint(`Nmm=`, evalf(Nmm), evalf(NM));
             if abs(Nmm) < NM then
                NM:=abs(Nmm);
                CC:=cc;
                CC1:=cc1;
             fi;
          fi;
        od;
     od;
  od;
if NM <1000 then 
   if kk=2 then 
      h:=f-subs(x=CC,f)/subs({c=CC,c1=CC1},q1)*subs({c=x,c1=CC1},q1);
   elif kk=3 then 
      h2:=f-subs(x=CC,f)/subs({c=CC,c1=CC1},q1)*subs({c=x,c1=CC1},q1);
      h:=h2-subs(x=CC,diff(h2,x))/subs({c=CC, c1=CC1},q2)*subs({c=CC,c1=CC1},q22);
   elif kk=4 then
      h2:=f-subs(x=CC,f)/subs({c=CC,c1=CC1},q1)*subs({c=x,c1=CC1},q1);
      #lprint(`q22=`,q22);
      h3:=h2-(subs(x=CC,diff(h2,x)))/subs({c=CC, c1=CC1},q2)*subs({c1=CC1,c=CC},q22);
      h:=h3-subs(x=CC,diff(diff(h3,x),x))
       *subs({c=CC,c1=CC1}, q2*diff(q22,c1)-diff(q2,c1)*q22)/subs({c=CC,c1=CC1},q1*q3);
   else
      ERROR(`failture`);
   fi:
RETURN([evalf(CC),evalf(NM),evalf(h)]);
else
ERROR(`failture`);
fi;

end: