"load utility routines";

hessian2(f):= block([f11,f12,f22],
        f11: diff(f,z1,2),
        f22: diff(f,z2,2),
        f12: diff(f,z1,1,z2,1),
        ratsimp(determinant(matrix(
                [f11,f12],[f12,f22]
        )))
)$

augmentedhessian2(f,g):= block([f11,f12,f22,g1,g2],
        f11: diff(f,z1,2),
        f22: diff(f,z2,2),
        f12: diff(f,z1,1,z2,1),
        g1: diff(g,z1),
        g2: diff(g,z2),
        ratsimp(determinant(matrix(
                [f11,f12,g1],[f12,f22,g2],[g1,g2,0]
        )))
)$

jacobian2(f,g):= block([f1,f2,g1,g2],
        f1: diff(f,z1),
        f2: diff(f,z2),
        g1: diff(g,z1),
        g2: diff(g,z2),
        ratsimp(determinant(matrix(
                [f1,f2],[g1,g2]
        )))
)$

cross2(f):= ratsimp(
	[-f[2],f[1]]
)$

d2(f):= ratsimp([diff(f,z1),diff(f,z2)])$

newton2(f):= ratsimp(cross2(d2(f)))$

critical2(f):= jacobian2(f[1],f[2])$




promote2(expr,hvar):= block([pret],
	if hvar[1]=z1 or hvar[1]=z2 or hvar[2]=z1 or hvar[2]=z2 then
		error("can't promote z1 or z2"),
	pret:expr,
	pret:subst(savez1,z1,pret),
	pret:subst(savez2,z2,pret),
	pret:subst(z1,hvar[1],pret),
	pret:subst(z2,hvar[2],pret),
	pret
)$

demote2(expr,hvar):= block([pret],
	if hvar[1]=z1 or hvar[1]=z2 or hvar[2]=z1 or hvar[2]=z2 then
		error("can't demote z1 or z2"),
	pret:expr,
	pret:subst(hvar[1],z1,pret),
	pret:subst(hvar[2],z2,pret),
	pret:subst(z1,savez1,pret),
	pret:subst(z2,savez2,pret),
	pret
)$

doit2(hvar,'doitfun,[doitexpr]):= block([doitret],
	doitret:doitexpr,
	doitret:promote2(doitret,hvar),
	doitret:apply(doitfun,doitret),
	doitret:demote2(doitret,hvar),
	doitret
)$



homog2(hvar,var,deg,expr):= ratsimp(
	hvar[2]^deg*subst(hvar[1]/hvar[2],var,expr)
)$

nhomog2(hvar,var,expr):= ratsimp(
	sublis([hvar[1]=var,hvar[2]=1],expr)
)$


discriminant(poly,var):= block([deg,lead],
	deg: hipow(poly,var),
	lead: ratcoeff(poly,var,deg),
	poly_discriminant(poly,var) / lead^(2*deg-2)
)$


prune(expr,var,deg):= block([ratweights,ratwtlvl],
	ratweight(var,1),
	ratwtlvl:deg,
	ratsimp(expr)
)$

ratcoeff2(expr,var1,deg1,var2,deg2):=
	ratcoeff(ratcoeff(expr,var1,deg1),var2,deg2)$


"load general definitions";


powerdisp:true$
algebraic:true$
tellrat(ee^4+ee^3+ee^2+ee+1)$

dw('dwfn,[dwargs]):= apply(doit2, append([[w1,w2],dwfn], dwargs))$
hz(deg,expr):= homog2([z1,z2],z,deg,expr)$
nhz(expr):= nhomog2([z1,z2],z,expr)$
hw(deg,expr):= homog2([w1,w2],w,deg,expr)$
nhw(expr):= nhomog2([w1,w2],w,expr)$

frac(expr):= ratsimp(expr[1]/expr[2])$
nfrac(expr):= [num(expr),denom(expr)]$
compose(f,g,hvar):=sublis([hvar[1]=g[1],hvar[2]=g[2]],f)$
compz(f,g):=compose(f,g,[z1,z2])$
compw(f,g):=compose(f,g,[w1,w2])$

lookatpole(form):= ratsimp(
	compw(compz(form,[z2,z1]),[w2,w1])
)$

lookatone(form):= ratsimp(
	compw(compz(form,[z1+z2,z2]),[w1+w2,w2])
)$



"compute and store icosahedral forms";


f12: ratsimp(
	z1*z2*(z1^10+11*z1^5*z2^5-z2^10)
);

h20: ratsimp(
	1/121 * hessian2(f12)
);

t30: ratsimp(
	1/20 * jacobian2(f12,h20)
);

if ratsimp(t30^2 + h20^3  - 1728 * f12^5) # 0 then
	error("basic identity failed")$


phi31_10: ratsimp(
	-t30 * [z1,z2]
);

phi31_01: ratsimp(
	h20 * newton2(f12)
);

phi31:
	w1*phi31_10 + w2*phi31_01
;


t6: ratsimp(
	z1^6+2*z1^5*z2-5*z1^4*z2^2-5*z1^2*z2^4-2*z1*z2^5+z2^6
);

w8: ratsimp(
	1/400 * hessian2(t6)
);

chi12: ratsimp(
	h20/w8
);



"compute sums of products of forms in z and w";


for i:1 thru 5 do 
	define (t[i](z1,z2),ratsimp(
			compz(t6,[ee^(3*i)*z1,ee^(2*i)*z2])
		))$

for i:1 thru 5 do 
	define (chi[i](z1,z2),ratsimp(
			compz(chi12,[ee^(3*i)*z1,ee^(2*i)*z2])
		))$

tchi: ratsimp(
	sum(t[i](z1,z2)*chi[i](w1,w2),i,1,5)
);

tt: ratsimp(
	sum(t[i](z1,z2)*t[i](w1,w2),i,1,5)
);

ttx2: ratsimp(
	sum(t[i](z1,z2)*t[i](w1,w2)^2,i,1,5)
);



"prepare to reexpress forms in terms of Z";



zexpress(form,deg):= express(form,deg,(-z1+z2)/1728,z2,f12^5,t30^2,5)$

rexpress(form,deg):= express(form,deg,r,1,t6^2,f12,1)$

express(form,deg,var1,var2,form1,form2,gap):=
	block([template,eqns,unkns,sols,g,i],
	local(a),
	template: sum(a[i]*form1^i*form2^(deg-i),i,0,deg),
	form: prune(ev(form,z2=1),z1,gap*deg),
	template: prune(ev(template,z2=1),z1,gap*deg),
	eqns: makelist(
		ratcoeff(form,z1,gap*i)=ratcoeff(template,z1,gap*i),
	      i,0,deg),
	unkns: makelist(a[i],i,0,deg),
	sols: solve(eqns,unkns),
	g: sum(a[i]*var1^i*var2^(deg-i),i,0,deg),
	ratsimp(ev(g,sols))
)$


"reexpress fhat,...";


f5z: zexpress(f12^5,1);
h3z: zexpress(h20^3,1);
t2z: zexpress(t30^2,1);

fhatcheck: ratsimp(hz(6,hw(12,
	( 91125 ) * z^6 +
	( - 133650*w^2 + 61560*w - 193536 ) * z^5 +
	( - 66825*w^4 + 142560*w^3 + 133056*w^2 - 61440*w + 102400 ) * z^4 +
	( 5940*w^6 + 4752*w^5 + 63360*w^4 - 140800*w^3 ) * z^3 +
	( - 1485*w^8 + 3168*w^7 - 10560*w^6 ) * z^2 +
	( 440*w^9 - 66*w^10 ) * z +
	( w^12 )
)))$

fhat: compz(f12, phi31) * f12^4 $
fhat: zexpress(fhat,7)$
fhat: ratsimp(fhat / f5z);

if ratsimp(fhat - fhatcheck # 0) then
	error("fhat is wrong")$

hhat: ratsimp(
	1/(11^2*12^2*h3z*t2z) * dw(hessian2, fhat)
)$

that: ratsimp(
	1/(240*t2z) * dw(jacobian2, fhat, hhat)
)$


if ratsimp(t2z*that^2 + h3z*hhat^3  - 1728 * f5z*fhat^5) # 0 then
	error("hat identity failed")$


"reexpress tchihat,...";


tchihatcheck: ratsimp(hz(7,hw(12,ratsimp( 100*z*(z-1)*(
	( 1215*w - 648 ) * z^4  + 
	( - 540*w^3 - 216*w^2 - 1152*w + 640) * z^3 + 
	( 378*w^5 - 504*w^4 + 960*w^3) * z^2 + 
	( 36*w^7 - 168*w^6) * z +
	( - w^9 )
)))));

tchihat: compw(tchi,phi31) * f12  * t30 $
tchihat: zexpress(tchihat,7)$

if ratsimp(tchihat - tchihatcheck # 0) then
	error("tchihat is wrong")$

tthat: compw(tt,phi31) * f12^4 $
tthat: zexpress(tthat,4)$

ttx2hat: compw(ttx2,phi31) * f12 * t30 $
ttx2hat: zexpress(ttx2hat,7)$

if ratsimp(tchihat - ttx2hat # 0) then
	error("tchihat and ttx2hat don't agree")$


"prepare to compute magic polynomials";

acoeff(i,j):= concat(a_,i,"_",j)$
aterm(i,j):= acoeff(i,j)*z^i*w^j$

killsoff(zap,gap,i,j):= if (zap-j)>(gap*i) then true else false$

livecoeff(i,j):=
	if killsoff(zaproot,gaproot,i,j)
		or killsoff(zappole,gappole,degz-i,degw-j)
		then false else true$

unkncoeffs(degz,degw,zaproot,gaproot,zappole,gappole):=
	block([i,j,l:[]],
	for i:0 thru degz do for j:0 thru degw do
		if livecoeff(i,j) then l:append(l,[acoeff(i,j)]),
	l
)$

unknform(degz,degw,zaproot,gaproot,zappole,gappole):=
	block([i,j,sum:0],
	for i:0 thru degz do for j:0 thru degw do
		if livecoeff(i,j) then sum:sum+aterm(i,j),
	rat(sum)
)$

magi(degz,degw,zaproot,gaproot,zappole,gappole,zapone,gapone):= block(
	[form,unkns,eqs,sols,i,j],
	form:unknform(degz,degw,zaproot,gaproot,zappole,gappole),
	unkns:unkncoeffs(degz,degw,zaproot,gaproot,zappole,gappole),
	oneform:ratsimp(subst([z=1+z,w=1+w],form)),
	eqs:[],
	for i:0 thru degz do for j:0 thru degw do
		if killsoff(zapone,gapone,i,j)
			then eqs:append(eqs,[ratcoeff2(oneform,z,i,w,j)=0]),
	sols:solve(eqs,unkns),
	form:ev(form,sols),
	form:content(form)[2],
	form:hz(degz,hw(degw,form)),
	form
)$
	

"make sure the magic works";

fmagic:magi(6,12,12,3,12,2,11,5);

if (ratsimp(fhat-fmagic)#0)
	then error("fmagic failed")$

tchimagic:magi(4,12,9,3,11,2,8,5);

if (ratsimp(tchihat-ratsimp(100*z1*z2*(z2-z1)*tchimagic))#0)
	then error("tchimagic failed")$

"test root finding";


top:factor(nhz(nhw(tchihat)));
bottom:nhz(nhw(fhat));

alpha:(3*sqrt(15)*( %i) + 5) / 2;
beta :(3*sqrt(15)*(-%i) + 5) / 2;
root1:ratsimp((alpha*w1-beta*w2)/(alpha^2-beta^2));
root2:ratsimp((alpha*w2-beta*w1)/(alpha^2-beta^2));
poly(r,zd):=r^5-10*zd*r^3+45*zd^2*r-zd^2;
zdtoz(zd):=1-1728*zd;

modnewt(f,v,d):=(v*diff(f,v)-d*f)/diff(f,v);
compute(expn,sub):=block([],return(rectform(ev(expn,sub))));
fpprec:32;
fpprintprec:16;

zdval:0.1b0;
"We are trying to solve...";
poly(r,zdval) = 0;

zval:zdtoz(zdval);
bottomval:compute(bottom,z=zval)$
topval:compute(top,z=zval)$

iter:modnewt(bottomval,w,12)$
iter:bfloat(ev(iter))$
expand(num(iter))/expand(denom(iter))$
define(f(w),%);

wval:1.0b0;
for i:1 thru 10 do display(wval:f(wval));

w1val:compute(topval/bottomval,w=wval);
wval:f(wval);
w2val:compute(topval/bottomval,w=wval);

root1:compute(bfloat(root1),[w1=w1val,w2=w2val]);
root2:compute(bfloat(root2),[w1=w1val,w2=w2val]);

compute(poly(root1,zdval),[]);
compute(poly(root2,zdval),[]);