Simo Kivelä pyysi laskinvalmistajia simuloimaan Buffonin neulaongelmaa omilla tuotteillaan. En ole laskinvalmistaja, mutta olen jossain yhteydessä suositellut Pythonia, joten tässä pientä testailua IPython Notebookilla.

In [1]:

%pylab inline
rcParams['figure.figsize'] = (8,8)

Populating the interactive namespace from numpy and matplotlib

Piirretään muutama yhdensuuntainen viiva. Sovitaan, että neulan toinen pää osuu yksikköneliön sisään, kuvassa keltaisella. Tehdään kuvan piirtämisestä saman tien funktio, jotta sitä voidaan hyödyntää myöhemmin.

In [2]:

def pohja():
    hlines([-1, 0, 1, 2], -1, 2)
    bar(0, 1, 1, color='yellow', alpha=0.2)
    xlim([-1.1, 2.1])
    ylim([-1.1, 2.1])
pohja()
show()

Piirretään kuvaan myös muutama mahdollinen neula. Jos alkupään koordinaatit ovat vaikka $(x, y)$ ja suuntakulma on $\theta$ , loppupään koordinaatit ovat $(x + \cos\theta, y + \sin\theta)$ .

In [3]:

def neula(x, y, theta):
    plot([x, x + cos(theta)],
         [y, y + sin(theta)],
         color='blue', lw=2)

pohja()
neula(0.5, 0.5, 0)
neula(0.1, 0.1, 1.1 * pi)
neula(0.2, 0.9, pi/6)
neula(0.8, 0.3, 3*pi/8)
show()

Ilmeisesti neula pysyy kahden viivan välissä, jos $0 < y + \sin\theta < 1$ . Tehdään tästä funktio:

In [4]:

def valissa(y, theta):
    return 0 < y + sin(theta) < 1

def ylittaa(y, theta):
    return not valissa(y, theta)

Sitten vain toistetaan monta kertaa. Tehdään tästäkin funktio, jotta voidaan kohta mitata aikaa.

In [5]:

def simulaatio(n):
    ylityksia = 0
    n = 1000000
    for i in range(n):
        y = random.uniform(0, 1)
        theta = random.uniform(0, 2*pi)
        if ylittaa(y, theta):
            ylityksia += 1
    return 1.0 * ylityksia / n

print simulaatio(1000)
print simulaatio(10000)
print simulaatio(1e6)

0.63641
0.637155
0.636206

Aikaa menee minun koneellani pari sekuntia:

In [6]:

%timeit simulaatio(1e6)

1 loops, best of 3: 1.84 s per loop

Nopeutta saa lisää vektoroimalla. Seuraava tuhlaa muistia mutta välttää silmukoiden ajamisen Python-tulkissa. Vielä vähän hienompi ratkaisu pilkkoisi tehtävän sopivan mittaisiin paloihin, joista kukin olisi erikseen vektoroitu, ja yhdistäisi tulokset Python-silmukalla.

In [7]:

def simulaatio2(n):
    y = random.uniform(0, 1, size=n)
    theta = random.uniform(0, 2*pi, size=n)
    sin(theta, out=theta)
    add(y, theta, out=y)
    return logical_or(greater(y, 1), less(y, 0)).mean()
print simulaatio2(1e8)
%timeit simulaatio2(1e8)

0.63663671
1 loops, best of 3: 6.37 s per loop

In [ ]: