1)  Derive a holonomic system for the internal triple integral of the expection of an Euler characteristic number

<< RISC`HolonomicFunctions` (* the pacakge is available in http://www.risc.jku.at/research/combinat/software/ergosum/installation.html#download *)

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

stheta=(2 s)/(1+s^2); (* make a universal transformation for Sin[theta] *)

ctheta=(1-s^2)/(1+s^2); (* make a universal transformation for Cos[theta] *)

sphi=(2 t)/(1+t^2); (* make a universal transformation for Sin[phi] *)

cphi=(1-t^2)/(1+t^2);(* make a universal transformation for Cos[phi] *)

R=s1 (b sphi stheta+d cphi ctheta-m11)^2+s2 (-b sphi ctheta+d cphi stheta-m21)^2+s2 (b cphi ctheta+d sphi stheta-m22)^2+s1 (-b cphi stheta+d sphi ctheta)^2; (* input R in the integrand *)

f=2/(1+t^2 )2/(1+s^2) (d^2-b^2) (s1 s2)/(2 Pi)^2 Exp[-1/2 R]; (*input the integrand of Example 1 *)

f1=Simplify[f/.{s1->2,s2->1,m11->1,m21->-1,m22->1}] ;(* specify the parameters as that in Example 1*)

ann=Annihilator[f1,{Der[d],Der[b],Der[s],Der[t]}];(* compute the differential equations ann of f1 *)

Timing[ann1=CreativeTelescoping[ann,Der[t]][[1]];] (* compute the differential equations ann1 of the inner single integral *)

supp1 = {1, Der[b], Der[d], Der[b]^2, Der[b] Der[d], Der[d]^2, Der[d]^3} (* a support set for a differential equation of the inner double integral *)

Timing[op1=FindCreativeTelescoping[ann1,Der[s], Support -> supp1][[1,1]];] (* compute a differential equation op1 of the inner double integral with support in supp1 *)

supp2 = {1, Der[b], Der[d], Der[b]^2, Der[b] Der[d], Der[d]^2, Der[b]^3, Der[b]^2 Der[d], Der[b] Der[d]^2, Der[d]^3};(* another support set for another differential equation of the inner double integral *)

Timing[op2=FindCreativeTelescoping[ann1,Der[s], Support -> supp2][[1,1]];](* compute another differential equation op2 of the inner double integral with support in supp2 *)

Timing[ann2 = OreGroebnerBasis[{op1, op2}];] (* compute a holonomic system ann2 of the inner double integral *)

UnderTheStaircase[ann2] (* ann2 has holonomic rank 6 *)

supp3 = Table[ Der[d]^i,{i, 0, 10}] (* a support set for a differential equation of the inner triple integral *)

Timing[ann3=FindCreativeTelescoping[ann2,Der[b], Support -> supp3][[1,1]];] (* compute a holonomic system ann3 of the inner triple integral with support in supp3 *)

2)  Use HGM and NIntegrate to compute the 4-fold  integral of the expection of an Euler characteristic number

LeadingCoefficient[ann3] // Factor (* this is the leading coefficient of ann3, it has factor d^3 and will cause problem when we do numerical solving of the differential equation ann3 around point d = 0 *)

ann4 = DFinitePlus[{ann3}, ToOrePolynomial[{Sum[RandomInteger[{1, 1000}] Der[d]^i, {i, 0, 3}]}, OreAlgebra[Der[d]]]][[1]]; (* Using desingularization technique to remove the factor d^3 in ann3 *)

LeadingCoefficient[ann4] // Factor (* we can reduce the multiplicity of d in ann3, but d is non-removable *)

uts3 = UnderTheStaircase[{ann3}] (* the parametric derivatives under the stair of ann3 *)

Timing[inits = NIntegrate[Together[ApplyOreOperator[uts3,1/2 f1] /. {d -> 1/2}], {b, -Infinity, Infinity}, {s, -Infinity, Infinity}, {t, -Infinity, Infinity},MaxRecursion->20,Method->{GlobalAdaptive,MaxErrorIncreases->10000}];] (* compute the initial value for the Pfaffian system of ann3 at d = 1/2. Since ann3 has a singularity at d = 0, it is not suitable to choose the initial point at d = 0 *)

eqnsd = ApplyOreOperator[ann3, y[d]] == 0; (* set up the differential equation for the inner triple integral *)

inits1 = Table[D[y[d], {d, i}] == inits[[i +1 ]], {i, 0, 9}] /. {d -> 1/2} (* set up the initial values for eqnsd at d = 1/2 *)

s1 = NDSolve[Join[{eqnsd},inits1], y, {d, 1/2, 30}]; (* numerical solving of the diffential equation eqnsd from 1/2 to 30 *)

Plot[Evaluate[y[d]/.s1],{d,1/2,30},PlotRange->All] (* y[d] is very small when d ≥ 5*)

HGM1 = Table[NIntegrate[y[d] /. s1, {d, i, 30}][[1]], {i, 1, 6}] (* Using numerical integration to compute the integral in Example1 for x from 1 to 6. The results are as accurate as that in Prof. Kuriki's mathematica notebook*)

3)  Use HGM directly to compute the 4-fold  integral  of the expection of an Euler characteristic number

ann4 = OreTimes[ann3, Der[d]] /. {d -> x} // Factor;(* ann4 is the annihilator for the 4-fold integral F[x] in Example1 *)

uts = UnderTheStaircase[{ann4}] (* the parametric derivatives of ann4 *)

Timing[init1 = NIntegrate[1/2 f1 , {d , 1/2, Infinity},{b, -Infinity, Infinity}, {s, -Infinity, Infinity}, {t, -Infinity, Infinity},MaxRecursion->20,Method->{GlobalAdaptive,MaxErrorIncreases->10000}];] (* compute the inital value F[1/2] *)

inits = Prepend[-inits, init1] (* the initial points F^(k)[1/2], k = 0, .., 10 *)

eqnsd = ApplyOreOperator[ann4, y[x]] == 0; (* set up the differential equation for the 4-fold integral in Example1 *)

inits1 = Table[D[y[x], {x, i}] == inits[[i +1 ]], {i, 0, 10}] /. {x -> 1/2} (* set up the initial values for eqnsd at d = 1/2 *)

s2 = NDSolve[Join[{eqnsd},inits1], y, {x, 1/2, 30}] (* numerical solving of the diffential equation eqnsd from 1/2 to 30 *)

Plot[Evaluate[y[x]/.s2],{x,1/2,30},PlotRange->All] (* y[x] is very small when x ≥ 5*)

HGM2 = Table[(y[x] /. s2)[[1]], {x, 1, 6}](* The results are as accurate as that in Prof. Kuriki's mathematica notebook*)

4) Monte Carlo study for the Expectation of Euler charateristic number in Example 1

EulerChar[A_,x_]:=Module[{b,d,a,c},{b,d}=SingularValueList[A];(*b and d are singular value of A*)If[TrueQ[b≥x],a=Sign[b^2-d^2],a=0];
If[TrueQ[d≥x],c=Sign[d^2-b^2],c=0];
Return[a+c];]; (* compute the Euler charateristic number of a random matrix A w.r.t. x*)

Clear[s1, s2];

sub={s1->2,s2->1,m11->1,m21->-1,m22->1}

S=DiagonalMatrix[{1/Sqrt[s1],1/Sqrt[s2]}]/.sub

M={{m11,0},{m21,m22}}/.sub

n =10000000 (* the number of iterations *)

ra=Table[S.{RandomVariate[NormalDistribution[],2,WorkingPrecision->10],RandomVariate[NormalDistribution[],2,WorkingPrecision->10]}+M,{i,1,n}]; (* generate n random matrix A*)

Timing[mc = Table[ N[Sum[EulerChar[ra[[i]],x],{i,1,n}]/n], {x, 1, 6}];] (*the results of expectation of Euler characteristic number by Monte Carlo study *)

mc (* The results are as accurate as that in Prof. Kuriki's mathematica notebook *)