O problema da solução ingênua da obtenção de valores aleatórios.

Considere que você quer fazer um programinha que “role” um dado de 6 faces. Geralmente, alguém faz algo assim:

int r;

srand( time( NULL ) );
r = rand() % 6 + 1;

A ideia é que obter o resto da divisão por 6 acarretará em valores menores que 6 (de 0 até 5) e somando 1, obtemos um dos 6 valores possíveis. Mas há um problema!

Considere que rand() retorna valores entre 0 e RAND_MAX que, em arquiteturas de 32 bits ou mais pode equivaler a 2^{31}-1, ou seja, o valor máximo é 2147483647 (31 bits setados). Isso representa um em 2147483648 possíveis valores e essa quantidade não é divisível por 6! Outra forma de pensar nisso é que temos 357913941 blocos de 6 valores possíveis e um bloco incompleto com os 2 primeiros valores, ou seja, no último bloco faltam 4 possibilidades para ter uma distribuição uniforme de valores possíveis.

Outra forma de pensar no assunto é essa: Considere que você tenha um dado de 6 faces e queira escolher, aleatóriamente, 1 entre 4 itens (a, b, c, d). As possibilidades são: (1=a, 2=b, 3=c, 4=d, 5=a e 6=b). Os itens a e b têm mais chance de serem escolhidos (2 chances em 6, cada) porque não temos uma distribuição de valores aleatórios uniforme! A solução é nos livrarmos dos 2 últimos valores, seja “jogando” o dado novamente, colhendo apenas um dos 4 valores possíveis, ou aplicando alguma regra maluca para nos livrarmos desses dois valores. Infelizmente, a segunda hipótese não é tão simples…

Considere agora o retorno de rand() como mostrado no fragmento de código lá em cima. Obteremos valores entre 1 e 6. Acontece que temos um espaço amostral de 2^{31}=2147483648 valores (RAND_MAX) dividido em \frac{2^{31}}{6}+1=357913942 blocos de 6 elementos, onde o último bloco só tem 2 (2^{31}\mod\,6=2). Assim, a chance de obtermos 1 ou 2 são de \frac{357913942}{2147483648} e de obter os outros 4 valores são de \frac{357913941}{2147483648}. Ora \frac{357913941}{2147483648} < \frac{1}{6} < \frac{357913942}{2147483648}. Não precisa ser um gênio para entender que 1 e 2 podem ocorrer mais vezes que os outros 4 valores.

A solução para esse problema não é trivial e, acredito, que nem mesmo seja possível num ambiente computacional binário. Em primeiro lugar, qualquer algorítmo que nos dê um valor aleatório conterá um valor entre 2^n outros. Ao obter o resto da divisão para limitar o escopo, criaremos, automaticamente, um último grupo incompleto, se o divisor não for, também, um valor 2^m menor ou igual a quantidade de amostras possíveis… Tudo o que podemos fazer é diminuir o erro relativo das distribuições e nisso rand() pode ser confiável, já que os erros relativos das distribuições, no nosso exemplo, são de, respectivamente, 1.86\cdot10^{-7}\,\% para 1 e 2; e 9,31\cdot10^{-8}\,\% para o resto:

\displaystyle \begin{matrix} E(\%)=100\cdot\left|\frac{X - X_{ideal}}{X_{ideal}}\right|\\ \\ E_{1,2}(\%) = 100\cdot\left|\frac{\frac{357913942}{2147483648} - \frac{1}{6}}{\frac{1}{6}}\right|\approx1.86\cdot10^{-7}\,\%\\ \\ E_{3..6}(\%) = 100\cdot\left|\frac{\frac{357913941}{2147483648} - \frac{1}{6}}{\frac{1}{6}}\right|\approx9.31\cdot10^{-8}\,\%\\ \end{matrix}

Ou seja, a cada 500 milhões (100\cdot\frac{1}{E_{1,2}(\%)}) de escolhas, aproximadamente, teremos mais 1’s e 2’s que o resto. Poderíamos melhorar esse erro de distribuição aumentando a quantidade de bits do valor aleatório. Se tivessemos 2^{64}=18446744073709551616 amostras (64 bits) teríamos 4 valores no último grupo e 3074457345618258603 grupos. As chances dos valores de 1 a 4 seriam de \frac{3074457345618258603}{18446744073709551616} e de 5 e 6 seria de \frac{3074457345618258602}{18446744073709551616}, Os erros de distribuição relativos da primeira faixa seriam de, aproximadamente, 1,8\cdot10^{-17}\,\% e da segunda, 1,2\cdot10^{-18}\,\%, ou seja, em 50 quatrilhões de “jogos” observaríamos o problema, que só se repetiria nos próximos 50 quatrilhões de jogadas, possivelmente.

Ahhh, mas eu sou um chato né? O erro é tão pequeno para aplicações práticas, não é? Acontece que RAND_MAX só tem 31 bits de tamanho em certos compiladores e certas arquiteturas. Outras podem usar valores muito menores, como 32767 (15 bits). De fato, a especificação ISO 9989 sugere que esse seja o valor mínimo para RAND_MAX, mas ele pode ser ainda menor!

Assim, neste caso, teríamos 5462 blocos e o último teria apenas 2 valores (como com 31 bits) e as distribuições seriam de \frac{5461}{32768} (para 3 até 6) e \frac{5462}{32768} (para 1 e 2), o que daria, para o maior erro relativo, um valor muito grande (0.012% para 1 e 2). A cada 10000 jogadas você perceberia o problema. Você já pode considerar esse dado “viciado”…

Ahhhh… note que não estou considerando que rand() geralmente é implementado usando LCG que oferece menos aleatoriedade nos bits inferiores. Assim, obter o resto de uma divisão por um valor pequeno tende a criar valores ainda menos aleatórios.