// Given a polynomium  P  and its derivative  Q
// Let P_c(z)=P(z)-P(c). 
//    Let R_c denote the roots of P_c,
//    BR_c the boundary of the convex hull of R_c
// Figure shows:
//   In blue : BR_{c_i} for the critical points  c_i
//   In red  : BR_c for c being the cursor position.

I=%i;
z=poly(0,'z');
//P0=z-z^4;
//P0=z^4 - (2-I)*z^2 +z;
P0=z^6 - z;
P0=z^5-z^4+4*z^3-z;
P0=z^6-z^3+z^2-z;
P0=z^6-3*z^3+z  // conjecture ne marche pas
// P0=z^6-z^3+3*z;  // conjecture ne marche pas
Q=derivat(P0);

function drawall(W,st);
  for i=2:dd; for j=1:i-1
      plot2d(real(W([i j])),imag(W([i j])),style=st,rect=B);
  end; end;
endfunction;

//****************************************************************
//   Given d points: W(1),...,W(d)  calculate all the segments of
//   the boundary of the convex hull   Conv{W(1),...,W(d)}
//        Method:  Run through all pairs  W(i),W(j).
//   They determine a halfplane H.
//   Check if all points are on one side of H.
//   If so save the segment [W(i);W(j)]

function L=convex(W)
  d=length(W);L=[];n=0;
  for i=2:d
     for j=1:i-1
	w=W(i)-W(j);       
        ps=imag((W-W(i))*conj(w)); 
	eps=max(ps)*min(ps);  
	if eps > -0.000001 then n=n+1;L(1:2,n)=W([i j]); end;
     end;
  end;
endfunction
//****************************************************************

//rseed=sum(getdate());
//rand("seed",rseed);

//****************************************************************
// calculate  BR_{c_i} for all critical points
Crit= roots(Q);
dd  = length(Crit);
LL  = [];

for i=1:dd;
     vc=horner(P0,Crit(i)); // value of P0 at Crit(i)
     PP=P0-vc;
     PPR=roots(PP);
     L=convex(PPR);
     LL=[LL L];
end;
//****************************************************************

B=2*[-1 -1 1 1];

clf;
xset("pixmap",1);

cx=0;cy=0;
ci=-1;

while modulo(ci,3)==-1,
  P = P0 -horner(P0,cx+%i*cy); RP= roots(P);

  clf;
  xset("wwpc");
  xtitle(sprintf(" cursor  =  (%8.5f)  +   ( %8.5f) i ",cx,cy));

  xset("thickness",1);
  plot2d(real(Crit),imag(Crit),style=-4,rect=B);
  plot2d(real(LL),imag(LL),style=2*ones(LL(1,:)),rect=B);
       
  LP=convex(RP);
  xset("thickness",1);
  plot2d(real(RP),imag(RP),style=-2,rect=B);
  plot2d(real(LP),imag(LP),style=5*ones(LP(1,:)),rect=B);

  xset("wshow");
  rep=xgetmouse();
  ci=rep(3);
  cx=rep(1);cy=rep(2);
end;
