Proper scaling of randoms

[oleg@pobox.com]

This message derives an integer interval mapping function that is not
only provably within the given range but also "uniform". This is the
function to use when mapping random integers from one interval to a
smaller one.

Problem: given an integer n within [0..maxn], design a scaling
function sc(n) that returns an integer within [s..t], t>=s.
The function sc(n) must be 'monotonic': 
	0<=a < b ==> sc(a) <= sc(b)
and it must map the ends of the interval:
	sc(0) -> s and sc(maxn) -> t.

One such function has been designed by Simon Thompson in his book:
	sct(n) = n `div` ((maxn+1) `div` (t-s+1)) + s

Another function has been derived in the previous message on this
list:
	sc1(n) = (n*(t-s)) `div` maxn + s

The sct() function has a serious drawback: it is not guaranteed that
sct(n) <= t (more about it below). The sc1() function is provably
within the given range, and monotonic. Alas, it is not "uniform."

When it comes to mapping of random integers, the mapping function must
possess an extra property: it must "map uniformly." That is, given
random numbers within [0..maxn] such that Prob(n==i) - const, we must
aim at Prob(sc(n)==i') - const forall i' in [s..t]. The mapping
function sc1(n) does not have this property.

Any mapping function sc(n) from [0..maxn] to [s..t] for (t-s+1) <
(maxn+1) maps several integers from the source interval to one integer
of the target interval. An integer of the target interval therefore denotes an
equivalence class of the integers from the source interval.
The uniformity property will be satisfied only if the equivalence
classes have the same size. This is only possible if (maxn+1) `rem`
(t-s+1) == 0.

Given this condition, the original Thompson's algorithm works.  If
(t-s+1) divides (maxn+1) exactly (i.e., exists denom
integer. denom*(t-s+1) == maxn+1), we have
    j*denom <= n <= j*denom+denom-1 ==> sct(n) = j+s, all j=0..(t-s) 
Thus each equivalence class has the size denom, and mapping is
perfectly uniform. As a consequence, s <= sct(n) <= t, given our
assumption (maxn+1) `rem` (t-s+1) == 0.

If this assumption does not hold, sct(n) may be greater than t (which
Dimitre Novatchev has demonstrated). Alas, in this case there is no
perfectly uniform mapping. The Thompson's algorithm misbehaves. The
worst case for the algorithm is when t = s + (maxn+1)/2. For example,
this occurs when we try to map [0..65535] to [0..32768]. In this case,
t-s+1 = 32769, and denom = 65536 `div` 32769 == 1, thus the mapping
becomes sct(n) = n. In almost half of the time, this mapping exceeds
the upper limit 32768.

OTH, the mapping sc1(n) = n*(maxn+1)/2 `div` maxn + s behaves in this
pathological case rather well. Not only sc1(n) <= t (this is always
guaranteed) but also the mapping is uniform: all equivalence classes
of the source interval [0..65535] except the last two classes have the
size of two.

Thus the best solution should be to "average" sc1(n) and sct(n).
This "averaging" is done by the formula
	sca = (n*(t-s) + (n*(denom-1) `div` denom)) `div` maxn + s
		where denom = (maxn+1) `div` (t-s+1)

It is easy to see that sca(0) -> s.
sca(maxn) = (maxn*(t-s) + (maxn*(denom-1) `div` denom)) `div` maxn + s
          = (maxn*(t-s) + p) `div` maxn + s
	    {where 0<= p = (maxn*(denom-1) `div` denom) <= maxn-1 }
	  = t-s+s
	  = t
It is easy to show that sca(n) is monotonic (as defined above).

If maxn+1 is evenly divisible by (t-s+1) (that is, denom*(t-s+1) ==
maxn+1), we have
    j*denom <= n <= j*denom+denom-1 ==> sca(n) = j+s, all j=0..(t-s) 

This requires a proof.
We have:
	sca(j*denom) 
      = (j*denom*(t-s) + (j*denom*(denom-1) `div` denom)) `div` maxn + s
      = (j*denom*(t-s) + j*(denom-1)) `div` maxn + s
      = j*(denom*(t-s) + denom-1) `div` maxn + s
      = j*maxn `div` maxn + s
      = j+s

	sca(j*denom+denom-1) 
      = ((j*denom+denom-1)*(t-s) + ((j*denom+denom-1)*(denom-1) `div` denom)) `div` maxn + s
      = ((j*denom+denom-1)*(t-s) + (j+1)*(denom-1)-1) `div` maxn + s
      = (j*maxn + (denom-1)*(t-s+1)) `div` maxn + s
      = j+s
The required result then follows by the monotonicity property.

In the most pathological case of maxn= 65535, s=0, t= 32768.
We have denom = 1 and sca(n) becomes sc1(n), which, in this case,
behaves very well.

Thus the sca mapping is the optimal one. It is guaranteed to keep the
result within the specified bounds at all times. When the perfect
uniformity of mapping is possible, it is always achieved.

The implementation isn't too complex either.

> scaleSequence :: Int -> Int -> [Int] -> [Int]
> scaleSequence s t
>   = map scale
>     where
>       scale n = (n*range + ((n*(denom-1)) `div` denom)) `div` maxn + s
>       range   = t - s
>       maxn    = modulus - 1
>       modulus = 65536
>	denom   = (maxn+1) `div` (t-s+1)

Main> any(>966) (scaleSequence 1 966 [0..65535])
False
Main> any(==966) (scaleSequence 1 966 [0..65535])
True 
Main> any(>32768) (scaleSequence 0 32768 [0..65535])
False 

> test = let sq = [32768-1..32768+1600] in
>   let result =
>        map (\x -> let t = length x in (t `seq` t)) $! groupBy (\(a,_) (b,_) -> a == b) $!
>        sortBy (\(a,_) (b,_) -> compare a b) $! 
>        zip (scaleSequence 0 32768 sq) $! sq
>   in (maximum $! result, minimum $! result)
==> (2,1) -- e.g., the mapping is roughly uniform