O algoritmo mais rápido do oeste

Dia desses um amigo me passou um probleminha bem divertido. O enunciado é simples: dado um vetor ordenado, construa uma binary search tree, balanceada, com os mesmos elementos do vetor. Esse problema tem uma solução trivial em O(n log n) e uma solução esperta em O(n), mas eu acabei achando uma terceira solução, com o algoritmo mais rápido do oeste: ele roda em O(0)!

Vamos revisar as respostas tradicionais antes. A solução trivial é usar uma árvore binária auto-balanceante, como a AVL ou a Red-Black, e inserir os elementos nela, um por um. Como você vai inserir n elementos e o custo de inserção é O(log n), então o total desse algoritmo é O(n log n).

Mas fazendo dessa maneira, você não está usando a informação de que o vetor estava ordenado. Com uma recursão simples você consegue aproveitar esse fato, e construir a árvore em apenas O(n). Um exemplo em python é assim:

def build_tree(vector):
  def partition(start, end):
    if (end < start):
      return None
    middle = start + (end-start) / 2
    return (vector[middle], 
            partition(start, middle-1), 
            partition(middle+1, end))
  return partition(0, len(vector) - 1)

Aparentemente, essa solução é ótima. Se você vai criar uma estrutura de dados nova usando os dados da antiga, o mínimo de trabalho que você precisa é copiar os elementos de uma pra outra, e isso limita em O(n) certo? Quase! Você pode contornar essa limitação se usar uma estrutura de dados implícita, e é isso que vamos tentar fazer.

Para usar o vetor ordenado original como uma estrutura implícita, nós precisamos criar um isomorfismo entre o vetor ordenado e a binary search tree. Se for possivel criar esse isomorfismo, então na verdade eu não preciso fazer nada pra converter o vetor, ele já é equivalente à árvore. E um algoritmo que não faz nada é um algoritmo O(0)! Eu poderia chamar essa estrutura de dados nova de Implicit Balanced Binary Search Tree, mas é muito comprido, então vou chamá-la de Cactus Binário (apropriado, daqui em diante as coisas vão ficar espinhosas).

Mas vamos ao isomorfismo então. Como o vetor está ordenado, então os valores dos elementos em si não importam muito, e nós podemos trabalhar só com os índices, sem perda de generalidade. Além disso, para as demonstrações ficarem mais simples, vamos dividir em dois casos. Digamos que o tamanho do vetor é N, então o primeiro caso que vamos tratar é quando N pode ser escrito como 2^h-1 para algum h>0. Nesse caso, o Cactus Binário é uma árvore perfeita (todas as folhas estão no mesmo nível e todos os nós internos tem dois filhos):

Para criar o isomorfismo, eu preciso fornecer duas funções que, para um dado índice, me retornem o índice do filho à esquerda e do filho à direita. Depois de pensar um pouco, acabei chegando nas funções abaixo:

left(x) = x - (x & (-x)) / 2
right(x) = x + (x & (-x)) / 2

Isso é mágica? Não, é matemática! Para entender porque as expressões funcionam, é só olhar para a árvore anterior usando a visão além do alcance:

Olha só que bacana, quando você olha para os índices da árvore em binário, parece que o nível em que está o índice é igual ao número de zeros no final da representação binária dele. Na verdade, nós podemos provar isso, usando indução na altura da árvore.

Vejamos: para uma árvore de altura unitária, o único elemento é o 1, que está no nível zero e tem 0 bits zero no final, então a base de indução está ok. Se a árvore tem altura h, então ela tem 2^h-1 elementos. O elemento central é o 2^h-1, que tem h-1 zeros no final e está no nível h-1. A sub-árvore da esquerda é uma árvore balanceada de tamanho 2^h-1-1, que funciona pela hipótese de indução. A sub-árvore da direita é igual à da esquerda, se você ligar o bit mais significativo de cada índice. Setar o bit mais à esquerda não muda a quantidade de bits zero à direita, então podemos usar a hipótese de indução novamente, e pronto: QED.

Agora ficou fácil achar a fórmula né? Se você tem um índice x no nivel h, então o filho à esquerda é o maior número com h-1 zeros à direita que seja menor que x, e o análogo vale para o filho à direita. Ou seja, pra descer do nível h para o nível h-1 pela esquerda, basta subtrair 2^h-1  (ou somar, se for pela direita).

E computacionalmente como eu faço isso? Basta achar o bit 1 mais à direita do número, dividir isso por dois e subtrair do número original. Para achar o primeiro bit 1, o truque é calcular x & (-x). Vamos provar: pela definição de complemento de dois, negar é o mesmo que inverter e somar um. Se o número em binário for algo do tipo xxxx10000, negado ele fica yyyy01111, e quando você soma um, o carry propaga até bater no zero, ficando yyyy10000. Quando você faz o AND, os primeiros bits dão zero porque o AND de um bit e seu complemento é sempre zero, e os últimos bits dão zero porque o AND de qualquer bit com zero dá zero. Sobra só o primeiro bit 1 à direita, QED.

Juntando tudo agora as duas fórmulas são evidentes. Mas pra completar o isomorfismo eu ainda preciso de uma fórmula que me indique onde está a raiz da árvore. A fórmula em si é simples: a raiz é a maior potência de dois menor ou igual ao tamanho do vetor, ou seja, basta achar o primeiro bit 1 à esquerda do tamanho.

Para achar o primeiro bit 1 à esquerda não tem nenhum truque tão simples quanto à direita, mas felizmente existe uma alternativa: os processadores x86 possuem um opcode meio desconhecido chamado BSR (bit scan reverse), que retorna exatamente onde está esse bit à esquerda. Em processadores bons, como o Core 2 Duo, o BSR é uma operação O(1), mas isso não é verdade em todos os processadores (por exemplo, não é verdade nos AMDs). Se você não tiver o BSR por hardware, pode fazer uma busca binária nos bits do número, que é uma operação O(log log n).

Isso conclui o primeiro caso. Vamos ver agora o que acontece quando o tamanho é um número qualquer. O pior caso do algoritmo é quando o tamanho é uma potência de dois, e o Cactus Binário fica assim:

Uepa! Essa árvore está balanceada? Bem, depende da sua definição de balanceada. Nessa árvore, o maior caminho entre a raiz e uma folha é 4, o que parece muito. Mas, na verdade, não dá criar uma árvore que tenha um caminho máximo menor. Se o caminho máximo fosse 3, então a árvore poderia ter no máximo 7 elementos, e pelo princípio da casa dos pombos a árvore de tamanho 8 com caminho máximo 3 é impossível. De fato, se você tentar balancear mais a árvore, o melhor que você consegue é algo assim:

Logo, o caminho máximo do Cactus Binário é o mesmo caminho máximo da árvore balanceada. Note, entretanto, que isso é melhor que as árvores auto-balanceantes! O caminho máximo da AVL, no pior caso, é 1.44log(n), e o da Red-Black é 2log(n), ambos maiores que o Cactus Binário.

Com isso, nós provamos o isomorfismo para os dois casos, e o algoritmo O(0) realmente funciona! Mas a grande dúvida agora é: ele serve pra alguma coisa? Na verdade, a grande vantagem da árvore binária é que ela faz inserções em O(log n), enquanto que uma inserção no Cactus Binário é apenas O(n). Por isso, nas aplicações práticas da vida real a árvore binária ganha. Aparentemente, a única utilidade do Cactus Binário é para escrever blog posts :)

Tatuagens em cadeia

Uma classe de artistas que eu admiro muito são os tatuadores, e por um motivo simples: eles não podem errar. Um desenhista pode usar borracha, um arte-finalista pode cobrir o erro com tinta branca, um artista digital pode apertar control-Z, mas o tatuador não pode fazer nada disso. Um escultor, se erra, tem a opção de jogar fora o bloco de pedra e começar de novo, mas o tatuador nem essa opção tem.

Exatamente por isso, tem um tatuador que eu conheço que se recusa a tatuar nomes de namorados. Tatuagens são permanentes, mas namoros são efêmeros: imagine namorar alguém que tem uma tatuagem com o nome do ex? Nesse caso, uma solução prática é reciclar a tatuagem e namorar alguém que tenha o mesmo nome do seu ex. Não é elegante, mas funciona.

Isso leva a conclusões curiosas. Veja o meu caso, por exemplo. Se eu tatuar ILA e minha esposa tatuar RICARDO, aparentemente ela estaria levando vantagem. É muito mais fácil arrumar outro Ricardo que arrumar outra Ila. Por outro lado, eu não precisaria procurar apenas Ilas, eu poderia arranjar uma Priscila ou uma Camila, bastando adicionar mais letras na tatuagem!

Isso naturalmente leva à pergunta: qual é o máximo de namoradas que um Serial Tatuator pode ter? Por exemplo, uma pessoa poderia ir de Ana para Iana e depois Fabiana. Será que dá pra achar uma cadeia de tatuagens de tamanho quatro? Cinco? Hora de fazer uma simulação!

Antes de começar, precisamos de um corpus de nomes. O jeito mais fácil de conseguir o corpus é escrevendo um crawler para um desses sites com nomes de bebês. Uma busca rápida e eu achei o Babyhold, que tem a vantagem de separar os nomes por sexo.

Script python para fazer crawling de nomes de meninas

O script retornou mais de 7100 nomes de meninas. Agora é só fazer um outro script que calcule a máxima cadeia de nomes. Como a quantidade de nomes não é tão grande assim, um algoritmo cúbico bobo dá conta do recado:

Script python para achar cadeias máximas

O script acabou achando várias soluções de tamanho seis, mas sem dúvida a seqüência mais curiosa foi a abaixo :)

Li -> Lin -> Elin -> Elina -> Angelina -> Angelina Jolie

Bônus!

Embora o algoritmo cúbico funcione bem para 7100 nomes, certamente há como otimizá-lo. Ao invés de fazer uma solução mais rápida, dessa vez eu vou abrir para os leitores enviarem suas soluções! Clicando no link abaixo você vai para o primeiro Ricbit Jam, feito com a engine do spoj. Ganha quem enviar a solução mais rápida nos próximos quinze dias. Não tem prêmio, mas o vencedor terá seu nome divulgado no próximo post :)

Ricbit Jam #1

Não adianta tentar mandar o programa em python acima, ele não acha todas as soluções e vai dar TLE no servidor do spoj. Mas você pode usar como base pra fazer a sua solução.

Como aprender computação

Dia desses o Lameiro me passou uma lista de livros para desenvolvedores, e eu achei a lista bastante curiosa. Se eu fosse criar uma similar, eu não colocaria nenhum dos livros citados naquela lista! Ao invés disso, a minha lista com os top 10 livros de computação seria a abaixo:


1 - Gödel, Escher, Bach (Hofstadter): Eu começaria com o GEB, por dois motivos. O primeiro é que ele é um ótimo reality check: se você não gostar do GEB, então mude de área, porque computação não é a sua praia :) O segundo motivo é que esse livro tem uma excelente introdução à lógica (proposicional e de predicados), que é a ferramenta básica onde você constrói a ciência da computação.

2 - Concrete Mathematics (Knuth): Você não precisa saber matemática para programar, mas se você quiser ser um bom programador, então matemática é essencial. O Concrete tem todo o básico que você precisa pra fazer análises de complexidade computacional, e tudo escrito de maneira extremamente bem-humorada.

3 - Algorithms in C++ (Sedgewick): Se você já sabe lógica e matemática, então agora pode partir pro estudo de algoritmos. O Sedgewick tem todos os algoritmos básicos, e é uma leitura bem leve: se você está começando com algoritmos agora, esse precisa ser o seu primeiro livro. Ele não vai muito fundo em nenhum tópico, mas isso é compensado pela extrema didática nos tópicos. O livro ainda tem código exemplo pra todos os algoritmos, e várias edições, uma pra cada linguagem (eu sei que tem, pelo menos, C, C++, Java e Pascal).

4 - Introduction to Algorithms (Cormen): Os algoritmos que você aprendeu no Sedgewick, você vai estudar em detalhes no Cormen. Esse livro é extremamente formal, e talvez por isso é o livro-texto usado nos cursos de computaçao do MIT. Ele também cobre algoritmos mais avançados, que o Sedgewick apenas cita (por exemplo, Fibonacci Heap)

5 - The Art of Computer Programming (Knuth): O tAoCP está para o Cormen assim como o Cormen está pro Segdewick, aqui você vai dissecar os algoritmos até o último bit deles. Além de ser uma coleção excelente, a encadernação é muito bonita (mas não se engane, ele fica ótimo na prateleira, mas melhor ainda na sua cabeça).

6 - Effective C++ (Meyers): Agora que você sabe os algoritmos, precisa de uma linguagem para programá-los. Pra falar a verdade, a escolha de linguagem nem importa tanto assim, mas se você escolher uma, aprenda-a tão bem quanto possível. Eu escolhi C++, e esse livro do Meyers é o que diferencia as crianças dos adultos (especialmente para aquelas horas quando você cria uma classe sem destrutor virtual e não sabe por que a memória está vazando).

7 - Effective STL (Meyers): Já teve uma época em que eu não gostava de C++, mas isso é porque eu sou velho o suficiente pra ter mexido em C++ antes que os compiladores tivessem templates. Com templates a linguagem fica muito mais atraente, e esse é o livro que vai te ensinar a dominar a STL.

8 - The Practice of Programming (Kernighan, Pike): Se você leu tudo até agora, então você já é um programador muito bom na teoria. Na prática, entretanto, tem um monte de skills que ainda faltam. Nesse livro você aprende sobre as coisas que usualmente não se aprende na escola: debugging, otimização, unit testing, documentação.

9 - Programming Pearls (Bentley): Com os livros lidos até agora, você já deve ser um excelente programador. O passo final é passar de programador para um true hacker, e esse é um passo que não requer só conhecimento, você também precisa de manha, criatividade e insight. Eu não sei se dá pra ensinar essas coisas, mas esse livro certamente é o que chega mais próximo disso.

10 - The Mythical Man-Month (Brooks): Depois de ler todos os livros acima você estará próximo do nirvana, mas pra atingir o zen da programação de verdade, é preciso lembrar que projetos não precisam só de computadores, precisam de pessoas também. O Mythical Man-Month é um livro de gerência de projetos de software escrito em 1975, mas é surpreendente como ele continua atual. A tecnologia avança, mas as pessoas continuam as mesmas :)

É claro que pra manter uma lista com só dez itens, muita coisa boa fica de fora. Mas a lista acima tem um mérito: foi com esses livros que eu aprendi computação de verdade (vale lembrar que eu sou autodidata, minha graduação foi em engenharia elétrica, e eu quase não tive computação em aulas). Se funcionou pra mim, pode ser que funcione pra você também :)

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 :)

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.

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.

Otimização com ábacos

Semana passada ficou pronta mais uma réplica da Máquina de Diferenças n°2 do Babbage. Essa é a segunda a ser construída a partir da especificação original, pesa 5 toneladas, e ficará em exposição no Computer History Museum de Mountain View (que eu visitei no ano passado).

Construir a partir da especificação original é um bocado caro, e reservado só pra museus e milionários mesmo. Porém, isso não impediu um hobbysta de criar sua própria máquina de diferenças usando LEGO, o que chega a ser até mais impressionante. Pra quem quiser simular a máquina de diferenças apenas em software, pode tentar o problema CMPLS no spoj, que é exatamente isso.

O que talvez não seja muito óbvio é que, apesar de usar apenas processamento mecânico, a máquina de diferenças é um computador digital (ela possui um clock, dado pela rotação de um eixo interno, e executa uma adição com carry a cada quatro rotações). Já os computadores analógicos podem ser bem mais bizarros, como o computador de sabão do post anterior.

Computadores analógicos possuem uma vantagem sobre os digitais: não estão presos aos mesmos limites computacionais. É bastante simples construir, por exemplo, um circuito, utilizando amplificadores operacionais, que resolva equações diferenciais complexas em tempo bem inferior a um computador digital (embora o processo prejudique a precisão).

Um exemplo mais curioso é o problema da ordenação de um conjunto de números. É possível demonstrar que um computador digital, usando como elemento básico de computação a comparação de dois valores, não pode ser melhor que O(n log n) (essa demostração está em qualquer livro básico de algoritmos, como o Cormen). Mas os computadores analógicos não tem essa restrição, sendo que existe até mesmo um algoritmo que resolve em O(1)!

Uma implementação simples desse algoritmo é com ábacos: primeiro você dispõe os números que você quer ordenar na base unária, como estão os números 2, 4, 1, 3, 3 na figura abaixo:


Depois, basta levantar o ábaco e deixar a gravidade fazer seu serviço:

Impressionante, não? Os números ficam perfeitamente ordenados. Apesar de ser um truque muito simples, esse método só foi inventado em 2002. No paper original, os autores chamam o método de Bead Sort e sugerem, além da implementação com ábacos, outras implementações (com redes resistivas, autômatos celulares e matrizes de flip-flops).

Como o Bead Sort usa operações analógicas, não dá pra analisar de verdade a complexidade computacional dele. Quando fazemos análises de big-oh, queremos saber quantos passos o algoritmo leva pra terminar, e o Bead Sort tem só um passo (levantar o ábaco), sendo que nesse aspecto ele é O(1) mesmo. Mas intuitivamente, o que queremos quando fazemos análise de complexidade é descobrir como o tempo de execução varia com o tamanho da entrada. Desse ponto de vista, a complexidade é O(sqrt(n)), se você considerar o tempo que uma bolinha leva pra ser puxada pela gravidade (no vácuo), ou então O(n), se você levar em conta que a bolinha vai atingir uma velocidade terminal devido à resistência do ar.

Otimização com água e sabão

Ano passado fui pela primeira vez à Estação Ciência, em São Paulo, e fiquei espantado com a qualidade. Eu já visitei o Natural History Museum de Londres, e o American Museum of Natural History em New York, então eu digo, com conhecimento de causa, que a Estação Ciência não fica nada a dever aos museus de ciência do primeiro mundo. Tem dinossauros, planetário, simulador de terremoto, e até simulador de tsunami (que não tinha nos museus lá de cima).

Na verdade, o museu brasileiro tem uma vantagem sobre os outros. Por lá, a média é de, mais ou menos, uns quarenta visitantes por guia, e aqui a média é de três guias por visitante! Se você tiver bastante tempo pra visitar, os guias irão lhe dar explicações extremamente detalhadas das exposições.

A minha predileta foi na seção de matemática experimental, onde fizemos a experiência do sabão. Você pega duas placas paralelas quaisquer, coloca quatro pregos nelas, e mergulha num balde de água com sabão. Eu estava esperando que o sabão grudasse nos pregos e fizesse um quadrilátero, mas na verdade o que ele faz é uma estrutura em duplo Y:

(Você pode fazer em casa esse experimento, mas para melhores resultados, coloque um pouco de glicerina na água. Experimente também com outras quantidades de pregos, e com figuras tridimensionais).

A explicação é que o sabão, como bom sistema físico, vai procurar a posição de energia mínima, que nesse caso é equivalente a minimizar a soma dos comprimentos da película. Anedoticamente, esse problema é conhecido como o problema da estrada (como ligar quatro cidades com estradas, gastando a menor quantidade possível de asfalto?), e a solução ótima é conhecida como Steiner Tree.

À primeira vista pode parecer que a Steiner Tree é equivalente à Minimum Spanning Tree, mas não é. Pra calcular a spanning tree, você só pode usar os pontos do grafo, já na Steiner tree você pode adicionar pontos novos. Isso faz uma diferença enorme na complexidade: existem algoritmos quase lineares para resolver a spanning tree, já a Steiner tree é NP-completa.

Na verdade, há quem diga que isso é uma prova de que P=NP, como os autores desse paper no Arxiv. Segundo eles, existe um computador analógico feito com sabão que resolve um problema NP-completo em tempo polinomial, logo P=NP :)

Pra quem quiser brincar com Steiner trees, o problema ELC do spoj pede pra resolver uma instância com 3k pontos. Como isso é demais para um brute force, você vai ter que aproximar: por exemplo, iterando a partir de uma spanning tree; ou então, calculando a triangulação de Delaunay e depois resolvendo a Steiner tree local de cada triângulo.

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