Code:
Clear[a, b, c, d, m, r, s, f, x0, x1];
a = 1/2;
b = 3/2;
c = 1/4;
d = 3/4;
m = 1;
r = 3/2;
s = 1/4;
(* set up data to calculate F(x) *)
g[aval_, bval_, cval_, dval_, mval_, fval_] :=
mval + aval*fval^bval - cval*(1 - fval)^dval;
x0 = N[g[a, b, c, d, m, 0]];
x1 = N[g[a, b, c, d, m, 1]];
myPoints = Table[{g[a, b, c, d, m, f], f}, {f, 0, 1, 0.01}];
myInverse = Interpolation[myPoints];
(* now construct integrand with interpolated value of F(x) *)
fval[x_, y_] := myInverse[x]^r*myInverse[y]^s*myInverse[x]*
(1 - myInverse[y]);
(* integrate in x-y plane *)
NIntegrate[fval[x, y], {y, x0, x1}, {x, x0, y}]
(* integrate in u-v plane *)
NIntegrate[u^(r + 1)*v^s*(1 - v)*(a*b*u^(b - 1) +
c*d*(1 - u)^(d - 1))*(a*b*v^(b - 1) + c*d*(1 - v)^(d - 1)),
{v, 0, 1}, {u, 0, v}]
(* compute symbolic value *)
i1 = (a*b)^2/((r + b + 1)*(r + s + 2*b + 1)*(r + s + 2*b + 2));
i2 = ((a*b*c*d)/(r + b + 1))*(Gamma[r + s + b + 2]/
Gamma[r + s + b + d + 3])*Gamma[1 + d];
i3 = a*b*c*d*Gamma[d]*(Gamma[2 + r]/((b + s)*(b + s + 1)*
Gamma[r + 2 + d]) + Gamma[b + s + r + 3]/
(Gamma[b + s + r + d + 3]*(b + s + 1)) -
Gamma[b + s + r + 2]/((b + s)*Gamma[b + s + r + 2 + d]));
i4 = (c*d)^2*((Gamma[1 + d]*Gamma[2 + r]*Gamma[3 + r + s])/
(Gamma[3 + r]*Gamma[4 + d + r + s]))*
HypergeometricPFQ[{1 - d, 2 + r, 3 + r + s},
{3 + r, 4 + d + r + s}, 1];
N[i1 + i2 + i3 + i4]