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.

Indo a extremos para evitar usar ponto flutuante. Você deveria fazer isso também!

Eis um problema que preciso resolver para o meu mais novo projeto de brinquedo: Preciso determinar se uma imagem (stream de vídeo) tem aspecto mais próximo do aspect ratio \frac{4}{3} ou do \frac{16}{9}. A solução mais “simples” é simplesmente dividir o comprimento (width) pela altura (height) e comparar o resultado com uma fração de ‘meio termo”. O “meio termo”, é claro, é a média aritmética simples dos dois aspectos:

\displaystyle \frac{\frac{16}{9} + \frac{4}{3}}{2}=\frac{16+12}{9}\cdot\frac{1}{2}=\frac{28}{18}=\frac{14}{9}

Por que a fração de “meio termo”? É que o vídeo pode ter resolução não padronizada e, neste caso, colocarei “barras pretas” horizontais ou verticais, dependendo do aspecto. Assim, temos:

int is_letterbox(unsigned int width, unsigned int height)
{
  double ratio = ( double )width / height;

  return ratio <= 14.0 / 9.0;
}

Nada mais simples, né? Se o aspecto for menor ou igual a \frac{14}{9}, então encaro o vídeo como candidato para ser convertido para letterbox, ao invés de widescreen. Ahhh… sim… note que o operador de comparação (<=) tem menor precedência que o de divisão (/)… Por isso a expressão intera faz o que queremos!

O problema é que nem a divisão de \frac{width}{height}, nem a divisão \frac{14}{9} são exatas. Bem… a primeira até pode resultar num valor exato, mas o quociente constante de \frac{14}{9} não! Isso poderia causar um “falso positivo” no caso de teste por “menor que” e um “falso negativo” no caso do teste de igualdade… Mas, aqui temos um outro problema: Imagine que o vídeo tenha resolução de 1400×901. Isso dará um aspecto de \frac{1400}{901}, necessariamente, que é muito próximo de \frac{14}{9}.

Falei “necessariamente” porque para obter a manor fração possível podemos usar o “máximo divisor comum” para fazer a divisão do dividendo e divisor (inteiras), obtendo os menores valores. Por exemplo:

#include <stdio.h>
#include <stdlib.h>

static unsigned int gcd( unsigned int, unsigned int );

int main( int argc, char *argv[] )
{
  unsigned int w, h, d;

  if ( argc != 3 )
  {
    fputs( "ERROR\n", stderr );
    return EXIT_FAILURE;
  }

  w = atoi( argv[1] );
  h = atoi( argv[2] );

  d = gcd( w, h );

  printf( "%u/%u -> %u/%u\n",
          w, h, w / d, h / d );

  return EXIT_SUCCESS;
}

unsigned int gcd( unsigned int a, unsigned int b )
{
  unsigned int t;

  if ( a < b )
  {
    t = a;
    a = b;
    b = t;
  }

  if ( ! b )
    return 1;

  while ( b )
  {
    t = b;
    b = a % b;
    a = t;
  }

  return a;
}

Compilando e testando:

$ cc -o test test.c
$ ./test 1280 720
1280/720 -> 16/9
$ ./test 1400 901
1400/901 -> 1400/901
$ ./test 1400 920
1400/920 -> 35/23

Se formos comparar o numerador e o denominador isoladamente a coisa vai complicar um bocado. Note que mesmo numa resolução maluca diferente (1400×920) o aspecto ainda tem o dividendo e divisor bem diferentes (\frac{35}{23}).

De volta ao ensino fundamental:

Como sabemos que, por exemplo, \frac{2}{3} é maior que \frac{3}{5}? A solução que você aprendeu lá no ensino fundamental é converter ambas as frações ao mesmo denominador. Daí podemos comparar diretamente os numeradores, certo? Basta multiplicar ambos o numerador e o denominador de ambas as frações pelo denominador da outra fração:

\displaystyle\frac{2}{3}\,\texttt{e}\,\frac{3}{5}\rightarrow\frac{2\cdot5}{3\cdot5}\,\texttt{e}\,\frac{3\cdot3}{5\cdot3}\rightarrow\frac{10}{15}\,\texttt{e}\,\frac{9}{15}

Agora, com ambas as frações “normalizadas”, podemos comparar diretamente os numeradores. No caso, 10 é maior que 9 e, por isso, \frac{2}{3} é maior que \frac{3}{5}

Vocẽ reparou que não precisamos mais nos preocupar com os denominadores se simplesmente multiplicarmos um numerador pelo denominador da outra fração, né? Isso facilita um bocado a nossa rotina is_letterbox(), evitando o uso de ponto flutuante:

int is_letterbox( unsigned int width, unsigned int height )
{
  unsigned int n1, n2;

  /* Não precisamos calcular os denominadores! */
  n1 = width * 9;   // 9 de 14/9
  n2 = 14 * height; // height de width/height

  return n1 <= n2;
}

Comparando as duas rotinas:

Compilando com otimizações e gerando o código em assembly:

$ cc -O2 -mtune=native -S -masm=intel test.c

Olhem só isso, com ambos os códigos lado a lado:

is_letterbox_fp:                  is_letterbox_int:      
  mov edi, edi                      imul  esi, esi, 14
  mov esi, esi                      lea eax, [rdi+rdi*8]
  pxor  xmm0, xmm0                  cmp eax, esi
  xor eax, eax                      setbe al
  pxor  xmm1, xmm1                  movzx eax, al
  cvtsi2sdq xmm0, rdi               ret
  cvtsi2sdq xmm1, rsi
  divsd xmm0, xmm1
  movsd xmm1, [rip+.LC0]
  ucomisd xmm1, xmm0
  setnb al
  ret

.LC0:
  dq 14.0 / 9.0  ; funciona no NASM.

Pra começar, a segunda rotina não tem divisões (que são LERDAS!). Ainda, a segunda rotina, além de ser menor, é melhor otimizada pelo processador (fusões, reordenações, etc). Outra vantagem é que, na primeira rotina, se height for zero, ela nos causará problemas porque a divisão resultará em NaN, que sempre retornará falso na comparação em ponto flutuante! Mais uma: A comparação, na segunda rotina, é sempre exata!

O único problema da segunda rotina poderia ser o caso de ocorrência de overflows, mas, como as resoluções gráficas que usaremos não resultam em frações com componentes muito grandes (pelo menos, não o suficiente para, quando multiplicados por 9 e 14 causem overflows), então a rotina é bem segura. Se, mesmo assim, você quer evitar a possibilidade de overflows, mude os tipos de n1 e n2 para unsigned long int ou unsigned long long. É prudente colocar um ‘UL‘ depois das constantes ou fazer um casting das variáveis nas multiplicações:

;  int is_letterbox( unsigned int width, 
;                    unsigned int height )
;  {
;    unsigned long n1, n2;
;  
;    n1 = width * 9UL;
;    n2 = 14UL * height;
;  
;    return n1 <= n2;
;  }
is_letterbox:
  mov esi, esi         ; Apenas para zeras os 32 bits
  mov edi, edi         ; superiores dos argumentos.

  lea rax, [rsi*8]     ; RAX = height*8
  lea rdx, [rdi+rdi*8] ; RDX = width*9
  sub rax, rsi         ; RAX = height*7
  add rax, rax         ; RAX = height*14

  cmp rdx, rax
  setbe al
  movzx eax, al
  ret

Mesmo assim o código ainda é melhor que o equivalente em ponto flutuante.

Siga o conselho do Tio Fred: Fuja de ponto flutuante sempre que puder.

Como NÃO tentar criar um gerador de valores aleatórios…

Já falei sobre o LCG (Linear Congruential Generator) e outros algoritmos como o Mersene Twster e o XorShift128+. Todos eles são “pseudo” geradores de valores aleatórios, já que não é possível gerar valores aleatórios com métodos determinísticos. Aliás, gerar valores aleatórios é um troço bem difícil porque a maioria dos fenômenos físicos tendem a oferecer um padrão, dada uma quantidade razoável de amostras. Por exemplo, se você jogar a mesma moeda, num jogo de “cara ou coroa”, digamos, 1 milhão de vezes, observará que a quantidade de “caras” ou “coroas” é maior que a outra. Ou seja, você pode prever, com maior segurança (embora não exatamente) qual será o resultado da próxima jogada…

De fato, até onde sei, os únicos fenômenos aleatórios são ao nível quântico (variação de fuga de carga num semnicondutor aquecido, por exemplo… Ou o decaimento radioativo num intervalo de tempo pequeno o suficiente)… De fato, também já falei por aqui, um desses fenômenos é usado, hoje em dia, para disponibilizar um DTRNG (Digital True Random Number Generator), presente dentro de todos os processadores da família Intel (pelo menos nos Ivy Bridge em diante) e acessíveis através da instrução RDRAND.

Mas, note: Mesmo que o LCG e outros métodos não criem valores verdadeiramente aleatórios, num escopo restrito e onde não há preocupação com a “aleatoriedade absoluta”, eles podem ser encacrados como geradores aleatórios reais. Abaixo, mostro um programinha que cria uma imagem de 1024×1024 pixels (RGB) com base na função rand() da biblioteca de C:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// Scale 'n' to 24 bits (0xBBGGRR).
#define rand2rgb(n) \
  ( ( int ) ( ( (n) * 16777215.0 ) / RAND_MAX ) )

int main ( void )
{
  int x, y, r;

  srand( time(NULL) );

  fputs( "P6\n1024 1024\n255\n", stdout );

  for ( y = 0; y < 1024; y++ )
    for ( x = 0; x < 1024; x++ ) 
    { 
      r = rand2rgb(rand());
      fwrite( &r, 3, 1, stdout );  // assume little endian. 
    }
}

Compilando e executando:

$ cc -O2 -o random random.c
$ ./random > rand1.ppm

Obtemos um arquivo rand1.ppm assim:

Ruído gerado por cores “aleatórias” dos pixels.

Recentemente vi alguém, num forum, tentando implementar seu próprio PRNG na rotina mcc_rnd(), abaixo:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <limits.h>

// Scale n to 24 bit value (0xBBGGRR).
// OBS: Uses long double because 'long' is 64 bits long in x86-64 systems.
#define rand2rgb(n) \
  ( ( int ) ( ( (n) * 16777215.0L ) / LONG_MAX ) )

// using size_t instead of ptrdiff_t (same thing)
typedef size_t mcc_rnd_t;

static long mcc_rnd ( mcc_rnd_t *seed, long min, long max );

int main ( void )
{
  int x, y;
  int r;
  long seed = 1;

  fputs ( "P6\n1024 1024\n255\n", stdout );

  for ( y = 0; y < 1024; y++ )
    for ( x = 0; x < 1024; x++ ) 
    { 
      r = rand2rgb ( mcc_rnd ( &seed, 0, LONG_MAX ) ); 
      fwrite ( &r, 3, 1, stdout );
    }
}

long mcc_rnd ( mcc_rnd_t *seed, long min, long max )
{
  // Most of the time tmp will point to the SAME memory block
  // for every process. Notice that, probably, tmp is aligned,
  // meaning the lower bits aren't 'random' at all... 
  // and the upper bits isn't random as well...
  void *tmp = malloc ( 1 );

  long val = time ( NULL ) + ( size_t ) tmp;

  free ( tmp );

  // This isn't 'random' seed. Taking the argument address means 
  // it will be copied to the stack and the stack pointer isn't
  // that 'random'... The value on stack isn't that unpredictable
  // as you think!
  if ( ! seed )
    seed = ( mcc_rnd_t * ) ( &seed );

  // Well... this isn't random, again, isn't it?
  if ( ! *seed )
    *seed = 1;

  val %= *seed;
  val *= clock();

  // The seed, again, isn't random!
  *seed <<= 1; 

  // Saturation get even less random values! 
  return ( val > max ) ? max : ( val < min ? min : val );
}

Os comentários na função mcc_rnd() são meus. A rotina tem várias falhas:

  1. malloc() retorna um endereço no virtual address space, que tende a ser o mesmo para todos os processos;
  2. Dependendo da velocidade com que a rotina é chamada, time() e clock() retornam o mesmo valor;
  3. clock() não mede a quantidade de ciclos de clock do processador, mas tem uma granulidade muito maior. Mesmo que medisse, teríamos valores praticamente sequenciais;
  4. A semente (*seed), usada como divisor, é tudo menos aleatória;
  5. A “saturação”, no caso de valores “fora da faixa” torna a rotina não aleatória!

De fato, compilando e executando, obtemos uma imagem assim:

Opa! tyemos padrões previsíveis.

Além da maioria dos valores serem zero, os pixels coloridos seguem padrões específicos. Notem os muitos pixels vermelhos e os gradientes do verde para o azul nos “blocos” acima…

Bem… minha dica: Não tente inventar um PRNG na base da tentativa e erro ou assumindo que certos “fenômenos” são aleatórios…

Ponto fixo para ponto flutuante…

Num outro artigo (este aqui) mostrei como obter o valor representado pela estrutura de um tipo ponto flutuante (float) e vice-versa, mas não mostrei um código fonte de exemplo. Eis aqui um, mas com falhas, porque não leva em consideração valores muito grandes, muito pequenos, subnormais, infinitos e NaNs. Ele converte um valor em ponto fixo (limitado), onde a parte inteira i e a parte fracionária f do valor é fornecida à função. Como está a parte fracionária jamais poderia ser 0.01 ou 0.00034, etc. Apenas valores como: 3.14, 100.4003, etc podem ser fornecidos.

A rotina converte o valor f na componente fracionária e junta as duas numa grande variável de 64 bits, deslocando os bits até isolar o bit de mais alta ordem no bit 32 (o 1.0. implícito num float) e daí obtemos o valor em ponto flutuante. Deixo o código para a análise de vocês:

#include <stdio.h>
#include <assert.h>

// União com a estrutura de um float, para faciliar
// as coisas...
union fp_u {
  float f;

  struct {
    unsigned int m:23;
    unsigned int e:8;
    unsigned int s:1;
  } __attribute__((packed));
};

// Calcula o log, na base 10, de x.
static int log10_( unsigned int x )
{
  int n = 0;

  while ( x >= 10 )
  {
    x /= 10;
    n++;
  }

  return n;
}

// Calcula 10^x (x precisa ser menor que 10).
static unsigned long long pow10_( unsigned int x )
{
  unsigned long long n = 1;

  assert( x < 10 );

  while ( x-- )
    n *= 10;

  return n;
}

// Converte um ponto fixo no formato (i,f) para float.
// Essa rotina tem multiplas falhas:
//  A parte fracionária não permite valores como 0.01, por exemplo.
//  Não verifico por subnormais, nans e infinitos.
//  Não trato valores negativos.
float fixed2float( unsigned int i, unsigned int f )
{
  unsigned int m;
  unsigned long long n;
  int e = 0;
  union fp_u fp = { .f = 0.0 }; // OK... 0.0f é apenas 0x00000000.

  // 0.0 é um caso especial.
  if ( i != 0 && f != 0 )
  {
    // Precisamos transformat a parte fracionária na represetação binária
    // do número fixo. Poderíamos fazer isso com ponto fluautuante, mas o objetivo
    // aqui é evitar ponto flutuante, não é?
    //
    // A mágica aqui é que na representação de número fixo (iiii.ffff) a parte fracionária
    // é (f / 2³²), mas temos apenas f inteiro. Então, coloco f na parte inteira (multiplicando por
    // 2³², e divido pela grandeza do valor (10^(log(f)+1)).
    //
    // Isso garante que m terá apeans 32 bits de tamanho final e todo o cálculo é inteiro (nada de
    // ponto flutuante aqui).
    m = ( ( unsigned long long )f << 32 ) / ( pow10_( log10_( f ) + 1 ) );

    // Agora que temos a representação fracionária correta, basta colocar os dois valores
    // nas posições corretas. A parte inteira são os 32 bits superiores, a parte fracionária, os
    // 32 bits inferiores.
    n = m + ( ( unsigned long long )i << 32 );

    // Um float normalizado sempre tem 1 implícito. Assim, se qualquer um dos 31 bits supeiores for 1, 
    // então temos que fazer deslocamentos para a direita, senão, para a esquerda, até que os 32 bits
    // superiores contenha apenas 1.
    //
    // OBS: ~0ULL << 33 == 0xfffffffe00000000ULL,
    //      ~0ULL << 32 == 0xffffffff00000000ULL e, é claro,
    //       1ULL << 32 = 0x0000000100000000ULL.
    if ( n & ( ~0ULL << 33 ) )
    {
      // Ao deslocar para a direita, incrementa o expoente.
      // OBS: Dá pra melhorar isso contanto a quantidade de 'leading zeros' e deslocando a
      //      quantidade necessária de uma vez só...
      while ( ( n & ( ~0ULL << 32 ) ) != ( 1ULL << 32 ) )
      {
        n >>= 1;
        e++;
      }
    }
    else
    {
      // Ao descolar para a esquerda, decrementa o expoente.
      // OBS: Dá pra melhorar isso contanto a quantidade de 'leading zeros' e deslocando a
      //      quantidade necessária de uma vez só...
      while ( ! ( n & ( 1ULL << 32 ) ) )
      {
        n <<= 1;
        e--;
      }
    }

    fp.m = n >> 9;  // Joga fora os 9 bits inferiores (precisams de apenas 23 bits).
    fp.e = 127+e;   // Calcula o expoente do float.
    fp.s = 0; // não lido com sinal, por enquanto!

    // Acerta o arredondamento. Precisamos arredondar a "mantissa" se o 9º bit for 1.
    if ( n & 0x100 )
      fp.m++;

    // Se a "mantissa" chegou a zero, temos que delocar mais um bit, decrementando o expoente.
    if ( ! fp.m )
      fp.e--;
  }

  return fp.f;
}

// Um pequeno teste.
int main( void )
{
  printf( "%.20f = %.20f\n", 3.1416f, fixed2float( 3, 1416 ) );
}

Compilando e rodando:

$ cc -o test test.c
$ ./test
3.14159989356994628906 = 3.14159989356994628906

.

Mifu…

Sabem aqueles grupos onde o pessoal costuma fazer perguntas de “problemas” de cursos de graduação de computação? Do tipo “faça um programa que some 5 números…”? Well… Eu acho isso enfadonho! Mas, recentemente, topei com um que, onde “me fudi” com uma resposta irônica de minha parte (e peço desculpas, porque a solução é mesmo interessante para o aprendiz!)…

O problema era esse ai mesmo que citei acima:

“Faça um programa que leia 5 valores inteiros e apresente a soma desses valores usando apenas uma variável

A parte em negrito é a que torna o problema interessante… Normalmente o estudante usa uma variável para ler os valores e outra para acumulá-los:

/* test1.c */
#include <stdio.h>

void main ( void )
{
  int sum, n;

  scanf ( "%d", &sum );
  scanf ( "%d", &n ); sum += n;
  scanf ( "%d", &n ); sum += n;
  scanf ( "%d", &n ); sum += n;
  scanf ( "%d", &n ); sum += n;}

  printf( "%d\n", sum );
}

Claro, isso usa duas variáveis, o que invalida a solução. Poderíamos usar um array:

/* test2.c */
#include <stdio.h>

void main ( void )
{
  int n[2];

  scanf ( "%d", &n[0] );
  scanf ( "%d", &n[1] ); n[0] += n[1];
  scanf ( "%d", &n[1] ); n[0] += n[1];
  scanf ( "%d", &n[1] ); n[0] += n[1];
  scanf ( "%d", &n[1] ); n[0] += n[1];

  printf( "%d\n", n[0] );
}

O programa acima é, essencialmente, a mesma coisa que o anterior e o estudante pode pensar que está usando apenas uma “variável”. No entanto, um arry é uma “sequência” de variáveis (a tradução direta da palavra “array” significa “sequência”), ou seja, está usando duas variáveis. Isso invalida a parte “usar apenas uma variável” do problema. Eis a solução engenhosa:

/* test3.c */
#include <stdio.h>

void main ( void )
{
  long long n;  // n tem 64 bits de tamanho!

  scanf ( "%d", ( int * )&n );

  scanf ( "%d", ( int *)&n + 1 ); 
  *( int * )&n += *( int * )(&n + 1);

  scanf ( "%d", ( int *)&n + 1 ); 
  *( int * )&n += *( int * )(&n + 1);

  scanf ( "%d", ( int *)&n + 1 ); 
  *( int * )&n += *( int * )(&n + 1);

  scanf ( "%d", ( int *)&n + 1 ); 
  *( int * )&n += *( int * )(&n + 1);

  printf( "%d\n", *( int * )&n );
}

Ao fazer o casting ( int * ) usamos a aritmética de ponteiros envolvendo o tipo int, que tem 32 bits de tamanho. Assim, ( int * )&n + 1 é uma expressão que é resolvida para os 32 bits superiores de n, enquanto ( int * )&n aponta para os 32 bits inferiores. Ao derreferenciar esse ponteiro, acessamos o seu valor, na pilha.

Em essência, estamos mexendo com um array, mas usando apenas uma variável!!! Isso atende a ambas os requisitos do problema, mas causa um problema de performance… Quando declaramos uma variável local qualquer o compilador, se tiver sido chamado com opções de otimização, tenta usar os registradores da CPU ao invés da pilha. Ao usar o operador & (endereço-de) forçamos o compilador a manter a variável na pilha. Como os scanf precisam do endereço da variável lida, os três programas, acima, geral praticamente o mesmíssimo código… Assim, cada chamda a scanf e cada acumulação usará um ou mais acessos à memória que, mesmo estando no cache L1d, poderá adicionar 1 ciclo de clock adicional às instruções (resolução do endereço efetivo).

A rotina abaixo, mais tradicional, tem o potencial de ser bem mais rápida e gastar menos espaço no cache L1i:

/* test4.c 
   Compile com: 
     gcc -O2 -fno-stack-protector -w -o test4 test4.c 
*/
#include <stdio.h>

void main ( void )
{
  int sum, n, count;

  count = 5;
  sum = 0;
  while (count--)
  {
    scanf ( "%d", &n );
    sum += n;
  }

  printf( "%d\n", sum );
}

Posso mostrar que count e sum não são mantidos na memória, mas em registradores. Apenas n não é otimizado dessa forma por causa do operador &:

; test4.asm - equivale a test4.c.
bits 64
section .rodata
scanf_fmt:  db `%d`
printf_fmt: db `%d\n`

section .text
global main
main:
  push rbp
  push rbx

  xor  ebp,ebp  ; sum=0
  mov  ebx,5    ; count=5

  sub  rsp,24   ; reserva espaço na pilha, alinhando-a

.loop:
  lea  rsi,[rsp+12]  ; n está alocado na pilha em rsp+12.
  xor  eax,eax       ; scanf() precisa disso (porquê?)
                     ; EAX, aqui, NÃO é stdin!
  mov  edi,scanf_fmt
  call __isoc99_scanf

  add  ebp,[rsp+12]  ; sum += n;
  sub  ebx,1         ; --count;
  jne  .loop         ; if (count != 0) goto loop.

  mov  edx,ebp
  mov  esi,printf_fmt
  mov  edi,1         ; stdout
  xor  eax,eax       ; printf() precisa disso (porquê?)
  call __printf_chk

  add  rsp,24   ; dealoca espaço previamente alocado na pilha.

  pop  rbx
  pop  rbp
  ret

Se não fosse pela necessidade do scanf, n também seria mantido em um registrador (R12D, por exemplo). Mas, note que temos apenas UM acesso à memória por iteração do loop (em add ebp,[rsp+12], mas, é claro, scanf também acessa memória!). No caso de test3.c temos 9 indireções, além das causadas pelo scanf, no caso de test4.c, 5.

Ahhh… lea rsi,[rsp+12] não é uma indireção, mas somente o cálculo do endereço efetivo…

Embora o exercício seja interessante, não passa de curiosidade acadêmica…

O bom, o mal e o feio

Neste artigo vou tentar mostrar 3 coisas: O algoritmo mais óbvio raramente é o melhor possível, um bom algoritmo pode ser sempre melhorado e que essa melhoria nem sempre vem sem custos. A ideia é usar um algoritmo custoso, em termos de processamento e, para isso, aqui vão os famosos testes para determinação se um valor inteiro positivo é ou não um número primo.

Nesses testes usarei o tipo unsigned int porque ele tem 32 bits de tamanho em ambas as plataformas Linux (e na SysV ABI) e Windows, tanto para o modo de operação do processador, x86-64 e i386. É bom reparar nisso, porque a determinação desse tipo de parentesco (se é primo!) em valores inteiros com mais bits pode ser muito demorada. Por exemplo, quando lidamos com o esquema de criptografia assimétrica, é comum o uso de chaves de 2048 bits de tamanho, o que torna uma das otimizações inviável (veremos adiante).

O feio

O algoritmo clássico é varrer toda a faixa de valores menores que o valor sob teste e verificar se existe algum divisor exato. Algo assim:

// Método tradicional: Tenta dividir todos os valores 
// inferiores a 'n' e conta quantos o dividem...
NOINLINE
int ugly_is_prime ( unsigned int n )
{
  unsigned int count, i;

  // caso especial.
  if ( !n ) 
    return 0;

  count = 1;
  for ( i = 2; i < n; i++ )
    if ( ! ( n % i ) )
      count++;

  return (count == 1);
}

Por que chamei essa rotina de “o feio”? Ela testa por todos os valores inferiores a n, começando por 2, se n for diferente de zero! Suponha que n=4294967291. Esse valor é primo, mas para a rotina determinar isso terá que percorrer todos os valores entre 2 e 4294967290. Veremos que, para isso, a rotina gasta cerca de 47 bilhões de ciclos de clock!

O mau

É claro que não precisamos verificar todos os valores! Por exemplo, 2 é o único valor par que também é primo. De cara podemos testar se n é par e diferente de 2. Se for, com certeza não é primo. Isso nos deixa os valores ímpares para teste, o que diminui pela metade o esforço computacional.

Outra melhoria que podemos fazer é a de que, de acordo com o Princípio fundamental da aritmética, todo número pode ser fatorado em números primos… Por exemplo, 9=3\cdot3, 10=2\cdot5, 76=2\cdot\cdot2\cdot19 e assim por diante (voltaremos a isso no algoritmo “bom”). Ainda, “lembre-se que a ordem dos fatores não altera o produto”. Tome 10 como exemplo… Ele pode ser escrito como 10=2\cdot5=5\cdot2.

Não é evidente, mas graças a isso só precisamos testar por divisores ímpares menores ou iguais a raiz quadrada de n:

// Método "ruim" (mas, melhor que o acima).
// Faz a mesma coisa, mas só verifica os valores inferiores a 'n' 
// que sejam ímpares e menores ou iquais a raiz quadrada de 'n'.
int bad_is_prime ( unsigned int n )
{
  unsigned int sr;
  unsigned int i;

  // casos especiais.
  switch ( n )
  {
    case 1:
    case 2:
      return 1;
  };

  // O único par primo é 2.
  if ( ! ( n & 1 ) )
    return 0;

  // Obtém o valor ímpar menor ou igual
  // à raíz quadrada de n.
  sr = sqrt(n);
  if ( ! ( sr & 1 ) ) 
    sr--;

  // Só precisamos testar os valores ímpares
  for ( i = 3; i <= sr; i += 2 )
    if ( ! ( n % i ) )
      return 0;

  return 1;
}

Aqui tomo um atalho para os 2 primeiros primos (1 e 2) para simplificar a rotina.

Mas, porque essa rotina é ? O problema aqui é que, mesmo eliminando testes, estamos testando muitos valores mais de uma vez… Por exemplo, se n não for divisível por 3 ele também não será para nenhum múltiplo de 3, como 9, 12, 15, 18, 21, 24, 27… No entanto, estamos testando esses valores também, se n for grande!

O bom

A solução, é claro, é eliminar todos os valores múltiplos dos anteriores, ou seja, usar apenas divisores primos. O problema é óbvio: Ora, bolas! Como é que eu vou fazer uma rotina que testa se um valor é primo usando como critério o teste se um divisor também o é?!

O método de Erastótenes (276 AEC ~ 194 AEC) mantém uma sequência onde marcamos os valores primos, percorrendo a lista do início ao fim, e verificando se os próximos valores são múltiplos. Se forem, marcamos como “não primo”. Na próxima interação, o próximo valor é, com certeza, primo e repetimos o processo, até que todo o array contenha apenas os primos “não marcados”. Isso exige que toda a faixa de valores seja colocada no array, mas os valores usados como base de testes podem ser limitados à valores menores ou iguais a raíz quadrada de n, como fizemos anteriormente.

O ponto negativo desse método, é claro, é o consumo excessivo de memória para conter todo o array, mas, e se mantivéssemos apenas um array com os valores primos ímpares inferiores ou iguais à raiz quadrada da faixa do tipo unsigned int? Isso não consumirá tanta memória quanto você pode imaginar: Em primeiro lugar, no pior caso, teríamos um array com \frac{\sqrt{2^{32}-1}}{2}\approx32767 itens. Se cada item tem 4 bytes de tamanho, aproximadamente 128 KiB serão usados para o array. No caso mais prático, vários valores desses 32767 são eliminados, porque não serão primos.

O programinha abaixo, que usa o bad_is_prime(), acima, monta esse array:

// gentable.c
//
// Cria tabela contendo os inteiros, primos, menores ou iguais ao
// valor ímpar mais pŕoximo (porém menor) que a raiz quadrada de 2³²-1.
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

extern int bad_is_prime ( unsigned int );

#define ITEMS_PER_LINE 16

void main ( int argc, char *argv[] )
{
  FILE *f;
  unsigned int i, sr;
  unsigned int count = 0, cnt2 = 0;

  if ( ! *++argv )
  {
    fputs ( "Usage: gentable <cfile>\n", stderr );
    exit ( 1 );
  }

  if ( ! ( f = fopen ( *argv, "wt" ) ) )
  {
    fprintf ( stderr, 
              "ERROR opening file %s for writing.\n", 
              *argv );
    exit ( 1 );
  }

  fputs ( "// *** criado por gentable.c ***\n\n"
          "// tabela começa de 3...\n"
          "const unsigned int primes_table[] = {\n", f );

  // ~0U é o máximo valor de unsigned int.
  // Obtém a raiz quadrada ímpar...
  sr = sqrt(~0U);
  if (!(sr & 1)) 
    sr--;

  // Só precisamos dos inteiros ímpares e primos.
  for ( i = 3; i <= sr; i += 2 )
    if ( bad_is_prime ( i ) )
    {
      if ( count ) 
      {
        fputc( ',', f );

        if (++cnt2 == ITEMS_PER_LINE)
        {
          fputc ( '\n', f );
          cnt2 = 0;
        }
        else
          fputc ( ' ', f );
      }

      fprintf ( f, "%uU", i );
      count++;
    }

  fprintf ( f, "\n};\n\n"
            "// tabela ocupa %zu bytes.\n"
            "const unsigned int primes_table_items = %u;\n",
            sizeof ( unsigned int ) * count, 
            count );

  fclose ( f );
}

Basta compilar e executar o programinha para gerar um arquivo texto contendo o código fonte de uma tabela com apenas as raízes quadradas desejadas! No caso, gerarei table.c que, como pode ser observado, contém 6541 valores, ocupando apenas 25.5 KiB:

// *** criado por gentable.c ***

// tabela começa de 3...
const unsigned int primes_table[] = {
3U, 5U, 7U, 11U, 13U, 17U, 19U, 23U,
29U, 31U, 37U, 41U, 43U, 47U, 53U, 59U, 
61U, 67U, 71U, 73U, 79U, 83U, 89U, 97U,
...
... 65447U, 65449U, 65479U, 65497U, 65519U, 65521U
}; 

// tabela ocupa 26164 bytes.
const unsigned int primes_table_items = 6541;

Para essa tabela de raízes o processamento é rápido. Considere o caso de valores de 64 bits. A máxima raiz quadrada é 4294967295. Como verificar se um valor é primo é relativamente lento, pré-calcular toda a tabela pode demorar um bocado… Existe ainda o fator espaço, que explico abaixo.

Usando essa tabela, a nossa rotina de teste good_is_prime fica assim:

extern const unsigned int primes_table[];
extern const unsigned int primes_table_items;

// Método "melhor": Faz quase a mesma coisa que acima, mas
// usa uma tabela pré calculada.
int good_is_prime ( unsigned int n )
{
  unsigned int i, sr, *p;

  switch ( n )
  {
    case 1:
    case 2:
      return 1;
  }

  if ( ! ( n & 1 ) )
    return 0;

  sr = sqrt(n);

  for ( i = 0, p = (unsigned int *)primes_table; 
        i < primes_table_items && *p <= sr; 
        i++, p++ )
  {
    if ( ! ( n % *p ) )
      return 0;
  }

  return 1;
}

Ela é bem parecida com bad_is_prime, mas usa uma tabela contendo apenas valores primos para os testes!

Nem sempre usar tabelas pré-calculadas é uma boa ideia

Em primeiro lugar, poderíamos usar uma tabela pré-calculada com todos os valores primos na faixa e a rotina final poderia ser uma simples busca binária, daí a performance da rotina seria proporcional a \log_2{n}, mas é claro que a tabela final seria enorme. Daí, escolhi armazenar apenas os valores que serão usados para verificar o divisor, eliminando os múltiplos…

A quantidade de valores primos, dado um valor máximo, pode ser aproximado (de longe!) por:

\displaystyle \Pi(x)\sim\frac{x}{\log{x}-1}

Assim, se quisermos montar uma tabela com todos os primos inferiores ou iguais à raiz quadrada do maior valor possível para um tipo, como no exemplo de unsigned int (que assume o valor máximo de 2^{32}-1), temos que o maior valor para a tabela é r_{max}=\lfloor\sqrt{2^{32}-1}\rfloor e, portanto, usando a aproximação acima, teríamos uma tabela com cerca de 17171 itens. Cada item cabe em 2 bytes (unsigned short), já que r_{max}=65535. O que nos daria um array de, no máximo, 33.5 KiB. Se usarmos unsigned int (como fiz) o array estimado ficaria com 67 KiB.

Observe o contraste com o que encontramos através de gentable.c… Temos menos da metade dessa quantidade! Mas isso não é garantido! Temos que assumir que a aproximação acima é uma quantidade média e usarmos como critério para estimar o custo da rotina.

Considere agora, 64 bits. Usando o mesmo critério, r_{max}=\lfloor\sqrt{2^{64}-1}\rfloor=4294967295, que cabe num unsigned int. Daí \Pi(r_{max})\approx497508081 itens com 4 bytes cada, ou seja, quase 2 GiB de array. Repare que dobramos o número de bits da faixa sob teste, mas o tamanho da tabela aumentou quase de 29 mil vezes.

Isso é esperado porque a estimativa é exponencial (proporcional a \frac{1}{\log{x}-1})… Então, não recomendo o método good_is_prime para valores com mais de 32 bits.

O quão bom é o “bom”?

No pior caso bad_is_prime é 10 milhões de vezes mais rápida que ugly_is_prime e good_is_prime é quase 4 vezes mais rápida que bad_is_prime (o que faz dela quase 40 milhões de vezes mais rápida que ugle_is_prime).

Estou informando esses valores, mas os medi, usando meu “contador de ciclos do homem pobre” em alguns casos:

  • Melhor caso: Um valor não primo pequeno;
  • Pior caso: Um valor primo muito grande: 4294967291
  • Um valor não primo grande cujos fatores sejam primos grandes: 144795733, por exemplo, pode ser fatorado em 12343 e 11731.

Em alguns casos, no “melhor caso” e somente às vezes, bad_is_prime ganha de good_is_prime, mas isso ocorre, provavelmente, por causa dos efeitos do cache L1 (que só tem 32 KiB!). Isso, talvez, possa ser compensado fazendo um prefetch da tabela antes de executar a rotina… Esse é especialmente o problema quando o valor não é primo! Quando a rotina precisa varrer grande parte do o array os efeitos do cache L1 entram em ação e ela fica um pouco mais lenta que a “má” rotina…

Prática horrível: Uso de MAX_PATH

Eis algumas coisas que considero muito erradas quando alguém tenta criar alguns tipos corriqueiros de rotinas: usar a pilha inescrupulosamente, e achar que o tamanho máximo de um pathname é de 260 caracteres. Ambas costumam levar a problemas sérios… Eis um exemplo:

void faz_algo(void)
{
  char fname[MAX_PATH+1];
  char *p;

  getcwd(fname, MAX_PATH);
  p = fname + strlen(fname) - 1;
  if (*p == '//') *p = '\0';
  strcat(fname, "/o-nome-maluco-de-algum-arquivo.txt");
  ...
}

Aqui, além do buffer apontado por fname ser alocado na pilha, gastando mais que 260 bytes, já que a pilha tem que estar, necessariamente, alinhada e ainda existe a red-zone, esse buffer tem tamanho máximo que parece bom… O que acontece se o diretório corrente for algo assim: /home/frederico/files/br/com/xpto/alpha/beta/theta/caralhos-me-levem/abobora/? Well… Ainda é provável que isso caiba num buffer de 260 caracteres, mas saiba que em certos ambientes o path máximo pode ter de 32 KiB! É o caso, por exemplo, do Windows (veja aqui).

A ideia do esquema acima pressupõe que tenhamos espaço suficiente e evita alocação dinâmica e, de bandeja, livra-se do buffer quando a função terminar. É interessante, mas não é assim que as linguagens que possuem o tipo primitivo para strings funciona. Eis uma maneira de fazer, com alocação, que não fica tão estranha assim:

char *fname, *cwd;

// get_current_directory() é uma extensão da glibc que
// aloca espaço da string pra você!
cwd = get_current_directory();

// asprintf() aloca espaço suficiente para todos os
// caracteres da string resultante!
asprintf(&fname,
         "%s/o-nome-maluco-de-algum-arquivo.txt",
         cwd);

// Não precisamos mais do bloco alocado em cwd.
free(cwd);
...
// Livra-se do buffer apontado por fname...
free(fname);

Ficou mais complicado? Que nada! O único detalhe é que você terá que se livrar dos buffers manualmente.

Os mais atentos verão que não faço nenhuma checagem de falha. Isso porque a alocação de pedaços pequenos de memória quase nunca falham. E free não liga se o ponteiro contém NULL… Mas, mesmo que isso seja importante pra você, alguns simples if‘s resolvem tudo.

Para lidar com strings dinâmicamente você pode, até mesmo, usar realocação, como, por exemplo:

char *_str_concat(char **pp, const char * const *s)
{
  char *tmp;
  size_t size;

  size = strlen(s);
  if (*pp) size += strlen(*pp);

  if (tmp = realloc(*pp, size+1))
  {
    *pp = tmp;
    strcat(*pp, s);
  }
  return *pp;
}

Neste caso, se o ponteiro para sua string for NULL, realloc alocará o espaço desejado (do mesmo tamanho da string apontada por s), caso contrário, os dois tamanhos são adicionados, o espaço da string original é aumentado e depois a concatenação é feita. O uso da dupla indireção pode te confundir, mas o uso é bem simples:

char *s = NULL;

_str_concat(&s, "Fred");
_str_concat(&s, " é ");
_str_concat(&s, "bom nisso!");
...
// Não se esqueça de livrar-se da string
free(s);

No fim das contas, o ponteiro s apontará para um buffer dinamicamente alocado contendo apenas a string “Fred é bom nisso\0“.

Você pode até mesmo garantir que sua string seja declarada corretamente usando uma macro:

#define DECLARE_STRING(name) char *name = NULL

Para usar, é claro, basta fazer DECLARE_STRING(s);, ao invés de escrever char *s = NULL;. E no caso do free, você pode fazer:

// Chame da mesma forma que chamou _str_concat.
void DELETE_STRING(char **pp)
{ if (*pp) free(*pp); *pp = NULL; }

Uma operação de assinalamento numa string já alocada não poderia ser mais simples:

char *_str_assign(char **pp, char *s)
{
  DELETE_STRING(pp);
  _str_concat(pp, s);
}

Todas as outras funções de string.h continuam funcionando do mesmo jeito e, ao contrário da crença “popular”, alocações de pequenos blocos são realmente rápidas. Você ainda pode criar rotinas de housekeeping, assim:

void alguma_funcao(void)
{
  jmp_buf jb;
  DECLARE_STRING(s);

  if (setjmp(jb))
  {
    // Coloque aqui seus DELETE_STRING()...
    DELETE_STRING(&s);
    return;
  }

  ...

  // Antes de sair chame longjmp()! :)
  // Não esqueça de fazer isso tb se tiver algum "return"
  // no código, fora daquele bloco de setjmp(), ai em cima...
  longjmp(jb, 1);
}

É claro, você pode criar MACROS para simplificar a coisa, tipo: DECLARE_HOUSEKEEPING, que declarará jb, HOUSEKEEPING que declarará o if (setjmp(jb)) e EXIT_FUNCTION, que chamará o longjmp(jb, 1):

// Coloque essas 3 num header!
#define DECLARE_HOUSEKEEPING jmp_buf _jb
#define HOUSEKEEPING if (setjmp(_jb))
#define EXIT_FUNCTION longjmp(_jb,1)

void alguma_funcao(void)
{
  DECLARE_HOUSEKEEPING;
  DECLARE_STRING(s);

  HOUSEKEEPING
  {
    // Coloque aqui seus DELETE_STRING()...
    DELETE_STRING(&s);
    return;
  }

  ...

  // Antes de sair chame EXIT_FUNCTION :)
  // Não esqueça de fazer isso tb se tiver algum "return"
  // no código, fora daquele bloco de setjmp(), ai em cima...
  EXIT_FUNCTION;
}

Ficou mais fácil de ler?