H := v * (1-v) * ( (1-v)*(1-u^2*v)*(1+u*v)^2 + x*v*(1-u*v)^2 );
# Move one branchpoint to infinity, one branchpoint to 0,  then
# the curve becomes Y^2 = v * Pol4 where Pol4 equals:
Pol4 := numer(factor(subs(v = 1/v, v = v+1, H)))/v;

# Plug in random values for u and x
Pol4 := subs(u = 87, x = 53, Pol4);


Humbert8 := 4*a1*a2*a3*a4*( (a1+a3)*(a2+a4) - 2*a1*a3 - 2*a2*a4)^2
	= (a2-a4)^2 * (a1-a3)^2 * (a1*a3+a2*a4)^2;

# Compute the 4 roots.  Note that if a1 is a root, then Q(a1) contains
# two roots, and we have to number those as a1, a3 in order for
# the Humbert8 equation to hold:

a1 := RootOf(Pol4,v);
a3 := evala(Roots( evala(Pol4/(v-a1)), v))[1,1];
a2 := RootOf(evala(Pol4/(v-a1)/(v-a3)),v);
a4 := RootOf(evala(Pol4/(v-a1)/(v-a2)/(v-a3)),v);

ShouldBeZero := evala( lhs(Humbert8) - rhs(Humbert8) );
# So our hyperelliptic curve, at least, for randomly selected values of u and x,
# satisfies the Humbert8 equation.

##########  Computing the Picard-Fuchs equation for I_{+}(u,x) in Theorem 1  ##################

Boven := 1 - u*v + v * 2 * t * u;
Onder := H;
u := 1/(1-2*t^2);
f := Boven / Onder^(1/2):  # This is I_{+}(u,x)  written in terms of t.

# With these next two lines, this file computes L_plus_x,  and without them,
# the file computes L_{+}^PF = the Picard-Fuchs equation for I_{+}(u,x)

	f := subs(x = -4*u^2*x, f); 
	Onder := subs(x = -4*u^2*x, Onder);

# We now want a differential equation for f.

# We integrate over v;  x is the parameter and t is a constant throughout this file.

#########################################################
# Computation of differential equation for the periods. #
#########################################################


N := 2;  # Order of the operator we're looking for
L := add( a[i]*Dx^i, i=0..N);

Lf := add( a[i] * diff(f, [x$i]), i=0..N):

# If L annihilates the period,  then integrating Lf over a loop in the v-plane should be 0.
# That means that Lf = diff(g, v)  for some element v of the function field.

# In order to find g,  we need to know its denominator, a degree-bound d for the numerator
# as well as all coefficients c[0] .. c[d] in the numerator.

# In the entire file, t is a constant, but from this point on, x is also viewed as a
# constant because we are differentiating with respect to v.
# So it is OK to get the denominator correct up to some factor in Q(x,t) since those
# are constant under d/dv.  But in that case, it's easy to see what denom(g) must be,
# just look at Lf.  Differentiating turns a denominator Onder^(..) to Onder^(..+1).
# Thus:

denom_g := Onder^(1/2 + N-1): # diff(g, v) below will make denom = Onder^(1/2 + N) = denom(Lf).

d := 9;
numer_g := add( c[i]*v^i, i = 0..d):

g := numer_g / denom_g :

# Note: increasing d is OK, decreasing it is not, because the pole order of diff(g, v)
# at v = infinity must be at least as high as that of Lf, so comparing them gives
# a lower bound for d.  Any d at least that bound is OK.

# If Lf = diff(g,v) then L annihilates the periods.

IsZero := numer(normal( Lf - diff(g,v) )):
IsZero := {coeffs(IsZero, v)}:
Vars := {seq(a[i], i=0..N), seq(c[i], i=0..d)}:

nops(IsZero), nops(Vars);

# Note: we have more equations than unknowns, so the only way there'll be a solution is
# if there is something very special about the function f = Boven / Onder^(1/2).

sol := solve(IsZero, Vars):

L := eval(L, sol);

L := collect(L/lcoeff(L,Dx), Dx, factor);
smp := proc(a) global x,t;
	if type(a,`*`) then map(procname,a) else sort(collect(a,x,factor),x) fi
end:

L := collect(collect(L,Dx,factor),Dx,smp);

# Smaller:

L := sort(L - coeff(L,Dx,1)*Dx + map(smp,convert(coeff(L,Dx,1),parfrac,x))*Dx,Dx);
lprint(%);