How can I speed up this Monte Carlo integration for the Stieltjes constant $\gamma_n$?

Question

Pre-requisite Definition

Stieltjes constants ($\gamma_n$) are the constants that occur in the Laurent series expansion of the Riemann zeta function. $${\displaystyle \zeta (s)={\frac {1}{s-1}}+\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{n!}}\gamma _{n}(s-1)^{n}.}$$

They have a finite integral representation as follows, $${\displaystyle \gamma _{n}={\frac {(-1)^{n}n!}{2\pi }}\int _{0}^{2\pi }e^{-nix}\zeta \left(e^{ix}+1\right)dx.}$$

Code

I have writen the following simple Monte Carlo algorithm for computing $\gamma_n$.

Note: I'm aware of MonteCarlo method in NIntegrate. The purpose of this code is more to understand/motivate the algorithm involved in crude MC Integration and I want to define it myself.

sint[n_, x_] := Exp[-n I x] Zeta[Exp[I x] + 1]
lowerlim = 0; upperlim = 2 Pi;
ParallelRepeatedStieltjesIntegral[n_, points_, repeat_] := 
((-1)^n n!)/(2 Pi) ParallelTable[(upperlim - lowerlim)/points Total[sint[n, RandomReal[{lowerlim, upperlim}, {points}]]], {repeat}]

My question is as follows, How can I speed up this code?

Extra Info

Time taken to execute this code on my PC

• ParallelRepeatedStieltjesIntegral[1, 10^3, 10^5] took 40 minutes!
• ParallelRepeatedStieltjesIntegral[2, 10^3, 10^5] took 32 minutes!

Also at the end of all this the values are quite poor (2 digit accuracy). I would prefer to run this code with larger points values.

I also want to repeat it, for the purpose of displaying the histogram (as this nicely demonstrates the law of large numbers).

PS: A not very relevant side question, how do you generate quasirandom numbers in Mathematica, instead of pseudorandom. (To demonstrate quasi-Monte Carlo integration in similar code as above.)

Answer 1

The slowest part of the code is clearly Zeta. You have to compute it only on a small parameter range and you also do not much very high accuracy (as Monte Carlo won't deliver much accuracy, too). So you can just sample Zeta over a fixed grid and use an interpolation function instead. Mathematica's interpolation functions are notoriously slow, so I implement my own compiled function here; it only uses piecewise-linear interpolation but nevertheless delivers 6 digits of precision.

sub = 2000;
tlist = Subdivide[0., 2. Pi, sub];
zeta = Zeta[1. + Exp[I tlist]];
cF = With[{zeta = zeta, tlist = tlist, L = sub/(2. Pi)},
   Compile[{{t, _Real}},
    Block[{i, \[Lambda]},
     i = Floor[Mod[t, 2. Pi]  L] + 1;
     \[Lambda] = (t - tlist[[i]]) L;
     (1. - \[Lambda]) zeta[[i]] + \[Lambda] zeta[[i + 1]]
     ],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];

I also made a couple of small modifictions in the remaining code:

sint[n_, x_] := Exp[-n I x] cF[x]
lowerlim = 0.; 
upperlim = 2. Pi;
ParallelRepeatedStieltjesIntegral[n_, points_, 
  repeat_] := (((-1)^n n!)/(2. Pi) (upperlim - lowerlim))  ParallelTable[
   Mean[sint[n, RandomReal[{lowerlim, upperlim}, {points}]]], {repeat},
   Method -> "CoarsestGrained"
   ]

Now this runs within a couple of seconds on my machine:

ParallelRepeatedStieltjesIntegral[2, 10^3, 10^5]; // AbsoluteTiming // First

2.83621

If you want to run the computations for various values of n anyways, you can save some time i) by using the same random points x for all n inequestion so that cF[x] has to be computed only once and ii) by computing only Exp[- I x] once and by computing Exp[-I n x] by multiplication. For example, the following does the computations for the first 100 coefficients at once, while it requires only about twelve times longer than for doing it for a single coefficient:

points = 10^3;
repeat = 10^5;
n1 = 1;
n2 = 100;

First@AbsoluteTiming[
  result = ParallelTable[
     Block[{x, Fx, z, factor, pow},
      x = RandomReal[{lowerlim, upperlim}, {points}];
      Fx = cF[x];
      z = Exp[(-I) x];
      factor = (upperlim - lowerlim)/(2. Pi)/points;
      pow = ConstantArray[1., Length[z]];
      Table[
       pow *= z;
       ((-1)^n n! factor) (pow.Fx)
       , {n, n1, n2}]
      ]
     , {repeat},
     Method -> "CoarsestGrained"
     ];
  ]

34.5558

For comparison, computing all the Exp[-I n x] directly lead to a runtime of 109.031 s.

No guarantees for correctness, though.