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

13 comentários:

  1. Apesar de "not delta" e "delta == 0" serem equivalentes neste caso, eu acho mais legível usar a igualdade, afinal delta é numérico e nunca será None. Prefiro deixar "not variavel" para quando a variável é booleana ou quando pode ser vazia ou None.

    Aliás, delta pode ser número em ponto flutuante, o que implica que comparações de igualdade raramente são uma boa ideia. Neste caso, parece ser inofensivo, mas ainda assim é uma questão a ser notada!

    ResponderExcluir
  2. Para reforçar que essas perdas de precisão não é algo artificial, lembro que quando paguei I.A. tive problemas com isso também. Nunca chegava nos valores adequados do bias no algoritmo de retro-alimentação. Para resolver, substitui o float pela tipo Decimal.

    --
    adorilson

    ResponderExcluir
    Respostas
    1. O problema do Decimal é que ele é bem mais lento que o float nativo. Dependendo do problema, compensa mais consertar a fórmula que apelar para o Decimal.

      Excluir
  3. Esse lance de evitar calcular com a equação tradicional por causa da precisão numérica é mesmo improtante. Eu encontrei esse problema na prática no meu estágio, e aprendi o macete no Numerical Recipes. Mas o que eles ensinam é um pouco diferente, vc primeiro calcula

    q = -1/2*(b+sgn(b)sqrt(b^2-4ac))

    e aí x_1 = q/a e x_2 = c/q. É interessante que nesse caso o x_1 vai pra infinito se a==0, mas o x_2 ainda vai ter a raiz certa.

    ResponderExcluir
    Respostas
    1. Eu também aprendi no Numerical Recipes, esse livro é bom ter em casa para referência :)

      Excluir
  4. ...Pq precisão é uma coisa suuuuper _improtante_!!

    ResponderExcluir
  5. Este artigo é interessante para quem estuda e/ou trabalha com Cálculo Numérico.

    ResponderExcluir
  6. Ricardo, muito bom a sua análise do código. Um "code review" bem comentado e explicado.

    Uma coisa que eu sempre tentei entender é esse problema de precisão com float, mesmo em operações simples. Olha esse exemplo em JS http://stackoverflow.com/questions/588004/is-javascripts-math-broken .

    ResponderExcluir
    Respostas
    1. O segredo do floating point é pensar em binário. 0.3 tem só uma casa decimal na base 10, mas em binário ele é uma dízima periódica. Como o computador não tem memória infinita, você precisa truncar em algum ponto, e aí o que era 0.3 vira alguma coisa do tipo 0.30000000000000004.

      Excluir
  7. Acho fantástico os seus posts! Parabens!

    ResponderExcluir
  8. Engraçado, eu pensava que esse bug da busca binária era uma feature. No meu post sobre busca binária, http://marathoncode.blogspot.com.br/2012/10/busca-binaria.html, eu resolvia o problema aumentando a precisão das variáveis :)

    ResponderExcluir