The following block of code creates a single static plot of all the roots for the family of polynomials given by z^n + z^k - 1, where 1 <= k <= n - 1. The user may input different values for the degree n. The plot includes the unit circle and the nth roots of unity (displayed in red).  The roots of unity may be deleted from the plot by deleting "g" from the "Show" command in the in the last line. Similarly, the unit circle may be deleted from the plot by deleting "b" from the "Show" command.  The plot may be further modified by changing the various plot options dealing with axes, fonts and colors.  See the documentation for plot options in the online help.

Clear["Global`*"];
n=20; (* Set the degree of the polynomial *)
rootdata=Table[{0.0, 0.0},{n*(n-1)}]; (* Create a table to hold the roots of the polynomial*)
i=1;
For[k=1, k≤n-1, k++,
For[root=1, root≤n, root++, z=N[Root[x^n+ x^k-1,root]];
rootdata[[i]]={Re[z], Im[z]};
i++]]
a=ListPlot[rootdata, AspectRatio->Automatic, AxesLabel->{x, y}, Frame->False, PlotStyle->{PointSize[.007],Blue}, Background->White,Axes->False, PlotRange->All];
b=Graphics[Circle[{0,0},1]];
g=ListPlot[Table[{Cos[2 Pi j/n], Sin[2 Pi j/n]}, {j, 0, n-1}], PlotStyle->{PointSize[.01], Red}];
c=Show[a,b, g]

The block of code below takes the rootdata table from above and creates an animiation.  It plots the complete set of roots for one polynomial in white.  It then plots the complete set of roots for the 2nd polynomial in white and changes the first set to magenta. Each subsequent polynomial’s roots is plotted in white with all “older” sets displayed in magenta.  It cycles through until all the roots from the entire family are displayed.  It then repeats.

table=Table[Show[ListPlot[Table[rootdata[[i]],{i, 1, m*n}],Axes->False,Background->Black,PlotStyle->{PointSize[.002], Magenta},AspectRatio->Automatic],
ListPlot[Table[rootdata[[i]],{i, (m-1)*n+1, m*n}],PlotStyle->{PointSize[.007],White}], PlotRange->{{-1.2,1.2}, {-1.2,1.2}}],{m, 1, n-1}];
ListAnimate[table, AnimationRunning->False]

The following block of code creates a single static plot of all the roots for the family of polynomials given by z^n +k^2 z^k - 1, where 1 <= k <= n - 1. While this family of polynomials was not addressed in the paper Locating Unimodular Roots, it provides a stunning example of the intricate patterns that may be obtained. The user may input different values for the degree n, as well as experimenting with the coefficient of the z^k term (here we used k^2 as it gives aestetically pleasing images).   The plot includes the unit circle and the nth roots of unity (displayed in red).  The roots of unity may be deleted from the plot by deleting "g" from the "Show" command in the in the last line. Similarly, unit circle may be deleted from the plot by deleting "b" from the "Show" command.  The plot may be further modified by changing the various plot options dealing with axes, fonts and colors.  See the documentation for plot options in the online help.

Clear["Global`*"];
n=60; (* Set the degree of the polynomial *)
rootdata=Table[{0.0, 0.0},{n*(n-1)}]; (* Create a table to hold the roots of the polynomial*)
i=1;
For[k=1, k≤n-1, k++,
For[root=1, root≤n, root++, z=N[Root[x^n+ k^2 x^k-1,root]];
rootdata[[i]]={Re[z], Im[z]};
i++]]
a=ListPlot[rootdata, AspectRatio->Automatic, AxesLabel->{x, y}, Frame->False, PlotStyle->{PointSize[.007],Blue}, Background->White,Axes->False, PlotRange->{{-1.5, 1.5}, {-1.5, 1.5}}];
b=Graphics[Circle[{0,0},1]];
g=ListPlot[Table[{Cos[2 Pi j/n], Sin[2 Pi j/n]}, {j, 0, n-1}], PlotStyle->{PointSize[.01], Red}];
c=Show[a,b, g]

The block of code below takes the rootdata table from above and creates an animiation.  It plots the complete set of roots for one polynomial in white.  It then plots the complete set of roots for the 2 nd polynomial in white and changes the first set to magenta. Each subsequent polynomial' s roots is plotted in white with all "older" sets displayed in magenta.  It cycles through until all the roots from the entire family are displayed.  It then repeats.

table=Table[Show[ListPlot[Table[rootdata[[i]],{i, 1, m*n}],Axes->False,Background->Black,PlotStyle->{PointSize[.002], Magenta},AspectRatio->Automatic],
ListPlot[Table[rootdata[[i]],{i, (m-1)*n+1, m*n}],PlotStyle->{PointSize[.007],White}], PlotRange->{{-1.5,1.5}, {-1.5,1.5}}],{m, 1, n-1}];
ListAnimate[table, AnimationRunning->False]

The program below counts the number of roots inside, outside and on the unite circle for the family of polynomials z^n+z^k-1 , where 1<=k<=n-1. The results are displayed in a table.  At the bottom of the table the program outputs the maximum number of roots that occur inside, outside , and on the unit circle. Finally it gives the maximum ratio of the number of roots occuring inside the unit circle compared to the degree of the polynomial. The user may change the value of n to experiment with different types of integers such as primes or numbers with certain factors.

Clear["Global`*"];
n=30;
jstart=n;
For[j=jstart, j ≤n, j++, 
maxin=0; maxout=0;maxon=0;
For[k=1,k≤j-1, k++, in=0; out=0; on=0;
For[root=1, root≤j, root++, z=N[Root[x^j+x^k-1, root]]; 
If[.9999999999<Abs[z]<1.0000000001, on++, If[Abs[z]>1.0000000001, out++, in++]]
]
If[in>maxin, maxin=in];
If[out>maxout, maxout=out];
If[on>maxon, maxon=on];
Print["n=", j, "  k=", k, "  On =", on,"  Inside=",  in, "     Outside=",  out]
]
Print["max in =", maxin, "   max out=", maxout, "   max on=", maxon, "   max ratio=", 1.0*maxin/n]
]