with(DEtools): _Envdiffopdomain := [Dx,x]: # Input: p4 is a square-free degree 4 polynomial in x. Let E: y^2 = p4 and let P1,P2 be the # two points at x=infinity. The divisor P1 - P2 must be torsion of order n = N/igcd(N,2), # otherwise the program will give an error "Divisor P1-P2 does not have order n". # If this torsion condition is met, then the program returns a polynomial F0 such # that F = F0 - sqrt(F0^2-1) = F0 - sqrt(p4) * polynomial # has precisely one pole of order n and one root of order n, at the points P1, P2, # and no other poles or roots. # TorsionCondition := proc(p4,x,n) local a,b,SQ,F,i,F0,ISZ,sol; options remember; if degree(p4,x) <> 4 then error "not implemented" fi; if lcoeff(p4,x) <> 1 then return procname( collect(p4/lcoeff(p4,x),x,evala), x, n) fi; SQ := p4sqrt(p4, x, 2*n); F0 := 1 + add(a[i]*x^i, i=1..n); F := F0 - SQ * add(b[i]*x^i, i=0..n-2); ISZ := convert(series(F, x=infinity, 2*n), polynom); sol := SolveTools:-Linear({seq(coeff(ISZ, x, i), i=-n+1..n)}, {seq(a[i], i=1..n), seq(b[i],i=0..n-2)}); F0 := collect(eval(F0, sol),x,factor); # Check if the expected conditions are met: if F0 = 0 then error "Divisor P1-P2 does not have order n" elif evala(p4 * eval(add(b[i]*x^i, i=0..n-2), sol)^2 - (F0^2-1)) <> 0 then error "unexpected" fi; F0 end: # This program was needed because series(sqrt(p4), x=infinity, 2*n) was too slow. p4sqrt := proc(p4,x,n) local p2,i,c; p2 := x^2; for i from 1 to -n by -1 do p2 := p2 + x^i * evala(solve(coeff((p2 + c*x^i)^2 - p4, x, 2+i),c)) od; end: # Input: Same p4 as in TorsionCondition. Such p4 is constructed from a point on X1(N). # Output: operator whose monodromy is expected to be O_zeta(N). # L_NP := proc(p4,x,N) local n,F0,qt,cc,p4m,P4; n := N/igcd(N,2); # Note: if we replace n with a multiple of n, then the output stays the same. F0 := TorsionCondition(p4,x,n); qt := numer(evala(diff(F0,x)/((F0-1)*(F0+1)))); if degree(qt,x) <> 1 then error "unexpected" fi; cc := evala(solve(qt,x)); p4m := p4/lcoeff(p4,x); P4 := sort(collect(subs(x=4*x, p4m) / 4, x, evala), x); # If we apply a change of variables x -> F0 to: # # L_Dn := (x^2-1)*Dx^2 + x*Dx - 1/(4*n^2) # Dihedral Galois group # # then we get this equation: # # p := p4m/(x-cc)^2; L := p * Dx^2 + diff(p,x)/2 * Dx - 1/4; # # with this solution: # # S := R^(1/(2*n)) + R^(-1/(2*n)) where R := F0 + sqrt(F0^2-1) # # If Z := S/sqrt(p4/x) = Sum a_i * x^i then Sum binomial(2*i,i) * a_i * x^i satisfies: P4 * Dx^2 + diff(P4,x) * Dx + diff(P4,x,x)/8 - (4*x - cc)^2 end: # Polynomials for X1(N) for N = 3..12 except 11, because 11-torsion requires an algebraic # extension (a "RootOf" in Maple notation). An example for N = 11 appears later in this file, # as well as an example for N = 37 and N = 21. # # The "subs" scales it to end up with no denominator in the series-solution. # P4[3] := subs(c1 = 1/3^3, x^4+c3*x^3+1/4*c3^2*x^2+c1*x); P4[4] := subs(c3 = 1/2, c2 = 1/64/c2+1/16, x^4+c3*x^3+c2*x^2+1/8*(-c3^2+4*c2)*c3*x); P4[5] := subs(c3 = (1-1/c)/5^3, x^4+c3*x^3+1/4*c3^2*(c^2-6*c+1)/(c-1)^2*x^2+1/2*c*c3^3/(c-1)^3*x); P4[6] := subs(c3 = (r-2)/(r^2-r)/6^3, x^4+c3*x^3-1/4*(3*r^2-4)*c3^2/(r-2)^2*x^2+1/2*c3^3*r*(r-1)/(r-2)^3*x); P4[7] := subs(c3 = 1/7^3 * 1/(r^3-r^2), x^4+c3*x^3+1/4*(r^4-6*r^3+3*r^2+2*r+1)*c3^2/(r^2-r-1)^2*x^2+1/2*r^2*c3^3*(r-1)/(r^2-r-1)^3*x); P4[8] := subs(c3 = 2^(-7)*(s^2-2)/( s*(s-1)*(s-2) ), x^4+c3*x^3+1/4*(s^4-8*s^2+4*s+4)*c3^2/(s^2-2)^2*x^2-1/2*s*c3^3*(s-1)*(s-2)/(s^2-2)^3*x); P4[9] := subs(c3 = 3^(-5) * (s^3-s^2-1)/( s^2 * (s-1) * (s^2-s+1)), x^4+c3*x^3+1/4*(s^6-6*s^5+9*s^4-10*s^3+6*s^2+1)*c3^2/(s^3-s^2-1)^2*x^2+1/2*s^2*c3^3*(s-1)*(s^2-s+1)/(s^3-s^2-1)^3*x); P4[10] := subs(c3 = 10^(-3) / s^5 / (s-1) / (s+1)^2 / (s^2+s-1), x^4+c3*s^2*(s^3+3*s^2-s-1)*x^3+1/4*s^4*(s^6+10*s^5+7*s^4-12*s^3-5*s^2+2*s+1)*c3^2*x^2+1/2*(s+1)^2*(s^2+s-1)*s^9*(s-1)*c3^3*x); P4[12] := subs(c3 = 2^(-5)*3^(-3) / (s^9-s^5) / (s-1) / (s^2-s+1), x^4+c3*(s-1)*(s^4+2*s^3-2*s^2+2*s-1)*x^3+1/4*(s^8+8*s^7-8*s^6+8*s^5-2*s^4-4*s^3+4*s^2-4*s+1)*(s-1)^2*c3^2*x^2 +1/2*s^5*(s-1)^4*(s+1)*(s^2-s+1)*(s^2+1)*c3^3*x ); for N from 3 to 10 do L[N] := L_NP(P4[N], x, N); od; L[12] := L_NP(P4[12], x, 12); # Remark: if N is odd, then L[N] and L[2*N] have the same monodromy O_zeta(N) = O_zeta(2*N). # So why don't these operators look the same? # When we construct p4 from X1(2*N) then the elliptic curve y^2 = p4 has 2-torsion, meaning # that p4 has another linear factor (in addition to the factor x). # So factor(lcoeff(L[2*N], Dx)) will have two linear factors in x, whereas factor(lcoeff(L[N], Dx)) # will only have one linear factor, namely x. # Indeed, if you take the factor of lcoeff(L[5],Dx) of degree 3 in x, parametrize and substitute to # make it reducible, then you get an equation that is equivalent to L[10]. # If we want to construct the hyperelliptic integral, then we have to parametrize # the equation: F0 = orthopoly:-T(n,T) # where F0 is the output of the program TorsionCondition. # # Since the algebraic relation that we have to parametrize is of the form: # # Polynomial in one variable = Polynomial in another variable # # it is not difficult to prove that this relation has genus 0. For instance, # the polynomial F0 branches above F0 = 1, F0 = -1, F0 = cc (from program L_NP) # and F0 = infinity. Also, F0^2 - 1 is a square times p4. After analyzing # orthopoly:-T(n,T) as well, we can find all the singularities, and we can also # find all the e-1's in the Riemann-Hurwitz formula. # # This proves that the genus is always 0, so there's a parametrization, and we # get a hyperelliptic integral just like in the paper. To avoid computing with sqrt's, # the program will return the square of the integrand: # IntegrandSquare := proc(p4, x, N, v) local Np, Pn, T, F0, Relation, dBas, NumGen, para, IS, i; Np := N/igcd(N,2); Pn := orthopoly:-T(Np,T); # Remark: with L and L_Dn from the comments in L_NP, a change of variables x -> Pn(x) # sends L_Dn to L_D1 = (x^2-1)*Dx^2 + x*Dx - 1/4, which is L with p = x^2-1. F0 := TorsionCondition(p4, x, Np); Relation := F0 - Pn; # Genus 0 by Riemann-Hurwitz, need to parametrize this. # Singular locus = denominator in algcurves[integral_basis](Relation, x, T) = dBas := gcd(diff(F0,x), (F0-1)*(F0+1)); # The following generates all functions with at most one finite pole, of order 1, # located at (x,T) = (0,1). We multiply by dBas*x to remove denominators: {Relation, diff(Pn,T) * (T^2-1) * x, dBas * x, normal((Pn-1)/(T-1)) * dBas}; NumGen := Groebner:-Basis(%, tdeg(T,x))[2]; # NumGen/(dBas*x) has only one pole, of order 1, at the point (x,T) = (0,1). # Hence it generates the function field. Now write x and T as functions of it: para := `algcurves/expr_in_p`(Relation, x, T, v, NumGen, dBas*x); # Note: the denominator in the integral in equation (19) in the paper is para[1] * (4*x - para[1]) # but to account for the denominator in (15), whose square is p4/x, we multiply by p4/x # (evaluated at x = para[1]) and this introduces a factor p4 while removing a factor x. IS := normal( (para[2]+1) / (subs(x=para[1],p4)*(4*x-para[1])) * diff(para[1],v)^2); # Size reduce: IS := subs(v = FindScaling({numer(IS), coeffs(collect(denom(IS),x),x)}, v), IS); mul(collect(i[1],x,factor)^i[2], i = sqrfree(IS,{x,v})[2]) end: # Input: a set S of polynomials in v. Output: a linear transformation v -> a*v + b # that is intended to reduce the size of the polynomials in S. FindScaling := proc(S, v) local V,m,Vm,n,Vmn,M,d,sb,sc,i; V := select(has,factors(convert(S,`*`))[2],v); m := min(seq(degree(i[1],v),i=V)); Vm := [seq(`if`(degree(i[1],v)=m,i,NULL),i=V)]; n := max(seq(i[2],i=Vm)); Vmn := sort([seq(`if`(i[2]=n,i,NULL),i=Vm)], length); M := Vmn[1]; d := degree(M[1],v); if nargs = 2 then sb := v - coeff(M[1], v, d-1)/lcoeff(M[1],v)/d; sc := procname(subs(v=sb, mul(i[1]^i[2], i = {op(V)} minus {M})), v, v); subs(v = sc*v, sb); else d := coeff(M[1],v,d-1)/lcoeff(M[1],v); if d = 0 then 1 else -d fi fi; end: N := 8; L[N]; IntegrandSquare( P4[N], x, N, v); # Next, lets look at larger N. A convenient way get a point on X1(N) is to use the # equation for X1(N) from Sutherland's website. If (x0, y0) is a solution of that # equation then the input polynomial for TorsionCondition is this: p4s := x^4 + c3*(x0^3*y0^2+x0^2*y0^2-2*x0^2*y0-2*x0*y0^2+3*x0*y0+y0^2-x0-2*y0+1)*x^3 + 1/4*(x0^6*y0^4+6*x0^5*y0^4-8*x0^5*y0^3-15*x0^4*y0^4+18*x0^4*y0^3+14*x0^3*y0^4-2*x0^4*y0^2-18*x0^3*y0^3-6*x0^2*y0^4 +8*x0^2*y0^3+4*x0^3*y0+3*x0^2*y0^2+2*x0*y0^3+y0^4-6*x0^2*y0-6*x0*y0^2-4*y0^3+x0^2+6*x0*y0+6*y0^2-2*x0-4*y0+1)*c3^2*x^2 + 1/2*y0^2*(x0*y0-1)*x0^3*(y0-1)*(x0-1)*(x0*y0-y0+1)*(x0^2*y0-x0*y0+y0-1)*c3^3*x; # This p4s is the polynomial you get by taking a point on X1(N) (encoded as a solution (x0,y0) of Sutherland's # equation for X1(N)), convert that point to an elliptic curve with a point of order N, and then move that point # to x = infinity. This doesn't necessarily mean that the torsion of P1-P2 is exactly N, because if you move a # point P of order N to P1, and you move -P to P2, then the divisor P1-P2 is only torsion of order n := N/igcd(N,2) # which explains the first line in the programs L_NP and IntegrandSquare. # Next, lets take some points (x0,y0) for various N's. We can download the equation F_N for X1(N) # from Sutherland's website, however, I constructed a fast algorithm that produces the same equation. ############# Intermezzo. Equation for X1(N) ############################ Fxy := proc(n::posint,x,y,s) options remember; local l; l := iquo(n,2); if n < 10 then 1 elif nargs=3 then procname(args, normal) elif irem(n,2)=0 then # Apply recursion relation for division polynomials P(n) # The omitted the factor P(l)/P(2) is handled in mul_diff. mul_diff(n, mul_list(Pxym(l+2), Pxym(l-1),2) , mul_list(Pxym(l-2), Pxym(l+1),2) ,x,y,s) else mul_diff(n, mul_list(Pxym(l+2), Pxym( l ),3) , mul_list(Pxym(l-1), Pxym(l+1),3) ,x,y,s) fi end: Pxy := proc(n::posint,x,y) local i; if nargs=4 then # Optional 4th argument gives product in list-form: [seq([nthpol(i[1],x,y),i[2]], i=Pxym(n))] else mul(i[1]^i[2], i=procname(n,x,y)) fi end: nthpol := proc(n::posint,x,y) if n > 7 then Fxy(n, x, y) else [x, y, x-1, y-1, x*y-1, x*y-y+1, x^2*y-x*y+y-1][n] fi end: Pxym := proc(n::posint) options remember; local i; [seq([i,predicted_mult(i,n)], i=1..7), seq(`if`(irem(n,i)=0,[i,1],NULL),i=10..n)] end: # For i in 1..7 this is the multiplicity of nthpol[i] # predicted_mult := proc(i::posint, n::posint) local v; v := [[-8/7, [2, 8/7, 11/7, 9/7, 9/7, 11/7, 8/7 ]], [-2/5, [1, 2/5, 3/5, 3/5, 2/5 ]], [ 2/5, [0, -2/5, -3/5, -3/5, -2/5 ]], [ 3/7, [0, -3/7, -5/7, -6/7, -6/7, -5/7, -3/7 ]], [-3/4, [1, 3/4 ]], [ 3/8, [0, -3/8, -1/2, -3/8 ]], [ 1/3, [0, -1/3, -1/3] ]][i]; v[1]*n^2 + v[2, 1 + irem(n,nops(v[2]))] end: mul_list := proc(f, g, m) local P, ml, p, i; # Computes f*g when f,g are products given in list-form. P := {seq(i[1],i=f)} union {seq(i[1],i=g)}; for p in P do ml[p] := 0 od; for p in f do ml[p[1]] := p[2] od; for p in g do ml[p[1]] := ml[p[1]] + m * p[2] od; [seq([p, ml[p]], p=P)] end: # Compute f - g and divide by predicted factors. mul_diff := proc(n, f, g, x, y, s) local i, P, c, d, f1, g1, ml, p; # Find the common factor c between f and g: P := {seq(i[1],i=f)} intersect {seq(i[1],i=g)}; for p in f do ml[p[1]] := p[2] od; for p in g do ml[p[1]] := min(ml[p[1]], p[2]) od; c := [seq(`if`(ml[p]=0, NULL, [p, ml[p]]), p=P)]; # Divide away the common factor c: f1 := mul_list(f, c, -1); g1 := mul_list(g, c, -1); d := mul(nthpol(i[1],x,y)^i[2], i=f1) - mul(nthpol(i[1],x,y)^i[2], i=g1); # The difference f - g # In the even case a factor was removed by omitting it in program "Fxy" # To compensate, add that to c if irem(n,2)=0 then c := mul_list(mul_list(Pxym(n/2),Pxym(2),-1), c, 1) fi; # Take expected multiplicity, here [1..-2] drops F_N. Then divide by c c := mul_list(Pxym(n)[1..-2], c, -1); # Now c is the predicted content of d, so divide it away: d := d / mul(nthpol(i[1],x,y)^i[2], i=c); # Apply the simplification specified by "s" which can be normal, evala, or a prime. if type(s,integer) then Normal(d) mod s else s(d) fi end: ############# Now back to differential equations ############################ ############# An example for N = 11 ############################ # To make examples, just take Fxy(N, x, y) for some N, and select a point, like this: N := 11; Fxy(N, x, y); # Quadratic in both variables. Lets take some random point, but don't take xs = 0 or 1 # because those correspond to cusps. Lets take this point: xs := -1; ys := RootOf( subs(x=xs, Fxy(N, x, y) )); p4 := collect(subs(x0=xs, y0 = ys, c3 = 1, p4s), x, evala); L11 := L_NP(p4, x, N); coeff(formal_sol(L11,T,x=0,terms=12)[1,2], T, 11) / binomial(2 * 11, 11); denom(%); ifactor(%); # There is still a denominator, to get rid of it, use the scaling factor c3. p4 := collect(subs(x0=xs, y0 = ys, c3 = 2^(-1) * 11^(-3), p4s), x, evala); L11 := L_NP(p4, x, N); coeff(formal_sol(L11,T,x=0,terms=12)[1,2], T, 11) / binomial(2 * 11, 11); denom(%); # This is the IntegrandSquare. IS := IntegrandSquare(p4, x, N, v); # The Picard-Fuchs equation of sqrt(IS) is L11. ############# An example for N = 37 ############################ zeta_n := proc(n) local x; RootOf(NumberTheory:-CyclotomicPolynomial(n,x),x) end: alpha_n := proc(n) local z,T; z := zeta_n(n); RootOf(sqrfree(evala(Norm(z + 1/z - T)))[2,1,1],T) end: N := 37; xs := RootOf(x^2-x-1); ys := solve(sort(evala(Factors(subs(x=xs, y^3+(-4*x-2)*y^2+(6*x+2)*y-x), alpha_n(7)))[2],length)[1,1],y); # Verify that (xs, ys) is a point on X1(37): ShouldBeZero := evala(subs(x = xs, y = ys, Fxy(N, x, y) )); p4 := collect(subs(x0=xs, y0 = ys, c3 = 1/N^2, p4s), x, evala); L37 := L_NP(p4, x, N); # Lets check a coefficient: coeff(formal_sol(L37,T,x=0,terms=12)[1,2], T, 11) / binomial(2 * 11, 11); denom(%); # OK. # To compute IntegrandSquare(p4, x, N, v) for N = 37 I may need to further speed # up the parametrization code because this involves parametrizing a degree 37 curve. ############# An example for N = 21 ############################ N := 21; xs := RootOf(x^3-3*x^2+3); ys := RootOf(subs(x=xs, y+x^2-2*x-1 )); # Verify that this is on X1(21): ShouldBeZero := evala(subs(x = xs, y = ys, Fxy(N, x, y) )); p4 := collect(subs(x0=xs, y0 = ys, c3 = N^(-2)/7, p4s), x, evala); L21 := L_NP(p4, x, N); coeff(formal_sol(L21,T,x=0,terms=12)[1,2], T, 11) / binomial(2 * 11, 11); denom(%); # IS := IntegrandSquare(p4, x, N, v); # # This last line takes about five minutes and produces a genus 10 hyperelliptic integral # whose Picard-Fuchs equation is L21.