"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),[]);