# StandardForm # # Phase 1 (default) # # It does various checks, and it selects a standard form that is unique up # to F-Mobius transformations (F -> 1-F and F -> 1/F). # It makes no x-Mobius transformations x -> (ax+b)/(cx+d). # # Phase 2 (called only if _Env_StandardForm = 2) # # If one sets _Env_StandardForm := 2 then the program will also standardize # the Belyi map F under x-Mobius transformations. The output is unique up # to a finite number of cases, except when j=0 or j=1728. The number of # cases depends on the "Singularity type of the Belyi map" discussed in # comments inside the StandardForm program below. # # This near-uniqueness also pushes the field of definition of F down to the moduli field # except if the "Singularity Type" is "2,2" or "4" (for notations see comments inside # the program). That in turn will be helpful for the Hensel lifting and reconstruction stage. # For the case "2,2" we will often need to rewrite F: conic -> P1 instead of P1 -> P1 # in order to get down to the moduli field. But that will be done in another file. # # Note that setting _Env_StandardForm := 2 may increase the size of the Belyi map, # but it often decreases the size of the 4-sing differential equation. The program DiffEq # computes this differential equation. # StandardForm := proc(f, x, p, L3, L4) local A,B,C, i,j,y,degF, SA,SB,SC,PA,PB,PC,dA,dB,dC,rA,rB,rC,sing,si,deg_ext,deg_j,l,OUT,RT, sel; global _a; if nargs = 1 or not type(x, name) then A := indets(f,name); if nops(A)<>1 then error "f is not a univariate rational function" fi; return procname(f, A[1], args[2..-1]) elif nargs=2 or not type(p, integer) then # We'll assume that the characteristic is 0 when not specified. return procname(f, x, 0, args[3..-1]) elif nargs > 5 and args[6] = "input is not sorted" then # L3 might not be sorted yet A := procname(f, x, p); if nops([A])>3 and [sort(A[3]),sort(A[4])] <> [sort(L3),sort(L4)] then lprint("WARNING: Type differs from what was specified") fi; return A elif nargs > 3 and nops(L3) = 3 then # Sort L3 by increasing denominator if denom(L3[1]) > denom(L3[2]) then return procname(1-f, x, p, [L3[2], L3[1], L3[3]], args[5..-1]) elif denom(L3[2]) > denom(L3[3]) then return procname( 1/f , x, p, [L3[1], L3[3], L3[2]], args[5..-1]) fi fi; C := EVAL( Normal(f), p ); A, B := numer(C), denom(C); C := EVAL( Expand(A-B), p); degF := max( degree(A,x), degree(B,x), degree(C,x) ); if nargs>4 and degF <> (2-convert(L4,`+`))/(1-convert(L3,`+`)) then # Not yet observed: error "Degree of F does not match the degree computed from", L3, L4 fi; SA,SB,SC := seq(SQ(i, x,p,degF,infinity), i=[A,B,C]); if nargs > 3 then dC, dA, dB := op(map(denom, L3)) elif max(seq(i[2], i=SC)) > max(seq(i[2], i=SA)) then return procname(1-f, args[2..-1]) elif max(seq(i[2], i=SA)) > max(seq(i[2], i=SB)) then return procname(1/f, args[2..-1]) fi; PC,dC,rC := CommonExponent(SC, x, dC, 1); PA,dA,rA := CommonExponent(SA, x, dA, 1-1/dC); PB,dB,rB := CommonExponent(SB, x, dB, 1-1/dA-1/dC); sing := [seq( (i[2]/dA)$degree(i[1],x), i=rA), seq( (i[2]/dB)$degree(i[1],x), i=rB), seq( (i[2]/dC)$degree(i[1],x), i=rC)]; sing := sort(sing, L4SORT); if nops(sing)<>4 then # This has happened a number of times: error "Number of singularities, degree, 2F1-expdiffs are:" , nops(sing), degF, [1/dC, 1/dA, 1/dB] elif nargs > 3 and [dC,dA,dB] <> map(denom,L3) then # Not yet observed: error "Type of F does not match", L3, [dA,dC,dB] elif nargs = 3 then return procname(args, [1/dC, 1/dA, 1/dB]) elif nargs > 4 and sing <> sort(L4,L4SORT) then # Not yet observed: error "Singularities of different types", sing, L4 fi; si := {op(sing)}; if nops(si) = 4 then si := "1,1,1,1" elif nops(si) = 3 then si := "1,1,2" elif nops(si) = 1 then si := "4" elif sing[1]=sing[2] and sing[3]=sing[4] then si := "2,2" else si := "1,3" fi; # The string si in { "1,1,1,1", "1,1,2", "4", "2,2", "1,3" } will be called # the "Singularity Type of the Differential Equation". It counts how many time # each exponent-difference occurs. # # We'll also define a "Singularity Type of the Belyi Map" which is similar except # that if two singularities with the same exponent-difference come from different # points in {0,1,infinity} then they be considered as # * the same in "Singularity Type of the Differential Equation" # * different in "Singularity Type of the Belyi Map". # # If L3 = [1/2, 1/3, 1/7] then the two Singularity Types correspond, but if # for instance L3 = [1/2, 1/3, 1/8] then the two Singularity Types need not # correspond because an exponent difference of 1/2 in L_Heun could come either # from 1 * (1/2) or from 4 * (1/8). The difference between these singularity # types could lead to a certain type of "field of definition" issue. j := mul(i[1], i = eval( [op(rA), op(rB), op(rC)], infinity=x-1)); j := algcurves[j_invariant](y^2-j,x,y); j := EVAL(Normal(j), p); deg_ext := `algcurves/degree_ext`(`algcurves/g_ext`(convert(f,RootOf)),{}); deg_j := 1; if indets(j,{RootOf, radical, nonreal}) <> {} then j := evala(Norm( x - convert(j, RootOf) )); j := sqrfree(j)[2][1][1]; deg_j := degree(j,x) fi; if deg_j <> deg_ext then lprint("WARNING: Algebraic degree of j differs from that of F") fi; l := EVAL(Normal(lcoeff(A,x)/lcoeff(B,x)), p); sel := proc() if _Env_StandardForm_Reg = true then args else # Default: don't include the part "Regular points are:" NULL fi end: OUT := l * mul(eval(i[1], infinity=x-1)^i[2], i=SA) / mul(eval(i[1], infinity=x-1)^i[2], i=SB), ["Degree of Belyi map and {Number Field}:", degF, {deg_ext, deg_j}], L3, sing, # ["Singularity type = ", si], # Not really needed [1/dC, "is at F = 1. Branchpoints are:", rC, sel("Regular points are:", PC) ], [1/dA, "is at F = 0. Branchpoints are:", rA, sel("Regular points are:", PA) ], [1/dB, "is at F = infinity. Branchpoints are:", rB, sel("Regular points are:", PB) ], # ["lcoeff =", l], # Not really needed ["j_invariant (or its minpoly) = ", j]; if _Env_StandardForm = 2 then Phase2(p, x, OUT) elif _Env_StandardForm_subs = true then OUT := [OUT]; RT := indets(OUT,RootOf); if nops(RT) = 1 then RT := RT[1]; OUT := subsop(2 = [op(OUT[2]), {_a = RT}], eval(OUT, RT = _a)) fi; op(OUT) else OUT fi end: EVAL := proc(a, p) if p=0 then evala(a) else a mod p fi end: L4SORT := proc(i,j) evalb( denom(i) < denom(j) or (denom(i) = denom(j) and i < j) ) end: SingSort := proc(i,j) evalb( i[2] > j[2] ) end: # subs inf = x-1 to turn (x-inf) into 1. # SQ := proc(A,x,p,degF,inf) local i,j,S; if type(A,list) then for i from 1 to nops(A)-1 do for j from i+1 to nops(A) do if A[i][2] = A[j][2] then return procname(subsop(j = NULL, subsop(i = [ A[i][1]*A[j][1], A[i][2] ], A)), args[2..-1]) fi od od; A else S := EVAL(Sqrfree(A), p)[2]; if degree(A,x) < degF then S := [op(S), [x-inf, degF-degree(A,x)] ] fi; procname(S, args[2..-1]) fi end: CommonExponent := proc(S, x, e, oth) local m,ma,i; m := max(seq(`if`(1/i[2] < oth and (type(e,name) or i[2]=e),degree(i[1],x),NULL), i=S)); ma := NULL; for i to nops(S) do if 1/S[i][2] < oth and degree(S[i][1])=m then ma := ma,i fi od; if nops([ma])>1 then m := sort( [seq(S[i][2], i=[ma]) ]); if oth - 1/m[-2] > 0 then lprint("WARNING: Type of f is not clear from f itself, choosing the max from this list:", m, e) fi; m := m[-1]; if type(e,integer) and m<>e then error "Made the wrong choice" fi; for i in [ma] do if S[i][2]=m then ma := i; break fi od; elif ma=NULL then lprint("WARNING: generic case: None of the points have multiplicity", e); return 1, e, sort(S, SingSort) fi; op(S[ma]), sort(subsop(ma=NULL, S), SingSort) end: Phase2 := proc(p, x, F, degF, L3, L4, C, A, B) local pts, POL, m; pts, POL := FindSingularPoints(p, x, C, A, B, F); # Number of elements of pts is: 0, 1, 2, or 4. # Move first element (if it exists) to infinity # Move second element (it it exists) to 0 # Move third element (if it exists) to 1. # POL = the remaining singular points. # Note that the SingType and the type of the FindPoints # output need not fully correspond. # E.g. 1/3 at x=0 and 1/9 at x=infinity then a pole of order 3 # and a root of order 1 give the same exponent (so you'd have a # number geq 2 in SingType) whereas they are distinct points 1,1 # in FindPoints. This discrepancy could result in the situation # that a Belyi map could potentially descend to a lower field of # definition if we don't let it branch above 0,1,infinity but # above three other points. Whether that actually happens remains # to be seen. # We'll treat the "2,2" and "4" cases as though they were "1,1,2" and "1,3" # cases. So for the "4" case (PointsType, not SingType), we just pick one # of the four points, send it to infinity, and send the rest to the roots # of x^3 - 3*t*x + 2*t. m := FindMobius(p, x, pts, eval(POL, infinity = x-1)); _Env_StandardForm := 1; StandardForm( EVAL( Normal(eval(F, x=m)), p), x, p, L3, L4) end: FindRoots := proc(p,x,V, F, B1) local r,i; if nargs=4 then procname(args, true) elif nops(V)=1 then procname(p,x,[V[1],1], F, false) elif has(V[1],infinity) then r := eval(V[1], infinity = x-1); if degree(r,x)=1 then [infinity, EVAL(solve(r,x), p)], V[2] elif degree(r,x)=3 then [infinity], r else error "this case should not occur" fi elif has(V[2],infinity) then procname(p,x,[V[2],V[1]], F, B1) elif coeff(V[1],x,0)=0 then r := normal(V[1]/x); if degree(r,x)=1 then [0, EVAL(solve(r,x), p)], V[2] elif degree(r,x)=3 then [0], r else error "this case should not occur" fi elif coeff(V[2],x,0)=0 then procname(p,x,[V[2],V[1]], F, B1) else r := EVAL( Roots(V[1], `if`(p=0, indets(F,{RootOf,radical,nonreal}),NULL)), p); r := sort([seq(i[1], i=r)],length); if r=[] then if B1 then procname(p,x,[V[2],V[1]], F, false) else error "Unable to normalize without a field extension" fi elif V[2]=1 then if p=0 then r, 1 else r[1], Normal(V[1]/(x-r[1])) mod p fi else r, V[2] fi fi end: FindSingularPoints := proc(p::nonnegint, x::name, C::list, A::list, B::list, F) local pts, POL, r; pts := [seq(`if`(degree(i[1],x) = 1, i[1], NULL), i = [op(B[3]),op(A[3]),op(C[3])])]; pts := [seq( EVAL( -coeff(i,x,0)/coeff(i,x,1), p), i = pts)]; POL := [seq(`if`(degree(i[1],x)>1, i[1],NULL), i = [op(B[3]),op(A[3]),op(C[3])])]; if pts = [] then lprint("WARNING: 2,2 or 4 ====> 1,1,2 or 1,3. Output not unique"); # By choosing one of the roots, a "4" or a "2,2" case can be turned # into something that looks to the rest of this program as a # "1,1,2" or "1,3" case. # Since there is no unique way to choose a root, we have to print a # message that the output of the program is not unique; the output # depends on an unpredictable choice. return FindRoots(p, x, POL, F) fi; pts, EVAL( Expand(eval(convert(POL,`*`), infinity=x-1)), p) end: # Has to send pts to infinity, 0, 1 and, if there are not that many points, send POL # to a standard form. FindMobius := proc(p, x, pts, POL) local P,m,M,c; if pts[1] <> infinity then if pts[1] = 0 then m := 1/x; else m := x + pts[1] fi; M := procname(p, x, ApplyMobius(m, p, x, pts, POL)); EVAL( Normal(eval(m, x=M)), p) elif nops(pts)>1 and pts[2] <> 0 then m := x + pts[2]; M := procname(p, x, ApplyMobius(m, p, x, pts, POL)); EVAL( Normal(eval(m, x=M)), p) elif nops(pts)>2 and pts[3] <> 1 then m := x * pts[3]; M := procname(p, x, ApplyMobius(m, p, x, pts, POL)); EVAL( Normal(eval(m, x=M)), p) elif nops(pts)>2 then x elif nops(pts) = 2 then # 1,1,2 # We have to map POL to x^2 + 2*t*x + t c := coeff(POL,x,1); if c=0 then lprint("WARNING: 1,1,2 case with (x^2-const)*x (so j=1728). Output not unique"); c := FindScaling( -coeff(POL,x,0)/coeff(POL,x,2), 2, p, x); c * x else c := EVAL( Normal(2*coeff(POL,x,0)/c), p); c * x fi else # 1,3 # We have to map POL to x^3 - 3*t*x + 2*t if degree(POL,x) <> 3 then error "ended up with wrong degree" fi; c := coeff(POL,x,2); if c=0 then if coeff(POL,x,0)=0 then # Reducing 1,3 to 1,1,2 case: return procname(p, x, [infinity,0], collect(POL/x,x)) fi; c := EVAL( -2/3 * coeff(POL,x,1)/coeff(POL,x,0), p); if c=0 then lprint("WARNING: 1,3 case x^3 - const (so j=0). Output not unique"); c := FindScaling( -coeff(POL,x,0)/coeff(POL,x,3), 3, p, x); c * x else EVAL(1/c, p) * x fi else c := EVAL(-c/3/lcoeff(POL,x), p); m := x + c; M := procname(p, x, ApplyMobius(m, p, x, pts, POL)); EVAL( Normal(eval(m, x=M)), p) fi fi end: ApplyMobius := proc(m, p, x, pts, POL) local d,i; d := degree(POL,x); if m=1/x then [seq(`if`(i=0, infinity, `if`(i=infinity, 0, EVAL(1/i , p))), i=pts)], add( coeff(POL,x,i)*x^(d-i), i=0..d) * `if`(nops(pts) + d < 4, x, 1) else d := m - x; if has(d,x) then d := m/x; if has(d,x) then error "unexpected input" fi; # m = d*x [seq( `if`(i=infinity,i,EVAL(i/d, p)), i=pts)], EVAL( Expand(subs(x=x*d, POL)), p) else # m = x + d [seq( `if`(i=infinity,i,EVAL(i-d, p)), i=pts)], EVAL( Expand(subs(x=x+d, POL)), p) fi fi end: DiffEq := proc(f, X, L3) local Dx,x, F, Le, e0, e1, ei, L, v,T,i,j; if not assigned(_Envdiffopdomain) then error "Should assign value to _Envdiffopdomain" elif nargs=1 then i := indets(f,name); if nops(i)<>1 then error "expecting univariate rational function" fi; return procname(f, i[1]) elif nargs=2 then _Env_StandardForm := 1; i := StandardForm( args ); return procname(i[1], X, i[3]) fi; Dx,x := op(_Envdiffopdomain); Le := (x-1)*x*Dx^2+(-1+2*x+e0-e0*x-e1*x)*Dx+1/4*(e0-1+e1+ei)*(-1+e0+e1-ei); # has exponents [0,e0] at x=0, [0,e1] at x=1, and exponent-difference # ei at x=infinity. L := eval(Le, {e0 = L3[2], e1 = L3[1], ei = L3[3]}); F := eval(f, X=x); L := add(DEtools[mult]( subs(x=F,coeff(L,Dx,i)) , (1/diff(F,x) * Dx)$i),i=0..2); L := collect(L/lcoeff(L,Dx), Dx, evala); L := DEtools[symmetric_product](L, Dx - coeff(L,Dx,1)/2); v := factors(denom(evala(coeff(L, Dx, 0))), indets(f,{RootOf,radical,nonreal}))[2]; v := add(min(seq(j[1], j=DEtools[gen_exp](L,T,x=RootOf(i[1],x)))) * diff(i[1],x)/i[1], i=v); L := evala(Primpart( DEtools[symmetric_product](L, Dx + v) , Dx)); L := collect(L,Dx,factor); sort(L,Dx) end: FindScaling := proc(c, n, p, x) local small,i,v; if n = 2 then small := [1,-1, 2,-2, 3,-3, 5,-5, 6,-6, 7,-7, 10,-10, 11,-11, 13,-13, 14,-14, 15,-15] elif n=3 then # small := [1, 2, 1/2, 3, 1/3, 5, 1/5, 6, 1/6, 7, 1/7, 10, 1/10, 11, 1/11, 12, 1/12] small := [1, 2, 3, 4, 9] else error "wrong input" fi; # Want to reduce x^n - c to x^n - i for some in in small. for i in small do v := EVAL( Roots(x^n - c/i), p); if v <> [] then return sort(v,length)[1][1] fi od; lprint("WARNING: FindScaling", x^n-c, "not scaled to", x^n-1); 1 end: