Estimando Pi com agulhas e linhas

Há alguns dias, assisti a um vídeo muito legal que mostrava que era possível estimar Pi pelo número de vezes que uma agulha cruzasse linhas verticais num tabuleiro no qual as linhas verticais fossem separadas por uma distância igual ao dobro do tamanho da agulha.

Para testar o que o sujeito fez no vídeo com fósforos, no lugar de agulhas, fiz um programinha em Javascript que deixo abaixo. Como o Javascript não é uma linguagem muito apropriada para isso, não esperem muita precisão nem muita performance. No final do post mostro uma solução mais elegante.

Atenção não coloque mais de 1.000.000 de agulhas ou seu browser pode travar.

Esse método de estimar Pi, é baseado no problema da agulha de Buffon, que foi proposto lá pelo século XVIII, pelo matemático francês Conde de Buffon.

A proposição diz o seguinte: Dada uma agulha com comprimento l e um assoalho com tábuas de largura t, a probabilidade de se jogar uma agulha no chão e ela cair sobre uma junção de tábuas é \frac{2l}{t\pi}.

Se l = \frac{t}{2}, então a probabilidade da agulha cair sobre uma junção é \frac{1}{\pi}

Portanto, se jogarmos N agulhas no chão e o número de agulhas que caírem sobre as junções for h, poderemos aproximar \pi pela fração \frac{N}{h}.

A demonstração dessa probabilidade você encontra aqui. Coloquei o link da Wikipedia em inglês porque essa página na Wikipedia lusófona está muito ruim.

Não satisfeito com a implementação porca que eu fiz em Javascript, o Pedro Paulo resolveu fazer uma muito mais eficiente e elegante em python, que deixo abaixo:

#!/usr/bin/env python
import math, random
 
def generateNeedles (n):
    cross = 0
    for i in xrange(n):
        x = random.random()
        theta = random.random()*(math.pi/2.0)
        if ((x + 0.25*math.cos(theta)) >= 1.0):
            cross += 1
        else:
            if ((x - 0.25*math.cos(theta)) <=0.0):
                cross += 1
    return cross
 
needleCount = 30000000
cross = generateNeedles(needleCount)
print cross,needleCount, needleCount / float(cross)

Ao executar o script acima obtive os seguintes resultados:

30.000 agulhas 3.13217790771
300.000 agulhas 3.14261171985
3.000.000 agulhas 3.14226609757
30.000.000 agulhas 3.1423348868

Esse Conde de Buffon deve ter perdido muita agulha deixando cair entre uma tábua e outra do piso de casa, só pode… 🙂

Comments on this entry are closed.

  • Carlo Pires

    Esse é o tipo de problema em python que pypy melhora muito o resultado final. Veja os tempos na minha máquina:
    carlo@carlo:/tmp$ time python2.7 bufon.py
    9547454 30000000 3.14219895692

    real 0m19.013s
    user 0m18.941s
    sys 0m0.004s

    Agora com pypy:
    carlo@carlo:/tmp$ time pypy teste.py
    9550457 30000000 3.14121093891

    real 0m3.457s
    user 0m3.424s
    sys 0m0.020s

    Belo ganho não?

  • Adoro estes posts do Zeletron que envolvem matemática e programação.

  • JB

    Fiz uma versão multiprocessada desse código e roda em menos de 2 segundos num Quad 2 core velho. 🙂

  • JB, cadê a versão multiprocessada?

  • JB

    Tá aqui: http://pastebin.com/rq89LH3y

    Usei a versão do PPJ como base e fiz algumas otimizações além da distribuição de carga entre vários processos. Escolhe o número de workers de acordo com número de núcleos da CPU (se tiver hiper threading, pode colocar o dobro).

    Se usar Pypy, fica ainda mais rápido.

  • Outra coisa: ainda acho que é roubo usar o pi no cálculo do pi 🙂

  • João Bernardo, foi só uma simulação.

    Se você quiser, pode jogar uma agulha no chão 1500 vezes para ver se dá certo. 🙂

  • JB

    Não quero espetar o pé.

    Fiz outra versão desse código usando numpy e numexpr, ficou 100x mais rápido que o do PPJ quando rodei num servidor 8-core HT: 1 bilhão de agulhas em 6.5 segundos contra 10 a 11 minutos

  • Faz com fósforos usados para não espetar o pé.

    Qual foi o resultado de 1 bilhão de agulhas?

    Deixa rodando 100.000.000.000 de agulhas de um dia para outro.

  • Pedro Paulo

    E cadê o código??

  • JB

    Aqui: http://pastebin.com/KX2fNAJz

    Roda assim:
    ./needle.py [numero de agulhas] [numero de processos]

    Coloquei 1 bilhão de agulhas e 32 processos (não muda muito após 16 nesta máquina 8-core HT)

    Só testei em linux e tem uma variável pra dividir o processamento em arrays que caibam na cache da CPU que eu estava usando: Em torno de 8MB.
    Eu faço um pequeno roubo para gerar metade dos números aleatórios, mas isso não influencia o resultado.

    Em torno de 1GB de RAM pra configuração acima.

  • JB

    Não posso deixar rodando que esse servidor faz parte do cluster que eu rodo meus testes de transcodificação de vídeo. Fiz só a noite, quando estava ocioso.

    Pra 1 bilhão ficava entre 3.1413… e 3.1418… as vezes ficava mais próximo, como 3.14160…

    root@c-1-0:~# time python nneedle.py 1000000000 32
    318286688 1000000000 3.14182162717

    real 0m6.403s
    user 1m36.746s
    sys 0m1.684s

  • Coloquei pra rodar 100 bilhões. Deu só 10 minutos 😀

    root@c-1-0:~# time python nneedle.py 100000000000 32
    31831158240 100000000000 3.1415759127

    real 10m3.982s
    user 158m1.473s
    sys 1m50.523s