A Fórmula de Bhaskara

Um tempo atrás o Karlisson teve uma idéia muito boa: criar uma coletânea com os 1001 algoritmos que você deve implementar antes de morrer. O objetivo é fazer uma implementação de referência, em python, para todos os algoritmos clássicos, desde os simples como o Bubblesort até os complexos como o Edmonds-Karp.

Se você também quiser participar, é só fazer um fork do repositório no github. É um ótimo exercício para quem é iniciante, e para os veteranos é uma boa oportunidade de praticar o hábito de fazer code reviews. Já que no post anterior eu falava de parábolas, vamos aproveitar pra revisar um dos códigos do projeto: uma função que acha as raízes reais da equação de segundo grau:

def bhaskara(a, b, c):
    delta = b ** 2 - 4 * a * c
    if delta < 0:
        return None
    else:
        raizes = []
        m1 = math.sqrt(delta)
        r1 =(-b + m1) / (2 * a)
        raizes.append(r1)
        r2 =(-b - m1) / (2 * a)
        raizes.append(r2)
        return raizes

A função é uma implementação direta da fórmula:

Conferindo com a fórmula, o código parece correto. Porém, olhando com cuidado, ele tem três problemas!

Primeiro problema: API inconsistente

O problema mais imediato do código é que ele faz uma divisão por zero quando a é igual a zero. Mas isso é um bug de verdade? A intenção do código era resolver equações do segundo grau, mas, se a é nulo, então a equação não é de segundo grau, é de primeiro!

Eu interpreto isso como um problema na API. Se a idéia do código é resolver somente equações de grau 2, então, antes de sair fazendo conta, ele deveria verificar se as entradas são válidas (ou seja, se a é não-nulo). Por outro lado, se a intenção é resolver equações de grau 2 ou inferior, então ele poderia retornar [-c/b] quando a for zero. Nos dois casos, faltou sinalizar exatamente qual o propósito da API.

A API ainda tem uma segunda falha, que é o valor de retorno. O que exatamente ele vai retornar? Eu poderia, por exemplo, dizer que a função sempre retorna uma lista com as soluções existentes. Nesse caso, quando o delta é negativo, o correto seria retornar uma lista vazia, ao invés de retornar None.

Além disso, o que devemos retornar quando o delta é exatamente zero? As duas possibilidades são retornar uma lista com uma única raiz, [-b/2a], ou então retornar uma lista com duas raízes iguais, ou seja, [-b/2a,-b/2a]. Qual resposta faz mais sentido?

A escolha correta vem do Teorema Fundamental da Álgebra, que diz que o número de raízes complexas de uma equação polinomial é igual ao grau da equação. O motivo de considerarmos que algumas equações tem múltiplas raízes iguais é pra fazer esse teorema ser válido em todas as situações.

Mas o teorema só vale para raízes complexas! Como estamos trabalhando com raízes reais, então eu acho que o mais correto é retornar uma lista com um único elemento mesmo.

Você também poderia argumentar que o nome da função não é apropriado, porque está descrevendo a implementação (fórmula de Bhaskara) ao invés de descrever a interface (solução da equação quadrática). Mas aí nós chegamos no segundo problema...

Segundo problema: Bhaskara quem?

Certamente você ouviu na escola que a solução da equação quadrática chama-se fórmula de Bhaskara. Mas, curiosamente, é só no Brasil que a fórmula tem esse nome!

Pode conferir: a wikipedia em inglês nem cita o Bhaskara. Talvez na wikipedia em hindi? Nada de Bhaskara lá também, na Índia eles chamam a equação de fórmula de Sridhar Acharya. Nos outros países a fórmula nem tem nome, é simplesmente "equação quadrática" ou "fórmula a-b-c".

Na verdade, os babilônicos já sabiam resolver alguns tipos de equações quadráticas há pelo menos 4000 anos. Os gregos conheciam soluções geométricas. E a fórmula em si, é do Bhaskara mesmo?

Bem, uma maneira de resolver a dúvida é procurando exatamente o que ele escreveu. O texto mais importante do Bhaskara foi o livro que dedicou à filha, Lilavati. O original era em sânscrito, mas tem uma tradução em inglês online para quem quiser ler. Eu procurei pela solução da quadrática, e o mais próximo que achei foi a proposição 63:

Uma quantidade, incrementada ou decrementada de sua raiz quadrada multiplicada por algum número, é dada. Então adicione o quadrado de metade do multiplicador da raiz ao número dado, e extraia a raiz quadrada da soma. Adicione metade do multiplicador, se a diferença foi dada; ou subtraia, se a soma for dada. O quadrado do resultado é a quantidade procurada.

Em notação moderna, o que o Bhaskara escreveu foi:

Eu não acho que seja correto creditar ao Bhaskara a solução da equação quadrática. A fórmula acima é bem parecida com a que usamos, mas note o detalhe: ela só fornece uma raiz! Para considerar a solução como correta, eu acho que ele deveria ter explicitado que algumas equações vão ter duas raízes.

Pois bem, podemos então mudar o nome da função, de bhaskara() para algo mais apropriado, como quadratic_roots(). Mas ainda resta o último problema...

Terceiro problema: Instabilidade numérica

Vamos testar o código original com a seguinte equação: x^2-2000000000x+4=0. As raízes são aproxidamadamente 2e9 e 2e-9. Porém, olhe só o resultado:

>>> bhaskara.bhaskara(1, -2e9, 4)
[2000000000.0, 0.0]

Como assim zero? Cadê a solução 2e-9? Essa raiz foi vítima de um problema que nem todo mundo presta atenção: muito cuidado com a perda de precisão em ponto flutuante. Em especial, quando você subtrai dois números de magnitude similar, a resposta pode não ter precisão suficiente para dar o resultado correto. Para a equação dada, o código original tem os valores b e m1 muito próximos, e a subtração perde a menor raiz.

A solução é evitar a subtração. Você inspeciona o sinal de b, e escolhe m1 com o mesmo sinal, de modo que os módulos são sempre somados, achando assim a primeira raiz x1. Tendo x1, você pode achar x2 usando a fórmula de Viète para o produto de raízes:

Dessa maneira a subtração problemática some, e chegamos na versão final do programa corrigido:

def quadratic_roots(a, b, c):
  """Find all real roots of a quadratic equation.

  Args:
    a, b, c: Coefficients of the quadratic, a*x*x+b*x+c=0
  Returns:
    A list with all real roots.
  Raises:
    ValueError if a==0 (the equation is not quadratic).
  """
  if not a:
    raise ValueError
  delta = b * b - 4 * a * c
  if delta < 0:
    return []
  if not delta:
    return [-b / (2 * a)]
  if b < 0:
    x1 = (-b + delta ** 0.5) / (2 * a)
  else:
    x1 = (-b - delta ** 0.5) / (2 * a)
  x2 = c / (a * x1)
  return [x1, x2]

Conferindo:

>>> roots.quadratic_roots(1, -2e9, 4)
[2000000000.0, 2.0000000000000001e-09]

Nice!

Você pode se perguntar se esse exemplo que eu dei não é artificial demais, e se coisas assim acontecem na prática. Mas acontecem sim! Um exemplo comum onde esse erro acontece é na implementação da busca binária.

Por exemplo, suponha que você quer resolver a equação ln(x)+sqrt(x)=5. A função ln() é monotônica, o sqrt() também, então dá pra achar a raiz por busca binária:

def search(epsilon):
  inf, sup = 0.0, 1e10
  while sup - inf > epsilon:
    med = (sup + inf) / 2
    x = math.log(med) + math.sqrt(med)
    if x > 5.0:
      sup = med
    else:
      inf = med
  return inf

Vamos testar essa rotina para várias precisões:

>>> binary.search(1e-4)
8.3093709690729156
>>> binary.search(1e-10)
8.3094326941933723
>>> binary.search(1e-15)

Oops! Quando epsilon vale 1e-15, a rotina trava! De novo, o problema é perda de precisão. Nesse caso, as variáveis sup e inf tem valores muito próximos, então a média med não sai do lugar, ficando presa no valor do inf para sempre.

A solução, novamente, é sumir com a subtração. Ao invés de controlar os extremos do intervalo, você usa uma variável para o valor inicial e outra para o tamanho do intervalo:

def smart_search(epsilon):
  inf, delta = 0.0, 1e10
  while delta > epsilon:
    delta /= 2
    med = inf + delta
    x = math.log(med) + math.sqrt(med)
    if x < 5.0:
      inf += delta
  return inf

Conferindo:

>>> binary.smart_search(1e-4)
8.3093709690729156
>>> binary.smart_search(1e-10)
8.309432694193374
>>> binary.smart_search(1e-15)
8.3094326942315693
>>> binary.smart_search(1e-30)
8.3094326942315693

Agora sim!

No fim das contas, é importante você ter a matemática sempre afiada, mas também é igualmente importante conhecer os limites impostos pela computação do mundo real.

(Agradecimentos ao Karlisson, ao Henrique, ao Wladimir, ao Chester e ao povo da lista eng-misc pela ajuda na pesquisa!)

A Terra Plana

Todo mundo já ouviu a história de que foi Colombo que confirmou que a Terra era redonda, e que antes dele as pessoas achavam que a Terra era plana. Porém, isso é um mito! Desde a época dos gregos já se sabia que a Terra é redonda.

Mesmo para os povos antigos as evidências eram fáceis de serem encontradas. Por indução, se o Sol é redondo e a Lua é redonda, por que a Terra não seria? Mais diretamente, durante um eclipse lunar, é possível ver claramente que a sombra projetada pela Terra é redonda.

Nem mesmo a Igreja disputava essa interpretação. Santo Agostinho, por exemplo, sabia que a Terra era esférica, embora não acreditasse que houvesse gente morando do outro lado do globo. São Tomás de Aquino começa a sua Summa Theologica citando duas diferentes demonstrações de que a Terra é redonda.

"Mas e aquela história de que os marinheiros da época tinham medo de cair pela borda da Terra?" Bem, para quem não entende como a gravidade funciona, é possível ter esse medo mesmo com a Terra sendo redonda! Se a gravidade puxa sempre para baixo, então você cai quando chega no Equador:

Aparentemente, todo esse mito começou com a publicação de "As Vida e as Viagens de Cristóvão Colombo", que você pode ler integralmente no Google Books. O livro realmente cita que o senso comum era de que a Terra era plana. O problema é que todo mundo achou que esse era um livro de história, quando na verdade era ficção! (Concluí-se, portanto, que o autor desse livro foi um precursor do The Onion).

Hoje em dia ainda há quem acredite que a Terra é plana, desde os ingênuos que nunca tiveram acesso à educação, até os malucos que formam seitas como a Sociedade da Terra Plana. Talvez você conheça alguém assim. Na verdade, é bem possível que você acredite na Terra Plana. Sim, você que está lendo esse texto!

Se você acha que não acredita na Terra plana, então tente responder o meu teste: Suponha que você está usando um estilingue para arremessar passarinhos em porquinhos. Na ausência de atrito com o ar e outras formas de perdas dissipativas, qual a curva que o passarinho descreve em seu trajeto?

Se você respondeu que a curva é uma parábola, surpresa! Você acabou de assumir que a Terra é plana!

"Mas peraí! Como assim?! Foi isso que meu professor ensinou!". Pois é, sempre é bom ter uma dose de ceticismo quando te ensinam alguma coisa. Na verdade, a prova de que a trajetória é uma parábola requer a existência de uma força gravitacional do tipo F=mg, onde g é uma constante de aproxidamente 9.81 m/s².

A pegadinha é que esse g não é uma constante! Ele varia com a altura, e, levando isso em conta, a curva não é mais uma parábola. Porém, se você assumir que a Terra é um plano infinito, aí sim a curva é uma parábola, porque nesse caso g é uma constante independente da altura.

E como provar que o g é realmente constante para uma Terra plana? Existem várias maneiras: se você souber um pouco de cálculo, pode usar o método do guerreiro, que envolve várias páginas de lutas contra integrais duplas. Se você for um mestre do cálculo, pode tentar o método do mago, que resolve com poucas contas, mas requer o uso de teoremas mágicos. Agora, se você não souber cálculo, o jeito é apelar para método do clérigo: acredite em mim e pule o quadro azul abaixo :)

O Método do Guerreiro

Antes de começar, convém lembrar a lei da gravidade de Newton, que é proporcional às massas e inversamente quadrática com a distância entre elas:

Como a lei tem simetria radial, as contas ficam mais fáceis se explorarmos essa simetria. Vamos calcular qual é a atração gravitacional causada por um disco de raio R e espessura infinitesimal, composto de materia uniforme com densidade σ, sobre um passarinho de massa m_P que esta a uma altura h.

Nós vamos picotar o disco em pedacinhos de massa dm e integrar sobre todos os pedacinhos. Como as dimensões de cada pedacinho são infinitesimais, então eles podem ser considerados pequenos retângulos, onde um dos lados vale dr, e o outro vale um pequeno comprimento de arco. Um comprimento de arco é igual ao raio vezes o ângulo, logo esse lado vale r.dϕ . A massa do pedacinho, então, é a densidade vezes a área: σ.r.dr.dϕ.

Além disso, se eu tenho o raio r e a altura h, por Pitágoras eu sei quanto vale d². A força total sobre o pedacinho é:

Antes de continuar, vamos decompor a força total em componentes radial (F_r) e vertical (F_h). Agora a simetria do problema se mostra útil, basta notar que para cada dm, tem um outro dm numa posição simétrica, do outro lado do disco; e nesse dm espelhado, o F_r vai ter o mesmo módulo e direção, mas sinal oposto. Então todos os F_r vão cancelar, e a força total sobre o passarinho vai ser apenas a soma dos F_h de cada dm.

E quanto vale o F_h? É a força total F vezes o cos θ, que por sua vez é o lado adjacente pela hipotenusa. Logo:

Para calcular o F_h total, é só integrar para o ϕ variando de 0 a 2π, e para o r variando de 0 a R:

Todo mundo é constante em relação a ϕ, então dá pra jogar pra fora, e a integral em dϕ fica simples:

Para continuar, uma mudança de variável:

O meio joga pra fora e cancela com o dois, o resto é um monômio simples, para integrar você soma um no expoente e passa dividindo:

A ultima expressão é a fórmula final da força gerada por um disco uniforme. Por fim, basta notar que o plano é o limite do disco quando o raio vai para o infinito:

Se o R vai para o infinito, então a primeira fração vai para zero, logo:

Nós estamos interessados no caso em que h é positivo, logo a fração simplifica e chegamos na forma final da força devida a um plano:

Como esperado, a força F não depende da altura h!

O Método do Mago

Ao invés de usar a lei da gravidade de Newton, podemos usar a lei da gravidade de Gauss, que diz que a integral do campo gravitacional sobre uma superfície fechada é proporcional a massa contida dentro dela:

Essas duas leis são equivalentes, pode-se usar uma para provar a outra. Mas em alguns casos a lei de Gauss é mais prática, especialmente se você puder bolar uma superfície que explore as simetrias do problema. No nosso caso, a superfície mágica é um cilindro centrado no plano, com as bases paralelas a ele:

Antes de mais nada precisamos da massa dentro da superfície. Vamos chamar a área da base do cilindro de A. Como o plano tem densidade σ, entao a massa dentro da superfície fechada é M = A.σ.

Nós já vimos antes que, pela simetria radial, todas as forças resultantes são perpendiculares ao plano. Na lateral do cilindro, todos os elementos de área dA são normais à lateral, logo são perpendiculares às forças, e portanto o produto escalar deles dá zero.

Além disso, também há uma simetria em relação à translações no plano, por isso todas as forças atuando a uma mesma altura são iguais entre si. Dessa maneira, a integral de superfície em cada base do cilindro vale esse valor constante (chamemos de F) vezes a área A.

Juntando tudo, temos:

Esse é o mesmo resultado do método anterior, mas com bem menos contas!

Mas se na Terra redonda a curva não é uma parábola, então que curva ela é? Bem, se você assumir que a Terra é uma esfera perfeita, homogênea, sem campos magnéticos, e assumir que não há atrito com ar e que ninguém está girando, então a curva correta é uma elipse.

Para quem leu o Principia, este é um resultado intuitivo. Na Proposição 71 do Livro 1, o Newton demonstra que a força resultante da atração por uma esfera homogênea é a mesma força de quando você concentra toda a massa da esfera num ponto material em seu centro. Daí o problema se reduz à atração entre dois pontos materiais, que sabemos ser uma elipse.

Esse raciocínio também serve para responder a uma dúvida comum da criançada: se a gravidade da Terra puxa tudo para o centro, por que os satélites entram em órbita ao invés de cair? A resposta mais imediata é que na verdade tudo que você joga pra cima entra em órbita, mas algumas órbitas são tão excêntricas que o objeto bate de volta na superfície da Terra.

Eu não gosto muito da maneira que ensinam Física no Brasil (o Feynman também não). Felizmente, hoje em dia você sempre pode remediar a sua formação usando os recursos online, como a Wikipedia e a Khan Academy!

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 R²=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(n²) 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 :)

Torto Reverso

No fim do ano passado o meu Nintendo DS quebrou, após anos de bons serviços prestados. Eu pensei em consertar, mas dado que o Nintendo 3DS ia sair em março, achei que valia a pena simplesmente esperar. Porém, o que fazer nesse meio tempo?

A minha solução foi ir à banca e comprar um punhado de cruzadinhas da Coquetel. No começo foi divertido, mas depois de algum tempo perdeu a graça: para quem está acostumado com a Denksel, os puzzles da Coquetel são fáceis demais.

Para aumentar a dificuldade dos puzzles, eu resolvi mudar as regras de alguns deles. A minha modificação predileta eu batizei de Torto Reverso. Antes de mostrar como funciona a minha variação, vamos relembrar as regras do Torto tradicional, tais como estão na revista. É dado um grid 3x6 de letras e a seguinte descrição:

"Deve-se formar as palavras seguindo em todas as direções, sempre ligando as letras em seqüência direta, sem cruzar, sem pular e sem repetir a mesma letra (para que uma palavra tenha a mesma letra repetida, é necessário que esta letra também esteja duplicada no diagrama). Damos como exemplo uma palavra encontrada. Só valem palavras de cinco letras ou mais. Na solução, constam 30 palavras formadas com este diagrama."

Se você tiver um bom vocabulário, esse puzzle não tem segredo. Muito mais difícil é o Torto Reverso: e se, ao invés de procurar palavras no grid, eu te der as palavras e pedir para achar um grid válido que contenha todas elas?

Depois de fazer alguns desses na mão, eu fiquei com vontade de implementar um solver. Eu poderia ter feito um solver em C++ ou Python, mas pensei aqui comigo: "eis uma boa oportunidade para aprender Prolog!"

Se você também não conhece Prolog, então acompanhe um pequeno tutorial antes de ver minha implementação do solver:

Introdução ao Prolog

Em linguagens de programação tradicionais, a metáfora que você usa para explicar é a da receita de bolo: programar é fazer uma seqüência de tarefas, como adicionar 200g de farinha e mexer até que a massa esteja consistente. Prolog usa um paradigma diferente, então a melhor metáfora é pensar que ele é um provador de teoremas: você entra com os axiomas e ele deriva os teoremas a partir deles.

Por exemplo, podemos usar o Prolog para definir os axiomas de Peano. Começamos com a definição de números naturais: zero é um número natural, e o sucessor de um natural também é um natural.

natural(0).
natural(s(X)) :- natural(X).

Com essas definições já podemos perguntar ao Prolog: será que 2, ou seja, o sucessor do sucessor do zero, é um número natural?

?- natural(s(s(0))).
true.

Bacana! E o sucessor de Togakure, é um número natural?

?- natural(s(togakure)).
false.

Está certo, afinal, se Togakure não é um número natural então o sucessor também não pode ser. Vamos agora definir a adição: zero é o elemento neutro, e se você somar qualquer coisa com o sucessor de outra coisa, o resultado é o sucessor da soma das duas coisas.

add(A, 0, A).
add(A, s(B), s(C)) :- add(A, B, C).

Podemos perguntar agora: será que 1+1=2? Será que 1+1=3?

?- add(s(0), s(0), s(s(0))).
true .
?- add(s(0), s(0), s(s(s(0)))).
false.

Até agora só usamos o Prolog para testar a veracidade de sentenças, mas a diversão começa quando usamos o Prolog para completar sentenças. Sempre que você deixar uma variável livre numa sentença, o Prolog tenta completar a sentença para torná-la verdadeira. Por exemplo, quanto é 2+2?

?- add(s(s(0)), s(s(0)), X).
X = s(s(s(s(0)))) .

Talvez a parte mais bacana do Prolog seja o fato da linguagem ser reversível. O mesmo programa que avalia a soma também é capaz de resolver uma equação. Por exemplo, qual é o número que somado com 2 dá 5?

?- add(X, s(s(0)), s(s(s(s(s(0)))))).
X = s(s(s(0))) .

E qual é o número que somado com 1 dá 0?

?- add(X, s(0), 0).
false.

O Prolog respondeu false, ou seja, no conjunto dos naturais não tem nenhum número que somado com 1 dá 0. Outra coisa bacana é que você pode deixar mais de uma variável livre. Por exemplo, quais são os pares de números cuja soma é 2?

?- add(A, B, s(s(0))).
A = s(s(0)),
B = 0 ;
A = s(0),
B = s(0) ;
A = 0,
B = s(s(0)) ;

Ou seja, as soluções possíveis são 2+0, 1+1 e 0+2.

Resolvendo o Torto Reverso

Agora que sabemos o básico do Prolog, já temos o material necessário para resolver o puzzle! Note que eu usei como interpretador o SWI-Prolog, porque ele tem a melhor biblioteca de built-ins dos Prologs que eu pesquisei. Tudo que estiver em verde são palavras reservadas do SWI-Prolog.

Em linguagens tradicionais você precisa descrever, passo a passo, um método para encontrar a solução do seu problema. Em Prolog é o contrário, você só precisa falar quais são as propriedades da sua solução, e o interpretador se vira para achar um método de solucionar.

Por exemplo, nós podemos criar um predicado que é verdadeiro se o nosso grid é uma lista com 18 elementos. Mas como Prolog é reversível, se você deixar a variável do grid livre, ele vai criar um grid para você!

validgrid(Grid) :-
  length(Grid, 18).

?- validgrid(X).
X = [_G282, _G285, _G288, _G291, _G294,|...].

O Prolog não tinha informação suficiente para deduzir quem eram os elementos do grid, então colocou placeholders no lugar (_G282, _G285, etc).

Note também que o Prolog não tem arrays bidimensionais, só listas. Por isso, vamos precisar de um predicado para checar se os índices do grid são válidos, e um para indexar o grid a partir dos índices:

valid(I,J) :-
  between(0, 2, I),
  between(0, 5, J).

validletter(Grid, Letter, X, Y) :-
  valid(X, Y),
  Pos is X + Y * 3,
  nth0(Pos, Grid, Letter).

Nesse trecho o importante é entender o comando is. Para o Prolog, 1+1*3 é diferente de 4, porque ele não faz a avaliação aritmética por default. Usando o comando is, você o força a avaliar a expressão.

Agora vamos definir a conectividade do grid. Dada uma posição válida qualquer, nós podemos nos mover para qualquer das oito direções, desde que você não saia do grid. Para isso, podemos definir as oito possibilidades para um predicado move:

move(I, J, left, X, J) :-
  X is I - 1,
  valid(X, J).

(...)

move(I, J, down_right, X, Y) :-
  X is I + 1,
  Y is J + 1,
  valid(X, Y).

Assim, já podemos perguntar ao Prolog: saindo da posição 0,0, quais são as posições atingíveis?

?- move(0, 0, Direction, X, Y).
Direction = right,
X = 1,
Y = 0 ;
Direction = down,
X = 0,
Y = 1 ;
Direction = down_right,
X = 1,
Y = 1.

Pronto, já temos as ferramentas básicas. Vamos descrever a nossa solução então. Nós queremos colocar uma palavra no grid, então quais as propriedades dessa palavra? Se a palavra tiver uma letra só, ela pode ser encaixada em qualquer posição do grid:

word(Grid, [Letter], X, Y) :-
  validletter(Grid, Letter, X, Y).

Agora, se ela tem várias letras, então você encaixa a primeira letra e faz a recursão para o resto da palavra, notando que o resto começa em um ponto que é um movimento válido a partir da primeira letra:

word(Grid, [Letter | Tail], X, Y) :-
  validletter(Grid, Letter, X, Y),
  move(X, Y, _, I, J),
  word(Grid, Tail, I, J).

Aqui temos mais dois operadores do Prolog: o operador pipe separa uma lista em head e tail, e o operador underscore significa "don't care", nós não ligamos para o que o Prolog vai encaixar ali.

Isso já é o suficiente para encaixar uma palavra no grid, por exemplo, a partir da posição 0,0:

solve(Word) :-
  word(Grid, Word, 0, 0),
  writegrid(Grid).

?- solve("torto").
tor
.ot
...
...
...
...

Eu pedi para imprimir só a primeira solução, mas nós poderíamos ter pedido para o Prolog imprimir todas as soluções. (Eu pulei a definição do writegrid, porque não é tão interessante.)

Antes de encaixar o resto das palavras, precisamos adicionar o restante das regras. Uma regra que não lidamos até agora é que uma palavra não pode repetir uma posição do grid durante o trajeto.

Para implementar essa restrição, podemos criar uma lista de posições já visitadas, e só permitir a inserção de uma letra se a posição dela ainda não estiver na lista. Vamos modificar nosso predicado word então:

word(Grid, [Letter], Visited, X, Y) :-
  validletter(Grid, Letter, X, Y),
  \+ member(pos(X, Y), Visited).

word(Grid, [Letter | Tail], Visited, X, Y) :-
  validletter(Grid, Letter, X, Y),
  move(X, Y, _, I, J),
  \+ member(pos(X, Y), Visited),
  word(Grid, Tail, [pos(X, Y) | Visited], I, J).

A novidade aqui é o operador \+, que é simplesmente um "not". Ou seja, insira a nova letra somente se a posição dela não estiver na lista Visited.

Mas ainda tem outra regra: a palavra não pode cruzar o próprio caminho. Na horizontal e na vertical isso não tem como acontecer, porque para isso ele teria que repetir uma posição. Mas pode acontecer na diagonal, como na figura abaixo:

Etna não é um vulcão válido.

Para evitar esse caso, podemos fazer uma nova lista. Sempre que passarmos na diagonal por um bloco delimitado por letras, adicionamos esse bloco na lista.

Agora, para checar se andamos na diagonal, poderíamos verificar a direção que o Prolog deduziu no predicado move. Ao invés disso, vamos fazer de um jeito diferente:

validBlock(X, _, X, _, Block, Block) :- !.
validBlock(_, Y, _, Y, Block, Block) :- !.
validBlock(X, Y, I, J, OldBlock, NewBlock) :-
  Bx is min(X, I),
  By is min(Y, J),
  \+ member(pos(Bx, By), OldBlock),
  NewBlock = [pos(Bx, By) | OldBlock].

Se um movimento repetir o X, ou se repetir o Y, então não é uma diagonal, e não precisamos adicionar nada na lista. Para isso, usamos o operador exclamação, que implementa um corte.

A idéia é simples: um predicado do tipo validBlock(0, 0, 0, 0, A, B) iria fazer match nas três instâncias declaradas, e por isso o Prolog iria deduzir três soluções para esse predicado. Mas, com o corte, ele desiste de procurar mais soluções após encontrar a primeira.

Note também que, para definir unicamente um bloco, independente da direção da diagonal, basta pegar o mínimo das posições X e Y.

Podemos agora incorporar esse predicado para chegar na versão final do word:

word(Grid, [Letter], Visited, _, X, Y) :-
  validletter(Grid, Letter, X, Y),
  \+ member(pos(X, Y), Visited).

word(Grid, [Letter | Tail], Visited, Block, X, Y) :-
  validletter(Grid, Letter, X, Y),
  move(X, Y, _, I, J),
  \+ member(pos(X, Y), Visited),
  validBlock(X, Y, I, J, Block, NewBlock),
  word(Grid,Tail,[pos(X, Y) | Visited],NewBlock,I,J).

Isso já encaixa qualquer palavra seguindo todas as regras! Agora é só falar para inserir uma lista de palavras, ao invés de uma palavra só, e terminamos nosso solver:

wordlist(_, []).
wordlist(Grid, [Word | WordList]) :-
  word(Grid, Word, [], [], _, _),
  wordlist(Grid, WordList).

solve(WordList) :-
  validgrid(Grid),
  wordlist(Grid, WordList),
  writegrid(Grid).

?- solve(["otica","dedal","tedio","cardeal","aldeia"]).
oti
dac
eda
tel
idr
oac

O solver em Prolog ficou curto, né? Note que nós descrevemos só a solução, e não o método de solução. Em outras linguagens você precisaria se preocupar com Backtracking, A*, Exact Cover e o escambau, enquanto que no Prolog isso é feito por baixo dos panos para você.

Se você quiser rodar o solver em casa, aqui está o fonte completo:

Solução do Torto Reverso em Prolog

A desvantagem de escrever em Prolog é que fica um pouco mais lento que as linguagens tradicionais. Com essa implementação mostrada, temos 30 palavras, onde cada uma tem uma posição inicial definida num grid de 18 posições, e no mínimo 5 letras, onde cada letra pode estar em 8 posições possíveis a partir da anterior.

Ou seja, o espaço de busca é de (18*8⁴)³⁰, o que dá um total de 106898987586074109109365096428568366376512477689679588996632532409238880450162648732507836110780231551768603122161241934850249808547434035647873024 grids possiveis. Daí, usar esse solver em Prolog para um grid com 30 palavras pode demorar um bocadinho.

Ricbit Jam #3

Você acha que consegue fazer melhor que o Prolog usando sua linguagem predileta? Então participe do Ricbit Jam! Quem conseguir fazer o solver mais rápido do Torto Reverso vai ganhar de brinde um Boliche de Dedo oficial do Google:

Boliche de dedo do Google

Para participar, é só ler as regras abaixo:

Regras do Ricbit Jam #3

Divirta-se, e boa sorte!

As Três Leis de Newton

Desde pequeno, eu sempre adorei ler textos sobre a história da ciência. Quando criança eu lia a Superinteressante, na época em que era boa, e mais recentemente eu lia os bons livros de divulgação, como o Drunkard's Walk do Mlodinow. Mas, em algum ponto no ano passado, caiu uma ficha aqui. Por que eu estava lendo informação de segunda mão sobre probabilidade, se podia ler diretamente os originais do Laplace?

O maior problema em ler textos originais é a língua. Afinal, o Laplace escrevia em francês, o Euclides em grego, o Einstein em alemão, e eu não domino nenhuma dessas línguas. A boa notícia é que praticamente todos os textos científicos relevantes já foram traduzidos para o inglês em algum ponto. A má notícia é que muitas vezes esses textos estão fora de catálogo e nunca foram digitalizados.

Mas como tipicamente acontece comigo, curiosidade se mistura com obsessão, e lá fui eu fui fazer arqueologia nos sebos da Califórnia, conseguindo no final uma lista bem considerável de textos originais! (Pra quem estiver curioso, eu coloquei online a lista com os textos que já consegui.)

Para dar um gostinho de como é divertido ler os originais, eu resolvi traduzir um trecho sobre um tema que todo mundo já viu no colégio: as Três Leis de Newton, tais como enunciadas em 1687. Mas antes de falar das leis em si, vale a pena falar um pouco sobre o livro em que elas aparecem, o Principia.

Ler o Principia é bastante diferente de ler um livro atual sobre Física. Em especial, tem duas diferenças que são bem gritantes. A primeira é que o Principia não tem equações. Nenhuma, nenhuma. Todos os cálculos são feitos com desenho geométrico.

Se você parar pra pensar, nem é tão diferente do que a gente faz normalmente. Antes de montar as equações em um problema de Física, você usualmente vai desenhar os vetores com as forças; o Newton simplesmente desenhava os vetores em escala e usava geometria para resolver o problema.

O motivo para usar desenho geométrico ao invés de equações só ele sabia; mas dado que ele era fã dos Elementos de Euclides, pode-se supor que ele queria deixar o Principia mais parecido com os tratados gregos.

A segunda diferença do Principia para os livros atuais é que ele não tem citações. Nas raras vezes que o Newton cita algum outro cientista, é sempre um cientista que já tinha morrido, como o Aristóteles ou o Galileu; e nunca um contemporâneo como o Leibniz ou o Hooke.

O motivo, claro, era para não dividir os créditos das descobertas com ninguém. Na esfera pessoal, não há dúvidas de que o Newton era um pilantra (mas, por outro lado, ele era o pilantra mais inteligente de sua época, o que faz dele uma espécie de Lex Luthor do século XVII).

Isaac Newton sonegou mais de 40 citações. E isso é terrível.

O Principia originalmente era dividido em três livros, mas para fins de exposição eu acho mais fácil dividi-lo nas seções abaixo. Preste atenção na ordem delas, porque tem uma pegadinha mais tarde:

1. Préfacios e Introduções
2. Definições
3. Axiomas, ou as Leis do Movimento
4. Introdução ao Cálculo
5. O Movimento dos Corpos
6. O Sistema do Mundo
7. Conclusão

A minha cópia do Principia é uma tradução recente, feita em 1999. Mas se você quiser acompanhar sem gastar dinheiro, no Google Books tem o scan da primeira edição em latim, e da primeira tradução para o inglês, feita em 1729 (se a fonte parecer esquisita demais, note que aquele símbolo que parece uma integral na verdade é um S longo).

Eu poderia falar horas sobre o Principia, mas já que o objetivo é mostrar as três leis eu vou cobrir só duas seções: as definições e os axiomas.

Definições

Um problema que o Newton tinha era evitar a ambigüidade que alguns termos de uso comum possuem. Espaço, tempo, reta, plano, essas coisas são todas bem definidas. Mas não havia um conceito correto para medir matéria e movimento.

Informalmente, media-se matéria através do peso, e movimento a partir da velocidade, mas o Newton sacou que essas medidas não são as melhores possíveis. Por isso, ele apresenta dois conceitos novos que serão mais úteis. (Tudo que estiver na caixa azul é texto traduzido do Principia):

Quantidade de matéria é uma medida da matéria que surge conjuntamente da sua densidade e de seu volume.

Se a densidade do ar dobrar em um espaço que também dobra, então existe quatro vezes mais ar, e seis vezes mais ar se o espaço triplicar. O caso é o mesmo tanto para neve e pó condensado por compressão ou liquefação, como para quaisquer outros corpos condensados de qualquer maneira por qualquer método. Por enquanto, eu não estou levando em conta o meio, se existir algum, que penetra pelos interstícios entre as partes do corpo. Além disso, eu sempre tenho essa quantidade em mente quando usar o termo "corpo" ou "massa" daqui em diante. Ela sempre pode ser deduzida a partir do peso de um corpo (por exemplo, fazendo experimentos precisos com pêndulos), e eu descobri que ela é sempre proporcional ao peso.

Nessa definição, o Newton apresenta pela primeira vez o conceito de massa. Ele já sabia que peso não é uma boa medida da matéria, afinal, ele muda de acordo com o ponto em que você está (Newton sabia que um corpo em Paris tinha um peso levemente diferente do que teria se estivesse em Londres). Além disso, faz sentido falar do peso da Lua ou de Júpiter?

Note também como o Newton apresenta suas equações de modo implícito. Primeiro ele fala que a massa é uma função da densidade e do volume, mas não fala que função é essa. Depois, pelos exemplos, fica claro que a função é uma multiplicação, então massa é igual a densidade vezes volume.

Quantidade de movimento é uma medida do movimento que surge conjuntamente da velocidade e da quantidade de matéria.

O movimento como um todo é a soma dos movimentos das partes individuais, e, portanto, se um corpo é duas vezes maior que outro, mas tem a mesma velocidade, então há duas vezes mais movimento; e se há o dobro da velocidade, então há quatro vezes o movimento.

É mais difícil empurrar um caminhão a 10km/h que uma bicicleta, por isso o Newton sacou que o movimento não é simplesmente a velocidade, mas sim a velocidade vezes a massa. Em terminologia moderna, o Newton nesse trecho está definindo o momento.

Axiomas, ou as Leis do Movimento

Esse é o trecho que interessa, onde ele define as três leis. Novamente, o formato é baseado nos Elementos de Euclides, você começa o livro com os axiomas e depois deriva todo o resto a partir deles.

Primeira lei: Todo corpo persevera em seu estado de repouso, ou de movimento uniforme e retilíneo, exceto se for compelido a mudar de estado por forças aplicadas.

Projéteis perseveram em seus movimentos, exceto se forem retardados pela resistência do ar, e impelidos para baixo pela força da gravidade. Um pião, que possui partes que impedem umas às outras de realizar movimento retilíneos por razão de sua coesão, não deixa de rodar, exceto se for retardado pelo ar. E corpos maiores, como planetas e cometas, preservam por mais tempo seus movimentos progressivos e circulares, pois eles acontecem em espaços com menor resistência.

A parte curiosa dessa explicação é que ele cria uma lei sobre movimento retílineo, mas todos os exemplos envolvem movimentos curvos :) Não sei se a parte do pião que gira é clara, o que ele quis dizer é que cada partícula individual do pião tem a tendência de sair pela tangente, e só não faz isso porque as partículas estão grudadas entre si.

Mas eu acho que a primeira lei, enquanto axioma, não é das melhores. O que essa lei diz é que, na ausência de forças, a velocidade não muda, ou seja, na ausência de forças, a aceleração é zero. Mas você não precisava de uma lei só pra isso, porque a segunda lei já vai dizer que F=ma. Se a massa é positiva e a força é zero, então a aceleração vai ser zero. A primeira lei é redundante.

Segunda lei: Uma mudança no movimento é proporcional à força motriz aplicada, e tem lugar na mesma linha reta em que a força é aplicada.

Se alguma força gera qualquer movimento, então o dobro da força gera o dobro do movimento, e o triplo da força gera o triplo do movimento, não importando se a força é aplicada de uma só vez ou sucessivamente em passos. E se o corpo já estava se movendo previamente, o novo movimento (já que o movimento é sempre na mesma direção da força motriz) é adicionado ao movimento original se ele está mesma direção, ou subtraído do movimento original se estiver na direção oposta, ou, se for numa direção oblíqua, é combinado obliquamente e composto de acordo com as direções de ambos os movimentos.

Na época de faculdade a gente sempre pregava peça nos calouros. A conversa começava com "Ei bixo, qual é a segunda lei de Newton?". O calouro respondia que era F=ma e aí você sacaneava ele: "Claro que não, F=ma só vale se a massa for constante! O que o Newton disse é que a força é a derivada do momento, e como momento é massa vezes velocidade, então a equação correta é assim:

Se você projetar um avião usando só F=ma, ele vai cair, porque você não levou em conta que ele perde massa enquanto queima o combustível!"

Esse argumento é bom pra sacanear calouro, mas é incorreto. O Newton nunca falou que a força é a derivada do momento. Lembra que eu disse que tinha uma pegadinha na ordem dos capítulos? Nesse capítulo o Newton não poderia ter falado isso, porque ele mostra o Cálculo pela primeira vez só no capítulo seguinte, e portanto não tem como definir força como sendo uma "derivada".

Na verdade, a segunda lei de Newton tal como está no Principia se aplica ao caso discreto, para diferenças de movimentos e não para as suas derivadas. O que engana aqui é o jargão: aquilo que o Newton chama de força não é o mesmo que nós chamamos de força. Ao invés disso, esse conceito apresentado na segunda lei é o que modernamente conhecemos como Impulso. Ou seja, a verdadeira segunda lei é:

Isso fica claro no trecho "não importando se a força é aplicada de uma só vez ou sucessivamente em passos", claramente ele está falando de uma aplicação impulsiva ao invés de uma aplicação contínua. Naturalmente, se você dividir dos dois lados por Δt e tomar o limite quando Δt tende a zero, você chega na definição que todo mundo conhece.

Mas eu acho que a segunda lei, enquanto axioma, não é das melhores. Imagine que eu quero medir o quão cansado eu fico ao empurrar um objeto. Eu fico mais cansado quando o objeto é pesado, e mesmo se o objeto for leve, eu fico cansado se o empurrar por muito tempo. Por isso, eu decido definir Cansaço como sendo o produto da massa pelo tempo: C=mΔt.

Mas tipo, e daí? Só porque eu dei um nome para um produto qualquer não quer dizer que isso seja uma lei da natureza. Da mesma maneira, só porque o Newton multiplicou aceleração por massa e chamou isso de força, não quer dizer que esse produto seja importante. A segunda lei não é uma lei propriamente dita: é só uma definição, e uma definição só faz sentido se tiver alguma propriedade que a acompanhe.

Terceira lei: Para cada ação, existe reação igual e oposta; em outras palavras, as ações de dois corpos um sobre o outro são sempre iguais e sempre opostas em direção.

Qualquer coisa que puxe ou empurre outra coisa, é igualmente puxado ou empurrado por ela. Se alguém aperta uma pedra com o dedo, o dedo também é apertado pela pedra. Se um cavalo puxa uma pedra amarrada por uma corda, o cavalo também será (por assim dizer) puxado igualmente na direção da pedra, já que a corda, esticada nas duas pontas, puxará o cavalo para a pedra e a pedra para o cavalo, no mesmo processo promovendo o movimento de um e impedindo o movimento do outro. Se um corpo, colidindo com outro corpo, muda a direção do segundo em razão de sua própria força, então, pela força do outro corpo (já que a pressão mútua é igual), também mudará seu movimento na direção oposta. Por meio dessas ações, mudanças iguais ocorrem no movimento, mas não na velocidade, exceto, é claro, se os corpos não forem impedidos por nada além. Por isso, mudanças na velocidade que ocorrem em direções opostas são inversamente proporcionais aos corpos, já que os movimentos são trocados igualmente. Essa lei também é válida para atrações.

Ahá, agora sim! A definição de força faz sentido porque força tem uma propriedade importante: ela é uma interação entre dois corpos que aparece em pares de igual intensidade e direções opostas. Nesse trecho eu também acho bacana os exemplos com pedras, cordas e cavalos: esse é um livro escrito em 1687, então cavalos eram muito mais comuns no dia-a-dia das pessoas que hoje em dia.

Mas eu acho que a terceira lei, enquanto axioma, não é das melhores. Imagine que você está em um carro que dá uma freada brusca. Você é acelerado para a frente, e se você multiplicar essa aceleração pela sua massa, descobre a força que te empurrou. Agora, onde está a reação dessa força, hein? hein?

A resposta é que essa força não tem reação. A terceira lei de Newton não funciona sempre, ela só funciona no caso especial em que o referencial não é acelerado. Essa força que você mediu é chamada de força inercial, ela não tem reação e só aparece em referenciais acelerados. Tem quem chame a força inercial de força fictícia, mas eu não gosto dessa nomenclatura: tente andar num carro sem cinto de segurança para ver o estrago fictício que o vidro vai fazer na sua cabeça!

"Mas peraí Ricbit, você não gosta de nenhuma das três leis!". Verdade, mas é porque dá pra fazer coisa melhor. Apresento a vocês a Lei Unificada de Newton:

Em um sistema isolado, o momento se conserva.

Pronto, uma única lei, bem mais simples que as anteriores, e você pode deduzir as três leis a partir dela. Na verdade, com esse enunciado você ainda ganha um monte de vantagens: ao contrário das três leis originais, essa versão funciona em referenciais não-inerciais, na relatividade especial e geral, e até na mecânica quântica!

Ler textos originais é um exercício muito divertido, mas para apreciar completamente você precisa se despojar do que aprendeu na escola. E a minha conclusão após ler bastante é que aprender é fácil, mas desaprender é muito mais difícil!