1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/security/nss/lib/freebl/mpi/utils/pi.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,164 @@
1.4 +/*
1.5 + * pi.c
1.6 + *
1.7 + * Compute pi to an arbitrary number of digits. Uses Machin's formula,
1.8 + * like everyone else on the planet:
1.9 + *
1.10 + * pi = 16 * arctan(1/5) - 4 * arctan(1/239)
1.11 + *
1.12 + * This is pretty effective for up to a few thousand digits, but it
1.13 + * gets pretty slow after that.
1.14 + *
1.15 + * This Source Code Form is subject to the terms of the Mozilla Public
1.16 + * License, v. 2.0. If a copy of the MPL was not distributed with this
1.17 + * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
1.18 +
1.19 +#include <stdio.h>
1.20 +#include <stdlib.h>
1.21 +#include <string.h>
1.22 +#include <limits.h>
1.23 +#include <time.h>
1.24 +
1.25 +#include "mpi.h"
1.26 +
1.27 +mp_err arctan(mp_digit mul, mp_digit x, mp_digit prec, mp_int *sum);
1.28 +
1.29 +int main(int argc, char *argv[])
1.30 +{
1.31 + mp_err res;
1.32 + mp_digit ndigits;
1.33 + mp_int sum1, sum2;
1.34 + clock_t start, stop;
1.35 + int out = 0;
1.36 +
1.37 + /* Make the user specify precision on the command line */
1.38 + if(argc < 2) {
1.39 + fprintf(stderr, "Usage: %s <num-digits>\n", argv[0]);
1.40 + return 1;
1.41 + }
1.42 +
1.43 + if((ndigits = abs(atoi(argv[1]))) == 0) {
1.44 + fprintf(stderr, "%s: you must request at least 1 digit\n", argv[0]);
1.45 + return 1;
1.46 + }
1.47 +
1.48 + start = clock();
1.49 + mp_init(&sum1); mp_init(&sum2);
1.50 +
1.51 + /* sum1 = 16 * arctan(1/5) */
1.52 + if((res = arctan(16, 5, ndigits, &sum1)) != MP_OKAY) {
1.53 + fprintf(stderr, "%s: arctan: %s\n", argv[0], mp_strerror(res));
1.54 + out = 1; goto CLEANUP;
1.55 + }
1.56 +
1.57 + /* sum2 = 4 * arctan(1/239) */
1.58 + if((res = arctan(4, 239, ndigits, &sum2)) != MP_OKAY) {
1.59 + fprintf(stderr, "%s: arctan: %s\n", argv[0], mp_strerror(res));
1.60 + out = 1; goto CLEANUP;
1.61 + }
1.62 +
1.63 + /* pi = sum1 - sum2 */
1.64 + if((res = mp_sub(&sum1, &sum2, &sum1)) != MP_OKAY) {
1.65 + fprintf(stderr, "%s: mp_sub: %s\n", argv[0], mp_strerror(res));
1.66 + out = 1; goto CLEANUP;
1.67 + }
1.68 + stop = clock();
1.69 +
1.70 + /* Write the output in decimal */
1.71 + {
1.72 + char *buf = malloc(mp_radix_size(&sum1, 10));
1.73 +
1.74 + if(buf == NULL) {
1.75 + fprintf(stderr, "%s: out of memory\n", argv[0]);
1.76 + out = 1; goto CLEANUP;
1.77 + }
1.78 + mp_todecimal(&sum1, buf);
1.79 + printf("%s\n", buf);
1.80 + free(buf);
1.81 + }
1.82 +
1.83 + fprintf(stderr, "Computation took %.2f sec.\n",
1.84 + (double)(stop - start) / CLOCKS_PER_SEC);
1.85 +
1.86 + CLEANUP:
1.87 + mp_clear(&sum1);
1.88 + mp_clear(&sum2);
1.89 +
1.90 + return out;
1.91 +
1.92 +}
1.93 +
1.94 +/* Compute sum := mul * arctan(1/x), to 'prec' digits of precision */
1.95 +mp_err arctan(mp_digit mul, mp_digit x, mp_digit prec, mp_int *sum)
1.96 +{
1.97 + mp_int t, v;
1.98 + mp_digit q = 1, rd;
1.99 + mp_err res;
1.100 + int sign = 1;
1.101 +
1.102 + prec += 3; /* push inaccuracies off the end */
1.103 +
1.104 + mp_init(&t); mp_set(&t, 10);
1.105 + mp_init(&v);
1.106 + if((res = mp_expt_d(&t, prec, &t)) != MP_OKAY || /* get 10^prec */
1.107 + (res = mp_mul_d(&t, mul, &t)) != MP_OKAY || /* ... times mul */
1.108 + (res = mp_mul_d(&t, x, &t)) != MP_OKAY) /* ... times x */
1.109 + goto CLEANUP;
1.110 +
1.111 + /*
1.112 + The extra multiplication by x in the above takes care of what
1.113 + would otherwise have to be a special case for 1 / x^1 during the
1.114 + first loop iteration. A little sneaky, but effective.
1.115 +
1.116 + We compute arctan(1/x) by the formula:
1.117 +
1.118 + 1 1 1 1
1.119 + - - ----- + ----- - ----- + ...
1.120 + x 3 x^3 5 x^5 7 x^7
1.121 +
1.122 + We multiply through by 'mul' beforehand, which gives us a couple
1.123 + more iterations and more precision
1.124 + */
1.125 +
1.126 + x *= x; /* works as long as x < sqrt(RADIX), which it is here */
1.127 +
1.128 + mp_zero(sum);
1.129 +
1.130 + do {
1.131 + if((res = mp_div_d(&t, x, &t, &rd)) != MP_OKAY)
1.132 + goto CLEANUP;
1.133 +
1.134 + if(sign < 0 && rd != 0)
1.135 + mp_add_d(&t, 1, &t);
1.136 +
1.137 + if((res = mp_div_d(&t, q, &v, &rd)) != MP_OKAY)
1.138 + goto CLEANUP;
1.139 +
1.140 + if(sign < 0 && rd != 0)
1.141 + mp_add_d(&v, 1, &v);
1.142 +
1.143 + if(sign > 0)
1.144 + res = mp_add(sum, &v, sum);
1.145 + else
1.146 + res = mp_sub(sum, &v, sum);
1.147 +
1.148 + if(res != MP_OKAY)
1.149 + goto CLEANUP;
1.150 +
1.151 + sign *= -1;
1.152 + q += 2;
1.153 +
1.154 + } while(mp_cmp_z(&t) != 0);
1.155 +
1.156 + /* Chop off inaccurate low-order digits */
1.157 + mp_div_d(sum, 1000, sum, NULL);
1.158 +
1.159 + CLEANUP:
1.160 + mp_clear(&v);
1.161 + mp_clear(&t);
1.162 +
1.163 + return res;
1.164 +}
1.165 +
1.166 +/*------------------------------------------------------------------------*/
1.167 +/* HERE THERE BE DRAGONS */
