Visualizing _PI challenge

[Jack002]

I haven't figured out how far my string math routines would go, I just know if it was accurate and infinite, it would take about 40 years to surpass the record pi calculation to date, of 31.4 trillion digits. Anyone know where I can find a computer that will last 40 years, without a %^&^$# Windows update?

Pete

That computer Deep Thought in Hitchhiker's Guide to the Galaxy never had these problems. :-)
 

[jack]

in the spigot routine, if the dimension in pi is greater than 9830 then I get an unhandled error #9, subscription out of range, anybody know why?
not that anybody would use that routine to calculate pi to that many places.

[codeguy]

Wasn't certain where this arbitrary precision Pi calculation algorithm should go.

Code: QB64: [Select]

1. DECLARE SUB atan239 (denom&)
2. DECLARE SUB atan5 (denom&)
3. DECLARE SUB PrintOut (words%)
4.  
5. 'Program to calculate pi, version 4.8
6. 'A major rewrite of version 4.2, this uses only two arrays instead of
7. 'three, and includes a host of speedups based on a similar C program.
8. 'A sampler: all the carries are reserved until the end, the divide and
9. 'add routines are combined, two terms are added at a time, and the number
10. 'of function calls is minimized. It's a big change for a small gain, since
11. 'the compiled version requires 28.6 seconds for 5000 digits on my 486 66MHz
12. 'computer, a 10% gain over version 4.2; like before, it's capable of about
13. '150,000 digits of pi.
14. '
15. 'This program has come a long way from version 1.0; thanks are due to
16. 'Larry Shultis, Randall Williams, Bob Farrington and Adrian Umpleby.
17. 'One final note for speed freaks: this program will run about 6 times faster
18. 'if written in C using an optimizing compiler. Likewise, if you can figure
19. 'out a way to do integer division and use both the quotient and remainder,
20. 'this program can easily be sped up by 25%.      -jasonp@isr.umd.edu
21.  
22. DEFINT A-Z
23. CLS
24. INPUT "how many digits"; digits&
25.  
26. words = digits& \ 4 + 3
27. DIM SHARED sum&(words + 1), term(words + 1)
28. start! = TIMER
29.                                        
30.                                          '---------------16*atan(1/5)
31. denom& = 3: firstword = 1: lastword = 2
32. sum&(1) = 3: term(1) = 3: sum&(2) = 2000: term(2) = 2000
33.  
34. DO UNTIL firstword >= words
35.    CALL atan5(denom&)
36.    denom& = denom& + 2
37. LOOP
38.  
39.                                          '------------   -4*atan(1/239)
40. denom& = 3: firstword = 2: remainder& = 4
41.  
42. FOR x = 2 TO words
43.    dividend& = remainder& * 10000              'crunch out 1st term
44.    term(x) = dividend& \ 239&
45.    remainder& = dividend& - term(x) * 239&
46.    sum&(x) = sum&(x) - term(x)
47. NEXT x
48.  
49. DO UNTIL firstword >= words
50.    CALL atan239(denom&)
51.    denom& = denom& + 4
52. LOOP
53.                                      
54.  
55. FOR x = words TO 2 STEP -1                        '-------finish up
56.    IF sum&(x) < 0 THEN                                  'release carries
57.       quotient& = sum&(x) \ 10000                       'and borrows
58.       sum&(x) = sum&(x) - (quotient& - 1) * 10000
59.       sum&(x - 1) = sum&(x - 1) + quotient& - 1
60.    END IF
61.    IF sum&(x) >= 10000 THEN
62.       quotient& = sum&(x) \ 10000
63.       sum&(x) = sum&(x) - quotient& * 10000
64.       sum&(x - 1) = sum&(x - 1) + quotient&
65.    END IF
66. NEXT x
67.  
68. CALL PrintOut(words)
69. PRINT "computation time: "; TIMER - start!; " seconds"
70. END
71.  
72. '------------------------------------------------------------------
73. SUB atan239 (denom&)
74. SHARED words, firstword
75.  
76. remainder1& = term(firstword)                        'first divide implicitly
77. remainder2& = 0: remainder3& = 0: remainder4& = 0
78. denom2& = denom& + 2: firstword = firstword + 1
79.  
80. FOR x = firstword TO words
81.    temp& = term(x)
82.    dividend& = remainder1& * 10000 + temp&
83.    temp& = dividend& \ 57121
84.    remainder1& = dividend& - temp& * 57121
85.  
86.    dividend& = remainder2& * 10000 + temp&
87.    temp2& = dividend& \ denom&
88.    remainder2& = dividend& - temp2& * denom&
89.    sum&(x) = sum&(x) + temp2&
90.  
91.    dividend& = remainder3& * 10000 + temp&
92.    temp& = dividend& \ 57121
93.    remainder3& = dividend& - temp& * 57121
94.  
95.    dividend& = remainder4& * 10000 + temp&
96.    temp2& = dividend& \ denom2&
97.    remainder4& = dividend& - temp2& * denom2&
98.    sum&(x) = sum&(x) - temp2&
99.    term(x) = temp&
100. NEXT x
101.  
102. firstword = firstword + 1
103. IF term(firstword) = 0 THEN firstword = firstword + 1
104.  
105. END SUB
106.  
107. '-------------------------------------------------------------------
108. SUB atan5 (denom&)
109. SHARED words, firstword, lastword
110.  
111. FOR x = firstword TO lastword + 1
112.    temp& = term(x)
113.    dividend& = remainder1& * 10000 + temp&
114.    temp& = dividend& \ 25
115.    remainder1& = dividend& - temp& * 25&
116.    term(x) = temp&
117.  
118.    dividend& = remainder2& * 10000 + temp&
119.    temp& = dividend& \ denom&
120.    remainder2& = dividend& - temp& * denom&
121.    sum&(x) = sum&(x) - temp&
122. NEXT x
123.  
124. FOR x = lastword + 2 TO words
125.    dividend& = remainder2& * 10000
126.    temp& = dividend& \ denom&
127.    remainder2& = dividend& - temp& * denom&
128.    sum&(x) = sum&(x) - temp&
129. NEXT x
130.  
131. IF term(lastword + 1) > 0 AND lastword < words THEN lastword = lastword + 1
132. IF term(firstword) = 0 THEN firstword = firstword + 1
133.  
134. denom& = denom& + 2
135. remainder1& = 0: remainder2& = 0
136.  
137. FOR x = firstword TO lastword + 1
138.    temp& = term(x)
139.    dividend& = remainder1& * 10000 + temp&
140.    temp& = dividend& \ 25
141.    remainder1& = dividend& - temp& * 25&
142.    term(x) = temp&
143.  
144.    dividend& = remainder2& * 10000 + temp&
145.    temp& = dividend& \ denom&
146.    remainder2& = dividend& - temp& * denom&
147.    sum&(x) = sum&(x) + temp&
148. NEXT x
149.  
150. FOR x = lastword + 2 TO words
151.    dividend& = remainder2& * 10000
152.    temp& = dividend& \ denom&
153.    remainder2& = dividend& - temp& * denom&
154.    sum&(x) = sum&(x) + temp&
155. NEXT x
156.  
157. IF term(lastword + 1) > 0 AND lastword < words THEN lastword = lastword + 1
158. IF term(firstword) = 0 THEN firstword = firstword + 1
159.  
160. END SUB
161.  
162. '------------------------------------------------------------------
163. SUB PrintOut (words)
164. PRINT "pi = 3+."
165. FOR i = 1 TO words \ 3
166.    PRINT " ";
167.    PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(3 * (i - 1) + 2))), 4);
168.    PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(3 * (i - 1) + 3))), 4);
169.    PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(3 * (i - 1) + 4))), 4);
170.    IF i MOD 5 = 0 THEN PRINT "  :"; 12 * i
171. NEXT i
172.  
173. PRINT " ";
174. FOR i = 3 * (words \ 3) + 2 TO digits
175.    PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(i))), 4);
176. NEXT i
177. PRINT : PRINT
178.  
179. END SUB
180. [code=qb64]
181.  

[Pete]

in the spigot routine, if the dimension in pi is greater than 9830 then I get an unhandled error #9, subscription out of range, anybody know why?
not that anybody would use that routine to calculate pi to that many places.

Sorry, never saw the movie, if there was one, or read the book, but great, now I miss Blondie. :(

Pete

[Pete]

in the spigot routine, if the dimension in pi is greater than 9830 then I get an unhandled error #9, subscription out of range, anybody know why?
not that anybody would use that routine to calculate pi to that many places.

Sure, i is dim as an integer and at that many iterations. i surpasses the 32767 limit. QB then returns a -32768 for i, which throws the subscript out of range error. Try DIM as _INTEGER64 instead of INTEGER for all dimension integers.

Pete

[doppler]

My PI challenge is

Apple
Blue berry
Cherry
and sometimes Peach.

Sorry I can't resist.

[Pete]

Q: What's the difference between a round pi and a square pie? A: $6.72.

A round pi is 3.14, and a square pi is 9.86

Ah, pi humor. And now back to our regularly scheduled program-ing...

Pete

[jack]

thanks Pete :)

[Bert22306]

So, pi humor aside, I thought searching for a fraction of two rational integers that approximates pi was kind of interesting. Of course, a fraction of two rational integers will always be a rational number, so it'll never be exact. So I wrote a search routine, and here are the first few results. I won't list the intermediate results except for the first two, just to show how the routine works. (Not terribly efficient, since I see the value 3.14 showing up repeatedly, as I increase the numerator or the denominator, in that last trial especially.)

Turns out, 355/113 is a pretty decent rational approximation! You have to go a lot ,longer to gain an extra decimal of accuracy.

Actual pi = 3.141592653589793238....

Iteration 1    1 / 1 = 1
Iteration 2    2 / 1 = 2
Iteration 3    3 / 1 = 3
Accuracy requested = 1
-------------------------------------
Iteration 1    1 / 1 = 1
Iteration 2    2 / 1 = 2
Iteration 3    3 / 1 = 3
Iteration 4    4 / 1 = 4
Iteration 5    4 / 2 = 2
Iteration 6    5 / 2 = 2.5
Iteration 7    6 / 2 = 3
Iteration 8    7 / 2 = 3.5
Iteration 9    7 / 3 = 2.333333333333333
Iteration 10   8 / 3 = 2.666666666666667
Iteration 11   9 / 3 = 3
Iteration 12   10 / 3 = 3.333333333333333
Iteration 13   10 / 4 = 2.5
Iteration 14   11 / 4 = 2.75
Iteration 15   12 / 4 = 3
Iteration 16   13 / 4 = 3.25
Iteration 17   13 / 5 = 2.6
Iteration 18   14 / 5 = 2.8
Iteration 19   15 / 5 = 3
Iteration 20   16 / 5 = 3.2
Accuracy requested = .1
----------------------------
Iteration 28   22 / 7 = 3.142857142857143
Accuracy requested = .01
----------------------------
Iteration 264   201 / 64 = 3.140625
Accuracy requested = .001
----------------------------
Iteration 438   333 / 106 = 3.141509433962264
Accuracy requested = .0001
----------------------------
Iteration 467   355 / 113 = 3.141592920353983
Accuracy requested = .00001
----------------------------
Iteration 467    355 / 113 = 3.141592920353983
Accuracy requested = .000001
----------------------------
Iteration 100122   75948 / 24175 = 3.141592554291623
Accuracy requested = .0000001