Python large random integers and the precision of the Mersenne Twister

Question

I am wondering what the "precision" of Python's random function means. It is described here:

Almost all module functions depend on the basic function random(), which generates a random float uniformly in the semi-open range [0.0, 1.0). Python uses the Mersenne Twister as the core generator. It produces 53-bit precision floats and has a period of 2**19937-1.

(http://docs.python.org/library/random.html?highlight=mersenne%20twister, accessed 20120727)

What interests me is that I can generate very large random integers (long integers) that appear to have considerably more than 2^53 precision. For instance (using Ipython):

In [1]: from math import factorial as F
In [2]: from random import randint as R
In [3]: R(1, F(900))
Out[3]: 55655511302846458744179265243566263049348396362730789786376014445325896599604354914431619960209388364677180234108513221468671377813842671874148746886513973171423907294544220953849330089822288697383171078250181973489187774341795574648920075697792011317798969959919449394758519496792725695600701199089972009688412593325291810024048811890509220571436407156566269358600296506017343255050788936280200352509087073097532486502694101150248815092174847010359868156616901409331336760344351058867833528749797221612169430654334458578364850198977511061993233818849689759090377347376020160658362459773356292085856906573553086825560047089834757501023094429371408722563891227474029563545206865055657504766128286451181119906678062837368414582707728324415466848186858173236300969443478496634754744888060794778485246692104851885847515244146665974598354436781340057667983223238998674622833320199904840957000014767293658171874973067958145430346745707636676061629278168015549755791407108399231392952706279787486238512258804098030513575025870504347283221015756832157863142353915612138589145084128778032995695113870365505775392647256056048691602676699581153972467494111720212363912926352356346807790816796784781384561736415741104584667536002819103176714157723039428367564698686945824679882523439229215035996634289075127375256728472056511244548311771570743103809147045947583819651257115044154025329883682429231394004470689760531056853018427649916035935302356382633012319775473728455377657692268855776796385819792347680100513177355101630543290088996770992548670273727988974570199179655691444984337837105283447276788151912408533352627494948390016029881755603243934955207024221452181883522004648595373130617729041347013155205217774450836687880723915563507108222768637840614647145898936109917167237397888104669458661404234553707323638883064861414284282190898741067404128885188113697448726481104763682489126524054241797759521120664366845719767486252884585742737830119890190213053751046461419643379561983590174574185268661318409035375114305279020423595250660644954841798619767985549553380200803904976806468796334648515423467654573415304912570341635682203261002606581817207689816015969520503052648773609840260050676394927780076948629298559638703440007364834579712680931643829764810072128419905903786966L

I am wondering in what sense 53 bits is the precision limit of the random function. And concretely, if I ask Python to return pseudo-random numbers between 1 and some very large upper bound, is it not true that all integers in that range have an equal likelihood of being returned?

Answer 1

What Python does is generate 64 bits of randomness from calling the 32-bit version of MT19937 twice, but since that number is constrained to [0.0, 1.0) the result is constrained to 53 bits of precision (limitation of floating-point format).

Python random.py is capable of pulling bits from /dev/urandom or Windows CryptGenRandom (if supported), if you use the random.SystemRandom class. Otherwise, it generates larger numbers by pulling successive bits from repeated calls to MT19937.