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…

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… ;)

Wallpaper animado usando VLC no Linux

Eis uma click-bait (e uma dica) proceis… Essa dica veio do “CLI magic”, mas eu a alterei para uso do VLC.

$ cvlc "file.mkv" --no-audio --loop --aspect-ratio '16/9' \
   --drawable-xid $(xwininfo -int -name 'Desktop' | \
                    awk '/Window id:/{print $4;}')

Repare que estou usando o CVLC para não usar a interface de janela do VLC… A linha acima retira o áudio, coloca o vídeo em loop infinito e ajusta o aspect-ratio para 16/9 (Widescreen). Outra coisa que ele faz é selecionar o “canvas” onde o vídeo é apresentado. Nesse caso, a janela do DESKTOP… Resultado (click-bait!):

É só colocar isso num script de inicialização e, voilà! :)

É conveniente que o vídeo esteja no mesmo aspect-ratio da tela, mesmo que tentemos forçar a barra para ajustar isso com a opção –aspect-ratio, o VLC tende a fazer um resize que não pode ficar muito agradável…

PS: Ainda tyem algumas coisas por fazer. Por exemplo, impedir que o título do arquivo seja mostrado por alguns segundos no início (adicione a opção --no-video-title-show); Zerar (colocar tudo preto) o fundo, caso o seu fundo não seja preto. Fundo não-preto causa blending do vídeo com o fundo, alterando as cores (não me preocupai em achar isso nas opções do VLC, ainda)… Deixo esses detalhes proceis.

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

.