See1 := false; # Set to true if you want to see the 4 animations with the moving branchpoints.

See2 := false; # Set to true of you want to see the paths from x = basePoint around each singularity.

ComputeMonodromy := true; # Set to true if you want to compute the four monodromy matrices.


L := Dx^4+4*(160*x^3+768*x^2+1170*x+729)/(4*x+9)/(16*x^2+56*x+81)/x*Dx^3+2*(12544*x^5+93120*x^4+271344*x^3
+395604*x^2+262440*x+59049)/x^2/(4*x+9)/(16*x^2+56*x+81)^2*Dx^2+72*(192*x^4+1280*x^3+2940*x^2+3282*x
+1215)/x^2/(4*x+9)/(16*x^2+56*x+81)^2*Dx+72*(8*x^3+50*x^2+80*x+93)/x^2/(4*x+9)/(16*x^2+56*x+81)^2;

L := primpart(L,Dx):   # L is the Picard-Fuchs operator for f = 1/sqrt(Onder) where:

Onder :=  subs( u = 1/2, v*(1-v)*(x*v*(-u*v+1)^2+(1-v)*(-u^2*v+1)*(u*v+1)^2) );

# The true singularities of L are S1, S2, S3, infinity, where S1, S2, S3 are the roots of discrim(Onder, v).

S1 := evalf( -7/4 + sqrt(-2) );
S2 := 0;
S3 := evalf( -7/4 - sqrt(-2) );

# For any x not in {S1, S2, S3, infinity},  the polynomial Onder has exactly 6 roots.
# We pick a basePoint x = 1 where the roots of Onder are:
#
#	r1 := -3-sqrt(17) ;
#	r2 := (3-sqrt(17))/2 ;
#	r3 := 0 ;
#	r4 := 1 ;
#	r5 := -3+sqrt(17) ;
#	r6 := (3+sqrt(17))/2 ;
#
# Then we loop, starting from x, around S1, S2, S3, and S4.  First some code:

basePoint := 1;

# Start at point p1. Then travel to p2 + 0.5*I in a straight line.
# Then travel to p2 - 0.5*I in a half-circle, on the left of p2.
# Then travel to p1 in a straight line.
Path_around_singularity := proc(p1, p2, t) local s,c;
	if type(t,name) then
		'procname'(args)
	elif t < 1/3 then
		s := t * 3;
		s * (p2 + 0.5*I) + (1-s) * p1
	elif t < 2/3 then
		s := t * 3 - 1;
		p2 + 0.5 * I * exp(Pi*I*s)
	else
		s := t * 3 - 2;
		s * p1 + (1-s) * (p2 - 0.5*I)
	fi
end:
SlowMiddle := proc(t)  # Makes the animation smoother in the middle:
	if t < 1/4 then
		3/2 * t  # Up to 3/8
	elif t < 3/4 then
		3/8 + (t-1/4)/2  # Up to 5/8
	else
		5/8 + 3/2 * (t-3/4) # Up to 1
	fi
end:
Path_around_infinity := proc(t)
	if type(t,name) then
		'procname'(args)
	else
		1/(  1/3 + 2/3 * exp(2*Pi*I*SlowMiddle(t)) )
	fi
end:

x_loop_1 := Path_around_singularity(basePoint, S1, t);
x_loop_2 := Path_around_singularity(basePoint, S2, t);
x_loop_3 := Path_around_singularity(basePoint, S3, t);
x_loop_4 := Path_around_infinity(t);

with(plots):
if See1 = true then
	hela := proc(x0) local i; global v,x,Onder;
		if indets(x0,name)<>{} then 'procname'(x0)
		else [seq([Re(i), Im(i)], i = fsolve(eval(Onder,x=x0), complex))]
		fi
	end:
	animate(pointplot, [hela(x_loop_1)], t=0..1);
	animate(pointplot, [hela(x_loop_2)], t=0..1);
	animate(pointplot, [hela(x_loop_3)], t=0..1);
	animate(pointplot, [hela(x_loop_4)], t=0..1);
fi;
if See2 = true then
	hola := proc(x0) local i; global v,x,Onder, S1, S2, S3, basePoint;
		if indets(x0,name)<>{} then 'procname'(x0)
		else [seq([Re(i), Im(i)], i = [x0, S1, S2, S3, basePoint])]
		fi
	end:
	animate(pointplot, [hola(x_loop_1)], t=0..1);
	animate(pointplot, [hola(x_loop_2)], t=0..1);
	animate(pointplot, [hola(x_loop_3)], t=0..1);
	animate(pointplot, [hola(x_loop_4)], t=0..1);
fi;
if ComputeMonodromy = true then

  Nsteps := 80;  Nterms := 80;  Digits := 60;

  # Want to compute   evalf(Int(1/sqrt(Onder), v=r1..r2))  but the numerical
  # integrator sometimes chokes on the denominator vanishing at r1 and r2.
  # So use a change-of-variables at the first and last 1/4 of the interval:
  #
  int_r1_r2 := proc(Onder, v, r1, r2) local I1,I2,I3,u;
	I1 := Int(2/sqrt(subs(v=u^2, quo(subs(v=v+r1,Onder),v,v) )), u=0..sqrt(r2-r1)/2);
	I2 := Int(1/sqrt(Onder), v = 3/4 * r1 + 1/4 * r2 .. 1/4 * r1 + 3/4 * r2);
	I3 := Int(2/sqrt(subs(v=-u^2, quo(subs(v=v+r2,-Onder),v,v) )), u=0..sqrt(r2-r1)/2);
	evalf(I1, 50) + evalf(I2,50) + evalf(I3,50)
  end;
  SolAtPoint := proc(L,x,Dx,T, xv, Init) local c,L0,A,isz,i;
	global Nterms;
	L0 := collect(subs(x = T + xv, L), Dx, expand);
	A := Init + add(c[i]*T^i, i=4..Nterms);
	isz := collect(add( coeff(L0,Dx,i)*diff(A,[T$i]), i=0..4 ), x);
	for i from 4 to Nterms do
		c[i] := solve(evalf(eval(coeff(isz,T,i-4))),c[i])
	od;
	eval(A)
  end:

  # Formal solutions at basepoint of the form T^i + O(T^4), i = 0..3
  FS := [seq(SolAtPoint(L,x,Dx,T, 1, T^i), i=0..3)];

  for i from 1 to 4 do
	xv := 1.0 + (i-2.5)*0.1 ;
	Onder :=  evalf(subs( u = 1/2, x=xv,  v*(1-v)*(x*v*(-u*v+1)^2+(1-v)*(-u^2*v+1)*(u*v+1)^2) ));
	r1,r2,r3,r4,r5,r6 := op(sort([fsolve(Onder)]));
	y1 := 2 * int_r1_r2(Onder, v, r1, r2);
	y2 := 2*I * int_r1_r2(-Onder, v, r2, r3);
	y3 := -2 * int_r1_r2(Onder, v, r3, r4);
	y4 := 2/I * int_r1_r2(-Onder, v, r4, r5);
	y5 := 2 * int_r1_r2(Onder, v, r5, r6);
	lprint(abs(y1 + y3 + y5));  # 0 up to working precision
	V[i] := [y1, y2, y4, y5];
	W[i] := subs(T = xv-1, FS);
  od:

  M := Matrix(4,4,[seq(W[i],i=1..4)])^(-1) . Matrix(4,4,[seq(V[i],i=1..4)]);

  # Now [y1 y2 y4 y5] = FS . M   and  FS = [y1 y2 y4 y5] . M^(-1)

  # Say that a loop sends FS to FS_loop = FS . A  for some matrix A.
  # Then [y1 y2 y4 y5] = FS . M goes to FS_loop . M
  #                                   =  FS . A . M = [y1 y2 y4 y5] . M^(-1) . A . M
  # so then the monodromy for basis y1 y2 y4 y5 is M^(-1) . A . M

  SolAtNextPoint := proc(L,x,Dx,T, xv, Delta_x, sol) local i;
	SolAtPoint(L,x,Dx,T,  xv+Delta_x, add(eval(diff(sol,[T$i]), T=Delta_x) * T^i / i!, i=0..3))
  end:
  SolAfterLoop := proc(L,x,Dx,T, path, t, sol) local S, i, old, new, Delta_x;
	global Nsteps;
	S := sol;
	for i from 1 to Nsteps do
		old := evalf(eval(path, t=(i-1)/Nsteps));
		new := evalf(eval(path, t=i/Nsteps));
		Delta_x := new-old;
		S := SolAtNextPoint(L,x,Dx,T, old, Delta_x, S);
	od;
	S
  end:
  for k from 1 to 4 do
	FS_loop := [seq(SolAfterLoop(L,x,Dx,T, x_loop_||k, t, i), i = FS)]:
	A := Matrix(4,4,[seq([seq(coeff(i,T,j), i=FS_loop)],j=0..3)]);
	Mon := M^(-1) . A^(-1) . M;   # See comment 1 below.
	E||k := map(round, 10^10 * Mon)/10^10;
	print("Matrix", E||k, "precision is:", evalf(max(map(abs,map(op,convert(Mon - E||k,listlist)))),4));
  od:
fi:

# Comment 1:
#
# In the line Mon := ...  we take A^(-1) instead of A, in order to ensure that
# the monodromy representation rho: pi_1 --> Gl4(C)
# is a homomorphism (instead of an anti-homomorphism).
# Suppose that A.1 is the matrix for loop.1, and A.2 is the matrix for loop.2.
# The notation   "loop.1 * loop.2"  usually denotes "First loop.1 and then loop.2"
# but the matrix A for that path will be  A.2 * A.1
# So if didn't invert A, then rho( loop_k ) = E||k = conjugate of A.k
# will produce an anti-homomorphism.   However, one expects rho to be a homomomorphism
# which we can accomplish by inverting A, as is done for example in equation (1.5)
# in  https://conservancy.umn.edu/handle/11299/1346