Computing PI -- a simple programming problem



#2

A couple of years ago there was a good discussion about calculating
PI
to many places
I'd like to offer a different PI programming
problem/challenge: What's the shortest program that you can write on a
programmable calculator to calculate PI to whatever the natural
accuracy of the calculator is.

Use of built-in transcendental functions is not allowed of course, nor
is input of the digits directly (e.g., 3.14...) nor indirectly (e.g.,
355/113) allowed. The challenge is to come up with a real algorithm
with reasonably fast convergence that terminates when the natural
accuracy of the calculator is reached. Without defining what
"reasonably fast convergence" means, say that it has to complete
within 1 minute on an original 12C. Assume that you know nothing about
the initial state of the calculator other than it's on, in run mode and
is set to show maximal decimal places.

The 12C is a good calculator to use for this exercise since it has no
built-in way to compute PI and has minimal programming capabilities
allowing one to concentrate on finding a good algorithm and not
"tricks". (Besides LBL instructions don't cost programming steps on
the 12C).

As a benchmark, I've found an algorithm that is 13 steps and runs in
about 12 seconds giving an answer of 3.141592649 on the original 12C.
Although the error is about 5e-9, that's due to the accumulated errors
of the 10-digit arithmetic in the calculator. On the latest 12C the
result is: 3.141592654 in 5 seconds, it does 12-digit arithmetic
"under the hood" -- it really comes up with 3.14159265365. (On the
32SII the same algorithm gives 3.14159265357 in about 2 seconds.)


#3

Katie, Longer and slower than the benchmark:

17 liner.
Based on Vieta's infinite product (1593)
A program that has escaped from forthcoming DataFile V25N5.
Runs for 16 seconds and joyfully ends with
3.141592654 in X, Y, Z and T!

Keystrokes |Display |
[f][P/R] | |
[f]CLEAR[PRGM] |00- |
0 |01- 0 |
[ENTER] |02- 36 |
2 |03- 2 |
[ENTER] |04- 36 |
[+] |05- 40 |
[x<>y] |06- 34 |
2 |07- 2 |
[g][LSTx] |08- 43 36 |
[RDN] |09- 33 |
[+] |10- 40 |
[g][SQRT] |11- 43 21 |
[/] |12- 10 |
[g][x<=y] |13- 43 34 |
[g][GTO]00 |14- 43,33 00 |
[g][LSTx] |15- 43 36 |
[x<>y] |16- 34 |
[g][GTO]04 |17- 43,33 04 |
[f][P/R] | |


#4

Dear Sir,
most respected Esquire!

You are a pro, I assume.

Ciao.....Mike

PS: one day I'll publish your PI routines here.


#5

Greetings Mike!
I first saw your pi program a year or so ago.
Do we have an eleven liner posted here yet?
Enjoy the Basle meeting!
Cheers,
Tony


#6

Hi Tony!

Last WE Thomas Rundel took me to Pforzheim where we met Michel Bel (!), Daniel, Gottfried, Gilles, Andrew, Peter, Stefan and some others. I was overwhelmed to meet them all after more than one year in my hermitage of depression.

Allschwil will be the next highlight within few days, I look forward to it since several month. I'll show a working HP Math Expander from a USB key chain, quite nice. Alas I have currently no way to publish it on my site for several reasons.

11 steps to compute PI? Reminds me the old saying about a fisherman's seasons what once was re-written in a PPC-Journal about programmers: first he tries to catch fish, then he tries to catch many fish, next he tries to catch big fish, after that he tries to catch rare fish, and one day he will be glad just to catch fish. With my error-prone brain I gave up all attempts to get the esoteric algorythms of your PI-routines <G>.

Have fun! (it's essential).....Mike


#7

Hullo Mike!

I recognise some famous names at last w/end's palmtop meeting!!!

The Vieta PI is a bit like a fisherman - it uses a transcendental hook (sqrt) to catch one (PI) :-)

Katie's solution has rhythm too, it goes: -STO-STO-STO-GTO-GTO-GTO

Wish I could afford the time and money to simply fly to Allschwil!
The HP Math Expander sounds interesting. Hope it can do PI<VBG>

Cheers,
Tony


#8

Tony,

Quote:
The Vieta PI is a bit like a fisherman - it uses a transcendental hook (sqrt) to catch one (PI) :-)

A cute phrase, but SQRT is the antithesis of a transcendental function.

-Katie


#9

Thanks for the correction Katie.

I was being irrational<G>.

P.S. Still no 11-liner 12C PI here.

Cheers,
Tony


#10

If you're going to be irrational, than I'll be denumerable (up to 11) (and silly) and give an 11-line solution for the 32SII:

[2]
[STO P]
{LBL X]
[SQRT]
[STO/ P]
[2]
[STOx P]
[x=y?]
[VIEW P]
[+]
[GTO X]

You need to [R/S] it with the PRGM at the TOP ([GTO ..] if you're not there. It will halt with PI on the display and you can put PI on the stack by hitting [ENTER]. It's not a total cheat.


#11

Katie,

You have indeed re-gained the lead, and even use a LBL!

So appropriate as the sto* and sto/ , the keys to the shortening, were discovered by you too :-)

Thanks for the challenge.

Cheers,
Tony

#12

Longer and somewhat faster: about 7 seconds on the good old HP-12C (output 3.141592655, FIX 9).

That's a straightforward implementation of the Brent-Salamin Algorithm, as detailed here:

http://mathforum.org/library/drmath/view/58283.html

Keystrokes      |Display       | 
[f][P/R] | |
[f]CLEAR[PRGM] |00- |
1 |01- 1 |
[STO]5 |02- 44 5 |
[STO]0 |03- 44 0 |
[.] |04- 48 |
5 |05- 5 |
[STO]3 |06- 44 3 |
[g][SQRT] |07- 43 21 |
[STO]1 |08- 44 1 |
[RCL]5 |09- 45 5 |
[RCL]1 |10- 45 1 |
[+] |11- 40 |
2 |12- 2 |
[/] |13- 10 |
[RCL]5 |14- 45 5 |
[RCL]1 |15- 45 1 |
[x] |16- 20 |
[g][SQRT] |17- 43 21 |
[x<>y] |18- 34 |
[STO]5 |19- 44 5 |
[RCL]1 |20- 45 1 |
[-] |21- 30 |
[ENTER] |22- 36 |
[x] |23- 20 |
[STO]2 |24- 44 2 |
[x<>y] |25- 34 |
[STO]1 |26- 44 1 |
2 |27- 2 |
[STO][x]0 |28- 44 20 0 |
[RCL]0 |29- 45 0 |
[RCL]2 |30- 45 2 |
[x] |31- 20 |
[STO][-]3 |32- 44 30 3 |
[RCL]2 |33- 45 2 |
[f][RND] |34- 42 14 |
[g][x=0] |35- 43 35 |
[g][GTO]38 |36- 43,33 38 |
[g][GTO]09 |37- 43,33 39 |
[RCL]5 |38- 45 5 |
[ENTER] |39- 36 |
[x] |40- 20 |
2 |41- 2 |
[x] |42- 20 |
[RCL]3 |43- 45 3 |
[/] |44- 10 |
[g][GTO]00 |45- 43,33 00 |
[f][P/R] | |

Thanks, Tony, for the listing utility :-)

I have tested the algorithm on the HP-50G as well (using LongFloat 3.89 Library, by Gjermund Skailand and Thomas Rast) :

1000 places: 04 min 09 sec
2037 places: 17 min 29 sec

Gerson.


Edited: 14 Oct 2006, 11:53 p.m.


#13

Gerson,

WOW! Indeed 7 seconds. Fast!

Cheers,
Tony


#14

It would be interesting to know the theory underlying the algorithm though. I must confess I have no idea why the method works... Other methods, like arctangent series for instance, are slower but they are easier to be understood.

There is a mistake in the key codes in line 37. These are the correct codes:


[g][GTO]09 |37- 43,33 09 |

For some reason the line number codes in lines 36 and 37 were missing in the output file. I made the mistake while inserting them manually.

I'd like to thank Katie and all contributors to this thread. Pi has always been a fascinating subject to me.

Best regards,

Gerson.


#15

This algorithm is actually a way to compute arctangent based on using arithmetic-geometric mean

Viete's infinite-product formula that Tony demonstrated is the same one that I used by chance and with his help we now have it down to just 12 instructions on the 12C. It's based on inscribed polygons and is nicely described in Petr Beckmann's A History of PI in chapter 9. A must have book if you're into PI.

I'll post the 12 instruction solution soon, but I'm hoping that someone will beat that.

Edited: 15 Oct 2006, 1:33 p.m.


#16

Thanks for the links and the book recommendation.

Gerson.

#17

I know you're looking for a shorter program, not necessarily a faster one. The following is even longer (47 steps) and takes about 9 seconds to give 3.141592654 (on the HP-12C simulated in Nonpareil). I think you'll find this interesting though. To my knowlegde, I am the first one to compute pi on a calculator using this series. :-)

The series is not so efficient but works for a few digits.

Gerson.


01 0
02 STO 2
03 1
04 STO 1
05 STO 0
06 1
07 CHS
08 STO* 1
09 2
10 STO+ 0
11 RCL 0
12 11
14 y^x
15 1/x
16 RCL 1
17 *
18 STO+ 2
19 RND
20 x=0
21 GTO 23
22 GTO 06
23 RCL 2
24 1
25 +
26 148635648
35 *
36 505.21
42 /
43 11
45 1/x
46 y^x
47 GTO 00


#18

I adapted your listing to run on an HP-41CX. I am impressed with the accuracy and speed. Using FIX 5, the program loops twice and produces a value for Pi that is equal to the built-in value!!! Usin FIX 4 or less, the program loops once and yields a value of Pi that has good accuracy (but does no equal to the built-in value).

Sooooooooooo ... now I am REALLY curious. What is the equation or p-code behind your method?

Namir

Edited: 16 Oct 2006, 2:32 p.m.


#19

Hello Namir,

I found this almost by chance. We know this very inefficient series:

1 -  1/3  +  1/5 -  ... =  pi/4
I changed the powers of the denominators and pi just to see what I would come up with:
1 - 1/3^3 + 1/5^3 - ... = k1 * pi^3

1 - 1/3^5 + 1/5^5 - ... = k2 * pi^5

1 - 1/3^7 + 1/5^7 - ... = k3 * pi^7

k1 = 1/32;

k2 = 5/1536;

The first two constants were obtained experimentally (using ->Q on the HP-48, or the HP-32SII). But I couldn't find a pattern... until I found this page:

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

What matters are the constants in formula (9). Follow the links (Sloane's A046976 and A053005) for more constants: 61/184320, 277/8257536, 50521/14863564800...

It appears the first few digits quickly obtained quickly, but then the rythm slows down. You can check it by yourself.

Gerson.


Edited: 16 Oct 2006, 3:43 p.m. after one or more responses were posted


#20

Your calculator code implements the following equation:

K1 = 148635648

K2 = 505.21

Pi = { K1/K2 * [ 1 + sum((-1)^i / (1 + 2i)^11) ]} ^(1/11) for i = 1,2,3,..., n

I checked it with a BASIC listing and it works.

Namir

Edited: 16 Oct 2006, 3:41 p.m.


#21

Quote:
Pi = { K1/K2 * [ 1 + sum((-1)^i / (1 + 2i)^11) ]} ^(1/11) for i = 1,2,3,..., n

You're right!

As I mentioned earlier, the table below demonstrates that fewer correct digits are obtained as the number of terms increases. The equation is in the HP-200LX Solver format. As you can see, the series is as slowly convergent as the Gregory-Leibniz series, despite the first three iterations giving twelve correct digits. That's why I said I had a nifty algorithm to better solve the first decimals. The first few decimals I meant :-)

P=(K1/K2*(1+SIGMA(I,1,N,1,(-1)^I/(2*I+1)^E)))^(1/E)

K1 = 14863564800
K2 = 50521
E = 11

N = 1 => P = 3.141592647876891
N = 2 => P = 3.141592653725998
N = 3 => P = 3.141592653581560
N = 4 => P = 3.141592653590661
N = 5 => P = 3.141592653589660
N = 6 => P = 3.141592653589820
N = 7 => P = 3.141592653589787
N = 8 => P = 3.141592653589795
N = 9 => P = 3.141592653589793

Anyway, these series are interesting as they can be expressed in terms of the Dirichlet Beta Function:

1 - 1/3n + 1/5n - 1/7n + ... = Beta(n);   n = 1, 3, 5...

Gerson.

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

I've just found this site:

http://home.att.net/~numericana/answer/analysis.htm

There we can read

Quote:
...

Every other time, such an integration gives an exact expression for the alternating sum of some new power of the reciprocals of odd integers.

...


which explains it all...


Edited: 16 Oct 2006, 10:28 p.m.

#22

Quote:
Usin FIX 4 or less, the program loops once and yields a value of Pi that has good accuracy (but does no equal to the built-in value)

Thanks! I hadn't noticed that. I don't have a 12C here, but on Nonpareil it takes less than 5 seconds to show 3.141592648, very close to the value obtained by Katie's program. Now all we have to do is squeezing the program to make it fit in 12 steps or less :-)

Gerson.


#23

Ok, here it is in 12 steps (11 seconds) on an original 12C, with a full precision accuracy 3.141592654:

Credit goes to Tony(nz) for writing this modified version of the one I started with. He got rid of one line by changing the order of the [STO/] and [STOx] thus avoiding the [1/x] that my code had, this switch also improved the accuracy. He also changed the initial 2 lines to shave off 1 second (I had used CLEAR SIGMA/SIGMA+ to accumlate the result in register 1).

The idea is the same however: Use a register to accumulate the product and keep doing [2][+][SQRT] on the stack until X equals 2.

[f][P/R]        |              | 
[f]CLEAR[PRGM] |00- |
2 |01- 2 |
[STO]0 |02- 44 0 |
[g][SQRT] |03- 43 21 |
[STO][/]0 |04- 44 10 0 |
2 |05- 2 |
[STO][x]0 |06- 44 20 0 |
[g][x<=y] |07- 43 34 |
[g][GTO]11 |08- 43,33 11 |
[+] |09- 40 |
[g][GTO]03 |10- 43,33 03 |
[RCL]0 |11- 45 0 |
[g][GTO]00 |12- 43,33 00 |
[f][P/R] | |

What's particularly neat about 12 lines is that it takes the same number of lines to write a program directly entering the digits:

[3][.][1][4][1][5][9][2][6][5][4][GTO 00]

Do we have a winner?

Edited: 16 Oct 2006, 9:15 p.m.


#24

Really nice!

I am not so good in writing short programs, as you can see. However, I have an obsessin for speed :-) (My second program runs in 10 seconds on the real thing. Six seconds using FIX 5, that would be cheating though.)

Too bad no one counted the letters in the words of my reply to Geir Isene. The phrases don't make much sense, but they aren't so bad for a non-native English speaker, are they? The one I like best is this:

"How I want a drink, alcoholic of course, after the heavy lectures involving quantum mechanics. All of thy geometry, Herr Planck, is fairly hard!"

I don't drink, but this always came to my mind during Quantum Mechanics classes... (Fortunately they were not so many.) :-)


Edited: 16 Oct 2006, 9:46 p.m.


#25

Gerson,

I think that everyone reading this thread figured to count the letters in your quotes. I've always found the phrases to remember things by harder than straight memorization of digits, spectrum colors, trigonometric formulas,... yadda, yadda, yadda. But they are often clever and fun to think up.

I didn't realize that you had developed this one yourself:

Quote:
Hey, I have a nifty algorithm to better solve the first decimals. Hopefully someone somewhere got to the solution also.

That's excellent for anyone, non-native English speaker or otherwise!

-Katie

#26

Katie,

Looks like you are using the following algorithm:

Q(1) = sqrt(2 + sqrt(2))

P(1) = 4 * sqrt(2) / Q(1)

For n = 2,3,4,…, until Q(m)->2

Q(n) = sqrt(2 + Q(n-1))

P(n) = 2 P(n-1) / Q(n)



Then Pi=P(m)



Namir

Edited: 17 Oct 2006, 11:15 p.m.


#27

Yes, that's it. Viète's formula is 400+ years old!

Edited: 17 Oct 2006, 11:39 p.m.


#28

Thanks!!! I needed a reference for the algorithm.

#29

Hi Namir,

For your algorithm Q(1)=0 and P(1)=2 would do. We used Q(1)=sqrt(2) and P(1)=4/Q(1) as it meant less lines :-)

Cheers,
Tony


#30

Thanks Tony for the simplification of the initial values.

The algorithm is a real gem.

Namir

#31

Thanks Katie and all for this wonderful thread. Nice to know an efficient way to calculate Pi.

Namir


#32

That would be interesting.

Gene


#33

100?? No use in just 100 decimal places. You can easily keep that in your head. You'll need a few hundred to help you on your way to pirvana.


#34

Hey, I have a nifty algorithm to better solve the first decimals. Hopefully someone somewhere got to the solution also.

(Actually the first nineteen decimals... :-)


--------

Thirty is better than the old nineteen now we realize amusingly!

Gerson.


Edited: 15 Oct 2006, 10:22 p.m.

#35

Quote:
You'll need a few hundred to help you on your way to pirvana.

A few hundred? Check out this guy...


#36

Quote:
TOKYO — A Japanese mental health counselor recited pi to 100,000 decimal places from memory on Wednesday, setting what he claims to be a new world record.

It's often said that people become psychologists/psychiatrists/etc. because they are in need of one themselves. Does this prove the point! :)


#37

Quote:
It's often said that people become psychologists/psychiatrists/etc. because they are in need of one themselves. Does this prove the point! :)

I know I should be more thick-skinned at times, and I have a very good sense of humour and try not to take myself too seriously, but I must say that I bristle at stuff like this. First of all, it is not generally true. Second, and most importantly for me, it demeans and stigmatizes those with bona fide psychiatric disorder (and having a love affair with Pi, though odd to some, does not count) and those of us who have the occasionally misunderstood and thankless job of being of assistance to them.

Just some food for thought.

Regards,

Leslie C. Wright, MD, MEd, FRCPC

Specialist in Adult General, Forensic, and Correctional Psychiatry


#38

Les,

I'm really sorry that my comment offended you, it was entirely in jest. I've been in therapy for many, many years and have gotten a lot out of it. I'm a huge proponent of psychiatric and related professions and thought about going into the field myself, many times! I'm more poking fun at myself with that comment than anyone else -- I've done a little bit of memorizing PI too!

-Katie

Edited: 18 Oct 2006, 11:13 a.m.


#39

Katie,

While studying for my master's n chem engineering in Ann Arbor, I remember memorizing the natural log of 100. One day, one of my professors was interactively solving a problem in class. One result ended in ln(100). He somewhat sarcasticaly said and the value of ln(100) is? The digits just came out of my mouth!!! He was taken back and stopped for few seconds. After that He smiled, since he realized he had "ok" students in his class.

And as for going to counseling, I highly recommend it myself. I have benefited from it greatly and my son has from the positive outlook it gave me on life.

I really have enjoyed the challenge. It has given me a chance to play with various calculators by keying in the different code versions.

Namir

#40

Les:

Now we learn that you are a software expert... albeit one specialized in a very particular platform; not the easiest to debug, indeed! :-)

Best regards from somone who very much appreciate your profession and the help it may give to us all (even those of us who are just slightly weird!!)

Closer to calculators, I may say that "setting your printer to TRACE" may help in finding what is wrong with your programs, even at a paper cost. Talking from a couch may have a similar effect.

#41

On an HP-33s try:

Z0001  LBL Z
Z0002 5419351
Z0003 ENTER
Z0004 1725033
Z0005 /
Z0006 R/S
That's just six steps on an HP-33s but sixteen steps on an HP-12.

It works on all my calculators ranging from the TI-59 and TI-85 to the HP-12C, HP-28C, HP-41C and HP-49+.


#42

Palmer,

You can save 3 steps on the 12C with 104348/33215 (ref: Datafile V24N1P11).

Cheers,
Tony


#43

You can find the same ratio in the following excerpt from the thread "CALCULATING MANY DIGITS OF PI" whhich was posted on March 19, 2004:

"The discussion in V10N4P26 [of TI PPC Notes] also proposed a simple-minded search which would compare the decimal equivalent of fractions against the value of pi (or any other number to be converted) loaded to the accuracy of the machine. Start with the fraction 1/1. If the value of the fraction is greater than the decimal to be converted then add a 1 to the denominator and try again. If the value of the fraction is less than the number to be converted then add a 1 to the numerator and try again. Sample programs for a CC-40 were included. Use of the simple minded search on a TI-95 yields two fractions of interest. 5419351/1725033 = 3.14159 26535 8985 ... which becomes 1E-12 less than the defined value on machines that truncate to thirteen digits such as the TI-59 and TI-66, and which becomes the defined value on machines that round to thirteen digits such as the TI-95 and 6565759/2089946 = 3.14159 26565 9009 ... which becomes the defined thirteen digit value whether truncated or rounded. The same technique on an HP-41 yields fractions such as 104348/33215 and 209051/66543 which will yield the ten digit value used for pi on that device."

V10N4 of TI PPC Notes was published in late 1985. The same article presented the ""Hicks" method which was to divide 2143 by 22 and take the square root two times. That is ten steps which will yield results on ten digit machines such as the HP-41 and HP-12C which are within one digit in the last place.

The 104348/33215 ratio does NOT yield the correct result on twelve digit machines such as the HP-28 and hp 33s. The result will be in error by 3.3E-10. Of course, there is no need for a means to find pi on machines which provide it.


Possibly Related Threads...
Thread Author Replies Views Last Post
  [OT] Mathematica free for Raspberry PI BruceH 32 6,000 11-23-2013, 05:24 AM
Last Post: Nick_S
  Simple Tetris. free for you to improve on cyrille de Brébisson 3 1,474 11-20-2013, 05:43 PM
Last Post: Erwin Ried
  Computing pi with the PC-1300S Kiyoshi Akima 0 772 11-17-2013, 12:24 AM
Last Post: Kiyoshi Akima
  [HP-Prime] Simple Game (Bugs) CompSystems 1 949 11-01-2013, 10:18 AM
Last Post: Han
  Calculating Pi LHH 9 2,127 09-27-2013, 10:50 PM
Last Post: Gerson W. Barbosa
  Visualization of pi Bruce Bergman 13 2,763 08-17-2013, 05:00 PM
Last Post: Howard Owen
  Simple Math Question Namir 2 1,029 08-09-2013, 06:13 PM
Last Post: Eddie W. Shore
  HP-42S with Electroluminescent screen and simple I/O port Jose Poyan 8 2,043 03-27-2013, 07:11 PM
Last Post: Jose Poyan
  Simple sample programs for the HP-41CX? Tom Lewis 5 1,481 03-25-2013, 07:11 PM
Last Post: Allen
  OT: Happy Pi Day! Eddie W. Shore 13 2,591 03-22-2013, 10:44 AM
Last Post: Les Koller

Forum Jump: