mt_rand([int $min, int $max])で $min > $max のときの挙動 (理由解説編)

前回の話の続き。引数の乱数範囲の値が min > max のときにおかしいという話ですが、これは普通のrand()でも同じマクロを使っているので同様です。
で、ソースコード眺めてたら、rand()関数の上にこんなコメントが。

/*
 * A bit of tricky math here.  We want to avoid using a modulus because
 * that simply tosses the high-order bits and might skew the distribution
 * of random values over the range.  Instead we map the range directly.
 *
 * We need to map the range from 0...M evenly to the range a...b
 * Let n = the random number and n' = the mapped random number
 *
 * Then we have: n' = a + n(b-a)/M
 *
 * We have a problem here in that only n==M will get mapped to b which
 # means the chances of getting b is much much less than getting any of
 # the other values in the range.  We can fix this by increasing our range
 # artifically and using:
 #
 #               n' = a + n(b-a+1)/M
 *
 # Now we only have a problem if n==M which would cause us to produce a
 # number of b+1 which would be bad.  So we bump M up by one to make sure
 # this will never happen, and the final algorithm looks like this:
 #
 #               n' = a + n(b-a+1)/(M+1) 
 *
 * -RL
 */

適当に訳すと、指定された範囲の乱数を作るのに、剰余を使うのを避けたかったということ。
0からRAND_MAX (unixだと、2147483647 = 0x7FFFFFFF)までの乱数をmin〜maxの間の範囲に収めるには、min + n % (max-min) と、余り(剰余)を求めるのがすぐ思いつく考え方ではあります。サイコロ代わりに1から6を出すなら、1 + n % 6 ですね。コンピュータがどうやって割り算をするのかは調べてて挫折したのですが、例えば、2という数で割る場合、デーモン小暮閣下の100,045という年齢だと余りは1。計算しなくても分かるのは、下1桁だけで判断できるから。2進数でも最下位ビットのみで判断でき、そうなるとそれ以外の30bitは無意味となってしまうのでいけてない。
なので、余りを計算するのではなく、31bitの乱数値を指定した範囲にぎゅっと縮めてマッピングするようにしたというわけです。

変な計算式の理由

上のコメントでは最終的に、n' = a + n(b-a+1)/(M+1)という謎の式が出てきています。これが、前回の問題の原因であったりするのですが、この式の意味について見ていきます。
以下の記号は、このような意味です

n: 乱数値。0 〜 Mの間の値 (int)
n': min <= n' <= maxの範囲にマップされた乱数値 (int)
min: 範囲の開始値 (int)
max: 範囲の終了値 (int)
M: RAND_MAX。nのとりうる最大値 (int)

指定した範囲に乱数をぎゅっと縮めるというのは、グラフにするとこういうイメージになります。

切片が min で、 n = M のときに n' = max となる直線ですね。但し、n'は整数値ですので、intに型キャストする必要があります。

すると、小数点を切捨ててしまうので、 n = M で確かに n' = max とはなるものの、範囲が均等に割り振られておらず、maxの数字が出てくる確率だけ、1/RAND_MAXという激レアになってしまいます。
しかたないので、

と、傾きに (1 + max - min)/M と1を加えます。そうすると、範囲が均等に割り振られますが、今度は1/RAND_MAXの確率で n' = max + 1 という値を出してしまいます。

この対策に、傾きを (1 + max - min)/(M + 1) としています。これで、n = M でもmaxの値を返すようになります。M >> 1 であるため、M ≒ M+1 でほぼ同じ値であり、影響はほぼ無視できます。

max < minのときの問題

では、今のが max < min であったときにどうなるかですが、(1 + max - min)/M と1を加えるのが、まず逆方向に作用してしまいます。

これで、max (という名前だけれど、範囲の最小値)が出現しなくなり、max+1 も1/RAND_MAXの激レアのままです。
そして、傾きを (1 + max - min)/(M + 1) として

max+1 も出現しなくなり、 min > max の場合、返却される乱数の範囲は、max + 2 <= n' <= min となります。