# My implementation

Are there some mistakes in my Mathematica code? Any help would be greatly appreciated.

### Attempt Version 1 (❌)

Clear["Global*"];

(*Function definitions*)
compMultiFoxH[params_, nsubdivisions_, boundaryTol_ : 0.0001] :=
Module[{boundaries, dim, signs, inputs, code, quad, volume},
boundaries = detBoundaries[params, boundaryTol];
Print["boundaries=", boundaries];
dim = Length[boundaries];
signs = Tuples[{-1, 1}, dim];
code = Tuples[Range[0, Floor[nsubdivisions/2] - 1], dim];
Do[Module[{points, integrandVals},
points = signs[[i]]  (code + 0.5)  (boundaries*2/nsubdivisions);
integrandVals = compMultiFoxHIntegrand[points, params];
quad += Total[integrandVals];], {i, 1, Length[signs]}];
volume = Times @@ (2*boundaries/nsubdivisions);

detBoundaries[params_, tol_] :=
Module[{boundaryRange, dims, boundaries, points, absIntegrand,
index}, boundaryRange = Range[0, 50, 0.05];
dims = Length[params[[1]]];
boundaries = ConstantArray[0, dims];
Do[points = ConstantArray[0, {Length[boundaryRange], dims}];
points[[All, dimL]] = boundaryRange;
absIntegrand = Abs[compMultiFoxHIntegrand[points, params]];
index =
Last[Flatten[
Position[absIntegrand, _?(# > tol*absIntegrand[[1]] &), {1},
boundaries[[dimL]] = boundaryRange[[index]];, {dimL, 1, dims}];
boundaries]

compMultiFoxHIntegrand[y_, params_] :=
Module[{z, mn, pq, c, d, a, b, m, n, p, q, npoints, dims, s, lower,
upper, mindist, sigs, num, cnorm, newdist, s1, prodGamNum,
prodGamDenom, zs}, {z, mn, pq, c, d, a, b} = params;
m = mn[[All, 1]];
n = mn[[All, 2]];
p = pq[[All, 1]];
q = pq[[All, 2]];
npoints = Length[y];
dims = Length[First[y]];
s = I  y;
lower = ConstantArray[0, dims];
upper = ConstantArray[0, dims];
Do[If[b[[dimL]] =!= {},
Module[{bj, Bj}, bj = b[[dimL, All, 1]][[;; m[[dimL + 1]]]];
Bj = b[[dimL, All, 2]][[;; m[[dimL + 1]]]];
lower[[dimL]] = -Min[bj/Bj];], lower[[dimL]] = -100];
If[a[[dimL]] =!= {},
Module[{aj, Aj}, aj = a[[dimL, All, 1]][[;; n[[dimL + 1]]]];
Aj = a[[dimL, All, 2]][[;; n[[dimL + 1]]]];
upper[[dimL]] = Min[(1 - aj)/Aj];], upper[[dimL]] = 1];, {dimL,
1, dims}];
mindist = Norm[upper - lower];
sigs = 0.5  (upper + lower);
Do[num = 1 - c[[j, 1]] - Total[c[[j, 2 ;;]]  lower];
cnorm = Norm[c[[j, 2 ;;]]];
newdist = Abs[num]/(cnorm + $MachineEpsilon); If[newdist < mindist, mindist = newdist; sigs = lower + 0.5 num c[[j, 2 ;;]]/(cnorm^2);], {j, 1, n[[1]]}]; s = s + sigs; s1 = Join[ConstantArray[1, {npoints, 1}], s, 2]; prodGamNum = 1 + 0 I; prodGamDenom = 1 + 0 I; Do[prodGamNum *= Gamma[(1 - s1 . c[[j]])];, {j, 1, n[[1]]}]; Do[prodGamDenom *= Gamma[(1 - s1 . d[[j]])];, {j, 1, q[[1]]}]; Do[prodGamDenom *= Gamma[(s1 . c[[j]])];, {j, n[[1]] + 1, p[[1]]}]; Do[Do[prodGamNum *= Gamma[(1 - a[[dimL, j, 1]] - a[[dimL, j, 2]] s[[All, dimL]])];, {j, 1, n[[dimL + 1]]}]; Do[prodGamNum *= Gamma[(b[[dimL, j, 1]] + b[[dimL, j, 2]] s[[All, dimL]])];, {j, 1, m[[dimL + 1]]}]; Do[prodGamDenom *= Gamma[(a[[dimL, j, 1]] + a[[dimL, j, 2]] s[[All, dimL]])];, {j, n[[dimL + 1]] + 1, p[[dimL + 1]]}]; Do[prodGamDenom *= Gamma[(1 - b[[dimL, j, 1]] - b[[dimL, j, 2]] s[[All, dimL]])];, {j, m[[dimL + 1]] + 1, q[[dimL + 1]]}];, {dimL, 1, dims}]; zs = z^-s; (prodGamNum/prodGamDenom) Product[zs, {2}]/(2 Pi)^dims // N] (*Example usage*) params1 = {{16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499}, {{0, 0}, {2, 1}, {2, 1}, {2, 1}, {2, 1}}, {{0, 1}, {1, 2}, {1, 2}, {1, 2}, {1, 2}}, {}, {0, 1, 1, 1, 1}, {{{1, 2}}, {{1, 2}}, {{1, 2}}, {{1, 2}}}, {{{1, 0.6666666666666666}, {3.5, 0.5}}, {{1, 0.6666666666666666}, {3.5, 0.5}}, {{1, 0.6666666666666666}, {3.5, 0.5}}, {{1, 0.6666666666666666}, {3.5, 0.5}}}}; result = compMultiFoxH[params1, 20]; Print[result];  ### Attempt Version 2 (❌) (*Define functions for gamma product \ computation*)ClearAll["Global*"] gammaProdNum[s1_, c_, n_] := Times @@ (Gamma[1 - c[[#]] . s1] & /@ Range[n[[1]]]) gammaProdDenom[s1_, c_, p_, q_] := Times @@ (Gamma[1 - c[[#]] . s1] & /@ Range[q[[1]]])* Times @@ (Gamma[c[[#]] . s1] & /@ Range[n[[1]] + 1, p[[1]]]) (*Compute boundaries*) detBoundaries[params_, tol_] := Module[{boundaryRange, dims, boundaries, points, absIntegrand, index}, boundaryRange = Range[0, 50, 0.05]; dims = Length[params[[1]]]; boundaries = ConstantArray[0, dims]; Do[points = ConstantArray[0, {Length[boundaryRange], dims}]; points[[All, dimL]] = boundaryRange; absIntegrand = Abs[compMultiFoxHIntegrand[points, params]]; index = Max[FirstPosition[ UnitStep[absIntegrand - tol*First[absIntegrand]], 1][[1]] - 1]; boundaries[[dimL]] = boundaryRange[[index]];, {dimL, dims}]; boundaries] (*Compute complex integrand of the multivariate Fox-H function*) compMultiFoxHIntegrand[y_, params_] := Module[{z, mn, pq, c, d, a, b, m, n, p, q, s, lower, upper, dims, s1, prodGamNum, prodGamDenom, zs, result}, {z, mn, pq, c, d, a, b} = params; {m, n} = Transpose[mn]; {p, q} = Transpose[pq]; {npoints, dims} = Dimensions[y]; s = I*y; lower = ConstantArray[0, dims]; upper = ConstantArray[0, dims]; Do[If[Length[b[[dimL]]] > 0, lower[[dimL]] = -Min[b[[dimL, All, 1]]/b[[dimL, All, 2]]], lower[[dimL]] = -100]; If[Length[a[[dimL]]] > 0, upper[[dimL]] = Min[(1 - a[[dimL, All, 1]])/a[[dimL, All, 2]]], upper[[dimL]] = 1];, {dimL, dims}]; mindist = Norm[upper - lower]; sigs = 0.5*(upper + lower); Do[num = 1 - c[[j, 1]] - Total[c[[j, 2 ;;]]*lower]; cnorm = Norm[c[[j, 2 ;;]]]; newdist = Abs[num]/(cnorm +$MachineEpsilon);
If[newdist < mindist, mindist = newdist;
sigs = lower + 0.5*num*c[[j, 2 ;;]]/(cnorm*cnorm)];, {j,
n[[1]]}];
s += sigs;
s1 = Transpose[Prepend[Transpose[s], ConstantArray[1, npoints]]];
prodGamNum = gammaProdNum[s1, c, n];
prodGamDenom = gammaProdDenom[s1, c, p, q];
Do[Do[prodGamNum *=
Gamma[1 - a[[dimL, j, 1]] -
a[[dimL, j, 2]]*s[[All, dimL]]];, {j, n[[dimL + 1]]}];
Do[prodGamNum *=
Gamma[b[[dimL, j, 1]] + b[[dimL, j, 2]]*s[[All, dimL]]];, {j,
m[[dimL + 1]]}];
Do[prodGamDenom *=
Gamma[a[[dimL, j, 1]] + a[[dimL, j, 2]]*s[[All, dimL]]];, {j,
n[[dimL + 1]] + 1, p[[dimL + 1]]}];
Do[prodGamDenom *=
Gamma[1 - b[[dimL, j, 1]] -
b[[dimL, j, 2]]*s[[All, dimL]]];, {j, m[[dimL + 1]] + 1,
q[[dimL + 1]]}];, {dimL, dims}];
zs = z^-s;
result = (prodGamNum/prodGamDenom)*Apply[Times, zs, {1}]/(2*Pi)^dims;
result]

(*Compute multivariate Fox-H function*)
compMultiFoxH[params_, nsubdivisions_, boundaryTol_ : 0.0001] :=
Module[{boundaries, dim, signs, code, quad, points, volume, result},
boundaries = detBoundaries[params, boundaryTol];
dim = Length[boundaries];
signs = Tuples[{1, -1}, dim];
code = Tuples[Range[0, Floor[nsubdivisions/2] - 1], dim];
Do[points =
DiagonalMatrix[sign]*((code + 0.5)*boundaries*2/nsubdivisions);
signs}];
volume = Apply[Times, 2*boundaries/nsubdivisions];
result]

(*Example usage*)
params1 = {{16.2982237081499, 16.2982237081499, 16.2982237081499,
16.2982237081499}, {{0, 0}, {2, 1}, {2, 1}, {2, 1}, {2, 1}}, {{0,
1}, {1, 2}, {1, 2}, {1, 2}, {1, 2}}, {}, {{0, 1, 1, 1,
1}}, {{{1, 2}}, {{1, 2}}, {{1, 2}}, {{1, 2}}}, {{{1,
0.6666666666666666}, {3.5, 0.5}}, {{1,
0.6666666666666666}, {3.5, 0.5}}, {{1,
0.6666666666666666}, {3.5, 0.5}}, {{1,
0.6666666666666666}, {3.5, 0.5}}}};
Print[compMultiFoxH[params1, 20]]

• There are several typos in your code. Have you translated your code from Python and Matlab to Mathematica? Commented Jul 21 at 7:27
• Version 2 has same typos as version 1. :) Commented Jul 21 at 9:46
• What is the definition of Multivariate FoxH? A example paper ? Commented Jul 21 at 11:36
• @MariuszIwaniuk The definition can be found at Eqn(1.1) of AN INTEGRAL INVOLVING GENERALIZED FUNCTION OF TWO VARIABLES and Eqn（A.1) of The H-Function Theory and Applications i.sstatic.net/D0lyUU4E.png Commented Jul 21 at 12:37
• thanks a lot :) Commented Jul 21 at 13:01

Now it works. The result is consistent with Python/MATLAB.

Any help for accelerating the Mathematica code would be greatly appreciated.

(*Define functions for gamma product computation*)ClearAll["Global*"]

(*Compute boundaries*)
detBoundaries[params_, tol_] :=
Module[{boundaryRange, dims, boundaries, points, absIntegrand,
index}, boundaryRange = Range[0, 50, 0.05];
dims = Length[params[[1]]];
boundaries = ConstantArray[0, dims];
Do[points = ConstantArray[0, {Length[boundaryRange], dims}];
points[[All, dimL]] = boundaryRange;
absIntegrand = Abs[compMultiFoxHIntegrand[points, params]];
index =
Max[Select[Range[Length[absIntegrand]],
absIntegrand[[#]] > tol*absIntegrand[[1]] &]];
boundaries[[dimL]] = boundaryRange[[index]];, {dimL, dims}];
boundaries]

(*Compute complex integrand of the multivariate Fox-H function*)
compMultiFoxHIntegrand[y_, params_] :=
Module[{z, mn, pq, c, d, a, b, m, n, p, q, npoints, dims, s, lower,
upper, mindist, sigs, num, cnorm, newdist, s1, prodGamNum,
prodGamDenom, zs}, {z, mn, pq, c, d, a, b} = params;
m = mn[[All, 1]];
n = mn[[All, 2]];
p = pq[[All, 1]];
q = pq[[All, 2]];
{npoints, dims} = Dimensions[y];
s = I   y;
lower = ConstantArray[0, dims];
upper = ConstantArray[0, dims];
Do[If[b[[dimL]] =!= {},
Module[{bj, Bj}, bj = b[[dimL, All, 1]][[;; m[[dimL + 1]]]];
Bj = b[[dimL, All, 2]][[;; m[[dimL + 1]]]];
lower[[dimL]] = -Min[bj/Bj];], lower[[dimL]] = -100];
If[a[[dimL]] =!= {},
Module[{aj, Aj}, aj = a[[dimL, All, 1]][[;; n[[dimL + 1]]]];
Aj = a[[dimL, All, 2]][[;; n[[dimL + 1]]]];
upper[[dimL]] = Min[(1 - aj)/Aj];], upper[[dimL]] = 1];, {dimL,
1, dims}];
mindist = Norm[upper - lower];
sigs = 0.5   (upper + lower);
Do[num = 1 - c[[j, 1]] - Total[c[[j, 2 ;;]]   lower];
cnorm = Norm[c[[j, 2 ;;]]];
newdist = Abs[num]/(cnorm + \$MachineEpsilon);
If[newdist < mindist, mindist = newdist;
sigs = lower + 0.5   num   c[[j, 2 ;;]]/(cnorm^2);], {j, 1,
n[[1]]}];
s = Table[s + sigEle, {sigEle, sigs}];
s = Flatten[s, 1];
s = ArrayReshape[s, {Length@s/dims, dims}];
s1 = Table[{1.00}~Join~sEle, {sEle, s}];
prodGamNum = 1 + 0   I;
prodGamDenom = 1 + 0   I;
Do[prodGamNum *= Gamma[(1 - s1 . c[[j]])];, {j, 1, n[[1]]}];
Do[prodGamDenom *= Gamma[(1 - s1 . d[[j]])];, {j, 1, q[[1]]}];
Do[prodGamDenom *= Gamma[(s1 . c[[j]])];, {j, n[[1]] + 1, p[[1]]}];
Do[Do[prodGamNum *=
Gamma[(1 - a[[dimL, j, 1]] -
a[[dimL, j, 2]]   s[[All, dimL]])];, {j, 1, n[[dimL + 1]]}];
Do[prodGamNum *=
Gamma[(b[[dimL, j, 1]] +
b[[dimL, j, 2]]   s[[All, dimL]])];, {j, 1, m[[dimL + 1]]}];
Do[prodGamDenom *=
Gamma[(a[[dimL, j, 1]] +
a[[dimL, j, 2]]   s[[All, dimL]])];, {j, n[[dimL + 1]] + 1,
p[[dimL + 1]]}];
Do[prodGamDenom *=
Gamma[(1 - b[[dimL, j, 1]] -
b[[dimL, j, 2]]   s[[All, dimL]])];, {j, m[[dimL + 1]] + 1,
q[[dimL + 1]]}];, {dimL, 1, dims}];
zs = Table[z^-sEle, {sEle, s}];
result = (prodGamNum/prodGamDenom)  *
Table[Times @@ zsEle, {zsEle, zs}]/(2   Pi)^dims // N;
result
]

(*Compute multivariate Fox-H function*)
compMultiFoxH[params_, nsubdivisions_, boundaryTol_ : 0.0001] :=
Module[{boundaries, dim, signs, code, quad, points, volume, result},
boundaries = detBoundaries[params, boundaryTol];
dim = Length[boundaries];
signs = Tuples[{1, -1}, dim];
code = Tuples[Range[0, Floor[nsubdivisions/2] - 1], dim];
Do[points =
Table[sign*((codeEle + 0.5)*boundaries*2/nsubdivisions), {codeEle,
code}];
signs}];
volume = Apply[Times, 2*boundaries/nsubdivisions];
result]

(*Example usage*)
params1 = {
{16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499}, (*z*)
{{0, 0}, {2, 1}, {2, 1}, {2, 1},{2, 1}}, (*mn*)
{{0, 1}, {1, 2}, {1, 2}, {1, 2}, {1, 2}}, (*pq*)
{}, (*c*)
{{0, 1, 1, 1, 1}}, (*d*)
{{{1, 2}}, {{1, 2}}, {{1, 2}}, {{1, 2}}},  (*a*)
{{{1, 0.6666666666666666}, {3.5, 0.5}},
{{1, 0.6666666666666666}, {3.5, 0.5}},
{{1, 0.6666666666666666}, {3.5, 0.5}},
{{1, 0.6666666666666666}, {3.5, 0.5}}} (*b*)
};
compMultiFoxH[params1, 20] // NumberForm[#, 20] & // Print
`
• This is good version (+1). Commented Jul 21 at 14:18