I have the formula:
    y = 2^n-1 + [rand() MOD (2ⁿ-1)]

[edit] or rather:
    y = 2^n-1 + [rand() MOD 2^n-1]