174 lines · c
1/*2 Name: pi.c3 Purpose: Computes digits of the physical constant pi.4 Author: M. J. Fromberger5 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 copy12 of this software and associated documentation files (the "Software"), to deal13 in the Software without restriction, including without limitation the rights14 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell15 copies of the Software, and to permit persons to whom the Software is16 furnished to do so, subject to the following conditions:17 18 The above copyright notice and this permission notice shall be included in19 all copies or substantial portions of the Software.20 21 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR22 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,23 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE24 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER25 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 THE27 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 37int g_radix = 10; /* use this radix for output */38 39mp_result arctan(mp_small radix, mp_small mul, mp_small x, mp_small prec,40 mp_int sum);41 42char g_buf[4096];43 44int 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 107CLEANUP: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 the116 result in sum.117 118 Computes atan(1/x) using the formula:119 120 1 1 1 1121 atan(1/x) = --- - ---- + ---- - ---- + ...122 x 3x^3 5x^5 7x^7123 124 */125mp_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 * x136 The initial multiplication by x saves a special case in the loop for137 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 else156 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 166CLEANUP:167 mp_int_clear(&v);168 mp_int_clear(&t);169 170 return res;171}172 173/* Here there be dragons */174