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.

Convertendo uma string para ponto flutuante

Não vou mostrar código aqui porque seria um bem grande e com alguns macetes e testes de condições especiais (como NaNs, overflows, underflows). Para dar uma olhada no código atual da glibc, obtenha o código-fonte (no Debian/Ubuntu: apt source libc6) e veja o arquivo stdlib/strtod_l.c.

O problema em questão aqui é como transformar uma string contendo um valor decimal que tenha componente fracionário para ponto flutuante. Para tanto, usarei o tipo float como base (porque tem menos bits significativos – me permitirá escrever menos!). Considere o valor de PI com algumas “casas decimais” de “precisão”: 3.14159

Converter a parte inteira é trivial: Basta realizar múltiplas divisões por 2 e os restos nos darão os bits à partir do MSB. Neste caso, 3 é 0b11. Mas, e quanto a parte fracionária? O motivo que dividirmos por 2 o valor inteiro até obtermos zero como quociente é que cada posição do componente inteiro é um múltiplo de 2, em binário. Ou seja, 3 é 1\cdot2^1+1\cdot2^0. No caso da parte fracionária, os fatores multiplicadores dos bits são da ordem de 2^{-n}. O expoente negativo indica que o valor de cada bit deverá ser dividido pela potência de dois. A operação inversa, então, é multiplicar o valor decimal fracionário por 2 e obter o valor inteiro que, necessariamente, será 0 ou 1. No caso, 0.14159\cdot2=0.28318. Obtivemos 0, na parte inteira, e 0.28318 na parte fracionária… Esse 0, inteiro, é o MSB da parte fracionária… No segundo passo, multiplicamos apenas a parte fracionária resultante (0.28318) por 2. A parte inteira do resultado é o próximo bit e repetimos tudo para os próximos bits (eis todos os cálculos para 25 bits):

0.14159 * 2 = 0.28318 (MSB fracionario = 0)
0.28318 * 2 = 0.56636 (0)
0.56636 * 2 = 1.13272 (1)
0.13272 * 2 = 0.26544 (0)
0.26544 * 2 = 0.53088 (0)
0.53088 * 2 = 1.06176 (1)
0.06176 * 2 = 0.12352 (0)
0.12352 * 2 = 0.24704 (0)
0.24704 * 2 = 0.49408 (0)
0.49408 * 2 = 0.98816 (0)
0.98816 * 2 = 1.97632 (1)
0.97632 * 2 = 1.95264 (1)
0.95264 * 2 = 1.90528 (1)
0.09528 * 2 = 1.81056 (1)
0.81056 * 2 = 1.62112 (1)
0.62112 * 2 = 1.24224 (1)
0.24224 * 2 = 0.48448 (0)
0.48448 * 2 = 0.96896 (0)
0.96896 * 2 = 1.93792 (1)
0.93792 * 2 = 1.87584 (1)
0.87584 * 2 = 1.75168 (1)
0.75168 * 2 = 1.50336 (1)
0.50336 * 2 = 1.00672 (1)
0.00672 * 2 = 0.01344 (0)
0.01344 * 2 = 0.02688 (0)
...

Assim o valor final de 3.14159 é, em binário, com 27 bits de precisão é 0b11.0010010000111111001111100.

Como estamos lidando com float e já temos 2 bits na parte inteira, precisamos fazer essas repetições apenas 23 vezes. Uma vez a mais para determinar se precisaremos arredondar o valor obtido. Se a precisão de um float é de 24 bits, é necessário obter 25 para usar esse bit extra para arredondamento.

Da sequência acima, obtemos o valor 0b11.00100100001111110011111. Agora, precisamos normalizar esse valor, deixando apenas 1 bit setado na parte inteira. Isso é feito deslocando o “ponto binário” 1 bit para a esquerda. O valor obtido é 0b1.10010010000111111001111_1. Como o bit inferior, além da precisão do float, ou seja, dos 23 bits da parte fracionária, é 1, arredonda-se o valor da “mantissa” (de 23 bits, em vermelho) para cima, acrescentando 1, obtendo: 0b1.10010010000111111010000.

Agora, note que, na equação que nos dá o valor contido na estrutura de um float:

\displaystyle v=\left(-1\right)^s\cdot\left(1+\frac{f}{2^{23}}\right)\cdot2^e

O valor inteiro f nada mais é do que a parte fracionária, binária, deslocada 23 bits para a esquerda, ignorando o valor inteiro normalizado, ou seja, f é 0b100_1001_0000_1111_1101_0000 (0x490fd0 ou 4788176). Agora temos f=4788176 e e=1:

\displaystyle v=\left(-1\right)^0\cdot\left(1+\frac{4788176}{8388608}\right)\cdot2^{1}=3.141590118408203125

Que é a representação mais próxima de 3.14159 que podemos conseguir com um float.

Neste ponto você deve estar pensando que é meio esquisito que uma rotina de conversão de string para ponto flutuante precise usar ponto flutuante. Isso é, mais ou menos, a história do ovo e da galinha (embora o “ovo” tenha surgido primeiro! A galinha evoluiu de outras espécies que botavam ovos!), mas a rotina geralmente usa aritmética de múltipla precisão e é até possível usar BCD (Binary Coded Decimal) para tanto… De novo, a rotina ____STR2F_INTERNAL() em strtod_l.c é bastante intrincada porque leva em contra aqueles underflows e overflows, mas também a notação xxx.xxx...E+/-eee, mas o princípio fundamental é esse ai em cima!

Transformando uma divisão inteira em multiplicação

É muito comum, na maioria das arquiteturas, que as instruções de divisão (quando existem) sejam bem mais lentas que as instruções de multiplicação (quando existem). Assim, é também comum que bons compiladores transformem divisões inteiras por divisores “invariantes” (constantes) em multiplicações:

\displaystyle \frac{n}{d} = \frac{1}{d}\cdot n

O problema aqui é que \frac{1}{d} resultará num valor em ponto flutuante e queremos evitar isso… Para resolver esse problema, vajamos o que acontece quanto calculamos \frac{1}{d} em ponto flutuante (usando precisão simples). Lembre-se que na estrutura de um float temos precisão de 24 bits e o valor final corresponde à equação:

\displaystyle v=\left(-1\right)^s\cdot\left(1+\frac{f}{2^{23}}\right)\cdot2^e

Vamos assumir, como exemplo, que d=13 e teremos:

\displaystyle \frac{1}{13}\approx\left(-1\right)^0\cdot\left(1+\frac{1935833}{2^{23}} \right )\cdot2^{-4}=\left(2^{23}+1935833\right )\cdot2^{-27}

Esse valor de f pode ser obtido manualmente ou através desse pequeno código:

#include <stdio.h>

struct fp_s {
  unsigned int f:23;
  int e:8;
  int s:1;
};

int main( void )
{
  float f = 1.0f / 13.0f;
  struct fp_s *p = ( struct fp_s *)&f;

  printf( "1 / 13 = %.27f = (-1)^%d * ( 1 + %u/2^23 ) * 2^%d\n",
    f, p->s, p->f, p->e - 127 );
}

Quando multiplicamos um valor n por \frac{1}{13} em ponto flutuante, na verdade estaremos multiplicando o valor 10324441 (ou seja: 2^{23}+1935833) e depois realizamos um shift para a direita de 27 bits. O único problema aqui é que float não tem precisão suficiente para a multiplicação de valores arbitrários para n, já que possui apenas 24 bits de precisão. Precisamos de 32 bits!

O GCC faz algo similar ao que fizemos acima, mas usa precisão de 32 bits (no caso de ‘int’s). O valor \frac{1}{13} é, em binário, 0b0.000_1001110110001001110110001001110_110…, onde apenas os 32 bits em vermelho nos interessam… Eles são arredondados para 0b1001110110001001110110001001111, ou 1321528399, porque o bit que segue o LSB em vermelho é 1. Como se pode observar, teremos que realizar o deslocamento de 34 bits para a direita no resultado da multiplicação desses valores inteiros (31 bits da precisão “fracionária” mais os 3 bits zerados “depois da vírgula” – exatamente como no caso do float, acima). Acontece que essa multiplicação resultatá num valor de 64 bits sempre. Então, ao invés de fazer um shift de 34 bits para a direita, só precisamos obter os 32 bits superiores e os deslocarmos para a direita em 2 bits. Ou seja:

; unsigned int f( unsigned int a ) { return a / 13; }
f:
  mov  eax,edi
  mov  edx,1321528399
  mul  edx             ; edx:eax = eax * edx
  mov  eax,edx         ; Obtemos os 32 bits superiores...
  shr  eax,2           ; ...e deslocamos em 2 bits para
                       ; direita.
  ret

A “divisão” com inteiros sinalizados é similar (só leva o sinal de a em consideração!).

; int f( int a ) { return a / 13; }
f:
  mov  eax,edi
  mov  edx,1321528399

  sar  edi,1           ; EDI = -1 se 'a' < 0 ou,...
                       ; EDI = 0, se 'a' >= 0.

  imul  edx            ; edx:eax = eax * edx
  mov  eax,edx         ; Obtemos os 32 bits superiores...
  sar  eax,2           ; ...e deslocamos em 2 bits para
                       ; direita, replicando o sinal.

  sub  eax,edi         ; Precismos subtrair 1, se o valor
                       ; de 'a' for negativo.
  ret

O motivo da subtração da unidade se a < 0 é que o quociente de divisão de valores positivos é q=\lfloor\frac{n}{d}\rfloor, enquanto para a <= 0 temos q=\lceil\frac{n}{d}\rceil.

Tenha em mente que aqui só expliquei como o compilador calcula esses valores “mágicos” ao substituir divisões com divisores invariantes por multiplicações. Se estiver interessado em toda a matemática por trás dessa “mágica”, recomendo a leitura deste paper (aqui). Cuidado que o troço é meio complicado… ;)

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

.

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…

O problema das pizzas – Ou, o problema da ignorância da matemática elementar!?

Bem… não foi surpresa nenhuma pra mim perceber que algumas pessoas não conseguem aplicar conhecimentos fundamentais no dia a dia. Eis um problema (bem simples, que um estudando ensino fundamental deveria saber resolver, e um adulto também!) que propus, recentemente:

“Você resolveu comprar pizza para a família toda. Uma pizzaria oferece pizzas “médias” com raio de 30 cm e “grandes” com raio de 45 cm, ao preço de R$ 40,00 e R$ 50,00 cada uma, respectivamente. Você compraria duas pizzas médias ou uma grande?”

O detalhe é que o preço é irrelevante no problema. O que se quer saber é qual das duas configurações vai te dar mais pizza para comer: Duas médias ou uma grande?

Como disse, não é surpresa ver que a maioria das pessoas escolhe as duas pizzas médias, achando que terá mais pizza do que uma grande. Infelizmente a geometria plana elementar desmente isso, já que a “quantidade” de pizza é medida pela área da circunferência… Quanto a isso, obviamente o problema considera que uma pizza seja isto:

ISSO é uma pizza!

Fazendo um cálculo rápido, qualquer um pode observar que uma pizza grande, do problema proposto, é 12,5% maior que duas pizzas médias (área da pizza média. 900\pi, área da pizza grande, 2025\pi)! Isso levanta a questão: Será que estou perdendo dinheiro toda vez que compro pizzas? Quando é vantajoso comprar duas médias ao invés de uma grande? Well…

Considere r como sendo o raio da pizza média e R como sendo o raio da pizza grande. Valerá a pena comprar duas pizzas médias ao invés de uma grande se:

\displaystyle 2\pi r^2 > \pi R^2

Ou seja, \frac{R^2}{r^2} < 2 \Rightarrow R < r\sqrt{2} (estou ignorando o valor negativo que satisfaz a inequação por motivos óbvios!) . Isso quer dizer que o raio R tem que exceder menos que 41% do radio r para ser vantajoso comprar duas médias… No exemplo do problema, se a pizza grande tiver menos que 42 cm (30\sqrt{2}\approx42,4), vale à pena comprar duas médias…

Ahhh… Esse tipo de racionalização também deveria ser óbvio para qualquer pessoa que passou por uma escola…

Arredondamento. Qual é o “correto”?

O último artigo gerou alguma controvérsia… Alguns leitores perguntam: “Qual é o jeito certo?” e não se satisfazem com a resposta de que ambos os métodos, de truncamento e de arredondamento para baixo, estão certos, dependendo da aplicação. Aqui quero provocá-los mostrando algo similar…

Já mostrei antes que operações em ponto flutuante dependem de algum método de arredondamento bem definido. O padrão IEEE 754 define quatro deles:

  • No sentido de +\infty;
  • No sentido de -\infty;
  • No sentido de 0;
  • O valor “mais próximo”;

Ao calcular f=1.0/10.0, o valor 0.1, que não pode ser representado numa variável do tipo ponto flutuante, já que, em binário, será uma dízima periódica:

\displaystyle (0.1)_{10} = (0.0110011001100110011001100...)_{2}

Portanto, precisa ser arredondado. Se estivermos lidando com o tipo float esse valor será (1.1001100110011001100110?)_{2}\cdot 2^{-2}, onde ‘?’ é o algarismo impreciso… Ou seja, se ele for zero, termos um valor levemente inferior a 0.1, se ele for 1, teremos um valor levemente superior a 0.1. Qual deles devemos escolher?

Por default, o padrão IEEE 754 escolhe o método do “mais próximo”. Repare:

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

extern float _mydiv(float, float);

void main(void)
{
  float f;
  unsigned int d = 0x3dcccccc;
  float *g = (float *)&d;
  unsigned int *p = (unsigned int *)&f;

  f = _mydiv(1.0f, 10.0f);

  printf("O valor '%.2f' é codificado como '0x%08X'.\n\n", 
    f, *p);

  printf("Mas, na verdade:\n"
         "0x%08X -> %.18f\n"
         "0x%08X -> %.18f\n",
    *p, f,
    d, *g);
}

O motivo de criar a função _mydiv é para evitar que o compilador otimize a operação de divisão entre duas constantes e não faça divisão alguma. Neste caso, o arredondamento será feito pelo compilador, não pelo nosso programa. A função _mydiv é, simplesmente:

/* mydiv.c */
float _mydiv(float a, float b) 
{ return a/b; }

Compilando e linkando tudo e depois executando, temos:

$ gcc -o test test.c mydiv.c
$ ./test
O valor '0.10' é codificado como '0x3DCCCCCD'.

Mas, na verdade:
0x3DCCCCCD -> 0.100000001490116119
0x3DCCCCCC -> 0.099999994039535522

É fácil perceber que o primeiro caso está mais próximo de 0.1 do que o segundo e é este que o processador escolherá (como isso é feito não interessa agora!). Mas, isso não significa que os outros métodos de arredondamento não sejam “corretos”. De fato, se mandarmos o processador arredondar de outra maneira ele o fará, Eis um exemplo:

/* test.c */
#include <stdio.h>
#include <fenv.h>

extern float _mydiv(float, float);

void main(void)
{
  float f;
  int oldround;

  oldround = fegetround();
  fesetround(FE_TOWARDZERO);
  f = _mydiv(1.0f, 10.0f);
  fesetround(oldround);

  printf("%.18f\n", f);
}

Ao compilar e executar esse código (será preciso especificar a libm com -lm, no gcc) você verá que a divisão de 1 por 10 será arredondada para o MENOR valor, ou melhor, em direção ao zero. Esse será o valor mais distante do verdadeiro, neste caso. Troque 1.0f por -1.0f e veja o que acontece…

Depois, modifique o método para FE_DOWNWARD (arredondamenteo no sentido de -\infty) e teste os dois casos (com 1.0f e -1.0f)… Note que o truncamento citado no artigo anterior é equivalente a FE_TOWARDZERO…

O que quero dizer aqui é que, tanto na programação quanto na matemática, às vezes existe mais de um jeito “correto” de fazer alguma coisa. Depende apenas da aplicação.