O aleatório não tão aleatório.

[fredericopissarra]

Suponha que você queira fazer um programinha que "jogue" um dado D6 n vezes e obtenha n valores, aleatoriamente. Nada mais simples:

srand( time( NULL ) );  // alimenta a semente com um valor "imprevisível".

n = NÚMERO_DE_JOGADAS;
while ( n-- )
  printf( "%d\n", rand() % 6 + 1 );

A lógica é de que rand() devolve um valor entre 0 e RAND_MAX que, ao obtermos o resto da divisão por 6, será limitado entre 0 e 5, somando 1 obteremos entre 1 e 6. Mas, eis o problema, que apresento de outra maneira: Suponha que você tenha 4 itens (a, b, c, d) e queira escolher entre eles de maneira aleatória usando um dado D6... É como se RAND_MAX fosse 5 e a sequência (a,b,c,d) fosse (0,1,2,3). Daí você faria algo assim:

// Suponha que RAND_MAX == 5.
r = rand() % 4;

Se rand() retornar 4 o resultado será 0, se rand() retornar 5, resultará em 1. Isso significa que as 6 possibilidades serão (0,1,2,3,0,1), significando também que 'a' e 'b' terão mais chances de serem escolhidas (2 chances em 6) que 'c' e 'd' (1 chance em 6)... O mesmo acontece com rand() da libc e a rotina do dado, lá em cima... No caso do exemplo dos 4 itens, podemos pensar que o espaço amostral (as 6 possibilidades de rand()) pode ser dividido em 2 grupos de 4 possibilidades, mas faltam duas no último grupo ou seja, o resto da divisão de 6 pelos 4 itens te dá quantos itens existem no último grupo (se for 0, temos um grupo final cheio).

No caso da rotina original, RAND_MAX tende a ser 2³¹-1, ou seja, temos 2³¹ (2147483648) valores possíveis, na maioria dos compiladores para arquiteturas de 32 ou 64 bits. Ao tomar o resto da divisão por 6 teremos (1+2147483648/6) grupos de 6 valores, onde o último grupo é incompleto e tem apenas 2 (2147483648 % 6). Ou seja, 1 e 2 têm mais chances de ocorrer do que 3, 4, 5 e 6.

Num cenário ideal cada um dos 6 valores deveriam ter 1 chance em 6 de ocorrerem (1/6), mas se temos 357913942  grupos (1+2147483648/6), com o último incompleto, 1 e 2 terão 357913942/2147483648 chances de ocorrer, enquanto 3, 4, 5 e 6, 357913941/2147483648. Se você fizer as contas verá que  

Mas ambos os valores são bem perto de 1/6.

Neste ponto você deve estar me achando um chato que adora certos "preciosismos", mas considere o seguinte: RAND_MAX não é sempre 2³¹-1. De fato, a especificação ISO 9989 sugere (mas não limita!) o valor mínimo de 32767 (15 bits). Com isso, a discrepância torna-se mais relevante, já que termos 5462/32768 chances para 1 e 2, e 5461/32768 chances para os outros valores.

Concluo que, usando o método de obtenção do resto por 6, terá um dado "viciado", mesmo que não perceba, depois de muitas jogadas...

Não assuma que lidar com "valores aleatórios" seja uma tarefa simples! Não é! Eis um bom livro sobre geradores de valores (pseudo) aleatórios de distribuição não-uniforme (aqui). Veja o tamanho do bicho (em páginas) e a quantidade de assuntos cobertos... É um excelente material para estudos...

[fredericopissarra]

Existe OUTRO problema com rand(), mas considerei a função um "perfeito" RNG (Random Number Generator), para fins do texto acima...

[fredericopissarra]

Se a limitação da faixa for um valor 2^n, então não temos o problema de distribuição não uniforme... Considere RAND_MAX como 2^m, onde m>n. O número de grupos de n itens é (m-n) e todos estarão cheios (2^m mod 2^n = 0). Neste caso, se rand() for um RNG perfeito, temos as mesmas chances para todos os n itens...

[fredericopissarra]

Outro detalhe importante para se lembrar sobre aleatoriedade é que ela não significa que você não possa obter um determinado valor mais vezes do que outros. Se fato, uma sequência do tipo (1,2,2,3,1,5,4) continua sendo aleatória, mesmo que 1 e 2 sejam repetidos... O problema todo está na distribuição das chances de obtenção dos valores...

Outro exemplo é a MEGA SENA... todo mundo faz jogos com valores "aparentemente" aleatórios, mas se esquecem que, por exemplo, a sequência (1,2,3,4,5,6) tem a mesma chance de ser sorteada que qualquer outra, se considerarmos a aleatoriedade perfeita (uniformemente distribuída).

[fredericopissarra]

Um exemplo de geração de ticket para MEGA SENA tendando corrigir a discrepância da não uniformidade de distribuição, por rejeição:

/* Sorteio MEGA SENA, corrigindo a distribuição não uniforme. */
#include <stdio.h>
#include <stdlib.h>
#if defined(__x86_64) || defined(__i386)
  #include <immintrin.h>
#endif

// Usar rand() não é uma solução ideal, mas pode ser preciso,
// se __RDRND__ não estiver definido.
#ifndef __RDRND__
  #include <unistd.h>
  #include <fcntl.h>
  #include <time.h>
#endif

/* Para obter aleatoriedade uniforme e garantir que
   o valor seja verdadeiramente aleatório, uso a instrução
   RDRAND, disponível nos processadores Intel/AMD. Usei
   _rdrand16_step() na esperança que a entropia não seja
   prejudicada entre chamadas.

   Note, também, que essa rotina retorna o resto da divisão por
   64 (6 bits). Isso garante que a distribuição seja uniforme
   e, mais adiante, simplesmente discarto os valores maiores que
   60. */
#ifdef __RDRND__
unsigned short get_random_64( void )
{
  unsigned short r;
  unsigned int retries;

  retries = 10;
  do
    if ( _rdrand16_step( &r ) )
      return r & 0x3f;  // o mesmo que 'r % 64'.
  while ( --retries );

  fputs( "\n\033[31;1mERROR\033[m: Could not get random value.\n", stderr );
  exit( EXIT_FAILURE );
}
#else
// Usa rand(), se RDRAND não estiver disponível.
unsigned short get_random_64( void )
{
  return rand() % 0x3f;
}

// Obtem uma semente mais aleatória que time(NULL)...
// Não funciona no Windows!
unsigned int get_random_seed( void )
{
  int fd;
  unsigned int r;

  fd = open( "/dev/urandom", O_RDONLY );
  if ( fd < 0 )
  {
    fputs( "\n\033[31;1mERROR\033[m: Cannot open random device!\n", stderr );
    exit( EXIT_FAILURE );
  }

  if ( read( fd, &r, sizeof r ) != sizeof r )
  {
    close( fd );
    fputs( "\n\033[31;1mERROR\033[m: Error reading from random device.\n", stderr );
    exit( EXIT_FAILURE );
  }

  close( fd );

  return r;
}
#endif

int main( void )
{
  /* LEMBRE-SE: O array começa no índice 0! */
  _Bool values[60] = {0};
  int count;
  unsigned int r;

  /* Alimenta a semente, se __RDRND__ não estiver definido. */
  #ifndef __RDRND__
    srand( get_random_seed() );
  #endif

  fputs( "Discarded: ", stdout );

  count = 6;
  do {
    r = get_random_64();

    // Somente entre 0 e 59 são aceitos.
    // NÃO é uma solução ideal o descarte dos valores
    // do último bloco, incompleto, mas é uma solução!
    if ( r > 59 )
    {
      printf( "%2u ", r );
      continue;
    }

    // Se o valor já foi sorteado, tenta de novo...
    if ( values[r] )
      continue;

    // Marca como sorteado.
    values[r] = 1;

    // Obtivemos o valor, continua enquanto ainda precisamos
    // de mais...
    count--;
  } while ( count );

  // Apresenta o ticket de maneira ordenada.
  fputs( "\nTicket:    ", stdout );
  r = 0;
  while ( r < sizeof values / sizeof values[0] )
  {
    if ( values[r] ) 
      printf( "%2u ", r + 1 );

    ++r;
  }
  putchar('\n');

  return EXIT_SUCCESS;
}

Funciona em Linux (qualquer plataforma), mas para melhores resultados na aleatoriedade, compile para Intel/AMD com:
 

$ gcc -O2 -march=native -o mega mega.c

 

[fredericopissarra]

Pode-se imaginar uma maneira, menos eficiente em termos de performance, para eliminar o problema da distribuição não uniforme trabalhando com ponto flutuante. Infelizmente isso tem alguns problemas...

Se RAND_MAX tem 31 bits de tamanho, usaremos o tipo double, que tem 53, para garantir a precisão e fazer algo assim:

#define RANGE 6

int r;

// r estará entre 0 e (RANGE - 1).
r = rand() / ( ( double )RAND_MAX / RANGE + 1 );

A distribuição dos valores resultantes estará em "clusters"... No caso, cada um dos 6 valores estará dentro de uma faixa de 0.1666 (1/6)... 0 está na faixa [0;0.1666...); 1, [0.1666...;0.33333...); ... Acontece que rand() pode devolver qualquer valor entre 0 e RAND_MAX, ou seja, valores entre [0; 0.1666...*RAND_MAX), neste exemplo, serão necessariamente 0... Ou seja, qualquer valor de rand() retornado entre 0 e 357913941 será zero!

Você acaba de obter um RNG bem menos aleatório que o original.

[fredericopissarra]

Yep... C++11 tem uma série de templates para "uniformização" de distribuições de inteiros e um monte de templates com "random engines". Mas, se você não souber usá-los, vai piorara a coisa... E mais: EU NÃO GOSTO DE C++!

[fredericopissarra]

O problema com rand() é o mesmo problema com a maioria dos PRNG (Pseudo Random Number Generators)... Eles não são, de fato, aleatórios. A não ser por efeitos físicos, como decaimento radioativo, flutuações quânticas, etc,  não há como obter valores realmente aleatórios - pode-se até mesmo usar meios físicos não tão especializados como contar a quantidade de batimentos de asas de um colibri se alimentando em 1 segundo ou alguma outra grandeza que pareça aleatória (algumas certificate authorities usam meios assim!)...

rand() geralmente é implementado usando-se o Linear Congruential Generator, que nada mais é do que uma equação linear:

Onde X0 é a "semente" informada via função srand(). Mas, esse esquema tem um problema... Se m for descartado ou for o tamanho da "palavra", então os bits inferiores tendem a ser menos aleatórios que os bits superiores (veja Knuth)... Felizmente a glibc tenta resolver isso com uma série de equações lineares e uma "mistura" binária "esperta":

// rand_r() é a versão therad-safe de rand() e a base do mesmo.
int rand_r (unsigned int *seed)
{
  unsigned int next = *seed;
  int result;

  next *= 1103515245;
  next += 12345;
  result = (unsigned int) (next / 65536) % 2048;

  next *= 1103515245;
  next += 12345;
  result <<= 10;
  result ^= (unsigned int) (next / 65536) % 1024;

  next *= 1103515245;
  next += 12345;
  result <<= 10;
  result ^= (unsigned int) (next / 65536) % 1024;

  *seed = next;

  return result;
}

Aqui temos 3 valores pseudo aleatórios obtidos com uma única semente que são "misturadas" com outro método PRNG chamado linar feedback shift registers. O quão aleatório isso é, não faço ideia... Mas, o problema permanece: Se eu sei qual é a semente inicial e a equação que é usada para obter os valores subsequentes, eles deixam de ser pseudo aleatórios para serem previsíveis!

Geralmente escolhemos uma semente que também pareça aleatória. Uma das formas é usar o epoch timestamp, obtidas em qualquer sistema baseado em UNIX com a função time():

// time() retorna o número de segundos ocorridos
// desde 1º de janeiro de 1970, às 0:00, UTC.
srand( time( NULL ) );

Mas um sujeito chato e metódico pode obter esse valor de antemão... Uma solução é confiar na "aleatoriedade" de certas APIs e dispositivos contidos no seu sistema operacional para obter uma semente aleatória e depois usar rand() para obter os demais valores mais rapidamente:

// Inicializa a semente com um valor "aleatório" no UNIX:
void setup_seed( void )
{
  unsigned int seed;
  int fd;
  ssize_t size;

  if ( ( fd = open( "/dev/urandom", O_RDONLY ) ) < 0 )
  {
    perror( "open" );
    exit( EXIT_FAILURE );
  }

  errno = 0;
  size = read( fd, &seed, sizeof seed );
  if ( errno || size < sizeof seed)
  {
    fputs( "read error\n", stderr );
    close( fd );
    exit( EXIT_FAILURE );
  }
  
  srand( seed );
  
  close( fd );
}

Ler /dev/urandom é, com toda certeza, mais lento do que usar rand(), mas resolve o problema da semente original. Não resolve o problema da aleatóriedade previsível do LCG ou qualquer outro PRNG, mas minimiza a coisa toda nos dando valores que "parecem" aleatórios.