Uniformly distributed random numbers

[albinahlback]

Arb doesn't seem to have this functionality (which is natural given its number theoretic orientation), so I added the functionality for it (and hopefully correctly).

Please have in mind that this is basically my first time doing something in C, so let me know if something is wrong or could be improved.

[fredrik-johansson]

I agree this is a very useful feature, but the implementation needs work:

• You can't generate "uniformly distributed random numbers in [0, 1) whose exponent is 0" because such numbers will always lie in [0.5, 1).
• It's not clear what the mag_bits parameter is for.
• MAG_MAN(x) = n_randlimb(state) will not create a valid mag instance.
• Doing fmpz_zero(ARF_EXPREF(x)); will have no effect when you assign to the arf_t later.

I think what you really want to do here is (in principle) sample a uniformly random real number from [0, 1), in the case of arf_t round to the nearest floating-point number in the prescribed rounding mode, and in the case of arb_t return an enclosure (which in practice always will be nonexact since exact dyadic numbers have measure zero). It's a bit tricky to do this in a mathematically rigorous way. You can have a look at how MPFR does it.

In practice, I'd perhaps just generate a uniformly random integer n with prec+128 bits and arf_set_round or arb_set_round the result to n*2^(-prec-128) with a precision of prec bits. This isn't theoretically perfectly uniform, but the deviation from ideal randomness would be completely negligible.

[albinahlback]

Thank you for your comments, Fredrik!

Yes, after some writing on paper, I now think that I understand how one can derive the proper distribution for such a random number (even though I did not write down an explicit formula for the mantissa). It seems like your suggested would be faster in either case, so I revised the code accordingly.

[fredrik-johansson]

I don't see the utility of arb_randtest_uniform returning a ball with random midpoint in [0,1) and random radius in [0,1). Do you have a use case for this? What I meant was to use exactly the same algorithm as arf_randtest_uniform but with arb_set_round_fmpz instead of arf_set_round_fmpz, so that you get an accurate ball representing a random real number in [0,1).

More comments:

• I would name these functions arf_urandom / arb_urandom, as randtest is meant for non-uniform testing random distributions.
• arf_urandom should probably take the rounding direction as input. Of course, with rounding up, you could get 1. Note that MPFR has a separate urandomb function which doesn't take a rounding mode and generates a uniformly distributed float, guaranteed to be < 1, while its urandom function takes a rounding mode and generates a uniformly random real, which gets rounded. We could have both versions.
• On that note, arf_urandom and arb_urandom should maybe explicitly sample real numbers from [0,1] rather than [0,1), since this is the more common definition of the uniform distribution.
• Instead of fmpz_sub_si(ARF_EXPREF(x), ARF_EXPREF(x), prec), use arf_mul_2exp_si (and likewise for the other types).
• You don't want fmpz_randtest_unsigned as this doesn't give a uniform distribution. Unfortunately there isn't a Flint function that gives you a uniform random integer in [0, 2^n-1], but you can use fmpz_randm. In fact, this will make more sense anyway assuming that you want to sample from [0,1], in which case you want to sample integers from [0, 2^n]. Alternatively (but it's more complex), you could use fmpz_randbits along with a randomly generated exponent following a geometric distribution (probability 1/2 for exponent 0, probability 1/4 for exponent -1, etc.).
• Make sure you put all variable declarations before code (compiler warnings).
• Some test code would be good. It's hard to test random numbers, but the test could be something very basic, like computing a few thousand random numbers and checking that the mean and variance are somewhere near what you expect.

[peteroupc]

My notes on uniformly distributed random numbers:

There is a closely related concept to arbitrary-precision uniform random numbers. Specifically, these numbers represent uniformly distributed random variates with an arbitrary number of fractional digits that are sampled on demand.

They were called "u-rands" by C.F.F. Karney in "Sampling exactly from a normal distribution", 2016, as well as the "geometric bag" by Flajolet et al., "On Buffon machines and numbers", 2010/2011. I call them "partially-sampled random numbers".

For example, take the following number: [2, [1, 3, 8]]. This represents a random number in the interval [2.138, 2.139]. As the need arises, we can sample additional digits uniformly at random, so that, for example, the number can become [2, [1, 3, 8, 7, 5]], bringing the number in the interval [2.13875, 2.13876]. Karney's paper, for example, uses this format in an algorithm to sample normally distributed variates to arbitrary precision and without calculating any transcendental functions.

[albinahlback]

Thanks for the patience, Fredrik.

I made a test for the arb one, which seems to pass. The current test should fail if it goes beyond one percent deviation from the "supposed" mean value 1/2 and variance 1/12.

[fredrik-johansson]

Thanks for keeping up the work :-)

• It's still not the arb_urandom I had in mind. What I had in mind is the same code as arf_urandom but with arb_set_round_fmpz instead of arf_set_round_fmpz.
• Then, arb_urandom does not need the rounding mode as input.
• In arb/test/t-urandom.c, make sure all arb_inits have matching arb_clears
• In arf_urandom, you need to clear n.
• Put your own name instead of mine in the copyright notices in the new files

[albinahlback]

* Put your own name instead of mine in the copyright notices in the new files

Are you sure? I have changed it in the new files, but unlike you I was not blessed with an ASCII-compatible Swedish name haha.

[fredrik-johansson]

At least your name is googleable. I share mine with death metal singers, bandy players, and at least two other scientists ;-)

All right, this looks good enough to merge. I will maybe make some minor changes later.