Looking for more algorithms for quasi-random numbers

[ttw]

This is one of the sequences from my paper in "Computational investigations of low-discrepancy point sets II" from the 1994 Las Vegas Conference on Monte Carlo and Quasi Monte Carlo Methods. I have made a single ad-hoc change (described below) that improves distribution for small numbers of points.

The Halton Sequence Phi(N,P) (for odd primes, 2 is a special case not considered here) can be described as:
1. Generate the digits of N in base P (for P an odd prime). Call these digits a(1) to a(k) where k is the maximum number of digits needed. (There should be lots of subscripts but I'll treat each prime separately to reduce index management.)
N=Sum from j=1 to k of a(j)*P^(j-1), that is: a(k)a(k-1)...a(2)(a(1).
2. Reverse the digits: a(1),a(2)....a(k-1),a(k) is resulting string.
3. Treat this string as a fraction with a decimal point (p-ary point?) in front.
Example: P=3, N=5: 5(3)=12. Reverse Phi(5,3)=.21(3) or 2/3+1/9 = 7/9
The sequence is very well distributed for example one starts with 0, 1/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9, 1/27, 10,27...
This sequence is uniformly in the unit d-cube using d different primes. For example in 3 dimensions using primes 3 and 5 gives the points: (skipping 0 which sits on the corner of the cube).
(1/3, 1/5, 1/7)
(2/3, 2/5, 2/7)
(1/9, 3/5, 3/7)
(4/9, 4/5, 4/7)
(7/9,1/25,5/7)
(2/9, 6/25, 6/7)
(5/9, 11/25, 1/49)
(8/9. 16/25, 8/49)
etc.
The process is sometimes termed a Kakatumi-von Neumann odometer.

There is a problem that I noticed about 1967 or so when I started working on quasi-Monte Carlo. For large base, the Halton Sequence produces strongly correlated points until enough points are generated. (This happens with all quasi-random sequences but not as severely.) Take the first few points using bases 101 and 103.
(1/101, 1/103)
(2/101, 2/103)
...
(100/101, 100/103)
(1/10201, 101/103)
(102/10201, 102/103)
(203/10203, 1/10609)
etc.
In 1993-1995 period, I figured out to multiply each numerator by a number (I called a spin) to break this up. Then I used the fact the fractional parts of square roots of primes are independent to do the following seemingly strange rule.

Give a prime P, to find a multiplier S, do the following.
1. Compute the nearest integers to the fractional part of the Sqrt(P), call these H and L (high and low, one is above and one below the number).
2. Compute the continued fraction of H/P and L/P; each generates a string of partial quotients. The multiplier S is the one of these satisfying the following:
3. A. Chose the one for which the sum of the partial quotients is smallest.
B. If tied, chose the one with the smallest maximum partial quotient.
C. If tied, chose whichever H/P or L/P is closest to the fractional part of the square root.
D. Ad Hoc Alert: if P=41, use 16. (To avoid 17/41 being close to 12/29. The only such case in all primes less than 2^32)

4. To generate the modified Phi sequence Phi(N,P,S): generate the digits of N base P as above and reverse. Multiply each digit by S modulo P and sum as above.
Examples: 3 dimensions: P=3, 5, 7 have S= 2, 2,and 5 respectively.

(2/3, 2/5, 5/7)
(1/3, 4/5, 3/7)
(2/9, 1/5, 1/7)
(8/9, 3/4, 6/7)
etc.
For the pathological case 101 and 103, the multipliers are 6 and 16 respectively (not the best but that's another post sometime).
(6/101, 16/103)
(12/101, 32/103)
(18/101, 48/103)
Clearly more spread out than the original Halton Sequence.

I've got some more but Mordechay Levin's paper ArXiv 1806 shows that even the original Halton Sequence hits the theoretical lower bound for dimensions 2 and up so great changes cannot be had by tinkering. I do have a bit better, but it's even longer to compute and I haven't tested the new ideas thoroughly.

HP50g code:(Number, Prime, Spin, Top, Bottom)

<< 0 1 -> N P S T B << WHILE N REPEAT
N S * P MOD T P * + 'T' STO
P 'B' STO* N P IQUOT 'N' STO END
T B />> >>

Not necessarily the fastest but eas to work with. (If I've done this right.)