1 /*
2 Name: pi.c
3 Purpose: Computes digits of the physical constant pi.
4 Author: M. J. Fromberger
5
6 Copyright (C) 2002-2008 Michael J. Fromberger, All Rights Reserved.
7
8 Notes:
9 Uses Machin's formula, which should be suitable for a few thousand digits.
10
11 Permission is hereby granted, free of charge, to any person obtaining a copy
12 of this software and associated documentation files (the "Software"), to deal
13 in the Software without restriction, including without limitation the rights
14 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
15 copies of the Software, and to permit persons to whom the Software is
16 furnished to do so, subject to the following conditions:
17
18 The above copyright notice and this permission notice shall be included in
19 all copies or substantial portions of the Software.
20
21 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
22 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
24 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
26 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
27 SOFTWARE.
28 */
29
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <string.h>
33 #include <time.h>
34
35 #include "imath.h"
36
37 int g_radix = 10; /* use this radix for output */
38
39 mp_result arctan(mp_small radix, mp_small mul, mp_small x, mp_small prec,
40 mp_int sum);
41
42 char g_buf[4096];
43
main(int argc,char * argv[])44 int main(int argc, char *argv[]) {
45 mp_result res;
46 mpz_t sum1, sum2;
47 int ndigits, out = 0;
48 clock_t start, end;
49
50 if (argc < 2) {
51 fprintf(stderr, "Usage: %s <num-digits> [<radix>]\n", argv[0]);
52 return 1;
53 }
54
55 if ((ndigits = abs(atoi(argv[1]))) == 0) {
56 fprintf(stderr, "%s: you must request at least 1 digit\n", argv[0]);
57 return 1;
58 } else if ((mp_word)ndigits > MP_DIGIT_MAX) {
59 fprintf(stderr, "%s: you may request at most %u digits\n", argv[0],
60 (unsigned int)MP_DIGIT_MAX);
61 return 1;
62 }
63
64 if (argc > 2) {
65 int radix = atoi(argv[2]);
66
67 if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX) {
68 fprintf(stderr, "%s: you may only specify a radix between %d and %d\n",
69 argv[0], MP_MIN_RADIX, MP_MAX_RADIX);
70 return 1;
71 }
72 g_radix = radix;
73 }
74
75 mp_int_init(&sum1);
76 mp_int_init(&sum2);
77 start = clock();
78
79 /* sum1 = 16 * arctan(1/5) */
80 if ((res = arctan(g_radix, 16, 5, ndigits, &sum1)) != MP_OK) {
81 fprintf(stderr, "%s: error computing arctan: %d\n", argv[0], res);
82 out = 1;
83 goto CLEANUP;
84 }
85
86 /* sum2 = 4 * arctan(1/239) */
87 if ((res = arctan(g_radix, 4, 239, ndigits, &sum2)) != MP_OK) {
88 fprintf(stderr, "%s: error computing arctan: %d\n", argv[0], res);
89 out = 1;
90 goto CLEANUP;
91 }
92
93 /* pi = sum1 - sum2 */
94 if ((res = mp_int_sub(&sum1, &sum2, &sum1)) != MP_OK) {
95 fprintf(stderr, "%s: error computing pi: %d\n", argv[0], res);
96 out = 1;
97 goto CLEANUP;
98 }
99 end = clock();
100
101 mp_int_to_string(&sum1, g_radix, g_buf, sizeof(g_buf));
102 printf("%c.%s\n", g_buf[0], g_buf + 1);
103
104 fprintf(stderr, "Computation took %.2f sec.\n",
105 (double)(end - start) / CLOCKS_PER_SEC);
106
107 CLEANUP:
108 mp_int_clear(&sum1);
109 mp_int_clear(&sum2);
110
111 return out;
112 }
113
114 /*
115 Compute mul * atan(1/x) to prec digits of precision, and store the
116 result in sum.
117
118 Computes atan(1/x) using the formula:
119
120 1 1 1 1
121 atan(1/x) = --- - ---- + ---- - ---- + ...
122 x 3x^3 5x^5 7x^7
123
124 */
arctan(mp_small radix,mp_small mul,mp_small x,mp_small prec,mp_int sum)125 mp_result arctan(mp_small radix, mp_small mul, mp_small x, mp_small prec,
126 mp_int sum) {
127 mpz_t t, v;
128 mp_result res;
129 mp_small rem, sign = 1, coeff = 1;
130
131 mp_int_init(&t);
132 mp_int_init(&v);
133 ++prec;
134
135 /* Compute mul * radix^prec * x
136 The initial multiplication by x saves a special case in the loop for
137 the first term of the series.
138 */
139 if ((res = mp_int_expt_value(radix, prec, &t)) != MP_OK ||
140 (res = mp_int_mul_value(&t, mul, &t)) != MP_OK ||
141 (res = mp_int_mul_value(&t, x, &t)) != MP_OK)
142 goto CLEANUP;
143
144 x *= x; /* assumes x <= sqrt(MP_SMALL_MAX) */
145 mp_int_zero(sum);
146
147 do {
148 if ((res = mp_int_div_value(&t, x, &t, &rem)) != MP_OK) goto CLEANUP;
149
150 if ((res = mp_int_div_value(&t, coeff, &v, &rem)) != MP_OK) goto CLEANUP;
151
152 /* Add or subtract the result depending on the current sign (1 = add) */
153 if (sign > 0)
154 res = mp_int_add(sum, &v, sum);
155 else
156 res = mp_int_sub(sum, &v, sum);
157
158 if (res != MP_OK) goto CLEANUP;
159 sign = -sign;
160 coeff += 2;
161
162 } while (mp_int_compare_zero(&t) != 0);
163
164 res = mp_int_div_value(sum, radix, sum, NULL);
165
166 CLEANUP:
167 mp_int_clear(&v);
168 mp_int_clear(&t);
169
170 return res;
171 }
172
173 /* Here there be dragons */
174