generating a normal distribution on gmp arbitrary precision

Question

So, I'm trying to use gmp to some calculations I'm doing, and at some point I need to generate a pseudo random number (prn) from a normal distribution.

Since gmp has a uniform random variable, that already helps a lot. However, I'm finding difficult to choose which method I should use generate the normal distribution from a uniform one. In practice, my problem is that gmp only has simple operations, and so for instance I cannot use cos or erf evaluations, since I would have to implement all by miself.

My question is to what extent can I generate prn from a normal distribution on gmp, and, if it is very difficult, if there is any arbitrary precision lib which already has normal distribution implemented.

As two examples of methods that do not work (retrieved from this question):

Ziggurat algorithm uses evaluation of f, which in this case is an non-integer exponential and thus not supported by gmp.

Box–Muller Transform uses cos and sin, which are not supported by gmp.

Related question: http://stackoverflow.com/questions/75677/converting-a-uniform-distribution-to-a-normal-distribution — gcbenison, Jun 11 '12 at 15:03
I've edited my question with more information about that question. I saw it before, now I've added a link to it. — Jorge Leitao, Jun 11 '12 at 15:10
gmp supports integer, rational and float numbers. What are you using? — salva, Jun 11 '12 at 15:36
sorry, I forgot, float numbers. Never heard of gaussian distributed rational numbers, but it is a valid question! xD — Jorge Leitao, Jun 11 '12 at 15:42
and even integers if your μ and σ parameters are big enough! — salva, Jun 11 '12 at 16:19

score 1 · Answer 1 · answered Jun 11 '12 at 15:26

1

The Marsaglia polar method would work, if your library has a ln.

answered Jun 11 '12 at 15:26

Thom Smith

13,916
6
45
91

1

And a square root, don't forget. Though I suppose one could implement Newton-Raphson to approximate a square root if needed. – Managu Jul 03 '12 at 15:03

score 1 · Answer 2 · answered Jun 11 '12 at 16:16

Combine a library able to generate a random numbers for a N(0,1) distribution as doubles with the uniform generator of GMP.

For instance, suppose your normal generator produced 0x8.F67E33Ap-1

Probably, just a few of those digits are really random, so truncate the number to a fixed number of binary digits (i.e. truncating to 16 bits, 0x8.F67E33Ap-1 => 0x8.F67p-1) and generate a number uniformly in the range [0x8.F67p-1, 0x8.F68p-1)

For a better approximation, instead of using a uniform distribution, you may like to calculate the values of the density function at the interval extremes (double precision is enough here) and generate a random number with the distribution associated to the trapezoid defined by those two values.

Another way to solve that problem is to just generate a table of 1000, 10000 or 100000 mpf values where N(x) becomes 1/n, 2/n, etc. then, use the uniform random generator to select one of these intervals and again, calculate a random number inside the selected interval using a uniform or linear distribution.

that is not the solution for this problem, but nevertheless is a good approximation. thanks. — Jorge Leitao, Jun 14 '12 at 22:08

score 1 · Answer 3 · answered Jul 03 '12 at 13:18

1

I ended up using mpfr which is essentially gmp with some more functionality. It already has a normal distribution implemented.

answered Jul 03 '12 at 13:18

Jorge Leitao

19,085
19
85
121

Indeed; the algorithm MPFR implements was given by [C.F.F. Karney](https://arxiv.org/pdf/1303.6257.pdf). Karney, C.F.F., 2016. Sampling exactly from the normal distribution. ACM Transactions on Mathematical Software (TOMS), 42(1), pp.1-14. – Peter O. Apr 03 '22 at 03:20

generating a normal distribution on gmp arbitrary precision

3 Answers3