FWIW, here are two approaches to getting an equi-distributed version of int(n*random()) where 0 < n <= 2**53.  The first mirrors the approach currently in the code.  The second approach makes fewer calls to random().

def rnd1(n):
    assert 0 < n <= 2**53
    N = 1 << (n-1).bit_length()
    r = int(N * random())
    while r >= n:
        r = int(N * random())
    return r

def rnd2(n, N=1<<53):
    assert 0 < n <= N
    NN = N - (N % n) + 0.0  # largest multiple of n <= N
    r = N * random()        # a float with an integral value
    while r >= NN:
        r = N * random()
    return int(r) % n