monte_carlo_integral looks deterministic

[edd8e884-f507-429a-b577-5d554626c0fe]

Tested on various computers, with various random seeds, the following invariably provides the same result:

var('x,y,z')
f = x*y*z
monte_carlo_integral(f, [0,0,0],[1,1,1], calls=1000)
(0.12001500488162796, 0.004393500801071119)

The result should depend on the current ranstate.

CC: @videlec

Component: algebra

Issue created by migration from https://trac.sagemath.org/ticket/32456

[nbruin]

comment:1

Well, it's certainly explicable: according to https://www.gnu.org/software/gsl/doc/html/rng.html the rng initialization sets up a default seed of 0 and no other seed is initialized in our code.

I don't know if the strict definition of monte carlo integration requires the use of true random numbers and/or the use of a pseudo-random generator with a non-deterministically initialized seed. The theory definitely needs it, because otherwise the probabilistic statements about accuracy could be engineered to fail for a particular integral, making use of the known sample sequence.

Fixing this is pretty simple: put a gsl_rng_set call in somewhere. The tricky bit is finding where to get the appropriate seed value from (since we still need to be able to get reproducible results)

[mkoeppe]

comment:2

sage.misc.randstate likely needs updating.