Brain Dump

domingo, 20 de julho de 2008

Eu podia estar roubando

É isso. Eu podia estar roubando, eu podia estar matando, mas estou aqui escrevendo no blog. Reclamações sobre a periodicidade do blog devem ser enviadas diretamente à Rockstar :)

Ao falar de GTA, sempre tem quem o associe ao terrorismo islâmico. Mas pra mim, o que incomoda mais é a associação imediata do islã com o terrorismo, em oposição à sociedade ocidental civilizada. Na verdade, já houve um período onde acontecia o oposto, enquanto a ciência florescia na cultura islâmica, os cristãos saqueavam e destruíam em nome da fé (mas na época chamavam isso de cruzadas, ao invés de terrorismo).

Foi nessa época em que viveu Al-Khwarizm (cujo nome deu origem às palavras algarismo e algoritmo), nessa época em que Bhaskara escreveu Lilavati e resolveu a equação de segundo grau, e também nesse período é que surgiram as primeiras demonstrações por indução finita. Mas a influência mais direta que a matemática islâmica teve em mim foi, sem dúvida, através das obras do Júlio César de Mello e Souza, que escrevia com o pseudônimo de Malba Tahan.

Malba Tahan se inspirava na cultura islâmica pra escrever seus livros de divulgação científica, sendo que o mais famoso deles é o Homem Que Calculava. Num post anterior eu tinha dito que o cálculo do Pi com método de monte carlo tinha sido o segundo puzzle que levei mais tempo pra resolver; pois o primeiro puzzle foi justamente tirado de um dos capítulos do livro: como escrever qualquer número usando apenas quatro quatros. No livro, Malba Tahan mostra as soluções até o 10, que são mais ou menos assim:

0 = 44 - 44
1 = 44 / 44
2 = 4/4 + 4/4
3 = (4 + 4 + 4) / 4
...

Eu tinha só dez anos quando li pela primeira vez, e consegui estender a seqüência até o 20. O posfácio do livro dizia que era possível escrever até o 100, mas eu levei mais de uma década até ver a solução completa!

Depois dos 20 primeiros, eu só consegui fazer algum progresso significativo quando estava no colégio técnico. Eu notei que o problema ficava mais simples se eu reduzisse o número de quatros, montando primeiro todas as soluções com um quatro, depois todas com dois quatros, e assim por diante. Eu não sabia na época, mas o que eu estava fazendo era um forward chaining manual. Mesmo assim, eu só consegui fazer os números pares até 100, minha solução tinha um monte de buracos nos ímpares.

Quando eu estava na faculdade, com os skillz mais afiados, eu larguei mão de fazer manualmente e escrevi um script pra fazer o serviço automaticamente. Eu coloquei no script não só as quatro operações básicas, mas também todas as outras que estão no Concrete Mathematics: raiz quadrada, função piso e função teto, binomiais, potências, fatorial, fatorial crescente e decrescente. A única regra é não usar nenhuma notação que envolva letras (como sin, cos e log; se você permitir log, o problema perde a graça).

Meu script conseguiu finalmente resolver todos os números até 100, e até encontrou múltiplas soluções pra eles! Atribuindo um peso a cada operação, eu mandei ele imprimir somente a soluçao mais simples (por exemplo, adição tem peso baixo, piso e teto tem peso alto). A tabela com as respostas está abaixo, junto com o programa que escrevi na época:

Solução gerada pelo script
Script original, escrito em C

O script pode resolver também os 5 cincos, os 6 seis, e qualquer outra combinação, desde que você esteja disposto a esperar o suficiente. O código original é bem ruinzinho, mas na verdade eu me orgulho de achar que meu código escrito há dez anos é ruim, significa que eu aprendi alguma coisa desde então :)

Marcadores: , , , ,

segunda-feira, 23 de junho de 2008

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

Marcadores: , , ,

sábado, 7 de junho de 2008

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(2n), 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, é pq2, na quarta pq3 e assim por diante. Então o valor esperado é:

X = 1p + 2pq + 3pq2 + 4pq3 + ...
X = p (1 + 2q + 3q2 + 4q3 + ....)

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 + 4q2 + ....))

Agora você tira da cartola y=1+q+q2+q3+... e substitui:

X = p (1 + q(2 + 3q + 4q2 + ....))
X = p (1 + q(y + 1 + 2q + 3q2 + ....))
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(107)=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.

Marcadores: , , , , , ,

quarta-feira, 21 de maio de 2008

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 n15. 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 é log2n:


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.

Marcadores: , , , , ,

sábado, 3 de maio de 2008

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.

Marcadores: , , , ,

sábado, 19 de abril de 2008

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

Marcadores: , , ,

terça-feira, 15 de abril de 2008

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

Marcadores: , ,

domingo, 13 de abril de 2008

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.

Marcadores: , ,