"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/PI_DEMO/pi_1.c" (21 Feb 2010, 5470 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_1.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * PI_1.C
4 *
5 * DEMO 'C' PROGRAM TO COMPUTE PI USING
6 * BORWEIN'S QUARTICALLY CONVERGENT ALGORITHM
7 *
8 *
9 * Output from 'pi-1 200' : (I line wrapped for this page).
10 *
11 * start Borwein-Quartic PI calculation
12 * PI now known to 2 digits of accuracy ...
13 * PI now known to 8 digits of accuracy ...
14 * PI now known to 18 digits of accuracy ...
15 * PI now known to 38 digits of accuracy ...
16 * PI now known to 82 digits of accuracy ...
17 * PI now known to 168 digits of accuracy ...
18 * Borwein-Quartic PI calculation complete
19 *
20 * Borwein PI = [3.1415926535897932384626433832795028841971693993751
21 * 058209749445923078164062862089986280348253421170679
22 * 821480865132823066470938446095505822317253594081284
23 * 8111745028410270193852110555964462294895493038196E+0]
24 */
25
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <string.h>
30 #include "m_apm.h"
31
32 #define FALSE 0
33 #define TRUE 1
34
35 extern void compute_PI(M_APM, int);
36
37
38 int main(argc, argv)
39 int argc; char *argv[];
40 {
41 int ii, dplaces;
42 char *buffer;
43 M_APM pi_mapm; /* declare the M_APM variable ... */
44
45
46 if (argc < 2)
47 {
48 fprintf(stdout,"\nUsage: pi_1 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 pi_mapm = m_apm_init(); /* now initialize the M_APM variable ... */
56
57 dplaces = atoi(argv[1]);
58
59 if (dplaces < 2)
60 {
61 fprintf(stdout,"\nUsage: pi-1 digits \n");
62 exit(6);
63 }
64
65 fprintf(stderr, "\nstart Borwein-Quartic PI calculation \n");
66
67 compute_PI(pi_mapm, dplaces);
68
69 fprintf(stderr, "Borwein-Quartic PI calculation complete \n\n");
70
71 /*
72 * find out how many digits so we can
73 * convert the entire value into a string.
74 * (plus some pad for decimal point, exponent, etc)
75 */
76
77 ii = m_apm_significant_digits(pi_mapm) + 12;
78
79 if ((buffer = (char *)malloc(ii)) == NULL)
80 {
81 fprintf(stdout,"PI demo: Out of memory \n");
82 exit(8);
83 }
84
85
86 /* now convert the MAPM number into a string */
87
88 m_apm_to_string(buffer, dplaces, pi_mapm);
89 printf("Borwein PI = [%s] \n",buffer);
90
91 m_apm_free(pi_mapm);
92 free(buffer);
93
94 m_apm_free_all_mem();
95
96 exit(0);
97 }
98
99 /****************************************************************************/
100
101 /*
102 * Calculate PI using the Borwein's Quartically Convergent Algorithm
103 *
104 * Init : U0 = sqrt(2)
105 * B0 = 0
106 * P0 = 2 + sqrt(2)
107 *
108 * Iterate:
109 *
110 *
111 * U = 0.5 * [ sqrt(U ) + 1 / sqrt(U ) ]
112 * k+1 k k
113 *
114 *
115 * sqrt(U ) * [ 1 + B ]
116 * k k
117 * B = -----------------------
118 * k+1 U + B
119 * k k
120 *
121 *
122 * P * B * [ 1 + U ]
123 * k k+1 k+1
124 * P = ----------------------------
125 * k+1 1 + B
126 * k+1
127 *
128 *
129 * P ---> PI
130 * k+1
131 *
132 */
133 void compute_PI(outv, places)
134 M_APM outv;
135 int places;
136 {
137 M_APM tmp1, tmp2, tmp3, c0_5, u0, u1, b0, b1, p0, p1;
138 int ct, nn, dplaces, dflag;
139
140 tmp1 = m_apm_init();
141 tmp2 = m_apm_init();
142 tmp3 = m_apm_init();
143 c0_5 = m_apm_init();
144 u0 = m_apm_init();
145 u1 = m_apm_init();
146 b0 = m_apm_init();
147 b1 = m_apm_init();
148 p0 = m_apm_init();
149 p1 = m_apm_init();
150
151 dplaces = places + 16; /* compute a few extra digits */
152 /* in the intermediate math */
153 dflag = FALSE;
154 ct = 0;
155
156 m_apm_set_string(c0_5, "0.5"); /* need this constant */
157
158 m_apm_sqrt(u0, dplaces, MM_Two);
159 m_apm_add(p0, u0, MM_Two);
160
161 /* b0 initialized to 0 from m_apm_init() */
162
163 while (TRUE)
164 {
165 m_apm_sqrt(tmp1, dplaces, u0);
166 m_apm_reciprocal(tmp2, dplaces, tmp1);
167 m_apm_add(tmp3, tmp1, tmp2);
168 m_apm_multiply(u1, c0_5, tmp3);
169
170 m_apm_add(tmp2, MM_One, b0);
171 m_apm_multiply(tmp3, tmp1, tmp2);
172 m_apm_add(tmp2, b0, u0);
173 m_apm_divide(b1, dplaces, tmp3, tmp2);
174
175 m_apm_add(tmp3, MM_One, u1);
176 m_apm_multiply(tmp2, b1, tmp3);
177 m_apm_multiply(tmp1, p0, tmp2);
178 m_apm_add(tmp3, MM_One, b1);
179 m_apm_divide(p1, dplaces, tmp1, tmp3);
180
181 if (dflag)
182 break;
183
184 /*
185 * compute the difference from this
186 * iteration to the last one.
187 */
188
189 m_apm_subtract(tmp2, p0, p1);
190
191 if (++ct >= 4)
192 {
193 /*
194 * if diff == 0, we're done.
195 */
196
197 if (m_apm_sign(tmp2) == 0)
198 break;
199 }
200
201 /*
202 * if the exponent of the error term (small error like 2.47...E-65)
203 * is small enough, break out after the *next* p1 is calculated.
204 *
205 * get the exponent of the error term, which will be negative.
206 */
207
208 nn = -m_apm_exponent(tmp2);
209
210 /*
211 * normally, this wouldn't be here. it's nice in the demo though.
212 */
213
214 fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(2 * nn));
215
216 if ((4 * nn) >= dplaces)
217 dflag = TRUE;
218
219 /* set up for the next iteration */
220
221 m_apm_copy(b0, b1);
222 m_apm_copy(u0, u1);
223 m_apm_copy(p0, p1);
224 }
225
226 /* round to the accuracy asked for */
227
228 m_apm_round(outv, places, p1);
229
230 m_apm_free(p1);
231 m_apm_free(p0);
232 m_apm_free(b1);
233 m_apm_free(b0);
234 m_apm_free(u1);
235 m_apm_free(u0);
236 m_apm_free(c0_5);
237 m_apm_free(tmp3);
238 m_apm_free(tmp2);
239 m_apm_free(tmp1);
240 }
241
242 /****************************************************************************/
243