(* Instructions: Make sure to input the functions in this first cell (use Shift-Enter).  They should turn black in color (previously blue);

The functions defined here are

roots[p]
poly[set]
double[set]

minplus[set1, set2]
maxplus[set1, set2]
randplus[set1, set2]

mintimes[set1, set2]
maxtimes[set1, set2]
randtimes[set1, set2]
 
You can type ?function to get a description of what the function does.  There are three examples below.
*) 


(* Go back and forth between sets and polynomials *)

roots::usage = "roots[p] finds the roots of a polynomial p.  The default variable is x.
If you want a different one (say y) use roots[p, y].";
roots[p_, x_:x]:=(x /. #1)&/@ Solve[p==0, x]

poly::usage = "poly[set] takes set (a list of numbers) and makes the polynomial with those numbers as roots.  The default variable is x.  If you want a different one (say y) use poly[set, y].";
poly[set_, x_:x]:=Times @@ (x - set) // Expand

double::usage = 
"double[set] doubles the elements in the set, so (for example) 
   double[{a, b, c}] = { a, a, b, b, c, c }.";
double[set_]:= Join[set, set]

(* Addition *)

minplus::usage = "minplus[set1, set2] does the minimum lemscourge addition of set1 and set2.  The sets need to be the same size.";
minplus[set1_, set2_] := Sort[set1] + Reverse[Sort[set2]]

maxplus::usage = "maxplus[set1, set2] does the maximum lemscourge addition of set1 and set2.  The sets need to be the same size.";
maxplus[set1_, set2_] := Sort[set1] + Sort[set2]

randplus::usage="randplus[set1, set2] does the random lemscourge addition of set1 and set2.  That is, if 
  X_A is a random symmetric matrix with herm(X_A) = A, 
  X_B is a random symmetric matrix with herm(X_B) = B, 
then this returns the roots of the polynomial 
  E Det[t I - X_A - X_B].
The sets need to be the same size.";
randplus[set1_, set2_] := Module[{p1, p2, d, p3, t}, 
p1 = poly[set1, t];
p2 = poly[set2, t];
d = Length[set1];
p3=Sum[D[p1, {t, i}]*(D[p2, {t, d-i}] /. t->0), {i, 0, d}]/d!;
roots[p3 // N, t] // Chop
] 

(* Multiplication *)

mintimes::usage = "mintimes[set1, set2] does the minimum lemscourge multiplication of set1 and set2.  The sets need to be the same size.";
mintimes[set1_, set2_] := Sort[set1] * Reverse[Sort[set2]]

maxtimes::usage = "maxtimes[set1, set2] does the maximum lemscourge multiplication of set1 and set2.  The sets need to be the same size.";
maxtimes[set1_, set2_] := Sort[set1] + Sort[set2]

randtimes::usage="randtimes[set1, set2] does the random lemcourge multiplication of set1 and set2.  That is, if 
  X_A is a random symmetric matrix with herm(X_A) = A, 
  X_B is a random symmetric matrix with herm(X_B) = B, 
then this returns the roots of the polynomial 
  E Det[t I - X_A X_B].
The sets need to be the same size.";
randtimes[set1_, set2_]:= Module[{p1, p2, d, p3, t},
p1 = CoefficientList[poly[set1, t],t];
p2 = CoefficientList[poly[set2, t], t];
d = Length[set1];
p3 = p1 * p2  * ((-1)^(d-#)* t^#  / Binomial[d, #] & /@ Range[0, d]);
roots[Total[p3] // N, t] //Chop
]

(* Example 1: 
Write the definition for randplus.
*)
?randplus

(* Example 2:
Use minplus, maxplus, and randplus on 2 sets, and then plot the resulting polynomials
*)
set1 = {1, 2, 3};
set2 = {1, 1, 2};
p1=poly[minplus[set1,set2]]
p2=poly[randplus[set1, set2]]
p3=poly[maxplus[set1, set2]]

Plot[{p1, p2, p3}, {x, 2, 5}]

(* 
Example 3: form 2 sets, and compares the two polynomials;
(1) doubles both sets and then randplus,
(2) randplus the sets and then double;

Proving that (1) ≤ (2) is an open conjecture.
*)
set1 = {1, 2, 3};
set2 = {1, 2, 4};
p1 = poly[randplus[double[set1], double[set2]]];
p2 = poly[double[randplus[set1, set2]]];
Plot[{p1, p2}, {x, 2, 7}]