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