BASIC, Házení jehlou (diskuse)

[Josef Kufner]

Což sinus, ale dal bych tam desetinou tečku místo čárky.

[fakenickname]

Ha! Tak chybu na řádku 60 jsem už vyřešil. Překladatel byl až příliš svědomitý a přeložil desetinou tečku na desetinou čárku. Teď mi straší řádka 70 a její IF.

> blassic monte_carlo 

Syntax horror in 70

[Josef Kufner]

70 IF X > U THEN 90

[fakenickname]

výborně! už to jede! díky!

[trekker.dk]

nemá tam být GOTO 90 ?

[Robo]

Tiez sa mi zda, ze tam ma byt:

IF X > U THEN GOTO 90

[michich]

V mnoha dialektech BASICu je možné GOTO za THEN vypustit. V některých za THEN ani nic jiného být nemůže.

[fakenickname]

sice to jede ale výstup je dost divný :( místo krásného pí vypadávají nějaké dvojky .. :(

[Josef Kufner]

Asi to bude rozbitý...

[Lu-Tze]

Že by proto, že se tam nikde nepřiřazuje do fí?

[Lu-Tze]

Tohle jsem teda nepochopil, počítá to pí a pracuje to se sin a pí? WTF? přece rovnoměrně [0-1] střelím x-ovou souřadnici a y-souřadnici a podívám se jestli mi to spadlo do toho jednotkovýho čtvrtkruhu nebo ne, takhle získám aproximaci obsahu toho čtvrkruhu (počet pokusů co se trefily/počet včech pokusů), tohle vydělim 4 a mám pí

[Lu-Tze]

na to "podívání se" mi stačí umocňování, žádnej sin (x^2 + y^2 <= 1)

[Käyttäjä 11133]

Taky se mi ta implementace zdá nějaká divná :)

[Käyttäjä 11133]

Mě se mnohem více líbí tento vztah pro výpočet. :)

[michich]

No to bude jiný pokus: Buffonova jehla

[fakenickname]

jj to je ono

[fakenickname]

dokonce je tam parádní javová app http://www.angelfire.com/wa/hurben/buff.html

[Andrej]

Z Archa jsem vydoloval blassic (a classic interpreter of basic), naťukal source do souboru a hodil mu to.

Možná stojí za úvahu použít skutečný kompilátor - FreeBasic. S různými těmi interprety nemám dobré zkušenosti.

[spang]

Chyba je nejspíš ve funkci RND (resp. její implementaci v blassicu, která není zdokumentovaná). To RND(X) a RND(F) totiž, zdá se, vrací pořád nulu. Kdyby to vracelo náhodný float od 0 do 1, tak to půjde. Zajímavé věci to dělá, pokud RND(X) nahradíš třeba RND(1) a RND(F) nahradíš RND(2), ale správně to pořád nevypadá

Tady je funční ekvivalent v pythonu:

import random, math

n = 0
print "No. OF THROWS ; PI"
for j in xrange(1,24):
  for k in xrange(1,500):
    x = random.random()
    u = math.sin(math.pi*random.random())
    if x < u:
      n += 1
  t = 500*j
  p = 2*float(t)/float(n)
  print t, ";", p

[fakenickname]

poklona pane dobrá práce, díky moc

[fakenickname]

Zajímavé. Někdy kolem roku 1970~1980 tento experiment zabral cca 53 sec počítačového času (On the Librascope of the Electr. Engin. Depth, Univ. of Colorado). Ještě dříve to pár nadšenců dělali jako rekreační sport :) A dneska..

real    0m0.053s
user    0m0.036s
sys     0m0.000s

[Käyttäjä 11133]

Zdá se ti to rychlé? CPU podtaktované na 800MHz.

marek@mantisha:~/pi> time python pi.py
10000 ; 3.11041990669

real    0m0.125s
user    0m0.108s
sys     0m0.016s
marek@mantisha:~/pi> gfortran pi-jehla.f95
marek@mantisha:~/pi> time ./a.out
 ***   3.122561

real    0m0.009s
user    0m0.004s
sys     0m0.004s
marek@mantisha:~/pi>

PROGRAM PI_JEHLA
integer                     ::    i = 1,pocet_iteraci = 10000,vnejsi_citac = 24, splneno = 0, foo = 0
double precision            ::    pi = 0, random = 0, x = 0, u = 0
double precision,parameter  ::    const_pi = 3.141592653589793238462643383279
!!print*,"Napiste pocet iteraci"
!read(*,*)pocet_iteraci
 vnitrni:  DO
  call random_number(random)
  x = random
  call random_number(random)
  u = random
  if(x < sin(u*const_pi)) splneno=splneno+1
  i=i+1
!  if(splneno>0) print*,i," -- ",(2*real(i)/real(splneno))
  if(i==pocet_iteraci) EXIT vnitrni
 END DO vnitrni
  print*,"***",(real(2*pocet_iteraci)/real(splneno))
END

Místní guruové by to jistě ještě nějak zoptimalizovali.

[Käyttäjä 11133]

Použitý pythoní kód

import random, math

n = 0
for k in xrange(1,10000):
  x = random.random()
  u = math.sin(math.pi*random.random())
  if x < u:
    n += 1
t = 10000
p = 2*float(t)/float(n)
print t, ";", p

[rad.ovan@seznam.cz]

    Předpokládám že ten program je původně psaný v GW-BASICu od Microsoftu, takže chybka by byla skutečně v té funkce RND. Ona se tam totiž chová dost divně, RND nebo RND(cokoliv) vrací následující číslo v posloupnosti, zatímco RND(0) vrací to poslední použité. No, co by člověk chtěl od Gatese, žejo.
    Chvilku jsem si s tím hrál v Chipmunk BASICu, a tam to fungovalo když se v obou případech použilo RND(1), které vrací číslo mezi 0 a 1. Tady je výsledek na 1200MHz P3:

radovan@PC1:~/basic$ time basic montecarlo.bas
No. OF THROWS   PI
500     3.125
1000    3.154574
1500    3.144654
2000    3.167063
2500    3.134796
3000    3.176284
3500    3.146067
4000    3.145891
4500    3.144654
5000    3.143666
5500    3.138374
6000    3.138896
6500    3.129514
7000    3.139717
7500    3.127606
8000    3.121342
8500    3.115265
9000    3.119584
9500    3.133245
10000   3.133323
10500   3.124535
11000   3.126332
11500   3.133942
12000   3.136025

real    0m0.041s
user    0m0.036s
sys     0m0.004s

    Jen mě tak napadlo, když se počítá poměr zásahů uvnitř a vně kruhu, nebylo by lepší, než používat umělou "náhodu", prostě udělat mřížku a vzít to bod po bodu? Zkusím na tom po obědě trochu zapracovat.

P.S. Syntax horror je fakt pěkné hlášení

[Radovan Garabík]

nebylo by lepší, než používat umělou "náhodu", prostě udělat mřížku a vzít to bod po bodu?

- bolo by to v tomto prípade jednoduchšie, ale:

• v prípade zložitých mnohodimenzionálnych problémov môže byť vytvorenie mriežky väčší problém, než náhodné strieľanie
• mriežka by spôsobovala systematickú chybu, buď by bol výsledok vždy podhodnotený alebo nadhodnotený, podľa toho, ako by si rátal body po obvode kruhu. Vyrovnať túto chybu si vyžaduje dodatočnú analýzu, napr. ako presne je zakrivená kružnica v danom bode mriežky - a zložitosť takej analýzy je asi porovnateľná s priamym výpočtom π "iným spôsobom" - tak to už rovno môžeme použiť ten "iný spôsob"

[Kyosuke]

Snad jen nedává moc smysl zvyšovat počet lineárně...zákon velkých čísel je neúprosný:

kyosuke@gondolin:~> ruby <<EOF
> require 'narray'
> def pi(n)
>         x, y = Array::new(2){NArray::float(n).random}
>         return 4 * ((x**2 + y**2) < 1).mean
> end
>   puts   "      N   |    Pi   \n  =================="
> 14.times do |x| 
>   x = (x%2==0 ? 1 : 3) * 10**(x/2)
>   printf " %8d | %f\n", x, pi(x)
> end
> EOF
      N   |    Pi   
  ==================
        1 | 4.000000
        3 | 4.000000
       10 | 2.800000
       30 | 2.800000
      100 | 3.040000
      300 | 3.200000
     1000 | 3.056000
     3000 | 3.150667
    10000 | 3.134400
    30000 | 3.150133
   100000 | 3.146520
   300000 | 3.137933
  1000000 | 3.141056
  3000000 | 3.141111
kyosuke@gondolin:~>

Možná právě proto se to takhle nepočítá...

[fakenickname]

jak dlouho to jelo?

[Kyosuke]

Nevím, asi něco pod dvě sekundy? Není to žádný zázrak, kdybych to chtěl rychlý, nepíšu to v Ruby, ne? A hlavně bych nepoužíval tenhle algoritmus...

[Käyttäjä 11133]

Kdyby si to chtěl rychlé tak by si to napsal ve fortranu IV?

[Kyosuke]

Ne, použil bych GMP...

[Käyttäjä 11133]

Není GMP náhodou knihovna pro práci s velkýma číslama?

[Kyosuke]

Jasně, a chytřejší než její autoři nebudeme. Takže až budu potřebovat hodně dlouhé pí, vezmu ten ukázkový zdroják ke GMP a nebudu si hrát se špendlíkama.

[Käyttäjä 11133]

Když budu potřebovat hodně dlouhé pí tak napíšu echo "scale=$kolik_des_cisel; 4*a(1)" | bc -l Upřímně, kolik lidí kdy v životě použilo z nějakého praktického důvodu více než 20 číslic z desetinného rozvoje čísla pí? :) Myslím, že jich moc nebude to by mělo stačit i pro mezihvězdné vzdálenosti.

[Lu-Tze]

Já myslím, že je dobré, že pí i é jsou 3, dobře se to pamatuje

[Käyttäjä 11133]

Jo, jo pí to je něco kolem tří.

[Kyosuke]

Mno, jednak to má pedagogickou hodnotu, jednak tam pořád ještě můžou být schovaný všechny ty MP3ky, co ještě nemám. Jen je najít, no.

[Michal Vyskočil]

Proč GMP? Mrknu do wikipedie a najdu si tam hezký odkaz na první 4 milióny čísel z Pi. To by mělo uspokojit většinu požadavků

[Käyttäjä 11133]

Michale, Michale ty mi kazíš celou radost.

[Kyosuke]

Ale to není ono, tam se nic nekompiluje ani nedebuguje, kam zmizela dětinská radost z úspěchu?

[Michal Vyskočil]

Inženýrský přístup, hledá se nejjednodušší řešení, které ještě funguje (v překladu, jsem líný jako veš a chci to mí rychle)

[Kyosuke]

(Hmm, když jsem to z hnusného Ruby přepsal do ještě hnusnějšího Lispu, prošlo to celé asi za čtvrt sekundy. Mno jo, kompilátor.)

[rad.ovan@seznam.cz]

    Matematika není zrovna moje silná stránka, přiznávám. Tu mřížku jsem myslel takhle nějak:

10 d=54: z=d^2
20 for x=-d to d
30   x2=x^2
40   for y=-d to d
50     if x2+y^2<=z then n=n+1
60   next y
70 next x
80 print n/z
90 exit

    Při přibližně stejném počtu pokusů to je dokonce rychlejší, s podobnou přesností:

radovan@PC1:~/basic$ time basic pi_mrizka.bas
3.136145

real    0m0.022s
user    0m0.020s
sys     0m0.004s

    Ale pravda je, že abych z toho dostal správně aspoň pět míst, musel jsem dát d=10000 a počkat si deset minut. On BASIC nikdy nebyl z nejrychlejších, už kvůli tomu počítání s reálnými čísly V tom Chipmunku jsem ještě nenastudoval jak to přepnout do celých, neznáte někdo obdobu DEFINT A-Z? Pak by ten čas mohl být zajímavější.

[spang]

http://en.wikipedia.org/wiki/Computing_pi -- To je dobrý článek. Je tam popsána i tato metoda.

Tady je python verze, jako bonus tu kružnici nakreslí do konzole (což má smysl jen dokud se vleze jeden řádek do konzole, tzn. pro r do 40):

import sys

r = 10
n = 0
for y in xrange(-r,r+1):
  for x in xrange(-r,r+1):
    if (x**2 + y**2)**0.5 <= r:
      sys.stdout.write(".")
      n += 1
    else:
      sys.stdout.write(" ")
    if x == r:
      sys.stdout.write("\n")
print "pi =", float(n)/float(r**2)

To mi připomíná jeden pěkný bash skript co jsem někde viděl -- do konzole takhle kreslil Mandelbrotův fraktál .

[fakenickname]

hezké :)