Keep It Simple, Stupid PerlMonks

### Re: Re: Pi calculator

by Fingo (Monk)
 on Feb 16, 2001 at 01:24 UTC ( #58721=note: print w/replies, xml ) Need Help??

in reply to Re: Pi calculator

Thank you. I did not come up with the method if that's what you are asking. Its based on a circle with a radius of 1 and a square dartboard with each sid being 2. The circle is oin the square. WHen you throw the dart, there is a propability of pi/4 that it will land.... Can you please help me with something? I modified the code to repeat the calculations a certain number of times adding 10 cycles each time. As it got farther than 50,000 it started going into 3.16... while it wa at 3.14 in the 50,000 area. Is this a nuence of the method, or is this have to do with srand not being completly random? Also after 50,000 or so, not once did it generate pi starting with 3.14 ?!?

Replies are listed 'Best First'.
Accuracy of Random Pi
by gryng (Hermit) on Feb 16, 2001 at 06:16 UTC
Running for 10,000,000 twice I got 3.1405 the first run, and 3.1424 the second. (for an average of 3.1415) It's quite likely your rand() isn't perfect. (neither is mine 10,000,000 runs should give me a digit or two more accuracy).

Anyway, here is my favorite approximation for pi, mainly because it only uses the number 2. Even though two is normally a computer friendly number, this algorithm isn't, because it also uses sqrt's. With my perl, 14 iterations is gives the maximum accuracy: 3.14159265480759

```#!/usr/bin/perl -w
use strict;

print "Enter how many iterations:\n";
chomp(my \$i = <>);

my \$x = \$i - 1;
my \$y = sqrt(2);
do { \$y = sqrt(2 + \$y) while (--\$x);
\$y = sqrt(2 - \$y);
} if \$x;

my \$z = \$y * (2 ** \$i);

print "Pi is close to: \$z\n";

So, while not the best, but I have some strange affinity to it. :)

Ciao,
Gryn

p.s. Sorry for the cryptic code for a quick decrypt its: (2**n)*sqrt(2-sqrt(2+sqrt(2+sqrt(2+sqrt(2))))) with the number of 2's inside the sqrt equaling n.

```#!/usr/bin/perl
# This is a quick program to calculate pi using the Monte Carlo method
+.
# I recomend inputting  a value for \$cycles greater that 1000.
# I am working on a detailed explanation of how and why this works.
# I will add it as soon as I'm done.
use strict;
open(PI, ">>pi.dat") || die "pi.dat";
my (\$i, \$j, \$yespi, \$pi) = 1;
my \$cycles = 1;
srand;
while (\$j <= 100000) {
\$cycles = \$j;
while (\$i <= \$cycles) {
my (\$x, \$y, \$cdnt) = 1;
\$x = rand;
\$y = rand;
\$cdnt = \$x**2 + \$y**2;
if (\$cdnt <= 1) {
++\$yespi;
}
\$i=\$i + 10;
\$pi = (\$yespi / (\$cycles / 10)) * 4; # since I add 10 every time.
print PI "\$cycles \$pi\n";
}
\$j = \$j + 10;
}
close(PI) || die "pi.dat";
Here is the modified script I was talking about. If you look at pi.dat it starts getting inacurate at one point. I will try using one of the ways to genereate seeds that are slightly more random than plain srand. Note: you need to make a file pi.dat for the script to run.
I'm not sure if setting srand would really help the situation. In order for the Monte Carlo method to be most effective the random numbers produced must be the most evenly spread out. In fact, you can be more accurate without random numbers and instead going through all four corners, then their mid-points, and their mid-points' mid-points, e.g. :
```.   .      . . .
->  . . .  ->  etc.
.   .      . . .

However, if you would like to use random numbers you need one of two types of random numbers. You either need true even-bias random numbers (pseudo-random can work but often don't :) ). Or you can use a particular kind of random numbers that are designed not to repeat and are garunteed to be spread evenly and more and more closely together (very similar to the systematic approach above).

Two examples are Hamilton and Sobol sequences. If you use them, you get a 1/N convergence, instead of a 1/N^2 (which you get for uniform numbers since there is nothing that keeps them from clumping up). I was going to give below "Numerical Recipe's in C"'s version of Antonov-Saleev's variant on Sobol's sequence. However it's simply too long. Also, I have to wake up for work tomorrow, and the first version I typed in (converting it to Perl from C) didn't work just right. Oh well.

Good luck, Hamilton sequences are much easier to do.

Ciao,
Gryn

Create A New User
Node Status?
node history
Node Type: note [id://58721]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others avoiding work at the Monastery: (4)
As of 2020-10-26 07:17 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My favourite web site is:

Results (250 votes). Check out past polls.

Notices?