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