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!

Sorteios e Aniversários

No fim de semana passado eu participei do IV ENSL: o Encontro Nordestino de Software Livre, que nessa última edição foi realizado em Natal, no Rio Grande do Norte. Eu achei bastante apropriado, dado que Natal sempre foi um pólo de inovação tecnólogica! Citando alguns exemplos:

• Natal é a sede do Instituto de Neurociências do Miguel Nicolelis, pioneiro das interfaces cérebro-máquina e um dos maiores cientistas vivos.
• Alex Kipman, diretor de incubação da Microsoft, se inspirou na cidade para criar o Xbox Kinect (que, vocês devem lembrar, tinha o nome-código de Project Natal).
• Na década de 80, era onde estudava Tadeu Curinga da Silva, pioneiro dos videogames nacionais e autor do jogo Em Busca dos Tesouros, o melhor jogo da plataforma ZX-81.
• Além disso, Natal também é onde mora o Karlisson Bezerra, autor das tirinhas do Nerdson!

Embora o turismo nerd tenha sido excelente, eu fiquei um pouco frustrado de não ter aproveitado tanto o turismo tradicional. Eu consegui ir à praia e provar a deliciosa manteiga de garrafa, mas não deu tempo de fazer o típico passeio de buggy pelas dunas. Por outro lado, posso não ter feito o passeio de buggy, mas presenciei um outro tipo de bug!

Tradicionalmente, todo encontro de software livre termina com um sorteio de brindes. Como é um evento de programadores, naturalmente não se usou um saco plástico cheio de papéis com os nomes de cada um. Ao invés disso, o povo programou uma solução na hora mesmo, usando só três linhas de python:

import random
x = open("participantes.txt").readlines()
print x[random.randint(0, len(x))]

Como isso foi feito no modo interativo, então basta repetir a última linha para realizar um novo sorteio. Mas enquanto eles digitavam isso no telão, eu pensei: putz, esse código tem dois bugs. Você consegue achá-los?

Primeiro Bug

O primeiro bug, na real, não é do programa e sim da API do Python. Quando você projeta uma API, precisa ter muito cuidado para que as interfaces sejam consistentes, e o randint() é um dos raros casos onde o Python pisa no tomate.

O problema é que quase todas as definições de intervalo do Python usam intervalos semi-abertos, mas o randint usa um intervalo fechado (se você não lembra o que é um intervalo aberto ou fechado, pergunte ao professor Joe). Na notação de list slices, ou mesmo no range(), o último número não faz parte da lista; mas no randint() o último número conta.

Por isso, se você tentar sortear com randint(0, len(x)), ele pode sortear um número além do fim da lista, e você vai ter um erro do tipo list index out of range. Murphy sendo Murphy, é claro que isso aconteceu ao vivo lá no telão :)

Depois que uma API já está publicada e sendo usada, é muito difícil consertar um bug assim. Num mundo ideal você consertaria a interface, mas isso poderia quebrar todos os códigos legacy que a usam. Por isso, o melhor que dá pra fazer é depreciar a função original e propor uma interface nova, que no nosso caso é a randrange():

print x[random.randrange(0, len(x))]

Nesse problema em específico tem uma outra solução até mais elegante, que é usar o choice():

print random.choice(x)

Mas nenhuma dessas soluções conserta o outro bug do programa...

Segundo bug

Esse é mais sutil. No começo o programa estava funcionando bem, mas em um certo ponto, eis que o programa sorteia uma pessoa que já tinha sido sorteada antes! Eram mais ou menos 700 participantes, então a chance disso acontecer deveria ser bem pequena né? Mas algum tempo depois dessa pessoa, teve mais uma que foi sorteada duas vezes, e depois dessa ainda mais uma. Como pode um evento tão raro acontecer tantas vezes?

Marmelada não era, porque todo mundo viu o código, então o erro deve estar no gerador de números aleatórios do Python, certo? Surpreendentemente, a resposta é não! Esse é um caso clássico onde a nossa intuição falha. A análise ingênua é que a chance de um número ser sorteado duas vezes é 1/700 vezes 1/700, ou seja, 0.0002%. Mas esta é a análise correta do problema errado.

Se eu tivesse feito só dois sorteios, entre 700 participantes, e perguntasse qual era a chance do Ricbit ganhar nos dois, aí sim a chance seria 0.0002%. Mas o nosso problema é diferente. Nós temos um número p=700 de participantes, e vamos fazer s=70 sorteios (de olho eu acho que foram uns 70 sorteios: tinha brinde pra caramba no evento, e alguns sorteados eram descartados porque não estavam na sala).

O enunciado correto então é: qual o valor esperado do número de pessoas que é sorteada duas vezes? Esse valor é bem mais alto do que a intuição imagina!

Vamos calcular passo a passo quanto dá isso. Primeiro, vamos calcular qual é a probabilidade do Ricbit ser sorteado no primeiro evento. Isso é fácil né, a chance é 1/p. A chance oposta, ou seja, do Ricbit não ser sorteado no primeiro evento, é de (1-1/p). Continuando, a chance do Ricbit não ser sorteado em nenhum evento, é a abaixo:



Qual é a chance do Ricbit ser sorteado no primeiro evento e não ser sorteado em nenhum outro? Pelo mesmo raciocínio acima, ela é:



E qual a chance do Ricbit ser sorteado uma única vez durante toda a série de sorteios? Ora, é a chance de ser sorteado só no primeiro sorteio, mais a chance de ser sorteado só no segundo, e assim por diante. Como essas chances são todas iguais, então a chance do Ricbit ser sorteado uma única vez é:



Uma outra maneira de pensar a fórmula acima é que ela representa a chance de ser sorteado uma única vez numa posição qualquer, vezes o número de posições possíveis.

Agora chegamos onde queríamos, qual a chance do Ricbit ser sorteado exatamente duas vezes na série toda? Analogamente, é a chance de ser sorteado duas vezes, multiplicada pelo número de combinações possíveis dos dois sorteios ao longo da série. Logo, ela vale:



Por fim, qual o valor esperado do número de repetições? Uma repetição é quando uma pessoa é sorteada duas vezes: então basta, para cada pessoa, somar a probabilidade de ser sorteado duas vezes. Daí, o valor esperado final é:



Substituindo os parâmetros, descobrimos que o valor esperado é 3.1, o que realmente é bem próximo das três repetições que eu presenciei (nos cálculos eu desprezei sorteios triplos, quádruplos, etc, porque esses são realmente raros). Eu fiz uma simulação de Monte Carlo para quem só acredita vendo:

Simulação de Monte Carlo em Python

Na verdade, esse resultado contra-intuitivo é bem conhecido na literatura, o povo chama de Paradoxo do Aniversário. Numa festa com 30 pessoas, qual a chance de haver duas pessoas com o mesmo aniversário? A intuição é que isso deve ser raro, mas a probabilidade, na verdade, é de 70%. A argumentação é exatamente a mesma anterior, e, de fato, se você usar a nossa fórmula com s=30 e p=365, vai descobrir que o valor esperado é 1.1. Ou seja, em média, festas com 30 pessoas costumam ter um par de pessoas com o mesmo aniversário.

Se você trabalha com computação, então isso tem uma aplicação prática. A chance de repetição em um sorteio é exatamente a mesma chance de uma colisão em uma hash table, assumindo que sua função de hash tem distribuição uniforme. Por isso, quando você for dimensionar o tamanho de uma hash, pode usar nossa fórmula para saber o número esperado de colisões (por exemplo, em 70 inserções numa hash de tamanho 700, você vai ter em média 3 colisões).

E como resolver o problema original do sorteio? Basta tirar da lista os números que já foram sorteados! Uma implementação curta em python é a abaixo:

import random
x = open("participantes.txt").readlines()
random.shuffle(x)

print x.pop()

Ao invés de sortear nomes de maneira independente, você simplesmente faz uma permutação aleatória da lista com o shuffle() e vai tirando os elementos de um em um usando o pop(). Esse programa é bem curto e não tem os problemas do original. Mas, no fim das contas, se o programa original não tivesse nenhum bug, o sorteio não teria sido engraçado como foi :)

Se você é veterano do FISL, considere no ano que vem ir também ao ENSL. O povo é divertido, a comida local é excelente, e nada supera ir a uma praia do Nordeste entre uma palestra e outra :)

O Altar de Apolo

Quando pensamos na influência da religião na ciência, usualmente o que vem à mente é o obscurantismo medieval, como o julgamento do Galileu, ou então as cruzadas modernas contra a evolução e a pesquisa com as células-tronco. Mas também houve episódios onde a religião ajudou a ciência, e deles o meu preferido é o do Altar de Apolo.

Consta que, por volta de 400 aC, uma grande peste se espalhou pela cidade grega de Delos. Ela foi tão intensa, que estima-se que um quarto da população local tenha morrido. Preocupados com a saúde da população, os governantes resolveram consultar a única pessoa que poderia saber a causa da peste: o Oráculo de Delphi.

No Oráculo, as perguntas eram feitas para a sacerdotisa de Apolo, que entrava em contato diretamente com o Deus para solicitar a resposta. Diz-se que o Oráculo era construído sobre uma fenda que exalava gases alucinógenos, e isso ajudava no contato com o divino.

A sacerdotisa disse que a peste era causada pelo descontentamento de Apolo, e o único jeito de deixá-lo alegrinho novamente era com uma oferenda: eles deveriam duplicar o altar de Apolo, que era um grande cubo de mármore.

Os governantes então mediram as dimensões do cubo, duplicaram as arestas e fizeram o novo altar. Porém, mesmo seguindo à risca a orientação, a oferenda não funcionou, e a peste continuou matando todo mundo! Revoltados, eles foram tirar satisfação com o Oráculo: como assim a gente segue sua recomendação e não dá certo?

O Oráculo, que de bobo não tinha nada, explicou o que aconteceu: quando eles duplicaram cada aresta, o volume do altar cresceu oito vezes! Na verdade, o que eles deveriam ter feito era duplicado o volume, e não a aresta.

Isso deixou os matemáticos da época perplexos. Para duplicar o volume do cubo mantendo as proporções, a aresta precisa ser maior que a original por um fator igual à raiz cúbica de dois. Agora, se o Apolo é tão chato a ponto de não aceitar duplicar a aresta, ele também não iria aceitar uma aproximação como 1,26. Teria que ser o valor exato.

Mas como construir o valor exato usando as ferramentas disponíveis na época, ou seja, régua e compasso? Esse foi um dos grandes problemas da antiguidade, e a resposta correta só foi aparecer quase 2200 anos depois: na verdade, é impossível construir esse número só com régua e compasso.

Mas não pense que os gregos desistiram! Quando viram que esse caminho da régua e compasso não estava dando frutos, resolveram apelar para outros métodos. O mais impressionante deles foi criado por Arquitas, um general grego que era amigo pessoal de Platão.

O raciocínio de Arquitas foi o seguinte: régua e compasso são ferramentas para resolver problemas no plano; como nosso problema é espacial, provavelmente precisamos de ferramentas que operem no espaço. Arquitas tinha uma visão espacial assombrosa, e conseguiu construir a raiz cúbica de dois usando a intersecção de três sólidos.

Para entender a solução de Arquitas, vamos precisar de computação gráfica e geometria analítica. O primeiro passo é construir o segmento com o lado do cubo (digamos que o comprimento seja a). Centrados em um dos vértices do segmento, nós desenhamos três círculos perpendiculares de raio a, paralelos a cada um dos eixos coordenados:

Nós vamos usar cada um dos círculos para construir um sólido diferente. Primeiro, pelo círculo perpendicular ao Ox nós traçamos um cone que parte da origem:

A equação desse cone é a abaixo (se não é óbvio pra você que isso é a equação do cone, eu coloquei no final do post um quadro azul com as demonstrações).

O próximo sólido é um cilindro, construído estendendo o círculo perpendicular ao eixo Oz:

A equação desse cilindro é a abaixo:

Por fim, nós pegamos o círculo restante e o rotacionamos em torno do eixo Oz, criando um toróide:

A equação desse toróide é a abaixo:

O insight do Arquitas é que essas três superfícies se encontram exatamente no ponto que queríamos!

Vamos verificar que isso é verdade, usando as equações das superfícies:

Tal como queríamos! A parte impressionante disso é pensar que o Arquitas concebeu essa solução sem usar computação gráfica e sem saber geometria analítica. Na verdade, ele nem sabia escrever equações: mesmo uma coisa simples como o sinal de igual só foi inventado 1900 anos depois!

Se você quiser renderizar em casa a solução do Arquitas, pode usar os meus scripts de povray abaixo (para renderizar equações implícitas é só usar o comando isosurface):

Scripts povray com a solução do Arquitas

Por fim, o quadro azul para quem não sabe como derivar as equações implícitas. Enjoy!

Cone

O cone é formado por um contínuo de círculos no plano yz. Portanto, a equação de cada círculo deve ser da forma:

Mas o círculo em x=0 tem raio 0, e em x=a tem raio a. Logo o raio é igual a x.

Cilindro

O cilindro é um contínuo de círculos no plano xy, todos com o mesmo raio a, e centrados no ponto (a,0,0). O z nesse caso é qualquer, então a equação é:

Toróide

Eu não achei jeito fácil de derivar o toróide, então vai o difícil mesmo. Vamos construir a superfície como a soma de dois vetores, um que gira no plano xy, e um que gira nos planos perpendiculares a esse. A equação paramétrica do círculo base do toróide (visto de cima), é:

Cada ponto desse círculo é o centro de um outro círculo, no plano perpendicular formado pelo y e por esse mesmo vetor que aponta pro centro. Daí, a equação desse segundo círculo é:

A equação paramétrica final é a soma dos dois vetores. Já separando em coordenadas:

Agora nós elevamos ao quadrado e somamos x e y:

Podemos elevar z ao quadrado e somar com o anterior:

Opa, agora é só elevar ao quadrado de novo e substituir:

QED. Se alguém souber de algum jeito mais fácil, me ensine :)