(* This file is a short supplement to "Discriminants of some Painleve 
polynomials".   Depending on your Mathematica system, you can either 

1)  read it in by <<(whatever you call this file) 
or
2)  just cut and paste the file in with the mouse.  

It might be best to keep two windows open: one to look at this file, as there
are comments interspersed throughout; one is your Mathematica session.  

Part I of this file implements the polynomials that the paper is about, namely 
p[m], h[m,n], and q[m,n].  Part II discusses discriminants, resultants and 
complex roots of these polynomials.  Part II also compares the behavior of
Painleve polynomials with "random polynomials" 

Part III is somewhat longer and I don't expect this Part to be of too much 
interest to others.  After writing the paper, I retyped in the main 
formulas to check for typos.   Part III shows these checks.   *)  



(* Part I.  Implementation of the polynomials.  *)

(* In all three cases we make use of the abbreviation b 
(which does not appear in the paper) *)

b[j_,f_] := D[f,{x,2}] f - D[f,x]^2 + j f^2

(* For m=0,1,2,... the Yablonsky-Vorobiev polynomial P_m is given as p[m]. 
Note: here and later, I've kept the code short by not making any attempt to
catch errors.  So if you type p[-1], you'll go into an infinite 
recursion which might crash your session.  Watch out! *)

p[0] = 1

p[1] = x

p[m_] := p[m] = -b[-x,p[m-1]]/p[m-2] //Factor

(* For m,n=0,1,2,... the biHermite polynomial H_m,n is given as h[m,n] *)

h[0,0]=1

h[1,0]=1

h[0,1]=1

h[1,1]= x

h[m_,0] := h[m,0] = b[m-1,h[m-1,0]]/(h[m-2,0]*(m-1)) //Factor

h[m_,1] := h[m,1] = b[m-1,h[m-1,1]]/(h[m-2,1]*(m-1)) //Factor

h[m_,n_] := h[m,n] = b[1-n,h[m,n-1]]/(h[m,n-2]*(1-n)) //Factor 

(* For m,n=0,1,2,... the Okamoto polynomial Q_m,n is given as q1[m,n]. 
For m,n= 0,-1,-2,-3,..., the Okamoto polynomial Q_m,n is given as q2[m,n].  
Both are then also called q[m,n].  
*)

q1[0,0] = 1

q1[1,0] = 1

q1[0,1] = 1

q1[1,1] = x 

q1[m_,0] := q1[m,0] = b[x^2+2 m - 3,q1[m-1,0]]/q1[m-2,0] //Factor

q1[m_,1] := q1[m,1] = b[x^2 + 2 m-2,q1[m-1,1]]/q1[m-2,1] //Factor

q1[m_,n_] := q1[m,n] = b[x^2 + 3 - m - 2 n,q1[m,n-1]]/q1[m,n-2] //Factor

q2[0,0] = 1

q2[1,0] = 1

q2[0,1] = 1

q2[1,1] = x 

q2[m_,0] := q2[m,0] = b[x^2+2 m +1,q2[m+1,0]]/q2[m+2,0] //Factor

q2[m_,1] := q2[m,1] = b[x^2 + 2 m+2,q2[m+1,1]]/q2[m+2,1] //Factor

q2[m_,n_] := q2[m,n] = b[x^2 -1 - m - 2 n,q2[m,n+1]]/q2[m,n+2] //Factor

q[m_,n_] := Which[
m>=0 && n >= 0,q1[m,n],
m<=0 && n<=0,q2[m,n]]




(* Part II. Discriminants, resultants and complex roots.  Comparison of 
Painleve polynomials and random polynomials. *)

(* randpol[d,c] gives a monic degree d polynomial with remaining 
coefficients integers randomly chosen from {-c,-c+1,...,c-1,c} *)

randpol[d_,c_] := x^d + Sum[Random[Integer,{-c,c}] x^k,{k,0,d-1}]

(* If f is a nonzero polynomial in x then disc[f] is its discriminant.  
disc[randpol[10,10]] //FactorInteger
should illustrate what happens with a "randomly chosen" polynomial: at 
least one large prime divides the discriminant.  On the other hand,
disc[p[10]] //FactorInteger
illustrates what happens with the Painleve polynomials here: all primes dividingthe discriminants are small.  
*)

disc[f_] := If[Exponent[f,x]==0,1, 
  If[Mod[Exponent[f, x], 4] <= 1, Resultant[f, D[f, x], x],
   -Resultant[f, D[f, x], x]]]

(* If f and g are polynomials in x then res[f,g] gives their resultant. 
res[randpol[10,10],randpol[10,10]] should illustrate that typically 
at least one large prime divides resultants.  
However, the resultants treated in the paper (and a few others in the 
biHermite and Okamoto cases) are unusual in that they are products of small
primes only .  Examples like res[p[4],p[8]] represent intermediate cases.  
They factor into small primes raised to large powers times large primes 
raised to small powers.  *)

res[f_,g_] := Resultant[f,g,x]

(* If f is a polynomial in x then rootplot[f] graphs its roots 
in the complex plane.  rootplot is designed to work at least somewhat on
systems with "terminal graphics".  So matter what your system, you should 
see a marked difference between e.g. rootplot[h[3,4]] and 
rootplot[randpol[12]] . To get pictures like Figure 6.1 
of the paper, one has to do a bit more formatting. *)    

rootplot[f_] := Module[{lis},
lis = x /. NSolve[f==0,x]; 
ListPlot[Transpose[{Re[lis], Im[lis]}], Axes -> False, Frame -> True,
FrameTicks->False]]

(* Part III.  Here are the checks.  Any command that ends D (for 
"difference") should return 0 for any allowable input.  Commands that end Q
(for "quotient") should return 1 or -1 for allowable inputs.   *)

(* First we check that the rational functions defined in Section 1 really do 
satisfy Painleve 2 or Painleve 4 *)

scalef[m_] := scalef[m] = 
(D[Log[p[m+1]/p[m]],x] /. x->2^(-2/3) x) 2^(-2/3) //Factor

pdiffD[m_] := 
D[scalef[m],{x,2}] - (2 scalef[m]^3 + x scalef[m] -m-1) //Factor

scaleh[m_,n_] := (D[Log[h[m,n+1]/h[m,n]],x] /. x->Sqrt[-2] x) *
Sqrt[-2]  //Factor

hdiffD[m_,n_] := D[scaleh[m,n],{x,2}] -
(D[scaleh[m,n],x]^2/(2 scaleh[m,n]) + 3 scaleh[m,n]^3/2 + 
4 x scaleh[m,n]^2 + 2(x^2-(m+2n+1)) scaleh[m,n] -2m^2/scaleh[m,n]) //Factor

scaleg[m_,n_] := (D[Log[q[m,n+1]/q[m,n]],x] + w x /.  x->Sqrt[-2/3] x) *
 Sqrt[-2/3]  //Factor

qdiffD[m_,n_] := D[scaleg[m,n],{x,2}] -
(D[scaleg[m,n],x]^2/(2 scaleg[m,n]) + 3 scaleg[m,n]^3/2 + 
4 x scaleg[m,n]^2 + 2(x^2-(m+2n)) scaleg[m,n] -
2/9 (3m-1)^2 /scaleg[m,n]) //Factor

(* Next we check that there are no typos in the quadratic relations 
appearing in the paper.  *)

pquad1D[m_] := p[m+1] p[m-1] - 
(D[p[m],x]^2 - D[p[m],{x,2}] p[m] + x p[m]^2) //Expand

pquad2D[m_] := p[m+2] p[m-2] -
((4 - (2 m + 1)^2) D[p[m],x]^2 +
(x^4 + (2m+1)^2 x) p[m]^2 -
4 x^2 p[m] D[p[m],x] + (2m+1)^2 p[m] D[p[m],{x,2}]) //Expand

pquad3D[m_] := p[m+2] p[m-1] - (-(2m+3) D[p[m],x] p[m+1] +
(2m+1) D[p[m+1],x] p[m] + x^2 p[m] p[m+1]) //Expand

hquad1D[m_,n_] := n h[m,n+1] h[m+1,n-1] - 
(h[m,n] D[h[m+1,n],x] - h[m+1,n] D[h[m,n],x]) //Expand

hquad2D[m_,n_] := m h[m+1,n] h[m-1,n+1] - 
(h[m,n] D[h[m,n+1],x] - h[m,n+1] D[h[m,n],x]) //Expand

hquad3D[m_,n_] := D[h[m,n],x]^2 - h[m,n] D[h[m,n],{x,2}] -
(-m(h[m+1,n] h[m-1,n] - h[m,n]^2)) //Expand

hquad4D[m_,n_] := D[h[m,n],x]^2 - h[m,n] D[h[m,n],{x,2}] -
(n(h[m,n+1] h[m,n-1] - h[m,n]^2)) //Expand

qquad1D[m_,n_] := q[m+1,n] q[m-1,n] - 
((x^2 + 2 m + n - 1) q[m,n]^2 + D[q[m,n],{x,2}] q[m,n] - 
D[q[m,n],x]^2) //Factor

qquad2D[m_,n_] := D[q[m,n],x] q[m+1,n] - 
(q[m,n] D[q[m+1,n],x] - x q[m,n] q[m+1,n] + q[m+1,n-1] q[m,n+1]) //Factor

qquad3D[m_,n_] := 
(3(1-m-n)-1) q[m+1,n-1] q[m-1,n+1] +
(3m-1) q[m,n+1] q[m,n-1] +
(3n-1) q[m-1,n] q[m+1,n] //Factor

(* Finally, we check that there are no typos in the discriminant and 
resultant formulas in the main theorems*)

pdiscQ[m_] := disc[p[m]]/Product[j^(j(2m+1-j)^2/4),{j,3,2m-1,2}] 

pres1Q[m_] := Resultant[p[m],p[m-1],x]/
Product[j^(j(2m+1-j)(2m-1-j)/4),{j,3,2m-3,2}]

pres2Q[m_] := Resultant[p[m],p[m-2],x]/
((2m-1)^((m-1)(m-2)/2) Product[j^(j(2m+1-j)(2m-3-j)/4),{j,3,2m-5,2}])

star[-1] = star[0] = starstar[-1] = starstar[0] = 1

star[n_] := star[n] = Product[j^j,{j,1,n}]

starstar[n_] := starstar[n] = Product[star[m],{m,1,n}]

e[m_,n_,j_] := Which[
j <= Min[m,n], j^2 - 2(m-j)(n-j),
Min[m,n] <= j <= Max[m,n],Min[m,n]^2,
Max[m,n] <= j,(m+n-j)^2]

hdisc[m_,n_] := Product[j^(j e[m,n,j]),{j,2,m+n-1}]

hdiscQ[m_,n_] := disc[h[m,n]]/hdisc[m,n] 

hres1Q[m_,n_] := res[h[m,n],h[m+1,n]]/
(hdisc[m,n] starstar[m+n-1]/(starstar[m-1] star[m]^n starstar[n-1]))

hres2Q[m_,n_] := res[h[m,n],h[m,n+1]]/
(hdisc[m,n] starstar[m+n-1]/(starstar[m-1] star[n]^m starstar[n-1]))

hres3Q[m_,n_] := res[h[m,n],h[m+1,n-1]]/
(hdisc[m,n] star[n-1]^m/star[m]^n)

f[x_,y_] := Which[
x >= 2 && y >= 2,Product[j^(j (x-(j+1)/3)(y - (j+1)/3)),{j,2,3 Min[x,y]-4,3}],
x <=-2 && y <= -2,
Product[j^(j(-x-(j-1)/3)(-y-(j-1)/3)),{j,4,3 Min[-x,-y]-2,3}],
True,1]

qdiscQ[m_,n_] := disc[q[m,n]]/(f[m,m] f[n,n] f[1-m-n,1-m-n])

qresQ[m_,n_] := res[q[m,n],q[m+1,n-1]]/(f[m,m+1] f[n,n-1] f[1-m-n,1-m-n])

