(* Choose the polynomial P and the exponent n. Nothing else needs to be changed *)
P = z(z-12)(z-2-8I)(z-8-7I);
n=15;


d = Exponent[P,z];
precision=10*d+n;

(* Finds the union of the roots of all orders of the derivative of P^n. Note that we remove the multiple roots at the zeros of P when m<n in order to speed up the computation of NSolve significantly *)
currentFunc = Expand[P^n];
roots = {};
Do[roots=Join[roots,{Re@#,Im@#}&/@(z/.NSolve[Simplify[currentFunc/(P^(n-m))]==0,WorkingPrecision->precision])];currentFunc=D[Expand[currentFunc],z],{m,1,n-1}]
Do[roots=Join[roots,{Re@#,Im@#}&/@(z/.NSolve[currentFunc==0,WorkingPrecision->precision])];currentFunc=D[Expand[currentFunc],z],{m,n,n*d-1}]

(* At this point, a simple call to ListPlot is enough to draw a basic shadow; the rest of the code below is just there to draw a more interesting plot, and can be omitted *)
(* ListPlot[roots] *) 



(* Determine some additional points for the final plot *)
rootsOfP = {Re@#,Im@#}&/@(z/.NSolve[Expand[P]==0,WorkingPrecision->precision]);
rootsOfPPrime = {Re@#,Im@#}&/@(z/.NSolve[D[Expand[P],z]==0,WorkingPrecision->precision]);
centerOfMassOfZP = Mean[rootsOfP];


(* Determine branch points of equation (1.4) (in the original paper) in the z-plane; a subset of these points serves as an appoximation of the boundary of the shadow *)
resStepLength = 1/100; (* A smaller step length yields more branch points *)
FF = Sum[((α-k)/k!)*D[P,{z,k}]*D[w^(d-k)],{k,0,d}];
RES = Expand[Resultant[FF,D[FF,w],w]];
resultantRoots = Table[({Re@#,Im@#}&/@(z/.NSolve[Simplify[RES]==0,WorkingPrecision->precision])),{α,0+resStepLength,d-resStepLength,resStepLength}];
allResultantRoots = {};
For[i=1,i<Length[resultantRoots],i++,allResultantRoots=Join[allResultantRoots,resultantRoots[[i]]]]; 


(* Additional properties of the final plot *)
rootRadius = Max[0.0015,(14054-Length[roots])/2265000]; (* The radius of the roots of (P^n)^(m) in the plot *)
convHullWidth = (Max[roots[[All,1]]]-Min[roots[[All,1]]]);
convHullHeight = (Max[roots[[All,2]]]-Min[roots[[All,2]]]);
rectSize=Max[convHullWidth,convHullHeight]/600; (* The size of the resultant zeros in the plot *)
critPointsRectSize=Max[convHullWidth,convHullHeight]/150; (* The size of the critical points in the plot *)
s=Max[convHullWidth,convHullHeight]/35;  (* The size of the center of mass of the zeros of P *)


(* Draw the plot, and export it to Directory[] as a PDF file *)
finalPlot = Show[Graphics[{{{Black,PointSize[.02],Point[rootsOfP]}},{{Opacity[0.7],Red,PointSize[rootRadius],Point@#}},{{Opacity[0.4],Blue,PointSize[.0015],Rectangle[{#[[1]]-rectSize,#[[2]]-rectSize},{#[[1]]+rectSize,#[[2]]+rectSize}]&/@allResultantRoots}},{{Opacity[0.6],RGBColor[0/255,175/255,0/255],PointSize[.001],Rectangle[{#[[1]]-critPointsRectSize,#[[2]]-critPointsRectSize},{#[[1]]+critPointsRectSize,#[[2]]+critPointsRectSize}]&/@rootsOfPPrime}},{{Opacity[0.7],RGBColor[0/255,255/255,0/255],PointSize[.0015],Triangle[{{#[[1]],#[[2]]+s/Sqrt[3]},{#[[1]]-s/2,#[[2]]-s/(2*Sqrt[3])},{#[[1]]+s/2,#[[2]]-s/(2*Sqrt[3])}}]&/@{centerOfMassOfZP}}},Black,Dashed},Axes->True,AxesLabel->{"Re","Im"},Frame->False,AspectRatio->Automatic,AxesOrigin->{Min[roots[[All,1]]]-convHullWidth/30,Min[roots[[All,2]]]-convHullHeight/30},PlotRange->{{Min[roots[[All,1]]]-convHullWidth/30,Max[roots[[All,1]]]+convHullWidth/30},{Min[roots[[All,2]]]-convHullHeight/30,Max[roots[[All,2]]]+convHullHeight/30}}]&@roots]

Export[StringJoin["Shadow-n",ToString[n],".pdf"],finalPlot]