Mostrando postagens com marcador cálculo numérico. Mostrar todas as postagens
Mostrando postagens com marcador cálculo numérico. Mostrar todas as postagens

sexta-feira, 16 de abril de 2021

As Figurinhas dos Vingadores são Honestas? (parte 2)

No último post do blog, eu comprei vários pacotinhos de figurinhas dos Vingadores, para tentar responder à pergunta: as figurinhas são honestas? Após fazer dois testes de hipótese, eu concluí que as figurinhas possuem distribuição uniforme, então parece que sim. Mas tem uma pegadinha: para serem honestas, garantir que são uniformes não é o suficiente, ainda tem outro teste que precisamos fazer!


A Distribuição dos Pacotinhos


Para fazer a análise das figurinhas, eu comprei 106 pacotinhos na banca. Mas enquanto eu abria os pacotinhos, notei uma propriedade curiosa: todos os pacotinhos tinham sempre 4 figurinhas diferentes. Mas olha só, se os pacotinhos não possuem repetidas entre si, então as figurinhas podem até ser uniformes, mas não são independentes! Cada pacotinho tem quatro figurinhas, e se eu tirei a figurinha 42 em um pacotinho, então eu sei que as outras três não são a figurinha 42.

Mas pode ter sido só o acaso né? Vai ver eu dei sorte de comprar exatamente 106 pacotinhos que não tinham repetidas. Mas felizmente nós sabemos como testar se foi sorte ou não, basta fazer o teste do ##\chi^2##, que vimos no último post!

Acho que o jeito mais fácil de testar é contar quantas figurinhas diferentes tem dentro do pacotinho. Quantos pacotinhos possíveis existem com 4 figurinhas diferentes? Esse é fácil, o álbum tem 180 figurinhas, então são ##180\choose 4## pacotinhos. Com 3 diferentes precisa de mais cuidado, existem ##180\choose 3## combinações de figurinhas, mas você precisa contar qual dessas três diferentes está repetida (os casos são AABC, ABBC, ABCC). Então para 3 diferentes são ##3 {180\choose 4}## pacotinhos.

Com 2 diferentes também tem três casos: AAAB, AABB, ABBB, logo temos ##3{180\choose 2}## pacotinhos. Para uma única figurinha diferente, são apenas ##180\choose 1## pacotinhos. Agora podemos montar a tabela com quantos pacotinhos eu deveria ter visto, e quantos pacotinhos eu vi:

1 diferente2 diferentes3 diferentes4 diferentes
Total de pacotinhos##180\choose 1####3{180\choose 2}####3{180\choose 3}####180\choose 4##
Total (numérico)18048330286758042296805
Total (normalizado)3.98e-61.07e-30.6340.936
Esperado (total*106)422e-60.1136.7299.163
Observado000106

Opa, nas duas primeiras categorias ##np\lt 5##, então eu não tenho amostras suficientes para usar o ##\chi^2##! Um jeito de consertar isso sem comprar mais pacotinhos é juntando as três primeiras categorias em uma só:

Menos de 4 diferentes4 diferentes
Esperado6.83799.163
Observado0106

Com isso conseguimos calcular o ##\chi^2## pela definição, chegando no valor de ##\chi^2=7.308##. Qual é o limite para rejeitar a hipótese? Se estivermos trabalhando com ##\alpha=5\%##, e um grau de liberdade, então o limite é:
InverseCDF[ChiSquareDistribution[1], 1 - 0.05]
>> 3.84146
Como ##7.308\gt 3.84146##, então temos que rejeitar a hipótese, e concluir que as figurinhas dos Vingadores não são independentes, e, como consequência, não são honestas, apesar de serem uniformes!


Felizmente para o colecionador, elas são desonestas mas estão jogando do seu lado. Se os pacotinhos diminuem a chance de ter repetidas, então o total esperado de pacotinhos que você precisa comprar é menor que no caso puramente honesto.

Dá para calcular o quanto você economiza com as figurinhas desonestas? Sure, vou mostrar dois métodos para fazer isso:

Cadeias de Markov


Esse primeiro método é mais complicado, mas vou começar por ele porque é mais genérico, e se aplica a mais casos. Imagine que eu já tenho 50 figurinhas no meu álbum, e eu abro mais um pacotinho. Se eu der azar, vão sair quatro figurinhas que eu já tenho. Se eu der sorte, pode sair até quatro que eu ainda não tinha. 

O método mais natural para modelar esse problema é com uma Cadeia de Markov: você considera o estado como sendo quantas figurinhas você tem no álbum, e as transições possíveis são quantas figurinhas você ainda não tinha no pacotinho.

Para montar a matriz de transição, nós precisamos da probabilidade de ir de um estado ##n## até um estado ##m##. Como cada pacotinho tem 4 figurinhas diferentes, então o ##m## precisa ser necessariamente ##n##, ##n+1##, ..., ##n+4##, e outras transições não são possíveis. 

Qual a chance de ir do estado ##n## para ##n+k##? Nós sabemos que as quatro figurinhas do pacotinho são diferentes entre si, então temos ##k## figurinhas novas e ##4-k## figurinhas que já tínhamos. A quantidade de pacotinhos que fazem essa transição é igual a ##{180-n\choose k}{n\choose 4-k}##, já que precisamos escolher ##k## figurinhas dentre que as não temos, e ##4-k## dentre as que temos. O total de pacotinhos é ##180\choose 4##, então a probabilidade total do estado ##n## para ##n+k## é ##{180-n\choose k}{n\choose 4-k}/{180\choose 4}## . Com isso podemos montar a matriz de transição ##M##:


$$ M= \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & 0.02 & 0.97 & 0 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & \approx0 & 0.04 & 0.95 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & \approx0 & \approx0 & 0.06 & 0.93 & \dots & 0 &  0 & 0 \\ \vdots &\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0.95 &  0.04 & \approx0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 &  0.97 & 0.02 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 &  0 & 1 \\ \end{bmatrix}$$


Pela matriz podemos conferir a intuição: se você não tem nenhuma figurinha, o primeiro pacotinho vai garantidamente te dar 4 figurinhas novas, então o estado 0->4 na matriz vale 1, e todos os outros nessa linha são zero. Por outro lado, se está faltando só uma figurinha para completar, comprar um pacotinho novo dá um chance de 2% de conseguir a faltante. E, claro, se você tem todas as figurinhas, então a chance de você não conseguir mais nenhuma é 100%.

O que queremos calcular é quanto economizamos devido às figurinhas desonestas. Para isso, vamos calcular qual o valor esperado do número de pacotinhos que precisamos comprar para completar o álbum. Tendo a matriz de transição em mãos, isso é equivalente ao número médio de transições para ir do estado 0 (álbum vazio) ao estado 180 (álbum completo). E como calcular isso? Usando um algoritmo secreto!

Really, um algoritmo secreto. Não sei o motivo, mas você não acha esse algoritmo em lugar nenhum; não tem na wikipedia e nem nos livros de processos estocásticos que eu tenho aqui em casa. Então aproveitem para anotar porque é bem útil:

  1. Dada a matriz de transição ##M##, calcule a matriz reduzida ##R##, retirando da matriz original todas as linhas e colunas associadas aos estados finais. (No nosso caso só a última linha e coluna são finais).
  2. Calcule a matriz fundamental ##N## do processo, através da fórmula ##N=(I-R)^{-1}##.
  3. O valor esperado do número de transições é a soma dos elementos da primeira linha da matriz ##N##.
Implementando o algoritmo acima com um scriptinho rápido no Mathematica, concluímos que o valor esperado do número de pacotinhos que você precisa comprar é ##257.98##. Como cada pacotinho custa R$ 1.20, então você precisa de R$309.58 em média para completar o álbum.

E quanto você gastaria com figurinhas 100% honestas? Esse cálculo eu já fiz aqui no blog anteriormente, se as figurinhas são uniformes e independentes então esse é o Coupon Collector Problem: para ##n## figurinhas a solução é ##n H(n)##, onde ##H(n)## é o número harmônico de ##n##. Essa fórmula dá o total de figurinhas, como elas são independentes, para o total de pacotinhos é só dividir por 4, chegando no total de ##259.78## pacotinhos, ou R$ 311.74, uma economia de 2 reais! (Eu falei que economizava, só não falei que era muito).

Embora o algoritmo secreto seja bacana e funcione em vários casos, ele é meio lento. Você precisa inverter uma matriz 181x181, o que não é muito fácil. O fato dela ser triangular e esparsa ajuda, mas tem um método melhor, que calcula o valor em O(n)!

Linearidade da Expectativa


Qual é o teorema mais surpreendente que você conhece? Só em probabilidade tem vários que são difíceis de acreditar à primeira vista, como o Monty Hall Problem e o Paradoxo do Aniversário. Mas um que sempre me espanta é a linearidade da expectativa (o valor esperado da soma é a soma dos valores esperados):

$$E[a+b]=E[a]+E[b]$$

Para mim a parte curiosa é que ele funciona mesmo se as variáveis ##a## e ##b## não forem independentes! Para o nosso problema isso ajuda bastante: o valor esperado para conseguir 180 figurinhas claramente é correlacionado com o valor esperado para conseguir 179 figurinhas, mas eu posso usar relações lineares entre eles sem medo.

Para essa análise eu acho mais fácil pensar invertido. Se estiver faltando 0 figurinhas para completar o album, quantos pacotinhos eu preciso? Nenhum né, ele já está completo. E se estiver faltando só uma figurinha? Aí depende do que sair no pacotinho, se sair a figurinha faltante, um pacotinho basta, se não sair, então você adiciona um pacotinho no seu valor esperado e repete o procedimento:

$$\begin{align*}\\E[1]&=p_0(1+E[0]) + p_1(1+E[1])\\ (1-p_1)E[1] &= p_1 + p_0(1+E[0]) \\ E[1] &=\frac{1}{1-p_1}\left(p_1+p_0(1+E[0])\right)\end{align*}$$

Olha só, eu sei que ##E[0]=0##; para calcular ##E[1]## eu preciso só do ##E[0]##; para o ##E[2]## eu preciso do ##E[1]## e ##E[0]##; e assim por diante. Como cada pacotinho tem no máximo 4 figurinhas que eu não tenho, então eu só preciso olhar quatro valores para trás na recorrência. Portanto, esse algoritmo é linear!

Preenchendo as probabilidades corretas, o algoritmo completo fica assim:

$$\begin{align*} \\ E[0] &= 0 \\ E[n] &= \frac{1}{1-{{180-n}\choose 4}/{180\choose 4}}\left(\frac{{180-n}\choose 4}{180\choose 4}+\sum_{1\le k \le 4}\frac{{n\choose k}{{180-n}\choose{4-k}}}{180\choose 4}E[n-k]\right) \\ \end{align*} $$

Essa implementação realmente é bem mais rápida que o anterior. A versão com matrizes no Mathematica leva 1s para rodar, enquanto que um scriptinho python com o algoritmo linear leva só 55ms!

Agora que temos os pacotinhos corretamente modelados, finalmente vamos conseguir bolar uma estratégia para minimizar o custo de completar o álbum. Mas vai ficar para o próximo post porque esse já ficou comprido demais :)


segunda-feira, 2 de junho de 2014

A Equação da Torre Eiffel

No mês de maio eu completei cinco anos de casado! Quando eu casei, aproveitei a lua-de-mel em Paris para falar sobre a ciência da Torre Eiffel, então nada mais apropriado que aproveitar o aniversário para falar sobre outro mistério da torre: a equação que descreve o seu formato!
  mini_ilafox_ricbit_paris  

A torre mais alta do mundo


O desafio proposto pelo Gustave Eiffel era grandioso: construir a torre mais alta do mundo. Na época, o recorde era o Monumento de Washington, com 169 metros. O Eiffel pretendia erguer sua torre com 324 metros, quase o dobro do recorde anterior.

Para entender o desafio, vale a pena comparar a Torre Eiffel com outras torres da antiguidade. O desenho abaixo está em escala, e mostra a Torre Eiffel comparada com a Grande Pirâmide de Giza no Egito, e com o Campanário de São Marcos em Veneza:

compare
Imagem: Wikipedia (CC-SA)

O primeiro desafio ao fazer uma torre é garantir que ela consiga suportar seu próprio peso. Os egípcios resolveram isso com força bruta: usaram blocos de pedra maciços e fizeram a base bem maior que o topo. O resultado final certamente funcionou, a pirâmide está aí de pé 4500 anos depois. Mas ela é um desperdício de material, com tecnologia moderna você não precisa de tanta matéria-prima para chegar nessa altura.

O Campanário de São Marcos é o extremo oposto. Ele foi construído em 1514, e nessa época a tecnologia já permitia criar uma torre com a largura da base igual à do topo. O problema é que a torre caiu! Em 1902 ela não resistiu ao próprio peso e colapsou (nenhuma vida humana foi perdida, mas o gato do zelador morreu). A torre que você visita hoje em Veneza é uma reconstrução.

Como fazer, então, uma torre alta que suporte o próprio peso? O Eiffel sabia a resposta porque tinha experiência no assunto. Antes de criar os ícones pelo qual é mais lembrado, ele construía pontes. Existem pontes do Eiffel por todo o mundo, desde a França e a Romênia, até o Vietnã e o Chile.

Para entender a tecnologia do Eiffel, vale a pena comparar com uma ponte antiga: a Pont du Gard, um aqueduto romano construído dois mil anos atrás. Ele cobre uma distância de 275 metros a uma altura de 48 metros. Para isso, ela precisa de cinco apoios no chão e três níveis de arcos para aguentar seu próprio peso:

Pont_du_gard_v1_082005
Imagem: Wikipedia (CC-SA)

Agora, compare com uma das mais belas pontes construídas pelo Eiffel: a Ponte Dona Maria, na Cidade do Porto, em Portugal. Ela cobre 353m a 60m de altura, em ambos os casos maior que a ponte romana. E faz isso usando só dois apoios no chão! O segredo do Eiffel para conseguir esse feito é a magia da treliça. Ao usar uma treliça, você consegue suportar mais peso com menos material.

Ponte_Maria_Pia_-_Porto Imagem: Wikipedia (CC-SA)

Todo mundo sabia que o Eiffel era expert no assunto, e que ele saberia construir uma torre que suportasse o próprio peso. Mas, ainda assim, na época ele ouviu muitas críticas. Uma coisa é você construir uma ponte na horizontal, outra é construir uma torre na vertical. Estruturas na vertical tem um problema adicional: será que o Eiffel sabia como fazer a torre resistir ao vento?

O método geométrico do Eiffel


Quando a torre ficou pronta, o Eiffel calou os críticos. Quem já visitou a torre sabe que no terceiro piso venta. Venta muito mesmo! Porém, a torre não se abala. Fica lá, imóvel, mesmo com os ventos mais fortes.

E qual foi o truque do Eiffel para deixar a torre tão estável com o vento? Até pouco tempo atrás, ninguém sabia. As plantas da torre Eiffel são conhecidas, você pode ver todas as especificações em um livrão publicado pela Taschen (todas mesmo, tem até o diâmetro e comprimento de cada parafuso). Mas as plantas não falam o truque da estabilidade, tudo que dizem é que o Eiffel determinou o formato da torre de forma a minimizar a ação do vento, e fez isso usando "métodos geométricos baseados em observações empíricas".

 DSC09174 

O formato da torre foi um mistério por mais de um século, até que Patrick Weidman, professor da Universidade de Colorado, matou o problema, e da maneira mais simples possível. Ele foi atrás dos diários do Eiffel e achou lá a dica que resolve o problema! O trecho do diário dizia:

"Suponha que as faces de uma treliça simples formam uma parede resistente às forças cisalhantes do vento, cujas componentes horizontais são P', P'', P''', P''''. Sabemos que, para calcular as forças atuantes no corte pelo plano MN, precisamos determinar a resultante P de todas as forças exteriores atuando acima da secção, e decompor essa resultante nas três forças passando pelas peças cortadas. Se a forma do sistema for tal que, para cada corte horizontal MN, a extensão das peças externas da treliça se intersectam na força exterior P, então as forças na barra interna serão zero e podemos excluir esse membro."

O diário tinha esse desenho acompanhando:

trelica_original

Lendo com cuidado dá para entender direitinho a motivação do formato da torre. Vamos recriar o raciocínio em câmera lenta. Nós começamos supondo que o vento é uma carga horizontal atuando sobre toda a torre (que ele representou como as forças P', P'', etc). O truque agora é pegar um corte horizontal arbitrário (que ele chamou de MN), e tirar toda a parte de cima da torre. As forças na estrutura não mudam se substituirmos a parte de cima por um ponto material de mesma massa, localizado no centróide do que foi tirado.

   diagrama 

Vamos então ver qual é o efeito do vento. Eu tenho uma força P horizontal, que é a resultante do vento. Tenho ainda três forças F1, F2, F3, que são as forças atuando sobre as vigas da treliça. Quando uma treliça está em equilíbrio, as forças sobre ela são sempre de tração ou compressão, ou seja, estão na mesma direção de suas vigas respectivas.

A torre está em equilíbrio estático, logo a soma das forças é zero. Além disso, ela também não está rotacionando; logo a soma dos torques também é zero. O torque é o produto da força pela distância; mas não é qualquer distância, é a distância na direção perpendicular à força.

Agora o truque do Eiffel é claro. Ele construiu as laterais da torre de modo que as tangentes às vigas mais externas sempre passam pelo centróide. Quando você faz isso, a direção das forças P, F1 e F3 apontam para o centróide, então o torque delas é zero. Só a força F2 está aplicando torque no centróide!

diagrama2

Podemos calcular a soma dos torques agora. A distância associada à força F2 é a distância d, em verde na figura. Note que ela é estritamente positiva.
$$\begin{align*} P\times 0 +F_1\times 0+F_2\times d+F_3\times 0&=0\\ F_2\times d&=0\\ F_2&=0 \end{align*}$$
Ahá! Se as tangentes externas apontam para o centróide, então a força na viga interna é zero! Se o vento horizontal fosse a única força atuando sobre a torre, então nós poderíamos até tirar a viga interna que não faria diferença!

É claro que na prática não podemos tirar a viga. Ela não ajuda a resistir ao vento, mas ajuda a sustentar o peso. Por outro lado, essa viga pode ser mais fraca que as outras. Com a idéia do Eiffel, só as vigas externas precisam ser fortes. Nessa foto que eu tirei do pilar oeste isso é bem visível: as vigas externas são largas e feitas de material rígido, enquanto que as internas são finas e feitas de mini-treliças, o que deixa a estrutura como um todo mais leve. E quanto mais leve, mais alto você pode subir.

ouest

A equação da Torre Eiffel


Agora que nós sabemos o truque do Eiffel, podemos calcular qual é o lugar geométrico dos pontos cujas tangentes passam pelo centróide, isso vai nos dar a equação que descreve o fomato da torre. Aparentemente o Eiffel nunca precisou fazer isso, já que ele usou "um método geométrico baseado em observações empíricas" (para mim isso é um jeito de dizer que ele calculou uma solução numérica ao invés de resolver a equação de fato).

Mas ele poderia ter resolvido a equação se quisesse. Vamos fazer as contas na caixa azul para conferir, pule se você tem medo de equações diferenciais:

Começaremos supondo que o formato da torre pode ser descrito por uma função f(x). Para simplificar as contas (e, acredite, simplifica mesmo), nós vamos rotacionar o gráfico de modo a colocar o eixo x na vertical, orientado para baixo, e o f(x) na horizontal. Vamos ainda supor que f(x) é real, contínua e pode ser diferenciada muitas vezes. grapheiffel
No gráfico acima, f(x) está desenhada em azul. Por um ponto qualquer x eu traço as retas tangentes (em vermelho), elas vão se encontrar em um ponto xC, que nós estamos supondo que sempre será o centróide da parte superior da torre (ou seja, do trecho que vai do ponto escolhido x, até o topo da torre xT). A posição do centróide é dado diretamente pela definição: $$x_C=\frac{\int_{x_T}^{x}t f(t)\,dt}{\int_{x_T}^{x} f(t)\,dt}$$ A tangente podemos calcular de duas maneiras. Nós sabemos que a tangente que passa por x é dada diretamente por f'(x). Por outro lado, ela também é a tangente do ângulo formado pelo eixo x e pela reta vermelha, que você calcula direto como cateto oposto pelo adjacente: $$f'(x)=\frac{\Delta y}{\Delta x}=\frac{f(x) -0}{x-x_C}=\frac{f(x)}{x-x_C}$$ Podemos isolar xC nessa última equação: $$\begin{align*} x-x_C&=\frac{f(x)}{f'(x)}\\ x_C&=\frac{x f'(x)-f(x)}{f'(x)} \end{align*}$$ E depois igualar as duas equações: $$\begin{align*} \frac{x f'(x)-f(x)}{f'(x)} &= \frac{\int_{x_T}^{x}t f(t)\, dt}{\int_{x_T}^{x} f(t)\, dt}\\ \left(x f'(x)-f(x)\right)\int_{x_T}^{x} f(t)\, dt &= f'(x)\int_{x_T}^{x}t f(t)\, dt\\ \left(f'(x)\int_{x_T}^{x}x f(t)\, dt \right)-\left(f(x)\int_{x_T}^{x} f(t)\, dt\right) &= f'(x)\int_{x_T}^{x}t f(t)\, dt\\ f(x)\int_{x_T}^{x} f(t)\, dt &= f'(x)\int_{x_T}^{x}(x-t) f(t)\, dt \end{align*}$$ Na última linha, veja que o x passou para dentro da integral, você pode fazer isso porque ele é constante em relação à variável de integração, que é t. Antes de prosseguir, vamos dar nomes às integrais, chamando-as de g(x) e h(x). Note que essas integrais são funções de x. $$\begin{align*} g(x)&=\int_{x_T}^{x} f(t)\, dt \\ h(x)&= \int_{x_T}^{x}(x-t) f(t)\, dt \end{align*}$$ Vamos substituir as integrais pelas funções que nomeamos, e depois derivar dos dois lados. Quando você deriva dos dois lados, você está comendo uma constante, então no final podem aparecer soluções espúrias. Depois de acabar tudo, temos que pegar as soluções e voltar para o começo, conferindo se elas são realmente válidas. $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ f'(x)g(x)+f(x)g'(x) &= f''(x)h(x)+f'(x)h'(x) \end{align*}$$ Como calcular g'(x) e h'(x)? Vamos começar por g'(x). Note, primeiro, que você pode usar o teorema fundamental do cálculo para escrever g(x) em função da primitiva de f(x), que iremos chamar de F(x): $$g(x)=\int_{x_T}^x f(t)\,dt=F(x) -F(x_T) $$ Já que F(xT) é uma constante, sua derivada é zero. F(x) é uma primitiva, e a derivada da primitiva de uma função é a própria função (desde que ela seja contínua): $$\begin{align*} g(x)&=F(x) -F(x_T)\\ g'(x)&=F'(x)-0=f(x) \end{align*} $$ Então g'(x)=f(x), essa foi fácil. Para calcular h'(x) é mais complicado. Lembremos a regra de integração por partes para integrais definidas: $$\begin{align*} \int_c^d a(t)b'(t)\,dt&=[a(t)b(t)]_c^d-\int_c^d a'(t)b(t)\,dt \end{align*} $$ Vamos escolher quem são a(t) e b(t): $$\begin{align*} a(t)=x-t&\implies a'(t)=-1\\ b(t)=F(t)&\implies b'(t)=f(t)\\ \end{align*} $$ Agora podemos quebrar h(x) em duas partes: $$\begin{align*} h(x)&=\int_{x_T}^x(x-t)f(t)\,dt\\ &=[(x-t)F(t)]_{x_T}^x-\int_{x_T}^x-F(t)\,dt\\ &=(x-x)F(x)-(x-x_T)F(x_T)+\int_{x_T}^xF(t)\,dt\\ &=-(x-x_T)F(x_T)+\int_{x_T}^xF(t)\,dt\\ &=-(x-x_T)F(x_T)+P(x)-P(x_T)\\ \end{align*} $$ No último passo, por falta de notação melhor, eu chamei de P(x) a primitiva de F(x). Agora podemos derivar h(x). Note que, se P(x) é a primitiva de F(x), então F(x) é a derivada de P(x). Além disso, como P(xT) é constante, sua derivada é zero: $$\begin{align*} h(x)&=-(x-x_T)F(x_T)+P(x)-P(x_T)\\ h'(x)&=-F(x_T)+F(x)-0\\ h'(x)&=\int_{x_T}^xf(t)\,dt\\ h'(x)&=g(x)\\ \end{align*} $$ Concluímos então que h'(x) é o mesmo g(x) que tínhamos visto antes! Antes de continuar, vale a pena resumir o que nós temos até agora: $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ f'(x)g(x)+f(x)g'(x) &= f''(x)h(x)+f'(x)h'(x)\\ f(x)&=g'(x)\\ f'(x)&=g''(x)\\ f''(x)&=g'''(x)\\ h'(x)&=g(x) \end{align*}$$ Note que dá para reescrever quase tudo usando só g(x) e suas derivadas. Fica faltando h(x), mas isso não é problema porque temos duas equações usando h(x), então é só isolar. Da primeira equação temos: $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ h(x)&=\frac{f(x)g(x)}{f'(x)}\\ h(x)&=\frac{g'(x)g(x)}{g''(x)}\\ \end{align*}$$ Da segunda equação, temos: $$\begin{align*} f'(x)g(x)+f(x)g'(x)&=f''(x)h(x)+f'(x)h'(x)\\ f''(x)h(x)&=f'(x)g(x)+f(x)g'(x)-f'(x)h'(x)\\ h(x)&=\frac{f'(x)g(x)+f(x)g'(x)-f'(x)h'(x)}{f''(x)}\\ h(x)&=\frac{g''(x)g(x)+g'(x)g'(x)-g''(x)g(x)}{g'''(x)}\\ h(x)&=\frac{g'(x)^2}{g'''(x)}\\ \end{align*}$$ Igualando as duas expressões para h(x): $$\begin{align*} \frac{g'(x)g(x)}{g''(x)} &= \frac{g'(x)^2}{g'''(x)} \\ \frac{g'(x)}{g(x)} &= \frac{g'''(x)}{g''(x)} \end{align*}$$ Tudo que falta agora é resolver a equação para g(x) e finalmente temos o formato da torre. Mas essa equação é uma equação diferencial não-linear de terceira ordem. E isso é os que os matemáticos usam para assustar os filhos: "se você não comer, vai aparecer o bicho-papão, o homem do saco e a equação diferencial não-linear de terceira ordem para te pegar!" Felizmente, essa equação pode ser resolvida se você tiver a visão além do alcance. O que acontece quando você tenta diferenciar log g(x)? $$\begin{align*} \frac{d}{dx}\log(g(x))=\frac{1}{g(x)}\times g'(x)=\frac{g'(x)}{g(x)} \end{align*}$$ Ahá! Essa é exatamente a forma da equação. Fazendo o análogo no outro lado, podemos reescrever a equação original: $$\begin{align*} \frac{g'(x)}{g(x)} &= \frac{g'''(x)}{g''(x)}\\ \frac{d}{dx}\log g(x) &= \frac{d}{dx}\log g''(x) \end{align*}$$ Agora nós podemos integrar dos dois lados (lembrando que vai aparecer uma constante), e depois tirar a exponencial dos dois lados: $$\begin{align*} \frac{d}{dx}\log g(x) &= \frac{d}{dx}\log g''(x)\\ \log g(x) &= \log g''(x) + C\\ g(x) &= k g''(x) \\ \end{align*}$$ Agora ficou fácil né? A equação medonha virou uma equação diferencial linear de segunda ordem, essa é bem-comportada, dá a patinha e finge de morta. A solução geral é da forma abaixo: $$g(x)=A e^{B x}+C e^{-B x}$$ Como f(x)=g'(x), então f(x) também tem a forma acima, só os coeficientes que serão diferentes. A natureza da função depende dos coeficientes; por exemplo, se A=0.5, C=0.5 e B=i, então f(x)=cos(x). Mas nem todas as combinações de coeficientes são soluções válidas! Para determinar os coeficientes corretos, precisamos voltar o f(x) lá na equação original: $$f(x)\int_{x_T}^{x} f(t)\, dt = f'(x)\int_{x_T}^{x}(x-t) f(t)\, dt $$ Daqui em diante é só fazer as contas. Se você substituir tudo, vai notar que o único jeito de conseguir uma solução real válida para qualquer x é fazendo C=0 e xT=-∞; nesse caso A pode ser um real qualquer, e B pode ser um real positivo qualquer. Estritamente falando, a torre teria que ser infinitamente alta para satisfazer a essa equação, mas na prática você pode aproximar "infinitamente alta" por "muito, muito alta". Concluindo, a solução final para o formato da torre Eiffel é: $$f(x)=A e^{B x} $$ Ou seja, a torre tem um perfil exponencial!

Matemática e Engenharia


Agora que sabemos qual o formato teórico da torre, podemos usar os dados reais da torre para estimar os coeficientes da exponencial. Eu rodei um best fit e cheguei em A=3.34323 e B=0.0102592:

bestfit

A aproximação é boa, mas... parece que não encaixa direitinho? É meio frustrante fazer esse monte de contas, e no fim a curva real não ser tão parecida assim com a equação deduzida.

Mas tem um motivo para isso, e o motivo é que o Eiffel era engenheiro, não matemático. Quando você plota uma exponencial em um gráfico semi-log, o resultado tem que ser uma linha. Mas olha só o que acontece com os dados reais da torre: eles formam duas linhas!

semilog 

Coisas que você aprende na sua primeira aula de engenharia: se você faz o projeto de um elevador, e determina que o cabo do elevador precisa aguentar no mínimo 10 toneladas, qual o cabo que você coloca no elevador? Um matemático diria que é um cabo de 10 toneladas, mas o engenheiro vai falar que o mínimo são 14 toneladas.

Isso é o resultado do fator de segurança. Por mais que o seu modelo matemático seja correto, sempre tem alguma coisa que você não levou em conta (pode ser que o aço que você usou tenha impurezas, ou que o solo não é tão firme quanto você achava, e assim por diante). Por isso, é costume sempre multiplicar o valor final por um fator de segurança para levar em contas esses imprevistos.

O Eiffel sabia disso, e no diário ele explica que imaginava que a base da torre estaria sujeita a torques maiores que o topo, e por isso usou um fator de segurança diferente na base e no topo. Quando você faz o best fit com duas exponenciais ao invés de uma só, o gráfico fica bem melhor!

twofit 

Agora finalmente podemos dar o mistério por resolvido: o formato da Torre Eiffel é descrito por duas exponenciais, escolhidas para minimizar o material necessário para resistir à carga do vento.

A Torre Eiffel é o literalmente o maior monumento à ciência construído, da próxima vez que passar por Paris aproveite para fazer uma reverência aos cientistas do passado :)

eiffeltower

terça-feira, 26 de junho de 2012

Estimule sua Criatividade #1

Um tempo atrás me perguntaram se existe algum método para ser mais criativo. Isso é bem difícil: para testar se um método funciona ou não, você precisaria medir a criatividade antes e depois de usar o método, e ninguém sabe como mensurar criatividade.

Dito isso, eu conheço dois truques que parecem funcionar bem. Neste post eu vou falar do primeiro método, e pra isso eu vou contar um história que aconteceu comigo lá nos antigos tempos de colégio técnico.


A Federal


Em 1991 eu comecei o curso técnico em eletrônica na ETFSP, popularmente conhecida como a Federal (e que hoje em dia é IFSP). A primeira coisa que nos passaram foi a lista de material para comprar. Além dos básicos livros e cadernos, a lista também incluía itens mais específicos, como caneta nanquim, normógrafo e calculadora científica.

A calculadora foi um problema. Naquela época, calculadoras científicas eram caras (ou melhor, até hoje em dia elas são caras). Como eu não tinha condições de comprar uma, resolvi encarar o desafio de fazer o curso usando só uma calculadora de quatro operações.

O primeiro ano foi fácil: no início do curso, o único componente usado era o resistor, e contas com resistores você consegue fazer só com as quatro operações. No segundo ano apareceram os capacitores, e aí a coisa complicou: agora as contas envolviam exponenciais e funções trigonométricas, que a minha calculadora simples não fazia.

O que fazer então? Dado que eu não conhecia ninguém para me emprestar uma calculadora científica, nem tinha dinheiro para comprar uma, o jeito foi apelar para a herança do meu avô.


A herança do meu avô não era dinheiro. Melhor que isso, era conhecimento! Mais especificamente, era uma estante enorme que ia de parede a parede, com um monte de livros sobre todo tipo de assunto. Os livros eram antigos, da década de 60, e iam desde o Tesouro da Juventude até a Gramática do Jânio Quadros.

O livro que me salvou foi uma curiosa coleção de Matemática para Seminaristas. Eu não sei exatamente por que um padre precisa saber matemática, mas o fato é que o livro era bem completo, começava com produtos notáveis e avançava até o Cálculo. E no meio do caminho, tinha o truque que eu precisava para usar minha calculadora!

A Exponencial


Os problemas que caíam na prova eram mais ou menos assim: dado o circuito RC-série abaixo, calcule a tensão no capacitor após 2 segundos.



No colégio técnico você não precisava resolver a equação diferencial. Ao invés disso, eles davam a fórmula final; você só precisava decorar e aplicar:


O meu problema era a exponencial, que só tem em calculadoras científicas. Mas, usando o livro dos seminaristas, eu descobri a fórmula que resolveu meu problema: a expansão em série de Taylor:


Parece complicado, mas é simples calcular essa fórmula numa calculadora de quatro operações, desde que ela tenha os botões de memória (M+, M-, MC e MR). No exemplo dado, em que eu preciso calcular exp(-0.2), a sequência de botões é curta:

C MC 1 M+ . 2 M- × . 2 ÷ 2 = M+ × . 2 ÷ 3 = M- × . 2 ÷ 4 = M+ MR

Com 29 toques eu tenho o valor de 0.818, correto com três casas decimais (e o professor só pedia duas). Se você tiver os dedos ágeis, praticamente não leva tempo!

Esse método também é bom para calcular senos e cossenos, que aparecem no cálculo de ângulos de fasores e no cálculo do fator de potência. Mas tem um caso onde o método não funciona...

O Logaritmo


Às vezes, o que caía na prova era o problema inverso: dado aquele mesmo circuito RC-série, quanto tempo leva pro capacitor ficar com metade da tensão da fonte? Nesse caso, a conta depende de log(0.5), e não tem logaritmo na calculadora de quatro operações.

Pior, nesse caso o livro também não ajudava. A única fórmula que tinha no livro era a seguinte:


Se você tentar usar essa fórmula na calculadora, logo vai perceber que ela não é prática: são necessários muitos termos para conseguir a precisão de duas casas decimais desejada. Olhando as fórmulas, a velocidade de convergência é clara: os coeficientes da exponencial caem com a velocidade do fatorial, enquanto que os coeficientes do logaritmo caem com a velocidade da série harmônica, que é muito mais lenta.

Hoje em dia, eu provavelmente usaria o método de Newton-Raphson para calcular o logaritmo. Mas eu não entendia de cálculo numérico, então nem sabia que esse método existia. Por isso, eu tive que apelar. Ao invés de decorar uma fórmula, eu decorei quatro números:


É fácil decorar quatro números de 7 dígitos né? Certamente todo mundo sabe pelo menos quatro números de telefones de cabeça, e a dificuldade é a mesma. Com esses quatro números, eu conseguia calcular tudo que queria! Precisa de log(0.5)?


Precisa de log(4.2)? Beleza:



E se for log(2.6)? Ops, aí não dá, 26 tem um fator primo 13 que não está na minha lista. Mas eu posso aproximar:


O resultado não é exato, mas tem precisão de duas casas. Quem precisa de calculadora científica, afinal?

Isso foi sorte, ou sempre é possível aproximar? Na época eu tinha uma intuição que sempre dava, mas hoje em dia eu consigo demonstrar que isso é verdade! A demonstração está na caixa azul, pule se você não sabe teoria dos números:

A demonstração usa diretamente o teorema da equidistribuição de Weyl: o conjunto das partes fracionárias dos múltiplos inteiros de um número irracional qualquer é denso em [0,1). Ou seja, para qualquer irracional α e real q, onde 0 ≤ q < 1, e para qualquer ε positivo dado, sempre existe um inteiro n tal que a fórmula abaixo é verdadeira:


Com um pouquinho de álgebra é fácil estender o domínio de [0,1) para todos os reais. Imagine que escolhemos um real x qualquer, tal que x=p+q, onde p é a parte inteira e q é a parte fracionária. Então podemos partir da equação anterior:


Note que o piso de  é um inteiro, e p é um inteiro. Vamos chamar de m a diferença dos dois:


Ou seja, eu posso escolher qualquer irracional α e qualquer real x, que sempre vai existir um par de inteiros m e n que satisfazem a inequação. Já que eu posso escolher qualquer número, então escolho α e x como abaixo:


Note que α é irracional, portanto podemos usar nossa inequação:


Ou seja, sempre podemos aproximar o log de um real y qualquer como a soma de múltiplos inteiros de log(2) e log(3), nem precisava ter decorado o log(5) e log(7)!

Infelizmente, essa demonstração é um exemplo de prova não-construtiva: ela garante a existência dos inteiros m e n, mas não fala como achá-los! No fim, você precisa descobrí-los pela intuição, exatamente como eu fazia :)

A Criatividade


Eu poderia terminar o post com uma moral do tipo "se a vida te deu limões, então arranje uns bastões de cobre e faça uma bateria", mas tem uma lição mais importante! A Marissa Mayer cantou a bola em uma entrevista de 2006: a criatividade adora restrições, especialmente se acompanhada de um desprezo saudável pelo impossível.

Tente se lembrar de alguma coisa que você viu e achou muito criativo. Que tal a artista que fazia retratos usando fita cassete? O cara que faz arte usando frutas? O doido que reescreveu The Raven, do Poe, de modo que o número de letras do poema batesse com os dígitos de pi?

Todos eles são exemplos onde o autor se auto-impôs uma restrição (de forma, de conteúdo, de matéria-prima). Se você quer expandir sua criatividade, impor limites costuma ser mais efetivo que retirá-los. Eu mesmo sempre tive problemas quando me pediam uma redação "com tema livre", era muito mais fácil quando o tema tinha alguma restrição.

E isso não vale apenas para arte, vale para computação também. Tente fazer aquele seu programa usar menos memória, menos CPU, tente programar a mesma coisa em menos tempo: tudo isso vai te forçar a usar soluções mais criativas.

Eu ainda tenho outro truque para estimular a criatividade, mas esse fica para o post seguinte :)

(Agradecimentos ao povo do Math Stack Exchange pela ajuda na demonstração)

quarta-feira, 13 de junho de 2012

O Bug mais Difícil

Essa é uma pergunta clássica em entrevistas: qual foi o bug mais difícil que você já consertou? No meu caso, o bug mais difícil que resolvi deu um trabalhão para encontrar, e por um motivo bastante curioso: o bug não era no programa, mas sim na pecinha entre a cadeira e o teclado!


A história


Essa história se passa em agosto de 1998. Nessa época, eu estava escrevendo um emulador de MSX, e tomando dois cuidados na implementação: o emulador tinha que ser rápido o suficiente para rodar em full speed nas máquinas da época, e tinha que ser preciso o suficiente para ser indistinguível de um MSX real.

O primeiro objetivo eu atingi escrevendo o emulador inteiramente em Assembly, o segundo através de uma metodologia de testes sólida. Na época eu ainda não fazia unit testing, mas, para cada feature adicionada, eu sempre fazia um programa de teste e comparava a execução dele no emulador e na máquina real.

As primeiras versões do emulador eram mudas, eu emulava só a CPU e o chip de vídeo. Quando esses dois ficaram estáveis, eu comecei a escrever a emulação do chip de audio: o AY-3-8912 da General Instrument, também conhecido como PSG. Esse chip tinha capacidade de produzir três canais de som e um canal de ruído, sendo que o ruído era usado para fazer efeitos sonoros (como tiros e explosões), e também para fazer sons de percursão (como baterias).

Para testar se a minha implementação do canal de ruído estava boa, eu fiz o seguinte programinha em MSX BASIC:

10 SOUND 7,183
20 SOUND 8,15
30 FOR I=8 TO 1 STEP -1
40 SOUND 6,I
50 PRINT I
60 FOR J=0 TO 300: NEXT J,I
70 BEEP

Os números mágicos no código assustam, mas o funcionamento é simples: o registro 7 é o mixer, com o valor de 183 eu ligo o gerador de ruído no canal A e desligo todo o resto. O registro 8 é o volume do canal A, que eu setei para 15, ou seja, volume máximo. Por fim, eu faço um loop alterando o valor do registro 6, que é a frequência do ruído.

Portanto, o que se espera desse programinha é que ele toque um barulhinho de volume constante, começando grave e ficando cada vez mais agudo. Isso na teoria. Na prática, aconteceu outra coisa!

O bug


Eu resolvi repetir o experimento da época e filmar. No video abaixo você pode ouvir esse programinha rodando em um MSX real, e depois rodando no meu emulador (eu usei o dosbox para compilar a versão de 1998, antes do bug ser corrigido):




Olha que curioso, embora eu tenha programado o volume para ficar constante, no MSX real o volume abaixa sensivelmente no final. E pior, no emulador isso não acontece!

A diferença é sutil, será que eu estava viajando? Um jeito de confirmar é capturando os dois sinais e jogando no Octave para comparar (quer dizer, na época ainda não tinha o Octave, então eu usava o Matlab mesmo):


De olho parece realmente que o turboR vai diminuindo, mas é difícil comparar. Um jeito mais eficiente é notando que a grandeza que perceptualmente nós sentimos como volume, fisicamente é causada pela potência do sinal. E potência a gente sabe calcular, é a integral do quadrado da forma de onda. Jogando isso no Octave:


Olha lá! Não era viagem não, nos dois últimos passos realmente o volume do turboR fica menor que o do emulador! Mas por que isso acontece? Talvez eu tivesse alguma pista analisando o hardware com cuidado.

O circuito


Para entender o funcionamento do PSG, não adianta usar o esquemático do MSX e nem pegar o datasheet do AY-3-8912, porque nenhum dos dois explica o que acontece dentro do chip. Em teoria daria para abrir o chip e tentar ler a pastilha com um microscópio eletrônico, mas como eu não tinha acesso a um, o jeito foi apelar para engenharia reversa mesmo.

Os canais de som do PSG foram fáceis de modelar. Ele tem três geradores de onda quadrada, e cada um funciona da seguinte forma:


O clock de entrada do contador é o clock do MSX dividido por 16, ou seja, 223506Hz. A cada pulso, o contador incrementa de um, e o resultado é comparado com o valor do registro que determina a frequência. Quando o valor é igual, o contador reseta e um flip-flop tipo T inverte de estado, gerando assim a onda quadrada desejada.

Mas e o gerador de ruído? Na época ninguém sabia como funcionava, uns chutavam que era ruído térmico de resistor, outros que era shot noise na junção de algum semicondutor.

Mas a verdade era muito mais simples, o PSG tinha um gerador pseudo-aleatório! Para ser mais exato, tinha um gerador de onda quadrada igual ao descrito acima, e a cada transição na saída ele fazia uma iteração em um Linear Feedback Shift Register.


Nós sabemos que todo LFSR gera um sinal periódico, então, capturando um período inteiro do gerador, dá para deduzir qual é o polinômio característico dele. Fazendo o experimento, o polinômio era o maximal de tamanho 17, ou seja, x17+x14+1.

Eu conhecia a forma de onda gerada, e conhecia até o exato LFSR que era usado. Ainda assim, isso não explicava por que o volume abaixava. Será que eu tinha implementado errado?

A implementação


Um dos meus requisitos era rodar o emulador a toda velocidade nas máquinas da época. Em 1998, isso significava emular esse circuito inteiro em um Pentium 1, e ainda sobrar tempo para emular a CPU e o VDP. A coisa era tão apertada que simplesmente escrever em Assembly não era suficiente, eu tinha que garantir que o inner loop rodasse inteiramente em registros.

O problema é que os meros 7 registros de uso geral do x86 não eram suficientes. Minha primeira solução foi apelar feio: durante a emulação do PSG, eu desligava as interrupções e usava o stack pointer como registro! Na minha casa isso funcionava bem, mas em casa eu usava DOS. A maioria dos meus usuários à época usava Windows 95, e aí essa abordagem de mexer no esp fazia a máquina travar.

Vamos para a solução dois então. Tinha um monte de variáveis que mudavam a cada nota tocada pelo PSG, mas eram constantes durante o inner loop. Sendo assim, era uma boa oportunidade para usar código auto-modificável: eu "recompilava" o inner loop cada vez que alguém mudasse um registro do PSG. Essa solução funcionava até no Windows e era rápida o suficiente.

Porém, ela era rápida, mas bastante complexa. Você pode ver esse código no github, não é a coisa mais legível do mundo, nem a mais fácil de entender. Então, seria de se esperar que o volume não estivesse abaixando por conta de algum erro de implementação.

Mas após pensar um pouco, eu concluí que a minha implementação deveria estar correta. Se ela estivesse errada, então esse erro deveria aparecer só no meu emulador. Mas na época exisitiam outros emuladores, e em todos os outros esse mesmo problema aparecia.

Se várias implementações diferentes tinham o mesmo problema, então o erro não era no código, tinha que ser conceitual. E se é conceitual, então talvez eu conseguisse alguma pista analisando matematicamente o sinal.

A transformada


Hora de fazer a modelagem do sinal então! Vamos considerar que os bits de saída do LFSR formam uma sequência discreta {di}, e que a cada bit 1 na entrada nós mandamos para saída uma forma de onda s(t). Dessa maneira, se a largura do pulso for T, então o ruído do PSG pode ser computado como:


Da teoria de processos estocásticos, nós sabemos que a densidade espectral de potência do sinal, nesse caso, vale:


...onde Gd(f) é a densidade espectral de potência da sequência {di}, e S(f) é a transformada de Fourier de s(t). Se admitirmos que o ruído em questão seja ruído branco, então o Gd(f) é uma constante. Já a forma de onda s(t) é um pulso retangular, cuja transformada é a função sinc. Portanto, a d.e.p. do sinal tem a forma de um sinc ao quadrado. 

O caso que dá problema no emulador é quando o registro 6 vale 1, ou seja, quando T=8.98µs. Assim, a densidade espectral de potência da saída é a abaixo:


Novamente, aquilo que perceptualmente nós percebemos como volume é a potência do sinal, e, pelo teorema de Parseval, a potência é a área abaixo do gráfico acima.

Depois de olhar um tempo para esse gráfico, finalmente caiu a ficha!


A solução


No fim, o problema todo estava entre a cadeira e o teclado! Mais especificamente, no fato de que tinha um ser humano entre a cadeira e o teclado! O sistema auditivo humano só consegue ouvir freqüências até 20kHz, e esse sinal estava espalhando potência para faixas inaudíveis! De fato, um humano só consegue perceber como volume essa parcela da potência:


A integral do sinc ao quadrado não dá para resolver analiticamente, mas dá para calcular numericamente. Eu posso, por exemplo, fazer uma tabela com a potência audível dividida pela potência total, para cada um dos valores de frequência que eu usei no programinha em BASIC:


Registro 6Porcentagem da potência audível
80.99
70.99
60.99
50.99
40.99
30.96
20.83
10.50


Se você comparar com o gráfico de potência do emulador versus turboR, dá pra notar que a porcentagem acima é exatamente o que falta para tornar a volume do emulador igual à máquina real! Adicionando essa atenuação no emulador, ele finalmente ficou fiel ao original como eu queria.

Leitores que estejam prestando atenção, a essa altura devem estar se perguntando: "Mas peraí, se o problema é puramente sensorial, então porque o gráfico de potência medido no Octave mostra a perda de volume?". A resposta é que isso também é um artefato do sistema auditivo humano, mas indiretamente.

Quando eu pluguei o turboR na entrada de microfone do meu notebook, a placa de som fez a captura do audio em qualidade de CD, ou seja, com amostragem de 44100Hz. E os projetistas da placa de som são espertos, eles sabem que antes de amostrar você precisa passar o sinal em um filtro passa-baixas com corte na frequência de Nyquist, para evitar o aliasing. E essa frequência foi calculada exatamente para jogar fora o que os humanos não ouvem, gerando o mesmo problema.

Se você quiser refazer o experimento, coloquei no github os scripts em Octave que usei para gerar as imagens e tabelas do post:

Scripts em Octave para análise de ruído

Por fim, esse bug deu um trabalhão para resolver, e no final desconfio que a diferença era tão sutil que ninguém além de mim notou a diferença :)

sábado, 2 de junho de 2012

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