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