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 nα é 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)

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, x¹⁷+x¹⁴+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 {d_i}, 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 G_d(f) é a densidade espectral de potência da sequência {d_i}, 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 G_d(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 6	Porcentagem da potência audível
8	0.99
7	0.99
6	0.99
5	0.99
4	0.99
3	0.96
2	0.83
1	0.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 :)

A Fórmula de Bhaskara

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

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

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

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

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

Primeiro problema: API inconsistente

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

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

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

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

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

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

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

Segundo problema: Bhaskara quem?

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

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

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

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

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

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

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

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

Terceiro problema: Instabilidade numérica

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

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

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

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

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

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

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

Conferindo:

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

Nice!

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

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

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

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

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

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

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

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

Conferindo:

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

Agora sim!

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

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

A Terra Plana

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

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

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

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

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

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

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

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

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

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

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

O Método do Guerreiro

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

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

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

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

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

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

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

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

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

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

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

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

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

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

O Método do Mago

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

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

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

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

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

Juntando tudo, temos:

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

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

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

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

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