Um cientista em minha vida

Eu li lá no blog do Kentaro que o meme da semana era "Um Cientista em minha vida", onde deveríamos falar sobre algum cientista que fez a diferença pra você. Como eu adoro constrained writing, resolvi participar (na verdade, eu adoro constrained anything, por isso que vivo criando programas em uma linha, programas que rodam em computadores de 8 bits, e assim por diante).

Eu já falo de cientistas aqui todo o tempo. Olhando no histórico, eu já falei do Knuth, do Erdös, do prof. Routo, do prof. Henrique, e de vários outros. Em comum, todos eles foram cientistas que eu conheci depois de adulto. Achei apropriado então que eu falasse de um cientista que fez a diferença quando eu era criança, e pra isso vamos ter que rebobinar até a década de 80.

Se você perguntar pra alguém sobre revistas de computador na década de 80, invariavelmente irá ouvir sobre a Micro Sistemas ("a primeira revista brasileira sobre microcomputadores"). A Micro Sistemas era muito legal, mas o que eu gostava mesmo era de outra revista, menos conhecida, chamada Microhobby.

A diferença da Micro Sistemas pra Microhobby era mais ou menos a diferença de Informática pra Computação. Na primeira, nós ficávamos encantados com as notícias da maravilhosa terra além da reserva de mercado (onde aprendíamos que a Apple planejava lançar um novo computador chamado McIntosh, que vinha com um periférico estranho e esquisito chamado mouse), enquanto que na segunda aprendíamos a calcular geodésicas e a usar o método de Bolzano para achar raízes de uma equação.

Mas o diferencial mesmo da Microhobby eram as colunas escritas pelo Renato da Silva Oliveira. Uma googlada rápida revela que o Renato é formado em Física, trabalhou nos planetários de São Paulo, Campinas, Vitória e Tatuí, e atualmente trabalha em uma empresa que vende planetários infláveis (how cool is that?!). Mas é claro que eu não sabia disso na época, o que eu sabia era que ele contava historinhas!

Foi lendo as historinhas do Renato que eu descobri que era possível escrever sobre ciência e computação, com clareza e bom humor. Pena que isso ainda não é muito difundido, a julgar pela quantidade de crianças que ainda acham que ciência é uma coisa chata :(

As historinhas que ele escrevia sempre tinham o mesmo formato: um certo sr. Nabor Rosenthal, em suas viagens pelo mundo, deparava-se com alguma situação que sugeria uma análise matemática (os tópicos eram os mais variados e iam de teoria dos grafos até contatos com extraterrestres). Depois de ponderar sobre o problema sem conseguir resolvê-lo, o Nabor tomava uma dose do raro Suco de Ramanujan, que o colocava num transe que ampliava suas capacidades analíticas, e conseguia solucionar o problema.

Mas a coluna sempre acabava antes que o Nabor mostrasse qual a solução! Ao invés disso, o leitor tinha um mês pra conseguir resolver o problema, e só no mês seguinte a solução era apresentada. Na década de 80 ainda não tinha spoj, então as colunas do Renato eram o que bombava pra quem gostava de puzzles computacionais.

Um dos puzzles apresentados foi o segundo puzzle mais difícil da minha vida, eu levei mais de dez anos pra conseguir resolver. Em uma das Microhobby, o Nabor entrou em transe após tomar o Suco de Ramanujan, e durante o transe ele sonhou "com um método para calcular o número pi, usando apenas o gerador de números aleatórios de seu micro" (essa era a época onde o micro mais avançado era o TK-82C, com 2kb de RAM).

Na época eu pensei muito e não consegui solucionar, achando que ia precisar de alguma matemática que eu ainda não tinha aprendido. Eu nunca consegui achar a revista seguinte com a solução, tive que passar pelo primário, pelo colégio técnico, e só no meio da faculdade é que caiu a ficha (e eu percebi que poderia ter solucionado ainda no primário, se tivesse insistido o suficiente :)

O truque é o seguinte: você vai fazer N experimentos, cada um consistindo no sorteio de dois números aleatórios escolhidos uniformemente entre 0 e 1. Se soma dos quadrados dos números for menor ou igual a 1, incremente um contador (digamos, M). Ao final dos experimentos, pi=4*M/N. O script abaixo implementa esse algoritmo:

Script em python para calcular pi usando números aleatórios

O funcionamento é bem simples e baseia-se na figura ao lado. Você começa inscrevendo um quarto de círculo num quadrado de lado 1. Os dois números que você sorteia a cada iteração podem ser interpretados como um ponto dentro do quadrado, e o teste feito é equivalente a testar se o ponto está dentro do círculo ou não. Como a distribuição dos pontos é uniforme, espera-se que a razão M/N seja igual à razão entre as áreas da figura. A área do quadrado é 1, a área do círculo é pi*r². Como o raio é unitário, então a área do quarto de círculo é pi/4. Isolando pi, chega-se em pi=4*M/N, QED.

A pergunta que deve ser feita ao encontrar qualquer algoritmo novo é: qual é sua complexidade? Infelizmente, esse método aleatório é bem ruim. No fundo, o que estamos fazendo é aproximar pi por uma fração, cujo denominador é N. Então a precisão máxima que podemos obter é 1/N, e se você quer calcular n dígitos de pi, esse método converge, no melhor caso, em O(10ⁿ), e na prática em bem menos que isso, porque os seus geradores de números aleatórios não são perfeitamente uniformes.

Eu nunca soube qual o método que o Nabor usou pra calcular o pi. Como ele tinha o Suco de Ramanujan e eu não, espero que tenha sido um método melhor que o meu :)

Programa de milhagem

Olhando no Google Analytics, eu descobri que alguém chegou aqui no blog procurando por "transformar kilometros em arquivo binarios". Se você é essa pessoa, desculpe, mas eu não entendi sua dúvida. Se você não é essa pessoa, puxe uma cadeira pra ver como até uma pergunta sem sentido pode ser desenvolvida num tema interessante :)

Uma coisa que me incomoda toda vez que venho pra Califórnia é ter que lidar com milhas e libras. Eu cresci com o sistema métrico: se você me disser que a distância de São Paulo pra Florianópolis é de 700km, eu sei o que isso significa. Agora, se você me disser que a distância de San Francisco pra Mountain View é de 100 milhas, minha intuição falha, e eu vou precisar converter pra km pra poder ter noção da distância.

Como esperado, converter mentalmente de milhas pra km é algo que faço todo o tempo por aqui. Existem várias maneiras de converter de cabeça, mas a filosofia Ricbit dita que, de várias maneiras equivalentes, o correto é escolher a mais bizarra! Sendo assim, vou mostrar a conversão utilizando números de Fibonacci.

Todo mundo conhece os números de Fibonacci, eles estão em todo lugar. Em minha época de estudante, uma das minhas diversões secretas era entrar escondido no andar superior da biblioteca do IME só pra ler a Fibonacci Quartely. Os leitores assíduos da revista conhecem um monte de fatos curiosos sobre os Fibonacci, e eu vou usar dois deles aqui.

O primeiro é que os números da seqüência de Fibonacci podem ser calculados diretamente, sem precisar fazer toda a recursão F(n+2)=F(n+1)+F(n). Pra isso, basta calcular qual o inteiro mais próximo da expressão abaixo:




Na expressão acima, o phi é a conhecida razão áurea. A demostração dessa fórmula é elementar e fácil de encontrar na web.

O segundo é que qualquer número pode ser escrito como uma soma de números distintos de Fibonacci. Por exemplo, 100 pode escrito como 89+8+3. Daí, se você enfileirar os números de Fibonacci, e atribuir a cada número 1 se ele é usado na soma, e 0 se não é, você pode atribuir a qualquer inteiro uma string de zeros e uns que funciona como uma espécie de base binária alternativa (o povo chama isso de base de Fibonacci).

Fazendo o processo com o número 100, chegamos em 10000101000. Essa base tem algumas propriedades curiosas, por exemplo, ela não é bijetora (de fato, você pode escrever 100 de outras maneiras, como 1100101000). Uma base numérica que tem representação múltipla tem utilidades bastante curiosas em design de circuitos elétricos (mas isso é uma história pra outro dia :).

Além disso, cada número possui pelo menos uma representação onde não há nenhuma seqüência com dois uns consecutivos (pela própria definição de Fibonacci, se houver dois uns em algum ponto, você pode apagá-los e trocar por um único 1 na posição seguinte). Esse fato é explorado em alguns tipos de sinalização, para fazer detecção de erro: se você receber dois uns seguidos, certamente recebeu um erro.

Mas qual a relação disso com milhas e quilômetros? É simples, para converter de milhas para km, basta fazer um shift de Fibonacci!

Uma milha equivale a 1.609344 km. A razão áurea é 1.61803399. Os dois números, apesar de não serem relacionados, são muito parecidos, e esse é o truque que vamos usar pra converter. Na base binária tradicional, um shift para a esquerda é equivalente a multiplicar por dois; na base de Fibonacci, um shift para a esquerda equivale a multiplicar pela razão áurea. Então, se você tiver um valor em milhas na base de Fibonacci, um shift irá transformar o valor para o equivalente em quilômetros.

Vamos conferir: 100 é 10000101000. Com um zero extra no final, fica 100001010000, que é 144+13+5=162. Se, ao invés disso, você converter diretamente, teria 160.9km, ou seja, o método realmente aproxima muito bem a conversão! O gráfico abaixo mostra a porcentagem do erro do método em relação ao ideal, e mais abaixo está o programa que converte a milhagem para quilometragem e gera o gráfico:

Script que plota o gráfico acima, em python

Como pode ser visto, o método esquenta bem rápido, pra distâncias superiores a 10 milhas o erro é inferior a 5%, e acima disso praticamente desaparece. Em comparação, o método naive de aproximar por 1.5 (multiplicar de cabeça por 3 e dividir por 2), tem um erro constante de mais ou menos 8%.

Primos aleatórios

Dia desses a Alice me perguntou se era possível criar um gerador de números aleatórios que só retornasse números primos. Eu respondi que sim, mas que provavelmente ela não iria gostar da resposta:

int random_prime(int n) {
 int x;
 do {
   x = random(n);
 } while (!is_prime(x));
 return x;
}

Eu sabia que o que ela queria na verdade era uma fórmula bonitinha; então, como esperado, ela não gostou :) Mas a verdade é que esse algoritmo é bem melhor que as alternativas!

Antes de mostrar porque isso é verdade, precisamos formalizar um pouco o problema. É claro que não existem algoritmos que geram números aleatórios: se você quiser aleatoriedade real, precisa pegar alguma fonte física, como o decaimento radiativo. Assumindo então que existe uma fonte física que gera uma distribuiçao uniforme sobre algum intervalo, para criar o algoritmo que retorna números primos aleatórios, basta criar uma função bijetora que leve naturais para primos. Ou seja, uma função que, para um dado um número n, retorne o n-ésimo primo.

O problema é que não existe nenhuma fórmula fechada que calcule isso de maneira eficiente. Você pode calcular alguma constante irracional que resolva o problema, no estilo da constante de Mills, só que mais cedo ou mais tarde a precisão vai te limitar. Você pode calcular o n-ésimo primo com base em alguma outra distribuição, como a função de Möbius, mas aí você só está empurrando o problema com a barriga, porque a outra função é tão difícil de calcular quanto a original.

Uma maneira sem as desvantagens acima é usar o teorema de Wilson pra chegar na seguinte fórmula:




Mas mesmo essa fórmula ainda está longe do ideal, primeiro porque você vai ter que lidar com números enormes nela (pra n=10 os valores intermediários ficam tão grandes que estouram o limite do que cabe num float), segundo porque, mesmo que você use uma lib para long floats, a complexidade é O(2ⁿ), ou seja, mais lento que os programadores do Duke Nukem Forever. Se ainda assim você quiser testar, minha implementação em python é a abaixo:

Implementação em python da fórmula acima

Sendo assim, quão melhor era a implementação original por tentativa e erro? Pra avaliar isso, precisamos calcular a complexidade daquele algoritmo. Não é difícil ver que a complexidade do algoritmo como um todo é a complexidade do is_prime() multiplicado pelo valor esperado do número de iterações do loop.

Se você estiver trabalhando numa faixa pequena de primos, pode tabelar todos os primos no intervalo e fazer um is_prime() que seja O(1), mas aí também não tem necessidade da tentativa e erro, você pode indexar seu número aleatório direto na tabela. O caso legal é quando você não pode tabelar, nesse caso você pode implementar o is_prime() usando, por exemplo, o algoritmo AKS, cuja complexidade é O((log n)^10.5).

O que resta então é calcular o valor esperado do loop. Lembrando que E[x]=sum(x*p(x)), o que precisamos é calcular qual é a probabilidade de ter uma iteração, duas iterações, e assim por diante. Ora, o teorema dos números primos nos garante que a quantidade de números primos menores que n é assintoticamente igual a n/log(n), então a chance de um número ser primo, num conjunto com n elementos, é 1/log(n). Vamos chamar isso de "p" só pra ficar mais fácil, e o complemento disso é q=1-p, ou seja, a chance de um número não ser primo.

Vejamos então: pra você acertar o primo de primeira, a chance é p. Se você acertar o primo na segunda, a chance é pq. Na terceira, é pq², na quarta pq³ e assim por diante. Então o valor esperado é:

X = 1p + 2pq + 3pq² + 4pq³ + ...
X = p (1 + 2q + 3q² + 4q³ + ....)

Quem tem prática com a transformada z sabe calcular isso de cabeça, mas dá pra calcular também só com matemática elementar. Se você isolar q na soma, fica com:

X = p (1 + q(2 + 3q + 4q² + ....))

Agora você tira da cartola y=1+q+q²+q³+... e substitui:

X = p (1 + q(2 + 3q + 4q² + ....))
X = p (1 + q(y + 1 + 2q + 3q² + ....))
X = p (1 + q(y + X/p)) = p + pqy + pXq/p = p(1+qy) + Xq
X - Xq = p (1 + qy)
X (1-q) = Xp = p (1 + qy)
X = 1 + qy

Mas y é só a soma de uma PG, e isso nós sabemos que vale y=1/(1-q)=1/p. Então:

X = 1 + q/p = (p+q)/p = 1/p

Como p=1/log(n), então o valor esperado que nós queríamos é tão somente X=log n (vocês também não se impressionam quando tudo simplifica no final?)

É claro que eu não iria resistir à tentação de implementar uma simulação pra ver se o valor bate mesmo. A nossa fórmula diz que, para a faixa de 10 milhões de números, o valor esperado tem que ser da ordem de log(10⁷)=16.1. A simulação abaixo retorna 15.2, bem próximo do valor que foi predito.

Simulação monte carlo do valor esperado, em C++

No fim das contas, a complexidade do algoritmo com tentativa e erro é apenas O(log n), se você tiver um tabelão de primos. Na prática, esse é o método usado por todos que precisam de primos aleatórios: a libgcrypt usada no gpg, por exemplo, utiliza esse método na função gen_prime(), com vários truques pra tornar o teste de primalidade bem rápido.

Python one-liners são Turing-complete

Quem programa em C há décadas normalmente não se dá conta de quão ilegíveis são as expressões mais comuns da linguagem. Quando eu era um garoto recém-saído do BASIC, eu lembro de ter me assustado com coisas básicas como for(i=0; i<10; i++). Mas isso é idiossincrasia do C, outras linguagens não sofrem disso, como o Python.

Python foi planejada para ser legível. Os programadores mais experientes citam o Zen of Python, que dita que "belo é melhor que feio", e "legibilidade conta". De fato, é até difícil escrever código ilegível em Python. Mas é claro que difícil não é impossível, e se o dr. Ian Malcolm fosse um programador, ele certamente diria que "obfuscation finds a way."

Aconteceu comigo semana passada: eu olhava alguns exercícios sobre listas para iniciantes, e notei que, embora eles fossem de fato muito simples, ficariam bem mais divertidos se eu tentasse resolvê-los usando apenas uma linha em cada. Abusando de programação funcional, eu consegui fazer os dez primeiros assim:

Soluções dos exercícios em uma linha de Python cada.

Depois de brincar algum tempo com one-liners, a pergunta que naturalmente se apresenta é: será que é possível fazer qualquer programa em uma linha de Python? A dificuldade vem do fato de que o Python diferencia statements de expressions, e você só pode ter um statement por linha. Em Python, statements incluem print, if, while, for e atribuições, ou seja, um one-liner só pode usar um único desses.

Então, colocando a pergunta de outra maneira: é possível demonstrar que um programa em Python com um único statement é Turing-complete? Existem dois caminhos pra demonstrar isso, o primeiro é construir um emulador para uma máquina de Turing universal em uma linha, o segundo é mostrar que é possível converter para uma linha de Python todos os programas possíveis de um sistema que seja Turing-complete, como o cálculo lambda, ou os tag-systems.

Eu resolvi abordar o problema com a filosofia Ricbit: se existem várias maneiras equivalentes de fazer alguma coisa, escolha a mais bizarra! Assim sendo, vou demonstrar que Python one-liners são Turing-complete através de redução ao Brainfuck (cuja universalidade já foi demonstrada várias vezes).

Vamos lá então: o estado de um programa em Brainfuck pode descrito em qualquer momento por uma quádrupla (mem, p, stdin, stdout), que são respectivamente a memória, o ponteiro, a entrada e saída. Vou implementar cada operação do Brainfuck como funções que recebem quádruplas e retornam quádruplas, descrevendo assim a transição de estado.

A operação mais simples é o ponto, que só adiciona o elemento apontado na saída:

dot = lambda mem, p, stdin, stdout: (mem, p, stdin, stdout+[mem[p]])

Para implementar a vírgula, eu preciso primeiro de alguma maneira de modificar um único elemento de uma lista. Se eu pudesse usar atribuições, bastaria algo do tipo mem[p]=value, mas como atribuições em Python são statements, preciso de uma função auxiliar. Além disso, eu preciso fazer o pop() do valor frontal da lista que guarda o stdin, o que me leva à outra auxiliar:

change = lambda mem, pos, value: [value if i==pos else a for i, a in enumerate(mem)]

get = lambda s: (s[0], s[1:]) if len(s) else (0,[])

comma = lambda mem, p, stdin, stdout: (lambda now, next: (change(mem, p, now), p, next, stdout))(*get(stdin))

Tendo a função change em mãos, fazer os comandos de mais e menos é simples:

plus = lambda mem, p, stdin, stdout: (change(mem, p, mem[p]+1), p, stdin, stdout)

minus = lambda mem, p, stdin, stdout: (change(mem, p, mem[p]-1), p, stdin, stdout)

Os comandos de esquerda e direita precisam tomar o cuidado de aumentar os limites da memória se necessário, a universalidade do Brainfuck requer uma fita infinita para os dois lados:

left = lambda mem, p, stdin, stdout: ([0]+mem if not p else mem, 0 if not p else p-1, stdin, stdout)

right = lambda mem, p, stdin, stdout: (mem+[0] if p==len(mem)-1 else mem, p+1, stdin, stdout)

Agora chegamos na parte complicada, que é o operador de loop. Como for e while são statements, e lambdas recursivos precisam de uma atribuição (fat = lambda x: 1 if x<=1 else x*fat(x-1)), então a única saída é apelar pra lazy evaluation, que no Python é implementada no módulo itertools. (Incluir o módulo itertools poderia tornar o programa um two-liner, mas felizmente é possível importar um módulo usando uma expression ao invés de um statement: a função __import__).

A solução para o operador de loop é criar uma lista infinita contendo [x, f(x), f(f(x)), f(f(f(x))), ...], onde cada f é uma aplicação do conteúdo do loop. Depois, para executar o loop, basta iterar nesta lista infinita, procurando o primeiro elemento onde o elemento apontado pelo ponteiro seja nulo. Precisamos então de uma função que calcule f^n e uma que gere a lista infinita:

composite = lambda f, n: lambda x: reduce(lambda a, b: b(*a), [f]*n, x)

infinite = lambda f, x: itertools.imap(lambda n: composite(f, n)(x), itertools.count(0))

Depois, basta criar um predicado que avalie quando o loop deve parar, e pegar o primeiro elemento da lista onde o predicado é verdadeiro:

predicate = lambda mem, p, stdin, stdout: mem[p] != 0

getfirst = lambda it: [i for i in itertools.islice(it, 1)][0]

loop = lambda f: lambda *x: getfirst(itertools.dropwhile(lambda x: predicate(*x), infinite(f,x)))

Tendo todos os comandos, só precisamos de uma função extra para encadeá-los, e depois, só para o programa não ficar grande demais, um shortcut que executa strings diretamente em Brainfuck:

chain = lambda f: lambda *x: reduce(lambda y, g: g(*y), f, x)

bf = {'+':plus, '-':minus, '.':dot, ',':comma, '<':left, '>':right}

run = lambda f: chain([bf[i] for i in f])

Feito! Agora é só fazer um script para parsear o Brainfuck original e gerar o one-liner. A título de ilustração, esse é o Hello World gerado pelo script:

print ''.join(chr(i) for i in ( (lambda itertools: (lambda change, get, chain, composite: (lambda comma, dot, plus, minus, left, right, infinite, predicate, getfirst: (lambda bf, loop: (lambda run: (chain ([run ("++++++++++"), loop (run ("<+++++++<" "++++++++++<" "+++<+>>>>-")), run ("<++.<" "+.+++++++..+++." "<++.>>" "+++++++++++++++" ".<.+++." "------.--------" ".<+.<.")])) ([0],0,[],[]) )( (lambda f: chain([bf[i] for i in f])) ) )( ({'+':plus, '-':minus, '.':dot, ',':comma, '<':left, '>':right}), (lambda f: lambda *x: getfirst(itertools.dropwhile(lambda x: predicate(*x), infinite(f,x)))) ) )( (lambda mem,p,stdin,stdout: (lambda now,next: (change(mem,p,now),p,next,stdout))(*get(stdin))), (lambda mem,p,stdin,stdout: (mem,p,stdin,stdout+[mem[p]])), (lambda mem,p,stdin,stdout: (change(mem,p,mem[p]+1),p,stdin,stdout)), (lambda mem,p,stdin,stdout: (change(mem,p,mem[p]-1),p,stdin,stdout)), (lambda mem,p,stdin,stdout: ([0]+mem if not p else mem, 0 if not p else p-1, stdin, stdout)), (lambda mem,p,stdin,stdout: (mem+[0] if p==len(mem)-1 else mem, p+1, stdin, stdout)), (lambda f,x: itertools.imap(lambda n: composite(f,n)(x), itertools.count(0))), (lambda mem,p,stdin,stdout: mem[p] != 0), (lambda it: [i for i in itertools.islice(it, 1)][0]) ) )( (lambda mem,pos,value: [value if i==pos else a for i,a in enumerate(mem)]), (lambda s: (s[0],s[1:]) if len(s) else (0,[])), (lambda f: lambda *x: reduce(lambda y,g: g(*y), f, x)), (lambda f,n: lambda x: reduce(lambda a,b:b(*a),[f]*n,x)) ) )(__import__("itertools")) )[3])

QED, em uma única linha, como prometido! (Eu não prometi que seria uma linha pequena :)

O script que converte de Brainfuck para Python one-liner está abaixo, para quem quiser brincar:

Conversor para Python one-liner.

Como nota final, vale lembrar que só porque você pode escrever qualquer coisa em uma linha, não significa que você deve fazer isso. Legibilidade conta :)

Potências ótimas

Olhando no Google Analytics, eu descobri que alguém chegou aqui no blog procurando por "como implementar em c++ potências". Se você é essa pessoa, a resposta está abaixo. Se você não é essa pessoa, puxe uma cadeira que o papo é divertido :)


Calcular potências aproximadas em ponto flutuante é trivial, basta incluir a biblioteca <cmath> e usar a função pow, que internamente é implementada como exp(y*log(x)). Mas existem várias aplicações onde você precisa do valor exato da potência, como, por exemplo, durante a criptografia RSA. Nesses casos, uma primeira abordagem pode ser como no código abaixo:

int natural(int x, int n) {
  int result = 1;
  for (int i = 1; i <= n; i++)
    result *= x;
  return result;
}


Esse código funciona, mas existem maneiras mais espertas. Nós, que temos dez dedos, não estamos acostumados a pensar em binário. Mas se fôssemos como os golfinhos, que tem um cérebro maior que o nosso e só duas barbatanas, poderíamos visualizar o expoente em binário e fazer um código assim:

int binary(int x, int n) {
  if (n == 0) return 1;
  if (n == 1) return x;

  int half = binary(x, n/2);
  if (n & 1)
    return half*half*x;
  else
    return half*half;
}

A versão original era O(n), essa é O(log n), uma melhoria significativa. Mas a melhor notícia sobre esse algoritmo é que você não precisa implementá-lo, ele já está pronto na biblioteca padrão do C++. Basta incluir <ext/numeric> e usar a função power. A função da STL ainda recebe como argumento opcional um functor de multiplicação, então você pode implementar com ela exponenciação modular, ou até mesmo potências de matrizes (ela funciona mesmo que sua multiplicação não seja comutativa).

Por outro lado, esse algoritmo é prático, mas não é ótimo. Em alguns casos, existem maneiras mais rápidas de calcular a potência, o primeiro caso onde isso acontece é para n¹⁵. O algoritmo binário precisa de seis multiplicações para resolver o problema, mas existem soluções com apenas cinco:



A melhor maneira de descrever as multiplicações necessárias para calcular uma potência é através de uma addition chain. Uma addition chain é uma espécie de generalização da seqüência de Fibonacci: enquanto no Fibonacci o próximo elemento é a soma dos dois imediatamente precedentes, numa addition chain o próximo elemento é a soma de dois anteriores quaisquer, ou até mesmo a soma de um elemento anterior com ele mesmo (com a restrição de que a seqüência precisa ser estritamente crescente).

Pela regra de formação dá pra perceber que, ao contrário da seqüência de Fibonacci, existem inúmeras addition chains. E melhor ainda, você pode associar uma addition chain com a seqüência de expoentes gerados no cálculo de uma potência. Para o exemplo de n=15 acima, as duas addition chains correspondentes são [1 2 3 6 7 14 15] e [1 2 4 5 10 15].

Achar uma seqüência ótima para calcular potências, então, é equivalente a achar uma addition chain de tamanho mínimo terminada no número que queremos. Infelizmente, eu tenho duas más notícias: o Erdös mostrou que a addition chain ótima não cresce mais lentamente que O(log n), então assintoticamente ela não é melhor que o método binário; e pior, o cálculo da addition chain ótima é NP-Completo. Abaixo eu implementei a addition chain ótima em C++ (usando brute force, então está bem lento):

Implementação da addition chain ótima em C++

É interessante também dar uma olhada nos casos onde o binário perde. No gráfico abaixo, a linha vermelha é o método binário, a linha azul é a addition chain ótima, e a linha verde é log₂n:


Olhando o gráfico, um golfinho certamente perceberia que o desvio do método binário é proporcional à quantidade de dígitos 1 na representação binária do expoente (os picos são em 63, 127, 255, e assim por diante). A demonstração disso é bem simples e está no Seminumerical Algorithms, junto com várias heurísticas para aproximar a addition chain ótima.

Paranóia x Matemática

No último post eu falei sobre Criptografia, então agora, pra balancear, o tópico é Criptanálise. Semana passada, a polícia prendeu uma gangue que estava instalando os mini-notebooks Eee PC dentro de caixas automáticos, para roubar senhas dos usuários. O vídeo com a matéria pode ser visto abaixo:

Eu tenho certeza que um monte de gente deve ter visto a matéria e pensado "omfg nunca mais vou usar caixas automáticos kthxbye", mas na verdade, mesmo com o notebook lá dentro, não é tão fácil conseguir roubar a senha!

Se o seu banco for como o meu, você não digita a sua senha diretamente, ao invés disso, a máquina associa dois dígitos para cada botão, e você aperta os botões correspondentes à sua senha. Assim, se algum xereta estiver atrás de você olhando, ele não vai conseguir descobrir sua senha, e isso vale também se tiver um notebook dentro da máquina registrando o que você digita.

Então esse sistema é à prova de sniffers? Não! Um jeito de quebrar esse sistema é fazendo várias observações. Se o xereta te olhar uma única vez, ele não consegue descobrir sua senha, mas reduz bastante as possibilidades. Se a senha tiver quatro dígitos, uma única observação reduz de dez mil possibilidades para apenas 16. Se ele olhar uma segunda vez, pode ser que consiga informação suficiente para reduzir ainda mais o espaço, e, eventualmente, repetindo o processo, ele pode conseguir deduzir a senha.

Para conseguir reconstruir computacionalmente a senha, tudo que ele precisa fazer é resolver uma matriz de exact cover (na verdade outros métodos podem ser usados, mas eu sou preguiçoso adepto da orientação à objeto e do reuso de código pronto). Assuma n observações: para cada um dos 4 dígitos da senha há dez possibilidades, então você tem 40 linhas. Além disso, para cada observação, você precisa garantir que os quatro dígitos são consistentes, o que dá 4n colunas, e ainda mais 4 colunas extras para garantir que um único número estará associado a cada dígito da senha. No total, a matriz terá 40x(4n+4) elementos.

E quantas observações ele precisa? Isso não dá pra saber a priori, depende de como os números foram sorteados. Se a máquina repetir sempre a mesma distribuição, ele nunca vai conseguir deduzir a senha (mas também nem precisaria, pois aí ele só precisa apertar os mesmos botões que você). Por outro lado, se você tiver azar, pode ser que só duas observações sejam suficientes, como no caso abaixo:

observação 1:
senha CADA, botões A=(6 2) B=(5 8) C=(4 3) D=(1 9) E=(7 0)
observação 2:
senha CAAC, botões A=(1 2) B=(8 4) C=(6 3) D=(5 0) E=(7 9)

Nesse exemplo, dá pra ver claramente que a única senha consistente com os dados é 3216. Se o ladrão for levemente mais esperto, ele pode até perceber que não precisa fazer observações suficientes para que a senha seja única, basta que ele reduza o espaço de possibilidades para 3 ou menos (já que ele pode chutar 3 senhas antes da máquina bloquear o cartão).

Embora não seja possível calcular a priori quantas observações são necessárias, é perfeitamente possível calcular qual é o valor esperado dessa distribuição. Como eu sou preguiçoso adepto das simulações computacionais, ao invés de calcular as probabilidades na unha, eu escrevi uma simulação de Monte Carlo para calcular esse valor. O resultado é que, para uma senha de 4 dígitos, um ladrão que queira achar a senha exata precisa de 2.3 observações, enquanto que o ladrão esperto, que se contenta em reduzir pra 3 ou menos possibilidades, precisa de apenas 2.1 observações.

Código do simulador de Monte Carlo

A lição prática dessa análise é que, se você estiver desconfiado que o caixa tem um sniffer, não precisa ficar preocupado, desde que digite sua senha uma única vez. Se você precisar fazer uma outra operação, é melhor fazer em outra máquina.

A Meta-Assinatura

Como eu já disse antes, eu sou uma criatura que se empolga fácil. Ainda não tinha feito nem duas semanas que eu e o Fábio tínhamos entrado na Poli, e nós já estávamos procurando iniciação científica pra fazer. Depois de alguma procura, achamos uma legal: o Routo Terada estava procurando alunos pra estudar Criptologia.

O nosso medo inicial era que o Routo não quisesse aceitar dois alunos de primeiro ano, mas isso foi mais simples que esperávamos: "Ah, eu posso passar uma tarefa simples pra vocês. O Schneier acabou de publicar um algoritmo novo chamado Blowfish, vocês tem seis meses pra quebrar". É claro que não conseguimos quebrar o Blowfish, mas aprendemos um bocado no processo :)

Assinaturas digitais, por exemplo. O Isaac Newton, quando queria provar que algum manuscrito era dele, podia simplesmente assiná-lo com uma pena; mas o Stephen Hawking não pode fazer isso! Pra ele, o ideal são as assinaturas digitais. Para assinar digitalmente, você precisa de algum tipo de problema que seja difícil de resolver, mas que seja fácil de checar se foi resolvido (como a fatoração de números, ou o problema da sacola).

Um exemplo simples de como isso funciona me veio à mente algum tempo atrás, enquanto eu lia um livro do Hofstadter (se você não conhece o Hofstadter, tem uma entrevista dele para a rede Globo disponível online). Suponha que eu fiz uma grande descoberta e quero divulgar isso para o mundo:

O Ricardo sabe onde está o Bin Laden.

Embora tenha meu nome ali, qualquer um pode alterar e trocar o nome, então não tem como garantir que fui eu que escrevi:

O Wilerson sabe onde está o Bin Laden.

O método que eu bolei, e que na falta de nome melhor eu chamo de Meta-Assinatura, consiste em adicionar informação auto-referente à sua sentença:

O Ricardo afirma que sabe onde esta o Bin Laden, nesta sentenca com dezessete letras a, vinte e sete letras e, seis letras i, sete letras o, quatro letras u e uma letra x.

Confira que a contagem de letras está certinha. Dessa maneira, o Wilerson não pode trocar o nome na frase, pois se ele trocar, a contagem de letras vai mudar. Assim, a frase com meta-assinatura garante quem é o autor. Nesse método, contar as letras é muito simples, mas consertar a frase para o número de letras bater, é bem difícil (quer dizer, só com seis letras e algum esforço, até dá pra consertar a frase, mas se você usar o alfabeto inteiro na sua contagem, aí fica realmente complexo).

Para criar a frase com meta-assinatura, você não pode tentar procurar a solução por força bruta, porque demora demais. Uma solução mais rápida é criar uma função que conte as letras da sentença e troque os números correspondentes, e depois cruzar os dedos e torcer pro ponto fixo dessa função ser um atrator. O script em python abaixo faz isso, tomando o cuidado de detectar loops para não ficar preso:

Meta-Assinatura em python

Eu ainda não consegui assinar uma sentença usando todas as letras do alfabeto (ie, gerando um pangram), porque esse método não garante convergência. Se você conseguir, me avise :)

Erdös e os logaritmos

Essa foi uma semana triste para a ciência, com a morte de dois grandes nomes: Edward Lorenz (criador dos atratores de Lorenz e criador da expressão "efeito borboleta"), e John Wheeler (orientador do Feynman e criador das expressões "buraco negro" e "buraco de minhoca"). Sempre que morre um cientista famoso, eu fico pensando que perdi a oportunidade de conhecê-lo pessoalmente pra agradecer pelo seu trabalho. Mas alguns cientistas eu consegui conhecer em vida, um deles foi o Paul Erdös.

Em novembro de 94, o Erdös deu uma palestra na USP, e, sabendo que ele estaria lá, fui correndo assistir. Pra quem não conhece, o Erdös foi o segundo matemático mais prolífico da história, só o Euler publicou mais que ele (embora anedoticamente ele seja mais conhecido pela brincadeira dos números de Erdös). É claro que um estudante de primeiro ano, como eu era, não tinha a menor chance de entender os detalhes da palestra que ele deu. Na verdade, o que me deixou impressionado foi que, em um dado momento, ele demonstrou que o upper bound de uma função era O(log log log log n), e eu pensei comigo mesmo que um dia ainda ia encadear logaritmos como ele fazia.

O tempo passou e eu ainda não consegui encadear quatro logaritmos, mas outro dia eu consegui pelo menos dois! Foi quando eu estava otimizando um código, e o seguinte problema apareceu no meio de um inner loop: achar a menor potência de 10 que seja maior ou igual a um inteiro dado. A implementação simples é a abaixo, vamos assumir que os inteiros em questão sejam de 64 bits, e que f(0)=1 por convenção:


unsigned long long simple_power10(unsigned long long i) {
  unsigned long long current = 10000000000000000000ULL;
  while (true) {
    if (current <= i)
      return current + !current;
    current /= 10;
  }
}

Esse código é razoavelmente rápido, roda em O(log n). O ideal seria rodar em O(1), fazendo uma tabela com os valores pré-calculados. Porém, uma tabela assim é inviável na faixa de valores de 64 bits. Um caminho mais esperto é usar uma busca binária para achar o valor correto:

if (i < 100) {
  if (i < 10)
    return 1;
  else
    return 10; 
} else {
  if (i < 1000)
    return 100;
  else
    return 1000;
}

Essa idéia é bem melhor, mas o problema agora é escrevê-la. Para a faixa de 64 bits, os ifs aninhados ficam muito longos, e um cara distraído como eu certamente iria errar alguma coisa na implementação. Felizmente, existe uma solução: template metaprogramming!

Usualmente pensamos em template metaprogramming para fazer cálculos em tempo de compilação, mas ele pode ser usado nesse caso também, pra gerar o código da busca binária. E ainda ganhamos uma vantagem, o código pode ser usado para qualquer tipo, não ficando preso ao unsigned long long, como no primeiro caso. Para implementar, começamos fazendo um template para gerar potências de dez:

template<class T, const int n>
struct p10 {
  const static T value = T(10) * p10<T, n-1>::value;
};

template<class T>
struct p10<T, 0> {
  const static T value = T(1);
};


Com ele em mãos, podemos fazer a busca binária propriamente dita:

template<class T, const int start, const int len>
struct compare10 {
  static T compare(const T x) {
    if (x >= p10<T, start + len/2>::value)
      return compare10<T, start + len/2, len/2>::compare(x);
    else
      return compare10<T, start, len/2>::compare(x);
  }
};

template<class T, const int start>
struct compare10<T, start, 1> {
  static T compare(const T x) {
    return p10<T, start>::value;
  }
};


E depois basta fazer o bootstrap, usando agora uma função pouco conhecida da biblioteca padrão do C++: o digits10, que volta a quantidade máxima de dígitos decimais que cabe num tipo qualquer.

template<class T>
T template_power10(T x) {
  return compare10<T, 0, numeric_limits<T>::digits10>::compare(x);
}


Abaixo, uma versão completa, já com benchmark, para comparar as duas versões. Na minha máquina, a versão com metaprogramming calcula um milhão de valores em metade do tempo da versão original. Isso é graças à complexidade reduzida da versão com metaprogramming, que é apenas O(log log n), com dois logaritmos, como eu queria demonstrar :)

Benchmark das duas versões

Variações sobre um tema

Pelos idos de 2006, um dos meus hobbies era resolver Sudokus. Como eu sou uma criatura que se empolga fácil, tinha sudokus em papel, jogos de sudoku no Nintendo DS, sudokus em tudo quanto é lugar. Depois de algum tempo eu tive que desistir dos sudokus, mas não foi porque eu enjoei, foi porque eu não achava mais puzzles no meu nível!

Os puzzles da Coquetel, por exemplo, podem ser resolvidos quase que exclusivamente com técnicas simples, e esses eu resolvia em poucos minutos. Só os puzzles marcados como "diabólicos" e "grande mestre" que precisavam de alguma técnica mais avançada, como X-Wing e unicidade de soluções.

Mas ainda assim eu não parei de comprar as revistinhas. O Hofstadter dizia, no Metamagical Themas, que "Making variations on a theme is really the crux of creativity", e o povo dos sudokus levou isso a sério. Os sudokus tradicionais já não tinham mais graça, mas a Coquetel também publica variações sobre o tema, como os Sudokus Irregulares:


Nesse irregular 6x6 valem as mesmas regras dos sudokus normais, cada linha, coluna ou região precisa ter todos os números de 1 a 6. Mas a maioria das técnicas avançadas não funciona, então esses puzzles ainda têm graça pra mim.

Computacionalmente, o sudoku irregular também pode ser resolvido com o exact cover. No caso desse irregular 6x6, a matriz resultante tem 216x144 elementos. Como eu já tinha a biblioteca de dancing links do post anterior, criar um solver para esse sudoku foi bem simples:

Solver do sudoku irregular 6x6
Entrada exemplo para o solver
Update da lib de exact cover

Dessa vez eu mudei um pouco a api do exactcover.h, a versão original só permitia contar o número de soluções possíveis, agora ela é um template que recebe um functor que é chamado a cada solução encontrada. Como o template ficou mais genérico, agora você pode implementar a variação que quiser sobre o processamento das soluções :)

Domino Dancing

Dois dos meus hobbies prediletos são resolver puzzles e desafios de programação. Por isso, é natural que eu me empolgue quando consigo fazer os dois ao mesmo tempo!

O caso em questão é o problema DOMINO do spojinho. O problema pede pra você resolver computacionalmente um puzzle clássico, que consiste de um tabuleiro 7x8 com números em cada casa, que deve ser completamente preenchido pelos 28 dominós. Esse puzzle é comum em coletâneas como a Denksel da Croco-Puzzle, e também tem disponível online em vários sites, como a applet java do Janko. Abaixo, um exemplo de puzzle no estado inicial e resolvido:


Como resolver computacionalmente esse puzzle? Quem tem alguma experiência logo nota que dá pra resolver com backtracking. Quem tem muita experiência, percebe que, na verdade, esse problema é uma instância do Exact Cover!

De fato, pra modelar o puzzle como um exact cover, basta verificar os constraints: você precisa colocar cada um dos 28 dominós no tabuleiro, e cada uma das 56 casas do tabuleiro deve estar ocupada por um único dominó. Assim, a matriz terá 84 colunas.

As linhas você obtém notando que cada bloco 2x1 ou 1x2 do tabuleiro só pode ser preenchido por um único dominó. Daí, você tem 49 jeitos diferentes de encaixar um dominó na horizontal, e 48 jeitos de encaixar na vertical, então o total de linhas será 97.

Com a matriz 97x84 em mãos o problema está praticamente resolvido, basta implementar um algoritmo de exact cover. As vantagens dessa abordagem são três: primeiro, reduzir a um problema mais genérico é elegante, segundo, você pode reusar código de algum outro exact cover que você já tenha feito na vida. Por fim, como o exact cover é um problema bem conhecido, é de se esperar que a literatura tenha soluções espertas para esse algoritmo, e há mesmo: o exact cover pode ser resolvido com um algoritmo popularizado pelo Knuth, os Dancing Links.


O exact cover em si é NP-Completo, e a solução naïve é exponencial. Usando Dancing Links, você pode implementar uma heurística muito rápida, que diminui consideravelmente o branching factor da busca pela solução. Quem quiser entender a fundo o algoritmo, só precisa ler o paper original do Knuth, a minha implementação em C++ é a abaixo:

Dancing Links em C++

(Warning: Se você tem dificuldade com listas ligadas e medo de listas circulares biligadas, passe longe do Dancing Links. A implementação requer uma lista bicircular pentaligada, o que pode causar graves coceiras em quem tem alergia a ponteiros).

O exact cover também pode ser usado pra resolver o Sudoku e o N-Queens, então é uma técnica que vale a pena conhecer.

Firefox e os Fractais

Números romanos são só um dos ritos de passagem que todo programador, mais cedo ou mais tarde, acaba fazendo. Certa vez eu notei que era uma vergonha nunca ter implementado o conjunto de Mandelbrot na vida. Resolvi isso rapidamente escrevendo uma versão em Actionscript, e acabei ficando impressionado com o resultado! Com um pouquinho de otimização, o arquivo swf resultante tinha menos de 512 bytes.

É claro que eu resolvi tomar como desafio fazer o mesmo em outras linguagens. Em javascript foi tranqüilo, em java eu tive que apelar: só consegui atingir a barreira de 512 bytes escrevendo o bytecode diretamente na unha (source). Em python foi tão tranqüilo que, com a ajuda dos amigos, eu consegui reduzir para menos de 256 bytes:

• Versão em flash
• Versão em javascript
• Versão em java
• Versão em python

De todas elas, a mais lenta certamente é a versão em javascript. Mas com todos falando bem do novo interpretador javascript do Firefox 3 Beta 5, eu resolvi usar esse fractal como benchmark. Fiz uma pequena modificação para imprimir o tempo gasto com o traçado, e eis os resultados:

1. Firefox 3: 4.0 s
2. Safari 3.1: 4.0 s
3. IE 6: 8.3 s
4. Firefox 2: 14.4 s
5. Opera: 21.3 s

Eu rodei todos os browsers na mesma máquina, um intel dual core com windows. O ganho foi como o esperado mesmo, o Firefox ficou mais ou menos 3x mais rápido. Mas em compensação ele não é tão mais rápido quanto dizem, só tem a mesma velocidade do Safari.

Ritos de passagem

Quase todas as comunidades possuem algum tipo de rito de passagem. Entre os cristãos, há o batismo, entre os judeus, o bar mitzva, entre as patricinhas, o baile de debutante. É claro que entre os programadores também existem ritos de passagem, sendo que o mais comum é o hello world. Nas universidades brasileiras da década de 90, um rito muito comum eram os números romanos.

Tanto a USP quanto a UFMG, nessa época, pediam aos alunos como primeiro exercício que escrevessem um programinha que fizesse a conversão de um inteiro para seu equivalente em números romanos. Parece uma tarefa extremamente simples, e é mesmo, embora seja apenas a primeira etapa de uma curva crescente de dificuldade (no meu ano, os exercícios seguintes foram o cálculo das forças em uma treliça, e um sistema para projeção de sólidos 3D).

Por outro lado, quanto mais simples a tarefa, maior a oportunidade para você exercer sua criatividade. Consta que um dos alunos da UFMG mandou o código abaixo como resposta para o exercício:

scanf("%d", &n);
if (n==1) printf("I");
if (n==2) printf("II");
...
if (n==3999) printf("MMMCMXCIX");


Quando eu vi essa solução, achei muito legal: é claro que não tem como o professor reclamar, afinal ela satisfaz ao enunciado proposto. Mas embora seja uma solução divertida, ela não é ótima. De fato, do jeito que está escrito, essa solução é apenas O(n).

Como fazer então uma solução que seja ao mesmo tempo sacana e ótima? Uma mudança simples seria criar um array de strings inicializado com os valores, assim você teria desempenho O(1). Por outro lado, simplesmente criar o array por extenso, como no caso acima, não tem muita graça. Mais legal seria criar o array em tempo de execução usando algum tipo de metaprogramming, como o template metaprogramming em C++.

Ainda assim, template metaprogramming já não é tão obscuro quanto costumava ser, muita gente já conhece o método. O que poucos conhecem é uma maneira de fazer metaprogramming usando apenas features presentes na linguagem C. A solução que eu acabei criando foi a abaixo, com um encantamento conhecido apenas pelos sacerdotes do ioccc, o preprocessor metaprogramming:

Números romanos usando preprocessor metaprogramming

O preprocessador do C não é turing-complete como os templates do C++, mas é suficiente para escrever rotinas mais simples. A minha solução usa um autômato finito, onde a transição de estado é feita através de um #include de si mesmo, e os estados são codificados bit a bit, em BCD.

Eu recomendo testar o código pra conferir que ele funciona, embora já avise de antemão que a compilação pode demorar um pouco. Pra ver o código gerado sem as diretivas do preprocessador, basta usar a flag -E do gcc.

Depois de terminado esse programinha, acabei ficando com mais uma idéia divertida de uso de metaprogramming, mas essa fica pra depois :)