"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/PI_DEMO/pi_3.c" (21 Feb 2010, 6281 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_3.c" see the
Fossies "Dox" file reference documentation.
1 /****************************************************************************/
2
3 /*
4 * PI_3.C
5 *
6 * DEMO 'C' PROGRAM TO COMPUTE PI USING
7 * J.BORWEIN-P.BORWEIN'S ALGORITHM (1995)
8 *
9 */
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include "m_apm.h"
15
16 #define FALSE 0
17 #define TRUE 1
18
19 #define Mult m_apm_multiply
20 #define Add m_apm_add
21 #define Subtract m_apm_subtract
22 #define Divide m_apm_divide
23 #define AssignString m_apm_set_string
24
25 extern void compute_PI_2(M_APM, int);
26
27
28 void m_apm_square(M_APM x, M_APM y) /* x = y^2 */
29 {
30 m_apm_multiply(x, y, y);
31 }
32
33
34 /*
35 * NOTE : The following 2 functions modify the
36 * input value ('y'). Use with caution
37 * in another application!
38 */
39
40 /* x = y^4 */
41
42 void m_apm_pow4(M_APM x, M_APM y)
43 {
44 m_apm_square(x, y);
45 m_apm_copy(y, x);
46 m_apm_square(x, y);
47 }
48
49 /* x = sqrt(sqrt(y)) */
50
51 void m_apm_root4(M_APM x, int dec_places, M_APM y)
52 {
53 m_apm_sqrt(x, dec_places, y);
54 m_apm_copy(y, x);
55 m_apm_sqrt(x, dec_places, y);
56 }
57
58
59
60 int main(argc, argv)
61 int argc; char *argv[];
62 {
63 int ii, dplaces;
64 char *buffer;
65 M_APM pi_mapm; /* declare the M_APM variable ... */
66
67 if (argc < 2)
68 {
69 fprintf(stdout,"\nUsage: pi_3 digits\t\t\t\t[Version 1.1]\n");
70 fprintf(stdout,
71 " Compute PI to \'digits\' decimal places of accuracy\n");
72
73 exit(4);
74 }
75
76 pi_mapm = m_apm_init(); /* now initialize the M_APM variable ... */
77
78 dplaces = atoi(argv[1]);
79
80 if (dplaces < 2)
81 {
82 fprintf(stdout,"\nUsage: pi digits \n");
83 exit(6);
84 }
85
86 fprintf(stderr, "\nstart Borwein-Borwein PI calculation \n");
87
88 compute_PI_2(pi_mapm, dplaces);
89
90 fprintf(stderr, "Borwein-Borwein PI calculation complete \n\n");
91
92 /*
93 * find out how many digits so we can
94 * convert the entire value into a string.
95 * (plus some pad for decimal point, exponent, etc)
96 */
97
98 ii = m_apm_significant_digits(pi_mapm) + 12;
99
100 if ((buffer = (char *)malloc(ii)) == NULL)
101 {
102 fprintf(stdout,"PI demo: Out of memory \n");
103 exit(8);
104 }
105
106
107 /* now convert the MAPM number into a string */
108
109 m_apm_to_string(buffer, dplaces, pi_mapm);
110 printf("Borwein PI = [%s] \n",buffer);
111
112 m_apm_free(pi_mapm);
113
114 free(buffer);
115
116 m_apm_free_all_mem();
117
118 exit(0);
119 }
120
121 /****************************************************************************/
122
123 /*
124 * Calculate PI using the J. Borwein - P. Borwein's Algorithm (1995)
125 *
126 *
127 * Init : A0 = 6 - 4 * sqrt(2)
128 * B0 = sqrt(2) - 1
129 *
130 * Iterate:
131 *
132 * 4
133 * 1 - quart( 1 - B )
134 * k
135 * B = --------------------
136 * k+1 4
137 * 1 + quart( 1 - B )
138 * k
139 *
140 *
141 * 4 2k+3 2
142 * A = A ( 1 + B ) - 2 B ( 1 + B + B )
143 * k+1 k k+1 k+1 k+1 k+1
144 *
145 *
146 *
147 *
148 * where quart(x) = sqrt(sqrt(x)) = 4th root of x
149 *
150 * 1
151 * A ---> -----
152 * k+1 PI
153 *
154 */
155
156 void compute_PI_2(outv, places)
157 M_APM outv;
158 int places;
159 {
160 M_APM tmp1, tmp2, tmp3, tmp4, a0, a1, b0, b1, MM_Six;
161 int kk, ct, nn, dplaces, dflag;
162
163 tmp1 = m_apm_init();
164 tmp2 = m_apm_init();
165 tmp3 = m_apm_init();
166 tmp4 = m_apm_init();
167 a0 = m_apm_init();
168 a1 = m_apm_init();
169 b0 = m_apm_init();
170 b1 = m_apm_init();
171 MM_Six = m_apm_init();
172
173 dplaces = places + 16; /* compute a few extra digits */
174 /* in the intermediate math */
175 dflag = FALSE;
176 ct = 0;
177 kk = 0;
178
179 AssignString(MM_Six, "6");
180
181 m_apm_sqrt(tmp1, dplaces, MM_Two);
182 Mult(tmp2, tmp1, MM_Four);
183 Subtract(a0, MM_Six, tmp2);
184 Subtract(b0, tmp1, MM_One);
185
186 while (TRUE)
187 {
188 m_apm_pow4(tmp1, b0); /* tmp1 = b0^4 */
189 Subtract(tmp2, MM_One, tmp1); /* tmp2 = 1 - b0^4 */
190 m_apm_root4(tmp3, dplaces, tmp2); /* tmp3 = sqrt[4](1 - b0^4) */
191 Subtract(tmp1, MM_One, tmp3); /* tmp1 = 1 - tmp3 */
192 Add(tmp2, MM_One, tmp3); /* tmp2 = 1 + tmp3 */
193 Divide(b1, dplaces, tmp1, tmp2); /* b1 = tmp1 / tmp2 */
194
195 Add(tmp1, MM_One, b1); /* tmp1 = 1 + b1 */
196 m_apm_pow4(tmp2, tmp1); /* tmp2 = (1 + b1)^4 */
197 Mult(tmp3, a0, tmp2); /* tmp3 = a0 * (1 + b1)^4 */
198
199 m_apm_square(tmp1, b1); /* tmp1 = b1^2 */
200 Add(tmp2, tmp1, b1); /* tmp2 = b1 + b1^2 */
201 Add(tmp1, tmp2, MM_One); /* tmp1 = 1 + b1 + b1^2 */
202 Mult(tmp2, tmp1, b1); /* tmp2 = b1(1 + b1 + b1^2) */
203
204 m_apm_integer_pow(tmp1, dplaces,
205 MM_Two, (2 * kk + 3)); /* tmp1 = 2^(2k+3) */
206
207 Mult(tmp4, tmp1, tmp2); /* tmp4 = 2^(2k+3) b1 (1 + b1 + b1^2) */
208
209 Subtract(tmp1, tmp3, tmp4); /* a1 = a0(1 + b1)^4 - tmp4 */
210 m_apm_round(a1, dplaces, tmp1);
211
212 if (dflag)
213 break;
214
215 /*
216 * compute the difference from this
217 * iteration to the last one.
218 */
219
220 m_apm_subtract(tmp2, a0, a1);
221
222 if (++ct >= 4)
223 {
224 /*
225 * if diff == 0, we're done.
226 */
227
228 if (m_apm_sign(tmp2) == 0)
229 break;
230 }
231
232 /*
233 * if the exponent of the error term (small error like 2.47...E-65)
234 * is small enough, break out after the *next* a1 is calculated.
235 *
236 * get the exponent of the error term, which will be negative.
237 */
238
239 nn = -m_apm_exponent(tmp2);
240
241 /*
242 * normally, this wouldn't be here. it's nice in the demo though.
243 */
244
245 fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(4 * nn));
246
247 if ((15 * nn) >= dplaces)
248 dflag = TRUE;
249
250 /* set up for the next iteration */
251
252 m_apm_copy(b0, b1);
253 m_apm_copy(a0, a1);
254 kk++;
255 }
256
257 /* round to the accuracy asked for after taking the reciprocal */
258
259 m_apm_reciprocal(tmp1, dplaces, a1);
260 m_apm_round(outv, places, tmp1);
261
262 m_apm_free(MM_Six);
263 m_apm_free(b1);
264 m_apm_free(b0);
265 m_apm_free(a1);
266 m_apm_free(a0);
267 m_apm_free(tmp4);
268 m_apm_free(tmp3);
269 m_apm_free(tmp2);
270 m_apm_free(tmp1);
271 }
272
273 /****************************************************************************/
274