"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/PI_DEMO/pi_2.cpp" (21 Feb 2010, 5316 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "pi_2.cpp" see the
Fossies "Dox" file reference documentation.
1 /****************************************************************************/
2
3 /*
4 * PI_2.CPP
5 *
6 * DEMO 'C++' PROGRAM TO COMPUTE PI USING
7 * BORWEIN'S QUARTICALLY CONVERGENT ALGORITHM
8 *
9 *
10 * Output from 'pi-2 200' : (I line wrapped for this page).
11 *
12 * start Borwein-Quartic PI calculation
13 * PI now known to 2 digits of accuracy ...
14 * PI now known to 8 digits of accuracy ...
15 * PI now known to 18 digits of accuracy ...
16 * PI now known to 38 digits of accuracy ...
17 * PI now known to 82 digits of accuracy ...
18 * PI now known to 168 digits of accuracy ...
19 * Borwein-Quartic PI calculation complete
20 *
21 * Borwein PI = [3.1415926535897932384626433832795028841971693993751
22 * 058209749445923078164062862089986280348253421170679
23 * 821480865132823066470938446095505822317253594081284
24 * 8111745028410270193852110555964462294895493038196E+0]
25 */
26
27
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31 #include "m_apm.h"
32
33 #define FALSE 0
34 #define TRUE 1
35
36 extern MAPM compute_PI(int);
37
38
39 int main(int argc, char **argv)
40 {
41 int ii, dplaces;
42 char *buffer;
43 MAPM pi_mapm; /* declare the MAPM variable ... */
44
45
46 if (argc < 2)
47 {
48 fprintf(stdout,"\nUsage: pi_2 digits\t\t\t\t[Version 1.1]\n");
49 fprintf(stdout,
50 " Compute PI to \'digits\' decimal places of accuracy\n");
51
52 exit(4);
53 }
54
55 dplaces = atoi(argv[1]);
56
57 if (dplaces < 2)
58 {
59 fprintf(stdout,"\nUsage: pi-2 digits \n");
60 exit(6);
61 }
62
63
64 fprintf(stderr, "\nstart Borwein-Quartic PI calculation \n");
65
66 pi_mapm = compute_PI(dplaces);
67
68 fprintf(stderr, "Borwein-Quartic PI calculation complete \n\n");
69
70 /*
71 * find out how many digits so we can
72 * convert the entire value into a string.
73 * (plus some pad for decimal point, exponent, etc)
74 */
75
76 ii = pi_mapm.significant_digits() + 12;
77
78 if ((buffer = (char *)malloc(ii)) == NULL)
79 {
80 fprintf(stdout,"PI demo: Out of memory \n");
81 exit(8);
82 }
83
84 /* now convert the MAPM number into a string */
85
86 pi_mapm.toString(buffer, dplaces);
87 printf("Borwein PI = [%s] \n",buffer);
88
89 free(buffer);
90
91 exit(0);
92 }
93
94 /****************************************************************************/
95
96
97 /*
98 * Calculate PI using the Borwein's Quartically Convergent Algorithm
99 *
100 * Init : U0 = sqrt(2)
101 * B0 = 0
102 * P0 = 2 + sqrt(2)
103 *
104 * Iterate:
105 *
106 *
107 * U = 0.5 * [ sqrt(U ) + 1 / sqrt(U ) ]
108 * k+1 k k
109 *
110 *
111 * sqrt(U ) * [ 1 + B ]
112 * k k
113 * B = -----------------------
114 * k+1 U + B
115 * k k
116 *
117 *
118 * P * B * [ 1 + U ]
119 * k k+1 k+1
120 * P = ----------------------------
121 * k+1 1 + B
122 * k+1
123 *
124 *
125 * P ---> PI
126 * k+1
127 *
128 */
129 MAPM compute_PI(int places)
130 {
131 MAPM t1, c0_5, u0, u1, b0, b1, p0, p1;
132 int ct, nn, dplaces, dflag;
133
134
135 dplaces = places + 16; // compute a few extra digits
136 // in the intermediate math
137
138 m_apm_cpp_precision(dplaces); // tell C++ about our precision
139
140 dflag = FALSE;
141 ct = 0;
142
143 c0_5 = "0.5";
144
145 u0 = sqrt("2"); // make sure we use the MAPM sqrt
146 b0 = 0;
147 p0 = 2 + u0;
148
149 while (TRUE)
150 {
151 /*
152 * U = 0.5 * [ sqrt(U ) + 1 / sqrt(U ) ]
153 * k+1 k k
154 */
155
156 t1 = sqrt(u0);
157 u1 = c0_5 * (t1 + 1 / t1);
158
159
160 /*
161 * sqrt(U ) * [ 1 + B ]
162 * k k
163 * B = -----------------------
164 * k+1 U + B
165 * k k
166 *
167 */
168
169 b1 = (t1 * (1 + b0)) / (u0 + b0);
170
171
172 /*
173 * P * B * [ 1 + U ]
174 * k k+1 k+1
175 * P = ----------------------------
176 * k+1 1 + B
177 * k+1
178 */
179
180 p1 = (p0 * b1 * (1 + u1)) / (1 + b1);
181
182
183 if (dflag)
184 break;
185
186 /*
187 * compute the difference from this
188 * iteration to the last one.
189 */
190
191 t1 = p1 - p0;
192
193 if (++ct >= 4)
194 {
195 /*
196 * if diff == 0, we're done.
197 */
198
199 if (t1 == 0)
200 break;
201 }
202
203 /*
204 * if the exponent of the error term (small error like 2.47...E-65)
205 * is small enough, break out after the *next* p1 is calculated.
206 *
207 * get the exponent of the error term, which will be negative.
208 */
209
210 nn = -(t1.exponent());
211
212 /*
213 * normally, this wouldn't be here. it's nice in the demo though.
214 */
215
216 fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(2 * nn));
217
218 if ((4 * nn) >= dplaces)
219 dflag = TRUE;
220
221
222 /* set up for the next iteration */
223 /* */
224 /* also, keep the number of significant */
225 /* digits from growing without bound */
226
227 b0 = b1.round(dplaces);
228 u0 = u1.round(dplaces);
229 p0 = p1.round(dplaces);
230 }
231
232 /* round to the accuracy asked for */
233
234 return (p1.round(places));
235 }
236
237 /****************************************************************************/
238