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.}$$


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.)

  • 1
    $\begingroup$ I would not use Monte-Carlo for this integral but rather FFT. $\endgroup$ Jun 25, 2021 at 13:04
  • 1
    $\begingroup$ Hi Henrik, my goal really is just a showcase of MC integration nothing else. Mathematica even has a built in function as StieltjesGamma[n]. FFT is without a doubt practically better, but that's not the reason I'm doing it. $\endgroup$
    – Boris
    Jun 25, 2021 at 13:13
  • 2
    $\begingroup$ I would not use Parallel* functions but rather use packed arrays and vectorized operations (which are automatically and very efficiently parallelized). Search for “[performance-tuning] packed” and sort by votes. The top two hits make good introductory reading. (Sorry on an iOS device w/o Mathematica at present and cannot show you. Several others on site know how, though.) $\endgroup$
    – Michael E2
    Jun 25, 2021 at 14:22
  • 1
    $\begingroup$ "Monte Carlo is an extremely bad method; it should only be used when all alternative methods are worse." - Alan Sokal. With that said, 1D integrals are usually a poor choice for demonstrating (quasi-)Monte Carlo; you might prefer using a 2D example instead, especially one where the region of integration is not necessarily amenable to a conversion to an iterated integral. $\endgroup$ Jun 26, 2021 at 15:16
  • 2
    $\begingroup$ With that said, the built-in "MKL" generator implements the Sobol' and Niederreiter sequences, which may be of interest. $\endgroup$ Jun 26, 2021 at 15:19

1 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


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;

  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]];
       pow *= z;
       ((-1)^n n! factor) (pow.Fx)
       , {n, n1, n2}]
     , {repeat},
     Method -> "CoarsestGrained"


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

No guarantees for correctness, though.

  • $\begingroup$ Marvelous, I have a long way to go in learning Mathematica. $\endgroup$
    – Boris
    Jun 29, 2021 at 15:35
  • 2
    $\begingroup$ Thank you! "Marvelous, I have a long way to go in learning Mathematica." -- But you actually did everything right! You used sint in the Listable way so that vectorization was enabled (in principle), you created packed arrays of random numbers with RandomReal and used parallelization only at the outer level. The only issue was that Zeta was implemented for correctness on (almost) the whole complex plane. Since you need only a simple, compact subregion of the complex plain and only quite low precision, this is a thing that we can exploit here. $\endgroup$ Jun 29, 2021 at 21:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.