Pi Day's ramblings



#32

Hi all, long time no see:

Today's Pi Day, so I feel like indulging in some quick Pi antics to commemorate the event, namely showing off two little-known, rarely seen, extremely short ways to compute this most ubiquitous constant.

I'll use the HP-71B (+ MathROM) to enter and run the algorithms described. If you'd like to key them in as shown and don't own an HP-71B, get the excellent freeware emulator Emu71 by J-F Garnier (Google for it) or any other capable emulator of your choice (Math ROM emulation required) or that failing try and adapt the code to your favorite HP calculator.

Let's see:

  1. An asymptotically exact, probabilistic way:

    This 50-byte one-liner will compute Pi using random numbers, you'll have to input how many tries (say 100,000) and it'll display the resulting Pi approximation at the end. The more tries the more correct Pi digits you'll get:

        10 INPUT K @ N=0 @ FOR I=1 TO K @ N=N-MOD(IROUND(RND/RND),2) @ NEXT I @ DISP 1-4*N/K
    To run it, key it in and then key in the following at the ">" prompt:
        >DESTROYALL @ RANDOMIZE 260 @ FIX 4    [Enter]
    >RUN [Enter]

    ? 100000 [Enter]

    3.1416

    As I said, the more tries the more correct digits you'll get. See if you can figure how and why it works, and if so check if you can prove that it actually produces the exact value of Pi in the limit.

  2. A faster, very precise, approximate way:

    Do the following:

    a) Compute the real root of the cubic equation x3 - 6*x2 + 4*x - 2:

        >DESTROYALL @ STD
    >FNROOT(-10,10,FVAR^3-6*FVAR^2+4*FVAR-2)

    5.31862821775

    b) raise it to the 24th power:
        >RES^24

    2.62537412641E17

    c) compute the natural logarithm of the result:
        >LN(RES)

    40.1091699911

    d) Divide the result by the square root of 163:
        >RES/SQR(163)

    3.14159265359

    You'll get Pi correct to all displayed digits and then some. Again, see if you can figure how and why this works and if so, check if you can prove that it won't produce an exact value of Pi but just a fairly good approximation.

That's all. Best regards from V.


#33

Hi Valentin. Good write-up. Give me a shout!


#34

Quote:
Hi Valentin. Good write-up. Give me a shout!

Thanks, Gene, I'll give you "a shout" this next weekend for sure, sooner if I can manage.

Best regards from V.

#35

Hi Valentín, welcome back!!!

I will try the two programs in my 71+MathROM, which I haven't used in ages. It should be fun.

Regards


#36

Hi, Fernando !! :D

Quote:
Hi Valentín, welcome back!!!

I will try the two programs in my 71+MathROM, which I haven't used in ages. It should be fun.


Yes, it should.

You might also want to try this Java version for the first one as well:

	private static void pi() {
double pi = 0.0;
for (long i = 0; i < Long.MAX_VALUE; i++) {
pi = pi - Math.round((Math.random()/Math.random()))%2;
if ((i&0x80000)==0x80000) {
System.out.printf("%d: %s%n", i, 1-4*pi/i);
}
}
}

I left it running yesterday on one of my servers, all night long (16 hours), and it produced the following:

...
487,446,282,239 tries -> 3.141592254139215
487,447,330,815 tries -> 3.141592246519439
487,448,379,391 tries -> 3.141592254556738
...

which is as expected, 7 correct digits for some 500 (US) billion tries. Only one core was used and it managed 22,000 tries per millisecond. Not bad.

The result also shows that Java's random number generator is working Ok, else there would be some bias which would be noticeable in the final result's accuracy.

Best regards from V.

#37

Regarding the second challenge: that's not too hard to solve once you got access to some symbolic math software - Wolfram Alpha will do. ;-)

Determine the real solution of the polynomial, raise the result to the 24th power, take the natural log and divide by sqrt(163), and you'll get a result, whose continued fraction representation looks quite familiar:

[3; 7, 1, 15, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, ...]

Right, these are the first sixteen terms of the continued fraction representation of Pi. Its first seventeen significant digits agree with the true value.

I'm sure someone else will provide a decutive proof. ;-)

Dieter

#38

Valentin wrote: "that it won't produce an exact value of Pi"

What is this exact value anyway? :-)


#39

Simple: the exact value of Pi is... Pi. :-)

Dieter

#40

Hello Valentin,

Welcome back!

I've used WolframAlpha to solve the equation ( solve x^3 - 6*x^2 + 4*x - 2 = 0 for x ) and asked for more digits and repeated the steps using the windows calculator:

5.31862821775018565910968015331802247^24

262537412640768767.99999999999922

40.109169991132519846765626111088

3.1415926535897932456228612777983

When I compared the result to pi
3.1415926535897932384626433832795 

I found the first 17 digits to match.

However, if 24 is subtracted from the result between steps b and c, the approximation is even more impressive:

ln(5.31862821775018565910968015331802246772191980883690026023^24-24)/sqrt(163) =

3.14159265358979323846264338327950319783820424737150850414

3.141592653589793238462643383279502884197169399375105820974945 (pi)

I have no idea why this works though, even after checking this MathWorld reference:

http://mathworld.wolfram.com/RamanujanConstant.html

Best regards,

Gerson.


#41

.

Hi, Gerson:

Quote:
Hello Valentin,

Welcome back!


Thanks.

Quote:
I've used Wolfram Alpha to solve the equation ( solve x^3 - 6*x^2 + 4*x - 2 = 0 for x ) and asked for more digits and repeated the steps using the windows calculator:

However, if 24 is subtracted from the result between steps b and c, the approximation is even more impressive:


Yes, I knew. However, I chose to omit the subtraction of 24 as the final result was close enough to Pi that it would nevertheless be indistinguishable in 10-digit and 12-digit calculators, and you wouldn't be able to notice it unless you resorted to multiprecision computations which would deny the simplicity, so I opted for one step less.

Best regards from V.

.

#42

The first one implies that the probability of one (evenly distributed between 0 and 1) random variable divided by another having a result between (.5 to 1.5) or (2.5 to 3.5) or (4.5 to 5.5) or (6.5 to 7.5) or ... is

(pi - 1) / 4.

Wikipedia talks about ratio distributions, and in fact that one in particular here.

Using the result there, you then have to integrate between the limits in the first paragraph, giving

1/4 + 1/2*((1 - 2/3) + (2/5 - 2/7) + (2/9 - 2/11) + (2/13 - 2/15) +...

(the first term is the hardest one: You have to break up the first integral into one from 0.5 to 1, of 1/2; and one from 1 to 1.5, of 1/(2x^2). The others are all of 1/(2x^2), for the whole range)

or -1/4 + ( 1-1/3+1/5-1/7+...)

That infinite series is the Leibniz formula for pi, so the probability is

-1/4 + pi / 4.

which is what was to be shown.


#43

By the way, I do think that is a more interesting example than the well-known Buffon needle problem. That is supposedly an example of pi coming up in a non-obvious situation, but I don't agree. It's obvious that when you toss a needle, the possible angles the needle could trace out when it lands will make a circle. And in fact you use circles and geometry to derive the probability.

In this example, however, it's really not obvious how pi is involved. And geometry is not used in the derivation.


#44

Quote:
By the way, I do think that is a more interesting example than the well-known Buffon needle problem.

Agreed. As far as obtaining Pi via Monte-Carlo methods (algorithms using random numbers) I think that this is by far the simplest and most elegant, apart from being utterly non-obvious.

Best regards from V.

.

#45

Quote:
The first one implies that the probability of one (evenly distributed between 0 and 1) random variable divided by another having a result between (.5 to 1.5) or (2.5 to 3.5) or (4.5 to 5.5) or (6.5 to 7.5) or ... is

(pi - 1) / 4.


Correct.

Best regards from V.

#46

How about using the following loop to get an estimate:

    RAN#
RAN#
R->P
IP
SIGMA+


Any takers to finish this?


Had an open day when I was studying mathematics at university. We did this dartboard PI estimate and the value was running a bit low. In stepped one of the older tutors who proceeded to correct the estimate a bit by putting a pile of darts on the same spot :-)


- Pauli


#47

HP-33s program

P0001 LBL P
P0002 STO N
P0003 CLSIGMA
R0001 LBL R
R0002 RANDOM
R0002 RANDOM
R0004 y,x->theta,r
R0005 IP
R0006 SIGMA+
R0007 RCL N
R0008 n
R0009 x<y?
R0010 GTO R
R0011 1/x
R0012 x<>y
R0013 SIGMAx
R0014 -
R0015 *
R0016 4
R0017 *
R0018 RTN

P: LN= 9 CK=B364
R: LN=66 CK=E2BB

Regards,

Gerson.


#48

Average of 14 runs of 7 XEQ P:

44/14 or 22/7 (YMMV)

Free42 might be an option for much faster results:

00 { 31-Byte Prgm }
01 LBL "PI"
02 STO 00
03 CLSIGMA
04 LBL 00
05 RAN
06 RAN
07 ->POL
08 IP
09 SIGMA+
10 RCL 00
11 RCL 16
12 x<y?
13 GTO 00
14 1/X
15 X<>Y
16 RCL 11
17 -
18 *
19 4
20 *
21 RTN
22 .END.

#49

May I humbly suggest some improvements? :-)

Since the number of loops is known I think using DSE is the obvious way to go here. When the program finishes the values for n and N are the same, so there is no need to distinguish between both. And since only the sum of all x-values is used, the relatively slow Sigma+ is not required either.

Here's a version for the HP-35s which also shows how to get along without the infamous R->P command. ;-)

P001 LBL P
P002 INPUT N ; get and save n
P003 STO K ; initialize loop counter
P004 CLX
P005 STO S ; initialize sum
P006 RANDOM
P007 RANDOM
P008 i
P009 x
P010 +
P011 ABS ; = sqrt(random1^2 + random2^2)
P012 IP
P013 STO+ S ; add 0 or 1
P014 DSE K ; decrement and check loop counter
P015 GTO P006
P016 RCL N
P017 RCL- S
P018 RCL/ N
P019 4
P020 x ; result = 4*(n-sum)/n
P021 RTN

My 35s runs 1000 loops in 2:45 minutes.

Edit: here's another version that requires just one memory register.

P001 LBL P
P002 INPUT N ; get and save n
P003 CLSTK ; clear sum
P004 RCL N
P005 R^ ; save n in stack
P006 RANDOM
P007 i
P008 x
P009 RANDOM
P010 +
P011 ABS ; = sqrt(random1^2 + random2^2)
P012 IP
P013 + ; add 0 or 1
P014 DSE N ; decrement and check loop counter
P015 GTO P006
P016 -
P017 X<>Y
P018 /
P019 4
P020 x ; result = 4*(n-sum)/n
P021 RTN

Dieter

Edited: 16 Mar 2011, 3:15 p.m.


#50

Quote:
May I humbly suggest some improvements? :-)

Of course, and don't be humble: your improvements are perfect!

I imagined Sigma+ was slow but I wanted to use all the instructionS Paul had provided. Prior to using Sigma+, I tested the idea in the third program below, which runs 1000 loops in 1 minute and 46 seconds (at first, I had used an HP-15C). As a comparison, the HP-33s program above runs in 2 minutes and 1 second, using Sigma+; I was expecting it to take more time.

The first program below is based in another (not so good) idea I had tried on the 15C. The indirect addressing slows it a bit. The second one is based on your first version, except the R->P instruction is still being used. It appears the newer machines don't benefit much from the faster DSE and ISG loop control instructions.

Regards,

Gerson.

T0001 LBL T
T0002 STO N
T0003 CLx
T0004 STO A
T0005 STO B
V0001 LBL V
V0002 RANDOM
V0003 RANDOM
V0004 y,x->theta,r
V0005 IP
V0006 1
V0007 +
V0008 STO i
V0009 STO+(i)
V0010 DSE N
V0011 GTO V
V0012 2
V0013 RCL* A
V0014 RCL+ B
V0015 RCL/ A
V0016 8
V0017 /
V0018 1/x
V0019 RTN

2' 13"


T0001 LBL T
T0002 STO N
T0003 STO A
T0004 CLx
T0005 STO S
V0001 LBL V
V0002 RANDOM
V0003 RANDOM
V0004 y,x->theta,r
V0005 IP
V0006 STO+ S
V0007 DSE A
V0008 GTO V
V0009 RCL N
V0010 RCL- S
V0011 RCL/ N
V0012 4
V0013 *
V0014 RTN


2'03"


T0001 LBL T
TO002 STO N
T0003 CLx
T0004 STO A
T0005 STO B
V0001 LBL V
V0002 RANDOM
V0003 RANDOM
V0004 y,x->theta,r
V0005 1
V0006 x>y?
V0007 STO+ B
V0008 STO+ A
V0009 RCL N
V0010 RCL A
V0011 x<y?
V0012 GTO V
V0013 1/x
V0014 RCL* B
V0015 4
V0016 *
V0017 RTN

1' 46"


#51

Hi Gerson,

Quote:
The first program below is based in another (not so good) idea I had tried on the 15C. The indirect addressing slows it a bit. The second one is based on your first version, except the R->P instruction is still being used. It appears the newer machines don't benefit much from the faster DSE and ISG loop control instructions.

This is an interesting point. Obviously every calculator has its specific "slow" and "fast" instructions. Regarding the 35s I can say that especially numeric constants are processed very slowly. So I often use workarounds like these:

   constant    workaround
(slow) (faster)
1 Clx ex (or simply SGN if x>0)
2 e IP
3 Pi IP
... ...

There also is another advantage: on the 35s the "fast" versions even use less memory.

Back to the various Pi-programs. Looking at the second 35s-version I wondered how fast or slow the 35s might process complex numbers, so I replaced lines P006 to P011 by a straightforward...

P006 RANDOM
P007 x2
P008 RANDOM
P009 x2
P010 +
P011 SQRT
...and guess what: it's much faster! 1000 loops are now done in
2:09.

Dieter


#52

Quote:
Looking at the second 35s-version I wondered how fast or slow the 35s might process complex numbers, so I replaced lines P006 to P011 by a straightforward...
P006 RANDOM
P007 x2
P008 RANDOM
P009 x2
P010 +
P011 SQRT
...and guess what: it's much faster! 1000 loops are now done in
2:09.

Likewise when these instructions replace lines V0002 through V0004 in the third HP-33s program above, the running time drops from 1min 46sec to 1min 04sec. However, I think Paul Dale was concerned about size rather than speed.

-------------

Actually SQRT is not necessary. If both terms of SQRT(random#12 + random#22) > 1 are squared, the inequality still holds. Thus, when SQRT is removed the running time for 1000 loops drops down to 53 seconds for the HP-33s program and to about 2 minutes for your 35s program.


Edited: 16 Mar 2011, 9:54 p.m.


#53

Quote:
However, I think Paul Dale was concerned about size rather than speed.

I wasn't. I just wanted to see what the folks here could do to the dart board approach to finding Pi :-)

I haven't been disappointed by the interesting programs.


- Pauli

#54

Right, the square root is obsolete - that's what I also realized after my last message yesterday.

The modified program now takes 1:55 minutes. So it's shorter and faster. ;-)

Dieter


#55

Hello Dieter,

DSE is really slower on the HP-33s, as we can see from the timings belows. Is it also slower on the HP-35s? Thanks!

HP-33s & HP-32SII

T0001 LBL T T0001 LBL T
TO002 STO N TO002 STO N
T0003 CLx T0003 CLx
T0004 STO A T0004 STO A
T0005 STO B T0005 STO B
V0001 LBL V V0001 LBL V
V0002 RANDOM V0002 RANDOM
V0003 x^2 V0003 x^2
V0004 RANDOM V0004 RANDOM
V0005 x^2 V0005 x^2
V0006 + V0006 +
V0007 1 V0007 1
V0008 x>y? V0008 x>y?
V0009 STO+ B V0009 STO+ B
V0010 STO+ A V0010 STO+ A
V0011 RCL N V0011 DSE N
V0012 RCL A V0012 GTO V
V0013 x<y? V0013 4
V0014 GTO V V0014 RCL* B
V0015 1/x V0015 RCL/ A
V0016 RCL* B V0016 RTN
V0017 4
V0018 *
V0019 RTN

1,000 loops timing

0' 53" (HP-33s) 1' 14" (HP-33s)
1' 28" (HP-32SII) 1' 26" (HP-32SII)


HP-42S

00 { 38-Byete Prgm }
01>LBL "PI"
02 STO 00
03 CLX
04 STO 01
05 STO 02
06>LBL 00
07 RAN
08 x^2
09 RAN
10 x^2
11 +
12 1
13 X>Y?
14 STO+ 02
15 STO+ 01
16 DSE 00
17 GTO 00
18 4
19 RCL* 02
20 RCL/ 01
21 RTN
22 .END.

30,000,000 loops timing

1' 10" ( Free42 Binary 1.4.67 @ 1.86 GHz )


#56

On the 35s I excpected a version without DSE to run even slower, as it requires a numeric constant ("1") within the loop. And so it was:

P001 LBL P
P002 INPUT N
P003 CLSTK
P004 RCL N
P005 R^
P006 ENTER
P007 Clx
P008 RANDOM
P009 x^2
P010 RANDOM
P011 x^2
P012 +
P013 IP
P014 +
P015 1
P016 STO- N
P017 Clx
P018 RCL N
P019 X>0?
P020 GTO P007
P021 RDN
P022 -
P023 X<>Y
P024 /
P025 4
P026 x
P027 RTN

Running time: 2:16

As mentioned before there are several ways to avoid numeric constants:

P015 RCL N                   P015 RCL N
P016 SGN P016 ENTER
P017 STO- N P017 SGN
P018 Clx P018 -
P019 RCL N P019 STO N
P020 X>0? P020 X>0?
P021 GTO P007 P021 GTO P007
... ...

Running time: 1:55 Running time: 1:59

So the two last solutions take about the same time as the original version with DSE.

Dieter

#57

Quote:
Average of 14 runs of 7 XEQ P:

44/14 or 22/7 (YMMV)




There's a good point in there. pi is known to be approximately 22/7, or, even better, approximately 355/113.

If you used a multiple of 113 as the number of trials, you might be able to get very close this way.

The same idea is suggested in the article on Buffon's needle I linked above.


#58

Crawl wrote:

Quote:
If you used a multiple of 113 as the number of trials, you might be able to get very close this way.

Since the experiment approximates Pi/4 I would suggest an integer multiple of 4*113 = 452.

This way n*452 trials will result in Pi ~= 355/113 if the number of hits is n*355.

Dieter

#59

Quote:
If you used a multiple of 113 as the number of trials, you might be able to get very close this way.

I had tried that earlier, but to no avail. I was trying single runs of 113 trials. I've finally made it with a single run of 452, as per Dieter's suggestion, the first time I tried. I little cheating though in order to ensure immediate success ;-)

On the real HP-42 and the 38-byte program above:

     94 SEED 452 XEQ PI  =>  3.14159292035  (after 1 min 45 sec)

On Free42 Binary 1.4.67:

66 SQRT SEED 452 XEQ PI  =>  3.14159292035  (instantaneously)

--------

00 { 34-Byte Prgm }
01>LBL "X"
02 STO 04
03>LBL 01
04 RCL 04
05 SQRT
06 SEED
07 452
08 XEQ "PI"
09 PI
10 -
11 ABS
12 1E-4
13 X>Y?
14 STOP
15 DSE 04
16 GTO 01
17 RTN
18 .END.

On the real 42S, line 05 is not necessary to find a "right" seed for the random numbers generator. That's because Free42 will always generate the same random number for seeds with same fractional parts, except when the seed is set to 0, when SEED behaves as it does on the real 42S, that is, always generates internally a different random number:

                         Free42               HP-42S
1 SEED RAN => 5.23054829358E-1 7.31362440213E-1
2 SEED RAN => 5.23054829358E-1 4.31362440213E-1
3 SEED RAN => 5.23054829358E-1 1.31362440213E-1
3.1 SEED RAN => 1.20349678274e-1 8.01362440213E-1
4.1 SEED RAN => 1.20349678274e-1 5.01362440213E-1
5.1 SEED RAN => 1.20349678274e-1 2.01362440213E-1
0 SEED RAN => undetermined undetermined
0 SEED RAN => undetermined undetermined

-----------------

On Free42 Binary 1.4.67 @ 1.86 GHz:

9236 SQRT SEED 132408 XEQ PI  =>  3.14159265301  (about half a second)

Edited: 18 Mar 2011, 10:48 p.m.

#60

I remember a "Spektrum der Wissenschaft" (the german release of "Scientific American") essay about "shooting at pi", using the random generator. If I remember correctly, I used my newly purchased Amiga to do the calculations and graphical representation. Lots of fun.
While googling for the essay, I stumbled across this:
Pi

Looks like someone invested quite a bit of time for pi day ;)


#61

A. K. Dewdney's column "Computer Recreations" in the April 1985 issue of Scientific American discussed this methodology. He used the analogy of a cannon firing into a square field which included an inscribed circular pond. Readers were asked to send the results for 1000 shots to him.

#62

Here's my approach for the second challenge.

The key is a constant found by Ramanujan who showed that

   R = e(Pi * sqrt(163))
is almost an integer. In fact, the exact value for R is
   R = 262 537 412 640 768 743,999 999 999 999 250...
and so
   R = ~= 262 537 412 640 768 744
In other words:
   Pi ~= sqrt(163) * ln R
The value of this expression agrees with Pi in its first 30 decimals.

Now, R also is the intermediate value 2,62537412641 E+17 mentioned in the challenge. At least is looks quite close - the given 12 digits are the same as in R.

So the remaining question is: Why is the real solution of the given cubic equation, raised to the 24th power, so close to R? Or, more precisely, indistinguishable in its first 12 digits.

The analytic solution to the given cubic equation is

   x  =  2 + (5 + sqrt(489)/9)1/3 + (5 - sqrt(489)/9)1/3
= 5,318 628 217 750 185 659 109 680 ...
Now, the 24th power of x is
   x24  =  262 537 412 640 768 767,999 999 999 999 251 ...
which is quite exactly R + 24. Which in turn is close enough to R to that also
   R ~= x24 ~= e(Pi * sqrt(163))
or finally
   Pi ~= ln(x24) * sqrt(163)
which is the value that was calculated in the challenge.

This also explains why the Pi-approximation gets much better if the
intermediate is corrected by 24, as Gerson mentioned.

However, there still is one thing left for the perfect solution:

Assume R resp. R1/24 is given. Now find a cubic equation with integer coefficents where its real solution is as close as possible to this value.

Dieter


Edited: 15 Mar 2011, 1:22 p.m.


Possibly Related Threads…
Thread Author Replies Views Last Post
  An amazing day: Giving a talk at HP about their calculators Geir Isene 9 4,647 12-16-2013, 06:14 PM
Last Post: aurelio
  [OT] Mathematica free for Raspberry PI BruceH 32 8,420 11-23-2013, 05:24 AM
Last Post: Nick_S
  Computing pi with the PC-1300S Kiyoshi Akima 0 1,129 11-17-2013, 12:24 AM
Last Post: Kiyoshi Akima
  Calculating Pi LHH 9 2,891 09-27-2013, 10:50 PM
Last Post: Gerson W. Barbosa
  HHC 2013 Day 2 Highlights Eddie W. Shore 6 2,634 09-23-2013, 04:03 PM
Last Post: Kimberly Thompson
  HHC 2013: Day 1 Highlights Eddie W. Shore 28 8,004 09-23-2013, 03:22 PM
Last Post: Brad Barton
  Visualization of pi Bruce Bergman 13 3,726 08-17-2013, 05:00 PM
Last Post: Howard Owen
  Happy Mother's Day! Eddie W. Shore 1 1,222 05-12-2013, 11:35 AM
Last Post: Walter B
  happy fibonacci day 5/8/13 Allen 8 2,914 05-09-2013, 01:48 AM
Last Post: Gerson W. Barbosa
  Stupid idea of the day: 41-compatible watch bhtooefr 0 1,033 03-31-2013, 08:23 AM
Last Post: bhtooefr

Forum Jump: