TAOCP 3.6 Exercise 11
Let the modulus be $MM=2^{52},$ and replace all arithmetic modulo $2^{30}$ by arithmetic modulo $2^{52}$.
Exercise 11. [**] [M25] Assuming that floating point arithmetic on numbers of type double is properly rounded in the sense of Section 4.2.2 (hence exact when the values are suitably restricted), convert the C routines run_array and run_start to similar programs that deliver double-precision random fractions in the range $[0,.,.,1)$, instead of 30-bit integers.
Verified: yes
Solve time: 1m19s
Let the modulus be
$MM=2^{52},$
and replace all arithmetic modulo $2^{30}$ by arithmetic modulo $2^{52}$. Since IEEE double-precision arithmetic represents all integers exactly up to $2^{53}$, the subtraction
$x-y$
is exact whenever $0\le x,y<2^{52}$. Hence the operation
#define MM (4503599627370496.0) /* 2^52 */
#define mod_diff(x,y) ((x)<(y)?(x)-(y)+MM:(x)-(y))
computes $(x-y)\bmod 2^{52}$ exactly in floating point arithmetic. The arrays ran_x and aa should therefore be declared as type double, and the recurrence becomes
$X_j=(X_{j-100}-X_{j-37})\bmod 2^{52}.$
The output fractions in $[0,1)$ are obtained by multiplying each generated value by
$2^{-52}=1/MM.$
The initialization routine ran_start is modified similarly. All occurrences of the old modulus $2^{30}$ are replaced by $2^{52}$, and the cyclic left shift by one bit is performed on 51-bit quantities:
ss+=ss; if (ss>=MM) ss-=MM-2;
This operation remains exact because all intermediate values are less than $2^{53}$. After the warm-up calls to ran_array, the routine returns numbers
aa[j]*(1.0/MM)
instead of integers. Thus the resulting program generates properly rounded double-precision random fractions with 52 random bits in the mantissa.