The Tor Browser / file comparison
comparison: security/nss/lib/freebl/mpi/utils/pi.c
security/nss/lib/freebl/mpi/utils/pi.c
- branch
- TOR_BUG_9701
- changeset 15
- b8a032363ba2
equal
deleted
inserted
replaced
| -1:000000000000 | 0:92ffbc5c3139 |
|---|---|
| 1 /* | |
| 2 * pi.c | |
| 3 * | |
| 4 * Compute pi to an arbitrary number of digits. Uses Machin's formula, | |
| 5 * like everyone else on the planet: | |
| 6 * | |
| 7 * pi = 16 * arctan(1/5) - 4 * arctan(1/239) | |
| 8 * | |
| 9 * This is pretty effective for up to a few thousand digits, but it | |
| 10 * gets pretty slow after that. | |
| 11 * | |
| 12 * This Source Code Form is subject to the terms of the Mozilla Public | |
| 13 * License, v. 2.0. If a copy of the MPL was not distributed with this | |
| 14 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ | |
| 15 | |
| 16 #include <stdio.h> | |
| 17 #include <stdlib.h> | |
| 18 #include <string.h> | |
| 19 #include <limits.h> | |
| 20 #include <time.h> | |
| 21 | |
| 22 #include "mpi.h" | |
| 23 | |
| 24 mp_err arctan(mp_digit mul, mp_digit x, mp_digit prec, mp_int *sum); | |
| 25 | |
| 26 int main(int argc, char *argv[]) | |
| 27 { | |
| 28 mp_err res; | |
| 29 mp_digit ndigits; | |
| 30 mp_int sum1, sum2; | |
| 31 clock_t start, stop; | |
| 32 int out = 0; | |
| 33 | |
| 34 /* Make the user specify precision on the command line */ | |
| 35 if(argc < 2) { | |
| 36 fprintf(stderr, "Usage: %s <num-digits>\n", argv[0]); | |
| 37 return 1; | |
| 38 } | |
| 39 | |
| 40 if((ndigits = abs(atoi(argv[1]))) == 0) { | |
| 41 fprintf(stderr, "%s: you must request at least 1 digit\n", argv[0]); | |
| 42 return 1; | |
| 43 } | |
| 44 | |
| 45 start = clock(); | |
| 46 mp_init(&sum1); mp_init(&sum2); | |
| 47 | |
| 48 /* sum1 = 16 * arctan(1/5) */ | |
| 49 if((res = arctan(16, 5, ndigits, &sum1)) != MP_OKAY) { | |
| 50 fprintf(stderr, "%s: arctan: %s\n", argv[0], mp_strerror(res)); | |
| 51 out = 1; goto CLEANUP; | |
| 52 } | |
| 53 | |
| 54 /* sum2 = 4 * arctan(1/239) */ | |
| 55 if((res = arctan(4, 239, ndigits, &sum2)) != MP_OKAY) { | |
| 56 fprintf(stderr, "%s: arctan: %s\n", argv[0], mp_strerror(res)); | |
| 57 out = 1; goto CLEANUP; | |
| 58 } | |
| 59 | |
| 60 /* pi = sum1 - sum2 */ | |
| 61 if((res = mp_sub(&sum1, &sum2, &sum1)) != MP_OKAY) { | |
| 62 fprintf(stderr, "%s: mp_sub: %s\n", argv[0], mp_strerror(res)); | |
| 63 out = 1; goto CLEANUP; | |
| 64 } | |
| 65 stop = clock(); | |
| 66 | |
| 67 /* Write the output in decimal */ | |
| 68 { | |
| 69 char *buf = malloc(mp_radix_size(&sum1, 10)); | |
| 70 | |
| 71 if(buf == NULL) { | |
| 72 fprintf(stderr, "%s: out of memory\n", argv[0]); | |
| 73 out = 1; goto CLEANUP; | |
| 74 } | |
| 75 mp_todecimal(&sum1, buf); | |
| 76 printf("%s\n", buf); | |
| 77 free(buf); | |
| 78 } | |
| 79 | |
| 80 fprintf(stderr, "Computation took %.2f sec.\n", | |
| 81 (double)(stop - start) / CLOCKS_PER_SEC); | |
| 82 | |
| 83 CLEANUP: | |
| 84 mp_clear(&sum1); | |
| 85 mp_clear(&sum2); | |
| 86 | |
| 87 return out; | |
| 88 | |
| 89 } | |
| 90 | |
| 91 /* Compute sum := mul * arctan(1/x), to 'prec' digits of precision */ | |
| 92 mp_err arctan(mp_digit mul, mp_digit x, mp_digit prec, mp_int *sum) | |
| 93 { | |
| 94 mp_int t, v; | |
| 95 mp_digit q = 1, rd; | |
| 96 mp_err res; | |
| 97 int sign = 1; | |
| 98 | |
| 99 prec += 3; /* push inaccuracies off the end */ | |
| 100 | |
| 101 mp_init(&t); mp_set(&t, 10); | |
| 102 mp_init(&v); | |
| 103 if((res = mp_expt_d(&t, prec, &t)) != MP_OKAY || /* get 10^prec */ | |
| 104 (res = mp_mul_d(&t, mul, &t)) != MP_OKAY || /* ... times mul */ | |
| 105 (res = mp_mul_d(&t, x, &t)) != MP_OKAY) /* ... times x */ | |
| 106 goto CLEANUP; | |
| 107 | |
| 108 /* | |
| 109 The extra multiplication by x in the above takes care of what | |
| 110 would otherwise have to be a special case for 1 / x^1 during the | |
| 111 first loop iteration. A little sneaky, but effective. | |
| 112 | |
| 113 We compute arctan(1/x) by the formula: | |
| 114 | |
| 115 1 1 1 1 | |
| 116 - - ----- + ----- - ----- + ... | |
| 117 x 3 x^3 5 x^5 7 x^7 | |
| 118 | |
| 119 We multiply through by 'mul' beforehand, which gives us a couple | |
| 120 more iterations and more precision | |
| 121 */ | |
| 122 | |
| 123 x *= x; /* works as long as x < sqrt(RADIX), which it is here */ | |
| 124 | |
| 125 mp_zero(sum); | |
| 126 | |
| 127 do { | |
| 128 if((res = mp_div_d(&t, x, &t, &rd)) != MP_OKAY) | |
| 129 goto CLEANUP; | |
| 130 | |
| 131 if(sign < 0 && rd != 0) | |
| 132 mp_add_d(&t, 1, &t); | |
| 133 | |
| 134 if((res = mp_div_d(&t, q, &v, &rd)) != MP_OKAY) | |
| 135 goto CLEANUP; | |
| 136 | |
| 137 if(sign < 0 && rd != 0) | |
| 138 mp_add_d(&v, 1, &v); | |
| 139 | |
| 140 if(sign > 0) | |
| 141 res = mp_add(sum, &v, sum); | |
| 142 else | |
| 143 res = mp_sub(sum, &v, sum); | |
| 144 | |
| 145 if(res != MP_OKAY) | |
| 146 goto CLEANUP; | |
| 147 | |
| 148 sign *= -1; | |
| 149 q += 2; | |
| 150 | |
| 151 } while(mp_cmp_z(&t) != 0); | |
| 152 | |
| 153 /* Chop off inaccurate low-order digits */ | |
| 154 mp_div_d(sum, 1000, sum, NULL); | |
| 155 | |
| 156 CLEANUP: | |
| 157 mp_clear(&v); | |
| 158 mp_clear(&t); | |
| 159 | |
| 160 return res; | |
| 161 } | |
| 162 | |
| 163 /*------------------------------------------------------------------------*/ | |
| 164 /* HERE THERE BE DRAGONS */ |
