# Compute the first 100 terms of the functions U(t;x) and U(t;x,b).
Utx  := add(binomial(2*n,n) * add(binomial(n,k) * binomial(n+2*k, 2*k) * binomial(2*k,k) * x^(n-k),   k=0..n)*t^n, n=0..100):
Utxb := t * diff(Utx, t) + b * Utx:

# Compute the first 100 terms of the functions V(t;x) and V(t;x,b).
Vtx  := add(binomial(2*n,n) * add(binomial(2*k,k)^2 * binomial(2*n-2*k, n-k)^2 / binomial(n,k) * x^k, k=0..n)*t^n, n=0..100):
Vtxb := t * diff(Vtx, t) + b * Vtx:

# Check equations (4.3)-(4.9)
Digits := 50;
evalf( subs(t=1/2160,  x=-324,  b=103/357, Utxb) - 30/(119*Pi));
evalf( subs(t=1/3645,  x=486,   b=0,       Utxb) - 10/(3*Pi));
evalf( subs(t=-1/1728, x=-324,  b=1/6,     Utxb) - 12*sqrt(375+120*sqrt(10))/(75*Pi));
evalf( subs(t=-1/160,  x=-20,   b=1/4,     Utxb) - sqrt(30)/(20*Pi) * (5+(145+30*sqrt(6))^(1/3)) / (145+30*sqrt(6))^(1/6));
evalf( subs(t=1/27648, x=-2160, b=289/1290,Utxb) - 16*sqrt(15)/(215*Pi));
evalf( subs(t=1/276480,x=12096, b=49/804,  Utxb) - 10*sqrt(15)/(67*Pi));
evalf( subs(t=2/135,   x=-27/8, b=5/24,    Utxb) - (5*sqrt(6)+4*sqrt(15))/(16*Pi));

# Check equations (4.10)-(4.11)
evalf( subs(t = 1/576, x=5,     b=5/28,    Vtxb) - 9*(2+sqrt(2))/(28*Pi));
evalf( subs(t = 1/576, x=-25/16,b=31/182,  Vtxb) - 189/(364*Pi));



# The following is a formula for V(t; -25/16) in terms of hyper-elliptic integrals:
#
#	(*)  V(t; -25/16) = ExtraFactor * Int(f1, v=0..TOP) * Int(f2, v=0..TOP)
#
# where
#
#	f1 = Boven_p / sqrt(Onder)
#	f2 = Boven_m / sqrt(Onder)

Boven_p :=  5/4*(2-128*t)^(1/2)*v + v + 1;
Boven_m := -5/4*(2-128*t)^(1/2)*v + v + 1;
Onder := (1600*t*v-9*v+16) * v * ( (1600*t*v-9*v+16) * (9-16*v-1600*t) * (1-v)^2
      + (   (12288*t^2-3)*(-(100*t+1)*(1600*t-9))^(1/2)-1433600*t^3-123136*t^2-350*t+9  ) * v * (1+v)^2 / (2*t) );
TOP := 1/1024*(1-((100*t+1)*(1-64*t))^(1/2))*(1600*t-9-3*((9-1600*t)*(1-64*t))^(1/2))/t/(9/16-100*t);
ExtraFactor := -8*(1-64*t)^(1/2)*((192*t-3)*(100*t+1)^(1/2)*(9-1600*t)^(1/2)-102400*t^2-1024*t+9)/Pi^2/t;


# One can prove (*) for sufficiently small t by computing a differential equation for both sides (they
# satisfy the same differential equation) and then comparing initial conditions.
# In the remainder of this file, we'll do a numerical check for equation (*).


############  Numerical verification of equation (*)  #################


# Integrate subs(t = tv, Boven/sqrt(Onder)) from v = 0 to v = TOP
#
NumericInt := proc(Boven, tv) global t, Boven_p, Boven_m, Onder, v, TOP;
        L := evalf( subs(t=tv, [Boven, Onder, TOP] ), 100);
	MyInt(L[1], L[2], v, 0, L[3]);
end:

# Assumption:  f = P/sqrt(Q) with a, b roots of Q.
# Goal: compute int(f, x=a..b) in a numerically stable way via change of variables.
MyInt := proc(P, Q, x, a, b)
	# x-a = t^2   # x = t^2 + a   # dx = 2*t*dt  # Q1 = Q/(x-a)
	# P/sqrt(Q) * dx = P/sqrt(Q1) * 1/t * dx = 2 * P/sqrt(Q1) * dt

	sq := sqrt( (b-a)/2 );
	I1 := Int( subs(x = t^2 + a,  2*P/sqrt(quo(Q,x-a,x))), t=0..sq) ;

	# b-x = t^2  # x = b-t^2
	I2 := Int( subs(x = b-t^2,  2*P/sqrt(quo(Q,b-x,x))), t=0..sq) ;
	evalf( I1 + I2 );
end:


# Compute:  ExtraFactor * Int(Boven_p / sqrt(Onder), v=0..TOP) * Int(Boven_m / sqrt(Onder), v=0..TOP)
# at the point t = tv.
#
Total := proc(tv) 
	global Boven_p, Boven_m, t, ExtraFactor;
	evalf(subs(t=tv, ExtraFactor)) * NumericInt(Boven_p, tv) * NumericInt(Boven_m, tv)
end:

Digits := 30;

Total(10^(-4));
Total(10^(-3));
Total( 1/576 );

# Check if the formula matches V(t; -25/16)

evalf( subs( t = 10^(-4), x=-25/16, Vtx) );
evalf( subs( t = 10^(-3), x=-25/16, Vtx) );
evalf( subs( t = 1/576,   x=-25/16, Vtx) );