Una vegada més, després d’una temporada d’exàmens i vacances m’he posat de nou amb la programació.
Parlant amb el meu professor de matemàtiques sobre la fórmula de Ramanujan (de la qual, ja vaig publicar alguna cosa), ell em va ensenyar una forma d’obtenir π que, una vegada vista, és molt més evident que la de Ramanujan o altres d’aquest estil. La idea de com obtenir π és bastant senzilla i evident (si caus en això, és clar), només necessites unes nocions bàsiques (gairebé per intuïció) de probabilitat i dues fórmules matemàtiques que saben fins i tot els nens de primària.
Investigant per la fantàstica e-enciclopèdia (proposo aquest nom a les enciclopèdies electròniques) Wikipedia i altres webs sobre π, em vaig informar que aquest procediment usa el mètode de Montecarlo, que curiosament també fou usat per a la construcció de la primera bomba nuclear; però deixem les anècdotes per a centrar-nos en el que avui, després d’un llarg període de sequera, ompli una pàgina de la meva bitàcora:
Si dibuixem un quadrat de costat 2 (l=2) i inscrivim una circumferència, aquesta tindria un radi d’1. Per tant, l’àrea del quadrat seria de 2^2=4 i la de la circumferència seria de π·1^2=π.
I és ara on entra la petita noció de probabilitat: Quina probabilitat hi ha que al llançar un dard a l’atzar en qualsevol posició del quadrat es clavi en la circumferència?
Doncs evidentment, la probabilitat serà de π/4 (Àrea circumferència/Àrea quadrat). Llavors, si llancéssim un nombre infinit de dards (això és el mètode de Montecarlo en si mateix) obtindríem la quarta part del valor exacte de π. Només quedaria multiplicar aquest valor per quatre per a obtenir el valor de π.
Una vegada plantejat el problema i una idea de com resoldre’l, em vaig posar mans a l’obra… i bo, això és el que he fet. No sé si serà un gran assoliment (no crec, ja que he vist molts programes d’aquest estil), però la sensació que produïx el resoldre un problema és suficient per a intentar resoldre’l.
A continuació, descriuré el codi que usa el programa escrit en Pascal. He de dir que el programa es podria haver fet per a consola i segur que es pot optimitzar, però em va semblar més bonic veure com el programa va acolorint la circumferència en cada interacció.
Primer de tot, definir el tipus TModPoint, que no és ni més ni menys que un punt en l’espai personalitzat per a admetre decimals, així:
TModPoint = record
X: Real;
Y: Real;
end;
Una altra part més interessant és la funció que utilitzo per a convertir un tipus TPoint a TModPoint. El tipus de Point no és sols que una estructura que conta amb una variable X i una Y, però cal saber, que en la informàtica les coordenades tenen diferents eixos que els tradicionals: En la pantalla, la component X d’un punt es conta del cantó superior esquerra (0) fins a la dreta de la pantalla; i la component I d’un punt es conta des del mateix cantó esquerra (0) cap avall. Així, no existeixen nombres negatius per a representar les coordenades, i a mi em feien falta per a esbrinar si un punt P(X,I) podia estar dintre de la circumferència o no. Per a això, abans cal fer-se altra pregunta: Com puc esbrinar si un punt està situat dintre d’una circumferència de radi X?
Partint de les idees bàsiques de trigonometria, si tenim una circumferència de radi 1, per a un valor X en aquesta component, el valor màxim d’Y perquè estigui dintre de la circumferència és de Sin(ArcCos(X)). Poso un dibuix (fet amb el Paint a última hora) per a explicar-me millor:
Llavors, per a aquest quadrant (i per al de la seva esquerra), la component Y tindrà un valor màxim de Sin(ArcCos(X)) per a ser vàlida perquè el punt P(X,Y) estigui dintre de la circumferència.
La idea està clara, no? Ara només hem de convertir les coordenades de la pantalla en coordenades vàlides per a usar les funcions trigonométricas.
Si el radi del cercle és de Z píxels, llavors, la component X=(Z-200)/200. Bueno, això es fa amb aquesta petita funció:
function PointToModPoint(Point: TPoint): TModPoint;
var
X, Y: Real;
ResPoint: TModPoint;
begin
X := (Point.X-200)/200;
Y := (200-Point.Y)/200;ResPoint.X := X;
ResPoint.Y := Y;Result := ResPoint;
end;
I una vegada tenim les components X i Y transformades per a treballar amb un cercle de radi, només hem de comprovar que la component Y sigui vàlida per a la component X donada. Atents a la funció següent:
function CheckPoint(ModPoint: TModPoint): Boolean;
begin
if ((ModPoint.Y >= 0) and (ModPoint.X >= 0)) or ((ModPoint.Y >= 0) and (ModPoint.X then
begin
if ModPoint.Y end else
if ((ModPoint.Y and (ModPoint.X or ((ModPoint.Y and (ModPoint.X >= 0)) then
begin
if ModPoint.Y >= -Sin(ArcCos(ModPoint.X)) then
Result := True else Result := False;
end;
end;
Perfecte, ja tenim tot el necessari, ara sols falta unir-ho tot, per exemple en l’event OnClick d’un botó:
procedure TForm1.Button1Click(Sender: TObject);
var
Point: TPoint;
I, Counter, Max: Integer;
begin
Randomize;
Image1.Canvas.Fillrect(Image1.Canvas.ClipRect);
Image1.Canvas.Ellipse(0,0,400,400);
Counter := 0;
Max := StrToInt(FloatToStr(JvSpinEdit1.Value));
for I := 1 to Max do
begin
Point.X := Random(400);
Point.Y := Random(400);
if CheckPoint(PointToModPoint(Point)) = True then
begin
Counter := Counter+1;
Image1.Canvas.Pixels[Point.X, Point.Y] := clRed;
end
else
Image1.Canvas.Pixels[Point.X, Point.Y] := clBlue;
Label2.Caption := Format(’Tirs: %d || Encerts: %d’, [I, Counter]);
Label2.Repaint;
Image1.Repaint;
end;
ShowMessage(FloatToStr((Counter/Max)*4));
end;
Bueno, aquí teniu una imatge del programa i a continuació un enllaç amb l’arxiu Pi_Montecarlo.zip que conté tot el codi i l’executable (Acabo de posar per separat l’arxiu .pas del programa principal).
![[FSF Associate Member]](http://www.elvoldelhomeocell.net/blog/wp-content/uploads/fsf_member.png){kind=small}