▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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.)
▼
Posts: 93
Threads: 2
Joined: Jul 2005
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] | |
▼
Posts: 858
Threads: 80
Joined: Feb 2009
Dear Sir,
most respected Esquire!
You are a pro, I assume.
Ciao.....Mike
PS: one day I'll publish your PI routines here.
▼
Posts: 93
Threads: 2
Joined: Jul 2005
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
▼
Posts: 858
Threads: 80
Joined: Feb 2009
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
▼
Posts: 93
Threads: 2
Joined: Jul 2005
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
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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
▼
Posts: 93
Threads: 2
Joined: Jul 2005
Thanks for the correction Katie.
I was being irrational<G>.
P.S. Still no 11-liner 12C PI here.
Cheers,
Tony
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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.
▼
Posts: 93
Threads: 2
Joined: Jul 2005
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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 93
Threads: 2
Joined: Jul 2005
Gerson,
WOW! Indeed 7 seconds. Fast!
Cheers,
Tony
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Thanks for the links and the book recommendation.
Gerson.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
Yes, that's it. Viète's formula is 400+ years old!
Edited: 17 Oct 2006, 11:39 p.m.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Thanks!!! I needed a reference for the algorithm.
Posts: 93
Threads: 2
Joined: Jul 2005
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
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Thanks Tony for the simplification of the initial values.
The algorithm is a real gem.
Namir
Posts: 2,247
Threads: 200
Joined: Jun 2005
Thanks Katie and all for this wonderful thread. Nice to know an efficient way to calculate Pi.
Namir
▼
Posts: 1,107
Threads: 159
Joined: Jan 1970
That would be interesting.
Gene
▼
Posts: 896
Threads: 183
Joined: Jul 2005
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 727
Threads: 43
Joined: Jul 2005
Quote:
You'll need a few hundred to help you on your way to pirvana.
A few hundred? Check out this guy...
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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! :)
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 1,193
Threads: 43
Joined: Jul 2005
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.
Posts: 901
Threads: 113
Joined: Jun 2007
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+.
▼
Posts: 93
Threads: 2
Joined: Jul 2005
Palmer,
You can save 3 steps on the 12C with 104348/33215 (ref: Datafile V24N1P11).
Cheers,
Tony
▼
Posts: 901
Threads: 113
Joined: Jun 2007
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.
|