Mostrando postagens com marcador programação funcional. Mostrar todas as postagens
Mostrando postagens com marcador programação funcional. Mostrar todas as postagens

segunda-feira, 19 de março de 2018

A Marcha dos Lemmings

Lemmings são pequenos roedores que vivem na tundra ártica. Eles possuem pelos longos e macios, rabo curto, são herbívoros e explodem.

Sim, eles explodem! Os lemmings possuem a capacidade de ligar um contador luminoso sobre suas cabeças. A cada passo que dão, o contador decrementa, até que finalmente explodem quando o contador chega a zero.


O que você faria se soubesse que vai explodir em seis passos? Parece que não tem muitos passeios possíveis com exatamente seis passos, mas a intuição é meio falha aqui. Na verdade, o número de passeios possíveis é exponencial no número de passos! Como fazer um algoritmo que calcula o número de caminhos possíveis que o lemming pode fazer antes de explodir?

As Regras do Jogo


Para modelar o problema, vamos supor que os passos que lemming pode fazer são descritos por um grafo orientado. Por exemplo, no grafo abaixo o lemming começa no nó A. Do A ele pode ir para B ou C, de B ele pode ir para C, e de C ele pode ir para A.


Se o nosso lemming pode dar seis passos, então existem sete caminhos diferentes que ele pode fazer:

ABCABCA
ABCACAB
ABCACAC
ACABCAB
ACABCAC
ACACABC
ACACACA

Logo ##p(6)=7##, e podemos até tabelar ##p(n)## para todos os ##n## de 0 a 6:

$$p(0..6)=1,2,2,3,4,5,7$$
Dado ##n##, como fazer um algoritmo que calcule ##p(n)##?

Eu vou mostrar três algoritmos possíveis usando a linguagem do Wolfram Mathematica. Para isso, o grafo precisa ser descrito de alguma maneira que os algoritmos entendam. Eu vou usar uma matriz de adjacências: para um grafo com ##k## nós, a matriz é ##k\times k## e o elemento na linha ##i## da coluna ##j## é 1 se existe uma ligação saindo do nó ##j## e indo em direção ao nó ##i##, e 0 em caso contrário. Para o nosso exemplo, a matriz é:

$$M=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right]$$

Força Bruta


Sempre que eu tenho que fazer um algoritmo do zero, prefiro começar com a força bruta. É fácil e rápido de implementar, e embora não seja eficiente, serve para conferir o resultado dos algoritmos mais avançados.

Nesse caso, para implementar a força bruta você passa o grafo, o nó em que o lemming está, quantos passos faltam para terminar, e o resto é uma recursão simples. Como toda recursão, não esqueça de fazer o caso base! Para esse problema, o caso base é ##p(0)=1##, ou seja, tem um único caminho possível sem dar passo nenhum, que é não sair do lugar.

No fim, a solução em força bruta fica assim:
count[graph_, pos_, 0] := 1
count[graph_, pos_, size_] := Sum[
    If[graph[[i, pos]] == 1, count[graph, i, size - 1], 0],
    {i, 1, Length[graph]}]
Vamos conferir o resultado:
graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}
Perfeito! Mas qual é a complexidade desse algoritmo? Bem, ele precisa visitar cada um dos caminhos possíveis, e o número total de caminhos é exponencial, logo a complexidade vai ser do tipo ##O(c^n)##, onde a constante ##c## depende do grafo. Ou seja, você não quer usar esse algoritmo quando o ##n## é grande.

Multiplicação de Matrizes


Como melhorar então? Tem uma maneira calcular ##p(n)## mais rápido usando multiplicação de matrizes. É mais fácil mostrar como funciona por exemplos: no estado inicial, nós temos um lemming no nó A, então vamos montar um vetor correspondente: ##V=[1 \;0\; 0]^T##. Se multiplicarmos ##M## por ##V##, qual o resultado?

$$M V=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right] \left[\begin{matrix}1\\0\\0\end{matrix}\right]=\left[\begin{matrix}0\\1\\1\end{matrix}\right]$$
Ou seja, um lemming que está em A vai parar em B ou C. Continuando:
$$M V=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right] \left[\begin{matrix}0\\1\\1\end{matrix}\right]=\left[\begin{matrix}1\\0\\1\end{matrix}\right]$$
Dado um lemming em B ou C, o resultado é um lemming em A ou C. Se o vetor de entrada mostra quantos lemmings tem em cada nó, multiplicar pela matriz M vai mostrar onde eles podem estar depois de dar um passo.

Agora é só estender o raciocínio. Se multiplicar por M correponde a um passo, então multiplicar ##n## vezes por M é o equivalente a ##n## passos. E se você começou o processo com um vetor unitário, então a soma do elementos de V é o número de caminhos possíveis partindo do nó inicial!

Uma implementação possível é a abaixo:
count[graph_, pos_, size_] :=
  Total[MatrixPower[graph, size] . UnitVector[Length[graph], pos]]
Conferindo:
graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}
Bateu certinho novamente! E qual é a complexidade? Se você implementou a multiplicação e a exponenciação de matrizes do jeito mais simples, vai ser da ordem de ##O(k^3 n)##, que certamente é melhor que exponencial mas ainda não é boa o suficiente para usar na prática.

No entanto, você pode melhorar isso com implementações melhores! A multiplicação de matrizes pode usar o algoritmo de Strassen, e a exponenciação pode usar o método binário, aí o total fica ##O(k^{2.8} \log_2 n)##, bem melhor; e se a matriz for esparsa dá para ficar ainda mais rápido.

Função Geradora


Mas tem um método que supera o anterior, que é usar funções geradoras! Nós começamos tudo com um grafo orientado, e todo grafo orientado é equivalente a um autômato finito. Por sua vez, todo autômato finito é equivalente a uma gramática regular, então deve ter um jeito de descrever os caminhos do lemming usando uma gramática. Para o exemplo dado, a gramática fica assim:

$$\begin{align*}\alpha &\to A \;|\; A \beta \;|\; A \gamma \\ \beta &\to B \;|\;  B \gamma \\ \gamma &\to C \;|\; C \alpha\end{align*}$$
Todos os caminhos possíveis produzem strings que são aceitas por essa gramática, e a gramática não aceita nenhuma string que não seja um caminho possível.

O truque agora é usar combinatória analítica, que para o problema de enumerar strings de uma gramática é bem simples: basta trocar cada token por um ##z##, e cada OR por uma adição. A gramática vira um sistema de equações:

$$\begin{align*}\alpha &= z + z\beta +z \gamma \\ \beta &= z+z \gamma \\ \gamma &= z+z \alpha\end{align*}$$
Resolvendo para ##\alpha##, chegamos na função geradora:

$$\alpha=\frac{z + 2 z^2 + z^3}{1 - z^2 - z^3}$$
E para que serve a função geradora? Olhe só o que acontece quando abrimos a função em série de potências:

$$\frac{z + 2 z^2 + z^3}{1 - z^2 - z^3}=z+2z^2+2z^3+3z^4+4z^5+5z^6+7z^7+\dots$$
Ahá! Os coeficientes de ##z## são exatamente o número de caminhos no grafo! Para calcular ##p(6)##, basta olhar o coeficiente de ##z^7## (você precisa somar um porque a função geradora está contando o número de tokens na strings, e nós queremos o número de passos).

Podemos ir além agora. Quando a função geradora é uma função racional, você sempre consegue inverter e conseguir uma fórmula explícita para ##p(n)##:

$$p(n)=0.957\times 1.325^n+0.426\times 0.869^n\cos(1.469+ 2.438 n)$$
Essa fórmula é exata (ou seria, se eu não tivesse truncado os números para caber na tela). Note que o resultado sempre vai ser um inteiro, e que a parte oscilante da fórmula sempre vai ser menor que 0.4, então dá para simplificar mais ainda:

$$p(n)=\left\lfloor 0.5+0.957\times 1.325^n\right\rfloor$$
Basta pegar a parte exponencial e arredondar para o inteiro mais próximo. Podemos observar duas coisas aqui. Primeiro, agora nós temos um algoritmo onde o número de operações independe do valor de ##n##, então efetivamente esse algoritmo é ##O(1)##. Na prática você vai implementar usando bignum, aí não fica ##O(1)## de verdade; mas os outros algoritmos mostrados também vão ficar mais lentos na mesma proporção.

Segundo, nós provamos que aquele primeiro algoritmo era de fato ##O(c^n)##! Conseguimos até calcular quem é o tal do ##c##: é aproximadamente 1.325; ou, se você quiser o valor exato com radicais:

$$c=\frac{\sqrt[3]{9+\sqrt{69}}+\sqrt[3]{9-\sqrt{69}}}{\sqrt[3]{18}}$$
Para implementar esse método computacionalmente, é só reescrever aquele sistema de equações em forma matricial:

$$\left[\begin{matrix}\alpha\\ \beta \\ \gamma\end{matrix}\right] = z \left[\begin{matrix}0 && 1 && 1\\0 && 0 && 1\\ 1 &&0 && 0\end{matrix}\right] \left[\begin{matrix}\alpha\\ \beta \\ \gamma\end{matrix}\right]  + z\left[\begin{matrix}1\\ 1\\1\end{matrix}\right] $$
A matriz ali no meio é a transposta da matriz de adjacências! Agora já podemos implementar:
count[graph_, pos_, size_] := With[{n = Length[graph]},
  SeriesCoefficient[LinearSolve[
      IdentityMatrix[n] - z Transpose[graph],
      ConstantArray[z, n]][[pos]],
    {z, 0, size + 1}]]
Como esperado, a resposta é correta:
graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}

Um Problema em Aberto


Nós chegamos em um algoritmo que é ##O(1)##, dá para otimizar mais? Não! Esse é um limite teórico.

Quer dizer, existem algoritmos que são ##O(0)##. Por exemplo, "crie um algoritmo que, dado uma circunferência e um diâmetro, calcule a razão entre os dois". Esse é algoritmo é ##O(0)## porque a resposta é sempre ##\pi##, independe do input. Mas, se a resposta não é constante, ##O(1)## é o melhor algoritmo possível. (A não ser que você crie um algoritmo que dê a resposta antes do input chegar. Nesse caso me ensine como faz, porque você inventou a viagem no tempo).

Mas existe uma pequena variação no problema dos lemmings que o torna muito, muito mais difícil. Suponha que agora o lemming não pode passar duas vezes pelo mesmo nó. A força bruta para essa variação ainda é ##O(c^n)##. Mas ninguém bolou um algoritmo que seja melhor que a força bruta nesse caso!

Quando o grafo é finito, você pode usar a força bruta para enumerar todas as soluções, já que eventualmente o grafo acaba. Mas se o grafo for infinito, aí danou-se. Para esse caso, tem gente procurando a solução desde a década de 50. Já são setenta anos procurando o algoritmo e ninguém achou nada para um grafo genérico, nem mesmo um assintótico! O melhor resultado até hoje é um de 2011 onde acharam o ##c## para um grafo específico, o lattice hexagonal. Para todos os outros grafos, o melhor que temos são heurísticas e aproximações numéricas.

Se você estiver sem nada o que fazer no próximo fim de semana, essa pode ser uma boa maneira de ocupar o tempo!

domingo, 31 de julho de 2011

Metaprogramming com constexpr

Até ano passado, quando eu queria fazer algum projetinho onde eu precisava da maior velocidade possível, acabava escolhendo escrevê-lo em C++. Mas a minha vida mudou quando eu conheci o C++0x. Ele é muito superior ao C++ tradicional, adicionando features que faziam muita falta, como lambdas, variáveis auto e range for.


Uma das features novas que mais gosto é o constexpr, uma maneira de avisar ao compilador que o resultado de função pode ser determinado em tempo de compilação. Combinando isso com funções recursivas, nós temos uma maneira nova de fazer metaprogramming em C++, bem melhor que o tradicional template metaprogramming.

Para mostrar a superioridade do constexpr, resolvi escrever de três maneiras diferentes uma função bem simples: o totiente de Euler, cuja definição é a quantidade de inteiros positivos menores que n, e que não possuem fatores em comum com n. Escrito em C++ tradicional, essa função fica assim:

int gcd(int a, int b) {
  if (b == 0)
    return a;
  else
    return gcd(b, a % b);
}

int totient(int x) {
  int sum = 0;
  for (int i = 1; i < x; i++) {
    if (gcd(x, i) == 1) sum++;
  }
  return sum;
}

int main() {
  cout << totient(1000) << endl;
}

O código em si é bem simples de entender e compila em qualquer lugar. Eu calculei o GCD usando o algoritmo de Euclides, que é O(log n), então essa implementação do totiente é O(n log n).

Compare agora com o mesmo algoritmo escrito em template metaprogramming:

template<int a, int b>
struct GCD {
  const static int result = GCD<b, a % b>::result;
};

template<int a>
struct GCD<a, 0> {
  const static int result = a;
};

template<int x, int i>
struct TotientHelper {
  const static int result =
      (GCD<x, i>::result == 1) +
      TotientHelper<x, i - 1>::result;
};

template<int x>
struct TotientHelper<x, 0> {
  const static int result = 0;
};

template<int x>
struct Totient {
  const static int result =
      TotientHelper<x, x - 1>::result;
};

int main() {
  const int x = Totient<1000>::result;
  cout << x << endl;
}

Essa versão é bem mais difícil de entender e é cheia de boilerplate. A vantagem dela é que o cálculo do totiente é feito em tempo de compilação; então, em runtime, o resultado é calculado em O(1).

Por fim, a versão em c++0x usando constexpr:

constexpr int gcd(int a, int b) {
  return b == 0 ? a : gcd(b, a % b);
}

constexpr int totient(int x, int i) {
  return i == 0 ? 0 : (gcd(x, i) == 1) + totient(x, i - 1);
}

constexpr int totient(int x) {
  return totient(x, x - 1);
}

int main() {
  constexpr auto x = totient(NUMBER);
  cout << x << endl;
}

Essa versão não é tão simples quanto a que calcula em runtime, mas também não é excessivamente complexa como a versão com templates.

Para fazer metaprogramming com constexpr, você só precisa seguir algumas regras simples:
  • O corpo da função precisa ter um único return.
  • Você não pode usar if. Se você precisar de condicionais, tem que usar o operador ternário.
  • Você não pode usar for. Se você precisar de loops, tem que usar recursão.
  • Você não pode fazer atribuições a variáveis.
  • Uma pegadinha: ele só faz o cálculo em compile time se você atribuir o resultado a uma variável constexpr. Se você tentar fazer diretamente cout << totient(100), ele vai calcular em run-time.
Além disso, obviamente, você precisa de um compilador que suporte constexpr. Até onde eu sei, o único que aceita constexpr é o gcc 4.7 ou superior.

Mas a grande vantagem do constexpr sobre os templates não é só a facilidade de programação: é o tempo de compilação. Um programa com constexpr compila muito mais rápido que uma versão com templates. Para comparar o tempo de compilação, eu fiz um script que mede o tempo para vários valores do totiente:

Script em C++ para calcular o totiente em run-time
Script em C++ para calcular o totiente com templates
Script em C++0x para calcular o totiente com constexpr
Script em Python para medir o tempo dos programas acima

Vamos primeiro comparar o tempo de compilação da versão em runtime com a versão em template:



Como esperado, o tempo de compilação da versão em runtime é constante e bem pequeno. Já a versão em template é bem ruinzinha. Eu fiz uma regressão para determinar que curva é essa, e o resultado é que é uma parábola (com coeficiente de determinação R2=0.999, ou seja, com bastante certeza).

Agora a fraqueza do método com templates é evidente: o algoritmo era O(n log n) em runtime, virou O(n2) em compile-time. Veja como isso não acontece com a versão com constexpr:



Ah, agora sim! A versão com constexpr continua sendo O(n log n), só trocou o tempo de run-time pelo tempo de compilação. Comparado o template nem tem graça: para n=30000, a versão com constexpr compila em 1 segundo, e versão em template compila em 2 minutos.

A essa altura já estamos convencidos de que o constexpr é o melhor jeito de implementar metaprogramming, mas ainda falta a pergunta crucial: na prática, pra que isso serve?

Para mim, a utilidade mais imediata é pra conseguir primeiro lugar nos problemas do spoj! Como o spoj só leva em conta o tempo de run-time para ordenar as soluções, se você jogar um pouco do tempo do run-time para a compilação, consegue subir algumas colocações no ranking :)

domingo, 18 de abril de 2010

ZenoReduce

"O movimento é uma ilusão", dizia o sábio grego Ζήνων (que translitera para Zeno, Zenon ou Zenão). Ele tinha várias demonstrações dessa afirmação, e uma delas começava com uma tartaruga desafiando o ligeiro Aquiles: "duvido que você chegue aqui onde eu estou". Ora, Aquiles só precisa de uma rápida corrida para chegar onde a tartaruga está. Ou não?


Para ajudar, a tartaruga colocou no caminho um monte de bandeirinhas. A primeira bandeirinha está na metade do caminho entre Aquiles e a chegada. A segunda bandeirinha está na metade do caminho entre a primeira bandeirinha e a chegada, e assim por diante.


Olha só que curioso: quando o Aquiles estiver na primeira bandeirinha, ele ainda vai precisar passar pela segunda, que está no meio do caminho que falta. Quando ele chegar na segunda, ainda precisa passar pela terceira, e assim por diante. Não importa em qual bandeirinha ele esteja, sempre tem mais uma no meio do caminho restante. E se sempre falta uma bandeirinha, então ele nunca vai chegar na tartaruga! Daí a conclusão do Zeno, de que o movimento do Aquiles é só uma ilusão.

Para os gregos isso era um paradoxo, mas isso é porque eles desconheciam limites e séries convergentes. Hoje em dia nós sabemos que o tempo entre uma bandeirinha e outra sempre diminui, e eventualmente tende a zero. Nesse caso o paradoxo se dissolve, e a nossa intuição de que o Aquiles consegue chegar na tartaruga volta a ser verdadeira.

Agora imagine como seria divertido se tivéssemos na computação um conceito análogo! Por exemplo, para percorrer uma lista, um computador leva sempre o mesmo tempo em cada elemento. Mas em algum universo paralelo poderíamos ter um computador que, a cada elemento da lista, leva só metade do tempo do passo anterior. Ao invés de serem máquinas de Turing, nesse universo paralelo os computadores seriam máquinas de Zeno.

Na verdade, eu fiquei tão curioso em saber como seriam esses computadores que escrevi um emulador deles pra brincar:

Emulador de máquina de Zeno, escrito em python
Unit tests do emulador, usando pyunit e mox

A api do emulador fornece uma única operação: o ZenoReduce. Ele funciona como se fosse um reduce normal do python, mas o tempo que ele leva em cada iteração da lista é metade do tempo da iteração anterior. A complexidade do ZenoReduce é curiosa: ela é igual à complexidade da operação binária que você usa, e independente do tamanho da lista. A demonstração é simples e está na caixa azul:

Livros clássicos de complexidade, como o do Papadimitriou, calculam complexidade contando o número de passos que o algoritmo executa. Mas isso só vale para máquinas de Turing, onde os passos tem sempre a mesma duração. Na máquina de Zeno, temos que medir complexidade usando o tempo de execução, ao invés do número de passos.

Suponha que o ZenoReduce vai executar uma operação binária qualquer em cada elemento da lista, e essa operação tem uma complexidade O(f(x)). Logo, o primeiro passo do ZenoReduce leva O(f(x)), o segundo leva O(f(x))/2 e assim por diante. Assim, se a lista tem n elementos, então a complexidade total é:



Mas 2.O(f(x)) é o próprio O(f(x)), logo provamos que o valor de n realmente não importa.

Escrever programas usando o ZenoReduce é divertido, porque os melhores algoritmos para ele usualmente são os piores na máquina de Turing. Por exemplo, para verificar se um número é primo, podemos usar a versão abaixo:
import zeno

def is_prime(n):
  return zeno.zenoreduce(
    lambda x,y: x and (n%y!=0), 
    xrange(2,n), True)
Para rodar esse ZenoReduce, você precisa usar a função run da api, que liga o módulo de emulação. O run requer uma função sem parâmetros que será executada na máquina de Zeno:
>>> zeno.run(lambda:is_prime(10001))
(False, 0.15625, 0.0)
A saída é uma tupla: o primeiro elemento é o resultado (False, ou seja, 10001 não é primo), o segundo elemento é o tempo que ele levou pra rodar na sua máquina (wall time, 0.15s), e o terceiro é quanto tempo ele teria levado se tivesse rodado em uma máquina de Zeno (zero, ou melhor, um tempo tão pequeno que é menor que a precisão do python).

A complexidade desse algoritmo é O(mod), ou seja, quem domina é o tempo para calcular o resto de uma divisão. Nossos computadores modernos fazem divisão por hardware, então esse tempo é O(1), porém só até um certo limite. Quando o número é grande, você precisa usar uma lib de bignum, e a complexidade aumenta. Mas aí é só usar o ZenoReduce novamente!
def mod(x, n):
  return zeno.zenoreduce(
    lambda x,y: x-n if x>=n else x, 
    xrange(x), x)
Agora nós reduzimos o O(mod) para O(subtração), e podemos continuar até bater em alguma operação que seja O(1) para qualquer entrada!

(Leitores atentos podem achar que usar um ZenoReduce dentro de outro é trapaça, mas não é. Sempre que você consegue achar um limite superior para o tamanho de todas as listas envolvidas, você também consegue reescrever os ZenoReduces encadeados como um único ZenoReduce rodando sobre o produto cartesiano das listas. Provar isso fica como exercício para o leitor atento.)

Quão poderosas são as máquinas de Zeno? Certamente mais que as máquinas de Turing! Na verdade, para elas é facinho provar que P=NP. Vamos fazer isso mostrando uma solução em O(1) para o problema do Subset Sum: dado uma lista qualquer de inteiros, queremos saber se existe um subconjunto dela cuja soma seja zero (esse problema é NP-completo).

Uma estratégia seria gerar uma lista com todos os subconjuntos possíveis, e testar cada um deles em um ZenoReduce. Mas isso não é muito bom porque você precisaria armazenar os 2n subconjuntos, ocupando assim espaço exponencial.

Ao invés disso, nós podemos simplesmente contar de 0 até 2n-1, onde, para cada número, os bits dele indicam se o elemento correspondente deve ser somado ou não. Assim nós podemos iterar pelos subconjuntos sem precisar armazená-los.

Antes de começar, duas funções auxiliares. Em python, o built-in len() é O(1), mas se não fosse poderíamos usar isso aqui:
def list_length(x):
  return zeno.zenoreduce(
    lambda x,y: x+1, x, 0)
Podemos também fazer um função O(1) para potenciação:
def power(n, x):
  return zeno.zenoreduce(
    lambda x,y: x*n, xrange(x), 1)
Dados uma lista e um número, cujos bits representam os elementos da lista a serem somados, podemos fazer essa soma condicional em O(1):
def conditional_sum(elements, n):
  return zeno.zenoreduce(
    lambda x,y: 
      x + (y[1] if (n & power(2, y[0]) > 0) else 0), 
    enumerate(elements), 0)
Por fim, é só testar todos os subconjuntos possíveis, novamente em O(1):
def subset_sum(elements):
  return zeno.zenoreduce(
    lambda x,y: x or (conditional_sum(elements, y) == 0), 
    xrange(1, power(2, list_length(elements))), False)
Vamos testar:
>>> zeno.run(lambda:subset_sum([-7, -3, -2, 5, 8]))
(True, 0.015625, 0.0)
>>> zeno.run(lambda:subset_sum(range(1,15)))
(False, 16.75, 1.89e-09)
Veja como esse algoritmo é horrível na máquina de Turing e quase instantâneo na máquina de Zeno!

Outra coisa que não seria possível nesse universo paralelo é a criptografia RSA. De fato, dá pra fatorar um número em O(1) usando ZenoReduces!
def greatest_factor(value, n):
  return zeno.zenoreduce(
    lambda x,y: 
      max(x, y if value % power(n, y)==0 else x), 
    xrange(value), 0)

def factorize(n):
  return zeno.zenoreduce(
    lambda x,y: (x + [(y, greatest_factor(n, y))]) 
      if is_prime(y) and n%y==0 else x, 
    xrange(2, 1+n), [])
Testando:
>>> zeno.run(lambda:factorize(300))
([(2, 2), (3, 1), (5, 2)], 12.65625, 1.77e-15)
Nesse universo paralelo onde facilmente se prova que P=NP e fatoração é O(1), certamente a computação seria bem diferente! Mas o que você me diria se eu te falasse que podemos fazer algo parecido aqui no nosso universo mesmo?

Imagine que você tem algum cálculo qualquer que demoraria um milhão de anos pra rodar, mas que numa máquina de Zeno demoraria apenas um dia. Você pode fazer esse cálculo usando uma máquina de Turing também, e o truque não é nenhuma magia: é relatividade!

Basta deixar seu computador executando o cálculo, e pegar uma nave espacial que viaje a uma velocidade próxima da luz (qualquer coisa mais veloz que 0.999999999999999997c serve). Para você, vai passar apenas um dia, mas para o computador que ficou na Terra, vai passar um milhão de anos, graças ao efeito de dilatação do tempo. Fica como exercício para o leitor fazer uma nave tão rápida, e garantir que o computador ainda vai estar lá funcionando depois de um milhão de anos :)



Bônus: Ricbit Jam #2

Nós vimos como usar o ZenoReduce para testar primalidade, resolver o subset sum e fatorar um número em O(1). Você consegue fazer algo semelhante para ordenar uma lista? Se você topa o desafio, é só participar do Ricbit Jam #2:

Regras do Ricbit Jam #2

Dessa vez o Ricbit Jam vai dar um prêmio para o primeiro colocado: um Cubo de Rubik oficial do Google!


Participe para ganhar!

quinta-feira, 29 de maio de 2008

Python one-liners são Turing-complete

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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


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


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

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

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

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

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

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

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

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

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

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


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

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

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


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

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

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

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

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

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

Conversor para Python one-liner.

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