% This program demonstrates how all complex roots of polynomial p1
% can be found by homotopy from polynomial p0 with known roots.
% This program trying to predict what happens when two roots z1 and z2 collide,
% and jump over the dangereous value T of t, just rotating two close roots
% around their midpoint zmid by 90 degrees 

function value=main()

% p0 is a polynomial with known roots; coefficients are stored in ascending order,
% the first one is the constant term, the last one correspond to the maximal degree
p0 = [0 1 0 1 0.2];
n = length(p0);
deg=n-1;
p=p0; % an intermediate polynomial of the same degree as p0
roots=[0 i -i -5]; % roots of p0
newroots=roots; 
Roots=roots;

% Here are our initial data: polynomial p1 and number of steps
p1 = [2 -3 -2 1 1]; % for this initial data there are two collisions
number_of_steps = 100;
dt=1./number_of_steps;

t=0;
while t<1-0.5*dt,
  t=t+dt; for x = 1:n, p(x) = (1-t)*p0(x) + t*p1(x); end
  for k=1:deg, z=roots(k); dz=P(p,z)/DP(p,z); newroots(k)=z-dz; end
  % ------- let's check if there is a collision and which roots are colliding
  distmin= 10^15; kmin=0; mmin=0; dzmin=0;
  for k=1:deg-1, for m=k+1:deg,
      dz=newroots(k)-newroots(m); dist=abs(dz); 
      if dist<distmin distmin=dist; kmin=k; mmin=m; dzmin=dz; end     
  end;           end
   newdz= dzmin;
   olddz=roots(kmin)-roots(mmin);
   if abs(newdz)>=0.75*abs(olddz), roots=newroots; Roots=[Roots; roots]; end
   if abs(newdz)< 0.75*abs(olddz),
     t=t-dt; for x = 1:n,p(x) = (1-t)*p0(x) + t*p1(x); end % one step back
     z1=roots(kmin)
     z2=roots(mmin)
     Z1=newroots(kmin)
     Z2=newroots(mmin)
     pause 
     for l=1:2,  
       z1=z1-P(p,z1)/DP(p,z1);
       z2=z2-P(p,z2)/DP(p,z2);
     end
     d1=(z1-z2)^2
     t1=t
     t=t+dt/2;
     for x = 1:n, p(x)=(1-t)*p0(x)+t*p1(x); end
     for l=1:2,  
       z1=z1-P(p,z1)/DP(p,z1);
       z2=z2-P(p,z2)/DP(p,z2);
     end
     d2=(z1-z2)^2     
     t2=t
     %     near to the moment of collision T, expression (z1-z2)^2 is 
     %     a linear function in T-t, let T = t2+x, then 
     x=(t2-t1)*d2/(d1-d2);
     x=abs(x);
     T=t2+x
     t3=t2+2*x
     zmid=(z1+z2)/2
     dz=(z1-z2)/2
     newz1=zmid+i*dz    
     newz2=zmid-i*dz
     pause
     t=t3;
     dt=1./(1-t3)/number_of_steps;
     for x = 1:n, p(x) = (1-t)*p0(x) + t*p1(x); end
     roots(kmin)=newz1;
     roots(mmin)=newz2;
     for l=1:2,
       for k=1:deg,z=roots(k);roots(k)=z-P(p,z)/DP(p,z);end
       Roots=[Roots; roots];
     end  
  end%if
end%while
% final improvement of roots using 5 steps of Newton
for l=1:5,  
  for k=1:deg,
  z=roots(k);  
  roots(k)=z-P(p1,z)/DP(p1,z);
  end
  Roots=[Roots; roots];
end
roots=roots % output   
%--------------------------------------------------------
[d,dd]=size(Roots);x=1:d;y=1:d;
for k=1:deg,
  for l=1:d,
    x(l)=real( Roots(l,k) );
    y(l)=imag( Roots(l,k) );
  end
  if k==2
     hold on 
  end
  plot(x,y,'o-');
end
r=2;
axis( [-r,r,-r,r] ); grid;
hold off;
%--------------------------------------------------------
function value = P(a,z)
n= length( a );
value = a(n);
for i = n-1:-1:1,
value = value*z;
value = value + a(i);
end

function value = DP(a,z)
n = length( a );
value = (n-1)*a(n);
for i = n-1:-1:2,
value = value*z;
value = value + (i-1)*a(i);
end