diff options
Diffstat (limited to 'plugins/supereq/nsfft-1.00/ooura/pi_fft.c')
-rw-r--r-- | plugins/supereq/nsfft-1.00/ooura/pi_fft.c | 1616 |
1 files changed, 1616 insertions, 0 deletions
diff --git a/plugins/supereq/nsfft-1.00/ooura/pi_fft.c b/plugins/supereq/nsfft-1.00/ooura/pi_fft.c new file mode 100644 index 00000000..c9a76bf8 --- /dev/null +++ b/plugins/supereq/nsfft-1.00/ooura/pi_fft.c @@ -0,0 +1,1616 @@ +/* +---- calculation of PI(= 3.14159...) using FFT ---- + by T.Ooura, ver. LG1.1.2-MP1.5a Sep. 2001. + +This is a test program to estimate the performance of +the FFT routines: fft*g.c. + +Example compilation: + GNU : gcc -O6 -ffast-math pi_fft.c fftsg.c -lm -o pi_fftsg + SUN : cc -fast -xO5 pi_fft.c fft8g.c -lm -o pi_fft8g + Microsoft: cl /O2 /G6 pi_fft.c fft4g.c /Fepi_fft4g.exe + ... + etc. +*/ + +/* Please check the following macros before compiling */ +#ifndef DBL_ERROR_MARGIN +#define DBL_ERROR_MARGIN 0.3 /* must be < 0.5 */ +#endif + + +#include <math.h> +#include <limits.h> +#include <float.h> +#include <stdio.h> +#include <stdlib.h> +#include <time.h> + + +void mp_load_0(int n, int radix, int out[]); +void mp_load_1(int n, int radix, int out[]); +void mp_copy(int n, int radix, int in[], int out[]); +void mp_round(int n, int radix, int m, int inout[]); +int mp_cmp(int n, int radix, int in1[], int in2[]); +void mp_add(int n, int radix, int in1[], int in2[], int out[]); +void mp_sub(int n, int radix, int in1[], int in2[], int out[]); +void mp_imul(int n, int radix, int in1[], int in2, int out[]); +int mp_idiv(int n, int radix, int in1[], int in2, int out[]); +void mp_idiv_2(int n, int radix, int in[], int out[]); +double mp_mul_radix_test(int n, int radix, int nfft, + double tmpfft[], int ip[], double w[]); +void mp_mul(int n, int radix, int in1[], int in2[], int out[], + int tmp[], int nfft, double tmp1fft[], double tmp2fft[], + double tmp3fft[], int ip[], double w[]); +void mp_squ(int n, int radix, int in[], int out[], int tmp[], + int nfft, double tmp1fft[], double tmp2fft[], + int ip[], double w[]); +void mp_mulh(int n, int radix, int in1[], int in2[], int out[], + int nfft, double in1fft[], double outfft[], + int ip[], double w[]); +void mp_squh(int n, int radix, int in[], int out[], + int nfft, double inoutfft[], int ip[], double w[]); +int mp_inv(int n, int radix, int in[], int out[], + int tmp1[], int tmp2[], int nfft, + double tmp1fft[], double tmp2fft[], int ip[], double w[]); +int mp_sqrt(int n, int radix, int in[], int out[], + int tmp1[], int tmp2[], int nfft, + double tmp1fft[], double tmp2fft[], int ip[], double w[]); +void mp_sprintf(int n, int log10_radix, int in[], char out[]); +void mp_sscanf(int n, int log10_radix, char in[], int out[]); +void mp_fprintf(int n, int log10_radix, int in[], FILE *fout); + + +int main() +{ + int nfft, log2_nfft, radix, log10_radix, n, npow, nprc; + double err, d_time, n_op; + int *a, *b, *c, *e, *i1, *i2, *ip; + double *d1, *d2, *d3, *w; + time_t t_1, t_2; + FILE *f_log, *f_out; + + f_log = fopen("pi.log", "w"); + printf("PI calculation to estimate the FFT benchmarks\n"); + fprintf(f_log, "PI calculation to estimate the FFT benchmarks\n"); + printf("length of FFT =?\n"); + scanf("%d", &nfft); + + printf("initializing...\n"); + for (log2_nfft = 1; (1 << log2_nfft) < nfft; log2_nfft++); + nfft = 1 << log2_nfft; + n = nfft + 2; + ip = (int *) malloc((3 + (int) sqrt(0.5 * nfft)) * sizeof(int)); + w = (double *) malloc(nfft / 2 * sizeof(double)); + a = (int *) malloc((n + 2) * sizeof(int)); + b = (int *) malloc((n + 2) * sizeof(int)); + c = (int *) malloc((n + 2) * sizeof(int)); + e = (int *) malloc((n + 2) * sizeof(int)); + i1 = (int *) malloc((n + 2) * sizeof(int)); + i2 = (int *) malloc((n + 2) * sizeof(int)); + d1 = (double *) malloc((nfft + 2) * sizeof(double)); + d2 = (double *) malloc((nfft + 2) * sizeof(double)); + d3 = (double *) malloc((nfft + 2) * sizeof(double)); + if (d3 == NULL) { + printf("Allocation Failure!\n"); + exit(1); + } + ip[0] = 0; + /* ---- radix test ---- */ + log10_radix = 1; + radix = 10; + err = mp_mul_radix_test(n, radix, nfft, d1, ip, w); + err += DBL_EPSILON * (n * radix * radix / 4); + while (100 * err < DBL_ERROR_MARGIN && radix <= INT_MAX / 20) { + err *= 100; + log10_radix++; + radix *= 10; + } + printf("nfft= %d\nradix= %d\nerror_margin= %g\n", nfft, radix, err); + fprintf(f_log, "nfft= %d\nradix= %d\nerror_margin= %g\n", nfft, radix, err); + printf("calculating %d digits of PI...\n", log10_radix * (n - 2)); + fprintf(f_log, "calculating %d digits of PI...\n", log10_radix * (n - 2)); + /* ---- time check ---- */ + time(&t_1); + /* + * ---- a formula based on the AGM (Arithmetic-Geometric Mean) ---- + * c = sqrt(0.125); + * a = 1 + 3 * c; + * b = sqrt(a); + * e = b - 0.625; + * b = 2 * b; + * c = e - c; + * a = a + e; + * npow = 4; + * do { + * npow = 2 * npow; + * e = (a + b) / 2; + * b = sqrt(a * b); + * e = e - b; + * b = 2 * b; + * c = c - e; + * a = e + b; + * } while (e > SQRT_SQRT_EPSILON); + * e = e * e / 4; + * a = a + b; + * pi = (a * a - e - e / 2) / (a * c - e) / npow; + * ---- modification ---- + * This is a modified version of Gauss-Legendre formula + * (by T.Ooura). It is faster than original version. + * ---- reference ---- + * 1. E.Salamin, + * Computation of PI Using Arithmetic-Geometric Mean, + * Mathematics of Computation, Vol.30 1976. + * 2. R.P.Brent, + * Fast Multiple-Precision Evaluation of Elementary Functions, + * J. ACM 23 1976. + * 3. D.Takahasi, Y.Kanada, + * Calculation of PI to 51.5 Billion Decimal Digits on + * Distributed Memoriy Parallel Processors, + * Transactions of Information Processing Society of Japan, + * Vol.39 No.7 1998. + * 4. T.Ooura, + * Improvement of the PI Calculation Algorithm and + * Implementation of Fast Multiple-Precision Computation, + * Information Processing Society of Japan SIG Notes, + * 98-HPC-74, 1998. + */ + /* ---- c = sqrt(0.125) ---- */ + mp_sscanf(n, log10_radix, "0.125", a); + mp_sqrt(n, radix, a, c, i1, i2, nfft, d1, d2, ip, w); + /* ---- a = 1 + 3 * c ---- */ + mp_imul(n, radix, c, 3, e); + mp_sscanf(n, log10_radix, "1", a); + mp_add(n, radix, a, e, a); + /* ---- b = sqrt(a) ---- */ + mp_sqrt(n, radix, a, b, i1, i2, nfft, d1, d2, ip, w); + /* ---- e = b - 0.625 ---- */ + mp_sscanf(n, log10_radix, "0.625", e); + mp_sub(n, radix, b, e, e); + /* ---- b = 2 * b ---- */ + mp_add(n, radix, b, b, b); + /* ---- c = e - c ---- */ + mp_sub(n, radix, e, c, c); + /* ---- a = a + e ---- */ + mp_add(n, radix, a, e, a); + printf("AGM iteration\n"); + fprintf(f_log, "AGM iteration\n"); + npow = 4; + do { + npow *= 2; + /* ---- e = (a + b) / 2 ---- */ + mp_add(n, radix, a, b, e); + mp_idiv_2(n, radix, e, e); + /* ---- b = sqrt(a * b) ---- */ + mp_mul(n, radix, a, b, a, i1, nfft, d1, d2, d3, ip, w); + mp_sqrt(n, radix, a, b, i1, i2, nfft, d1, d2, ip, w); + /* ---- e = e - b ---- */ + mp_sub(n, radix, e, b, e); + /* ---- b = 2 * b ---- */ + mp_add(n, radix, b, b, b); + /* ---- c = c - e ---- */ + mp_sub(n, radix, c, e, c); + /* ---- a = e + b ---- */ + mp_add(n, radix, e, b, a); + /* ---- convergence check ---- */ + nprc = -e[1]; + if (e[0] == 0) { + nprc = n; + } + printf("precision= %d\n", 4 * nprc * log10_radix); + fprintf(f_log, "precision= %d\n", 4 * nprc * log10_radix); + } while (4 * nprc <= n); + /* ---- e = e * e / 4 (half precision) ---- */ + mp_idiv_2(n, radix, e, e); + mp_squh(n, radix, e, e, nfft, d1, ip, w); + /* ---- a = a + b ---- */ + mp_add(n, radix, a, b, a); + /* ---- a = (a * a - e - e / 2) / (a * c - e) / npow ---- */ + mp_mul(n, radix, a, c, c, i1, nfft, d1, d2, d3, ip, w); + mp_sub(n, radix, c, e, c); + mp_inv(n, radix, c, b, i1, i2, nfft, d1, d2, ip, w); + mp_squ(n, radix, a, a, i1, nfft, d1, d2, ip, w); + mp_sub(n, radix, a, e, a); + mp_idiv_2(n, radix, e, e); + mp_sub(n, radix, a, e, a); + mp_mul(n, radix, a, b, a, i1, nfft, d1, d2, d3, ip, w); + mp_idiv(n, radix, a, npow, a); + /* ---- time check ---- */ + time(&t_2); + /* ---- output ---- */ + f_out = fopen("pi.dat", "w"); + printf("writing pi.dat...\n"); + mp_fprintf(n - 1, log10_radix, a, f_out); + fclose(f_out); + free(d3); + free(d2); + free(d1); + free(i2); + free(i1); + free(e); + free(c); + free(b); + free(a); + free(w); + free(ip); + /* ---- benchmark ---- */ + n_op = 50.0 * nfft * log2_nfft * log2_nfft; + printf("floating point operation: %g op.\n", n_op); + fprintf(f_log, "floating point operation: %g op.\n", n_op); + /* ---- difftime ---- */ + d_time = difftime(t_2, t_1); + printf("execution time: %g sec. (real time)\n", d_time); + fprintf(f_log, "execution time: %g sec. (real time)\n", d_time); + fclose(f_log); + return 0; +} + + +/* -------- multiple precision routines -------- */ + + +#include <math.h> +#include <float.h> +#include <stdio.h> + +/* ---- floating point format ---- + data := data[0] * pow(radix, data[1]) * + (data[2] + data[3]/radix + data[4]/radix/radix + ...), + data[0] : sign (1;data>0, -1;data<0, 0;data==0) + data[1] : exponent (0;data==0) + data[2...n+1] : digits + ---- function prototypes ---- + void mp_load_0(int n, int radix, int out[]); + void mp_load_1(int n, int radix, int out[]); + void mp_copy(int n, int radix, int in[], int out[]); + void mp_round(int n, int radix, int m, int inout[]); + int mp_cmp(int n, int radix, int in1[], int in2[]); + void mp_add(int n, int radix, int in1[], int in2[], int out[]); + void mp_sub(int n, int radix, int in1[], int in2[], int out[]); + void mp_imul(int n, int radix, int in1[], int in2, int out[]); + int mp_idiv(int n, int radix, int in1[], int in2, int out[]); + void mp_idiv_2(int n, int radix, int in[], int out[]); + double mp_mul_radix_test(int n, int radix, int nfft, + double tmpfft[], int ip[], double w[]); + void mp_mul(int n, int radix, int in1[], int in2[], int out[], + int tmp[], int nfft, double tmp1fft[], double tmp2fft[], + double tmp3fft[], int ip[], double w[]); + void mp_squ(int n, int radix, int in[], int out[], int tmp[], + int nfft, double tmp1fft[], double tmp2fft[], + int ip[], double w[]); + void mp_mulh(int n, int radix, int in1[], int in2[], int out[], + int nfft, double in1fft[], double outfft[], + int ip[], double w[]); + void mp_squh(int n, int radix, int in[], int out[], + int nfft, double inoutfft[], int ip[], double w[]); + int mp_inv(int n, int radix, int in[], int out[], + int tmp1[], int tmp2[], int nfft, + double tmp1fft[], double tmp2fft[], int ip[], double w[]); + int mp_sqrt(int n, int radix, int in[], int out[], + int tmp1[], int tmp2[], int nfft, + double tmp1fft[], double tmp2fft[], int ip[], double w[]); + void mp_sprintf(int n, int log10_radix, int in[], char out[]); + void mp_sscanf(int n, int log10_radix, char in[], int out[]); + void mp_fprintf(int n, int log10_radix, int in[], FILE *fout); + ---- +*/ + + +/* -------- mp_load routines -------- */ + + +void mp_load_0(int n, int radix, int out[]) +{ + int j; + + for (j = 0; j <= n + 1; j++) { + out[j] = 0; + } +} + + +void mp_load_1(int n, int radix, int out[]) +{ + int j; + + out[0] = 1; + out[1] = 0; + out[2] = 1; + for (j = 3; j <= n + 1; j++) { + out[j] = 0; + } +} + + +void mp_copy(int n, int radix, int in[], int out[]) +{ + int j; + + for (j = 0; j <= n + 1; j++) { + out[j] = in[j]; + } +} + + +void mp_round(int n, int radix, int m, int inout[]) +{ + int j, x; + + if (m < n) { + for (j = n + 1; j > m + 2; j--) { + inout[j] = 0; + } + x = 2 * inout[m + 2]; + inout[m + 2] = 0; + if (x >= radix) { + for (j = m + 1; j >= 2; j--) { + x = inout[j] + 1; + if (x < radix) { + inout[j] = x; + break; + } + inout[j] = 0; + } + if (x >= radix) { + inout[2] = 1; + inout[1]++; + } + } + } +} + + +/* -------- mp_add routines -------- */ + + +int mp_cmp(int n, int radix, int in1[], int in2[]) +{ + int mp_unsgn_cmp(int n, int in1[], int in2[]); + + if (in1[0] > in2[0]) { + return 1; + } else if (in1[0] < in2[0]) { + return -1; + } + return in1[0] * mp_unsgn_cmp(n, &in1[1], &in2[1]); +} + + +void mp_add(int n, int radix, int in1[], int in2[], int out[]) +{ + int mp_unsgn_cmp(int n, int in1[], int in2[]); + int mp_unexp_add(int n, int radix, int expdif, + int in1[], int in2[], int out[]); + int mp_unexp_sub(int n, int radix, int expdif, + int in1[], int in2[], int out[]); + int outsgn, outexp, expdif; + + expdif = in1[1] - in2[1]; + outexp = in1[1]; + if (expdif < 0) { + outexp = in2[1]; + } + outsgn = in1[0] * in2[0]; + if (outsgn >= 0) { + if (outsgn > 0) { + outsgn = in1[0]; + } else { + outsgn = in1[0] + in2[0]; + outexp = in1[1] + in2[1]; + expdif = 0; + } + if (expdif >= 0) { + outexp += mp_unexp_add(n, radix, expdif, + &in1[2], &in2[2], &out[2]); + } else { + outexp += mp_unexp_add(n, radix, -expdif, + &in2[2], &in1[2], &out[2]); + } + } else { + outsgn = mp_unsgn_cmp(n, &in1[1], &in2[1]); + if (outsgn >= 0) { + expdif = mp_unexp_sub(n, radix, expdif, + &in1[2], &in2[2], &out[2]); + } else { + expdif = mp_unexp_sub(n, radix, -expdif, + &in2[2], &in1[2], &out[2]); + } + outexp -= expdif; + outsgn *= in1[0]; + if (expdif == n) { + outsgn = 0; + } + } + if (outsgn == 0) { + outexp = 0; + } + out[0] = outsgn; + out[1] = outexp; +} + + +void mp_sub(int n, int radix, int in1[], int in2[], int out[]) +{ + int mp_unsgn_cmp(int n, int in1[], int in2[]); + int mp_unexp_add(int n, int radix, int expdif, + int in1[], int in2[], int out[]); + int mp_unexp_sub(int n, int radix, int expdif, + int in1[], int in2[], int out[]); + int outsgn, outexp, expdif; + + expdif = in1[1] - in2[1]; + outexp = in1[1]; + if (expdif < 0) { + outexp = in2[1]; + } + outsgn = in1[0] * in2[0]; + if (outsgn <= 0) { + if (outsgn < 0) { + outsgn = in1[0]; + } else { + outsgn = in1[0] - in2[0]; + outexp = in1[1] + in2[1]; + expdif = 0; + } + if (expdif >= 0) { + outexp += mp_unexp_add(n, radix, expdif, + &in1[2], &in2[2], &out[2]); + } else { + outexp += mp_unexp_add(n, radix, -expdif, + &in2[2], &in1[2], &out[2]); + } + } else { + outsgn = mp_unsgn_cmp(n, &in1[1], &in2[1]); + if (outsgn >= 0) { + expdif = mp_unexp_sub(n, radix, expdif, + &in1[2], &in2[2], &out[2]); + } else { + expdif = mp_unexp_sub(n, radix, -expdif, + &in2[2], &in1[2], &out[2]); + } + outexp -= expdif; + outsgn *= in1[0]; + if (expdif == n) { + outsgn = 0; + } + } + if (outsgn == 0) { + outexp = 0; + } + out[0] = outsgn; + out[1] = outexp; +} + + +/* -------- mp_add child routines -------- */ + + +int mp_unsgn_cmp(int n, int in1[], int in2[]) +{ + int j, cmp; + + cmp = 0; + for (j = 0; j <= n && cmp == 0; j++) { + cmp = in1[j] - in2[j]; + } + if (cmp > 0) { + cmp = 1; + } else if (cmp < 0) { + cmp = -1; + } + return cmp; +} + + +int mp_unexp_add(int n, int radix, int expdif, + int in1[], int in2[], int out[]) +{ + int j, x, carry; + + carry = 0; + if (expdif == 0 && in1[0] + in2[0] >= radix) { + x = in1[n - 1] + in2[n - 1]; + carry = x >= radix ? -1 : 0; + for (j = n - 1; j > 0; j--) { + x = in1[j - 1] + in2[j - 1] - carry; + carry = x >= radix ? -1 : 0; + out[j] = x - (radix & carry); + } + out[0] = -carry; + } else { + if (expdif > n) { + expdif = n; + } + for (j = n - 1; j >= expdif; j--) { + x = in1[j] + in2[j - expdif] - carry; + carry = x >= radix ? -1 : 0; + out[j] = x - (radix & carry); + } + for (j = expdif - 1; j >= 0; j--) { + x = in1[j] - carry; + carry = x >= radix ? -1 : 0; + out[j] = x - (radix & carry); + } + if (carry != 0) { + for (j = n - 1; j > 0; j--) { + out[j] = out[j - 1]; + } + out[0] = -carry; + } + } + return -carry; +} + + +int mp_unexp_sub(int n, int radix, int expdif, + int in1[], int in2[], int out[]) +{ + int j, x, borrow, ncancel; + + if (expdif > n) { + expdif = n; + } + borrow = 0; + for (j = n - 1; j >= expdif; j--) { + x = in1[j] - in2[j - expdif] + borrow; + borrow = x < 0 ? -1 : 0; + out[j] = x + (radix & borrow); + } + for (j = expdif - 1; j >= 0; j--) { + x = in1[j] + borrow; + borrow = x < 0 ? -1 : 0; + out[j] = x + (radix & borrow); + } + ncancel = 0; + for (j = 0; j < n && out[j] == 0; j++) { + ncancel = j + 1; + } + if (ncancel > 0 && ncancel < n) { + for (j = 0; j < n - ncancel; j++) { + out[j] = out[j + ncancel]; + } + for (j = n - ncancel; j < n; j++) { + out[j] = 0; + } + } + return ncancel; +} + + +/* -------- mp_imul routines -------- */ + + +void mp_imul(int n, int radix, int in1[], int in2, int out[]) +{ + void mp_unsgn_imul(int n, double dradix, int in1[], double din2, + int out[]); + + if (in2 > 0) { + out[0] = in1[0]; + } else if (in2 < 0) { + out[0] = -in1[0]; + in2 = -in2; + } else { + out[0] = 0; + } + mp_unsgn_imul(n, radix, &in1[1], in2, &out[1]); + if (out[0] == 0) { + out[1] = 0; + } +} + + +int mp_idiv(int n, int radix, int in1[], int in2, int out[]) +{ + void mp_load_0(int n, int radix, int out[]); + void mp_unsgn_idiv(int n, double dradix, int in1[], double din2, + int out[]); + + if (in2 == 0) { + return -1; + } + if (in2 > 0) { + out[0] = in1[0]; + } else { + out[0] = -in1[0]; + in2 = -in2; + } + if (in1[0] == 0) { + mp_load_0(n, radix, out); + return 0; + } + mp_unsgn_idiv(n, radix, &in1[1], in2, &out[1]); + return 0; +} + + +void mp_idiv_2(int n, int radix, int in[], int out[]) +{ + int j, ix, carry, shift; + + out[0] = in[0]; + shift = 0; + if (in[2] == 1) { + shift = 1; + } + out[1] = in[1] - shift; + carry = -shift; + for (j = 2; j <= n + 1 - shift; j++) { + ix = in[j + shift] + (radix & carry); + carry = -(ix & 1); + out[j] = ix >> 1; + } + if (shift > 0) { + out[n + 1] = (radix & carry) >> 1; + } +} + + +/* -------- mp_imul child routines -------- */ + + +void mp_unsgn_imul(int n, double dradix, int in1[], double din2, + int out[]) +{ + int j, carry, shift; + double x, d1_radix; + + d1_radix = 1.0 / dradix; + carry = 0; + for (j = n; j >= 1; j--) { + x = din2 * in1[j] + carry + 0.5; + carry = (int) (d1_radix * x); + out[j] = (int) (x - dradix * carry); + } + shift = 0; + x = carry + 0.5; + while (x > 1) { + x *= d1_radix; + shift++; + } + out[0] = in1[0] + shift; + if (shift > 0) { + while (shift > n) { + carry = (int) (d1_radix * carry + 0.5); + shift--; + } + for (j = n; j >= shift + 1; j--) { + out[j] = out[j - shift]; + } + for (j = shift; j >= 1; j--) { + x = carry + 0.5; + carry = (int) (d1_radix * x); + out[j] = (int) (x - dradix * carry); + } + } +} + + +void mp_unsgn_idiv(int n, double dradix, int in1[], double din2, + int out[]) +{ + int j, ix, carry, shift; + double x, d1_in2; + + d1_in2 = 1.0 / din2; + shift = 0; + x = 0; + do { + shift++; + x *= dradix; + if (shift <= n) { + x += in1[shift]; + } + } while (x < din2 - 0.5); + x += 0.5; + ix = (int) (d1_in2 * x); + carry = (int) (x - din2 * ix); + out[1] = ix; + shift--; + out[0] = in1[0] - shift; + if (shift >= n) { + shift = n - 1; + } + for (j = 2; j <= n - shift; j++) { + x = in1[j + shift] + dradix * carry + 0.5; + ix = (int) (d1_in2 * x); + carry = (int) (x - din2 * ix); + out[j] = ix; + } + for (j = n - shift + 1; j <= n; j++) { + x = dradix * carry + 0.5; + ix = (int) (d1_in2 * x); + carry = (int) (x - din2 * ix); + out[j] = ix; + } +} + + +/* -------- mp_mul routines -------- */ + + +double mp_mul_radix_test(int n, int radix, int nfft, + double tmpfft[], int ip[], double w[]) +{ + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_csqu(int nfft, double dinout[]); + double mp_mul_d2i_test(int radix, int nfft, double din[]); + int j, ndata, radix_2; + + ndata = (nfft >> 1) + 1; + if (ndata > n) { + ndata = n; + } + tmpfft[nfft + 1] = radix - 1; + for (j = nfft; j > ndata; j--) { + tmpfft[j] = 0; + } + radix_2 = (radix + 1) / 2; + for (j = ndata; j > 2; j--) { + tmpfft[j] = radix_2; + } + tmpfft[2] = radix; + tmpfft[1] = radix - 1; + tmpfft[0] = 0; + rdft(nfft, 1, &tmpfft[1], ip, w); + mp_mul_csqu(nfft, tmpfft); + rdft(nfft, -1, &tmpfft[1], ip, w); + return 2 * mp_mul_d2i_test(radix, nfft, tmpfft); +} + + +void mp_mul(int n, int radix, int in1[], int in2[], int out[], + int tmp[], int nfft, double tmp1fft[], double tmp2fft[], + double tmp3fft[], int ip[], double w[]) +{ + void mp_copy(int n, int radix, int in[], int out[]); + void mp_add(int n, int radix, int in1[], int in2[], int out[]); + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_i2d(int n, int radix, int nfft, int shift, + int in[], double dout[]); + void mp_mul_cmul(int nfft, double din[], double dinout[]); + void mp_mul_cmuladd(int nfft, double din1[], double din2[], + double dinout[]); + void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]); + int n_h, shift; + + shift = (nfft >> 1) + 1; + while (n > shift) { + if (in1[shift + 2] + in2[shift + 2] != 0) { + break; + } + shift++; + } + n_h = n / 2 + 1; + if (n_h < n - shift) { + n_h = n - shift; + } + /* ---- tmp3fft = (upper) in1 * (lower) in2 ---- */ + mp_mul_i2d(n, radix, nfft, 0, in1, tmp1fft); + rdft(nfft, 1, &tmp1fft[1], ip, w); + mp_mul_i2d(n, radix, nfft, shift, in2, tmp3fft); + rdft(nfft, 1, &tmp3fft[1], ip, w); + mp_mul_cmul(nfft, tmp1fft, tmp3fft); + /* ---- tmp = (upper) in1 * (upper) in2 ---- */ + mp_mul_i2d(n, radix, nfft, 0, in2, tmp2fft); + rdft(nfft, 1, &tmp2fft[1], ip, w); + mp_mul_cmul(nfft, tmp2fft, tmp1fft); + rdft(nfft, -1, &tmp1fft[1], ip, w); + mp_mul_d2i(n, radix, nfft, tmp1fft, tmp); + /* ---- tmp3fft += (upper) in2 * (lower) in1 ---- */ + mp_mul_i2d(n, radix, nfft, shift, in1, tmp1fft); + rdft(nfft, 1, &tmp1fft[1], ip, w); + mp_mul_cmuladd(nfft, tmp1fft, tmp2fft, tmp3fft); + /* ---- out = tmp + tmp3fft ---- */ + rdft(nfft, -1, &tmp3fft[1], ip, w); + mp_mul_d2i(n_h, radix, nfft, tmp3fft, out); + if (out[0] != 0) { + mp_add(n, radix, out, tmp, out); + } else { + mp_copy(n, radix, tmp, out); + } +} + + +void mp_squ(int n, int radix, int in[], int out[], int tmp[], + int nfft, double tmp1fft[], double tmp2fft[], + int ip[], double w[]) +{ + void mp_add(int n, int radix, int in1[], int in2[], int out[]); + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_i2d(int n, int radix, int nfft, int shift, + int in[], double dout[]); + void mp_mul_cmul(int nfft, double din[], double dinout[]); + void mp_mul_csqu(int nfft, double dinout[]); + void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]); + int n_h, shift; + + shift = (nfft >> 1) + 1; + while (n > shift) { + if (in[shift + 2] != 0) { + break; + } + shift++; + } + n_h = n / 2 + 1; + if (n_h < n - shift) { + n_h = n - shift; + } + /* ---- tmp = (upper) in * (lower) in ---- */ + mp_mul_i2d(n, radix, nfft, 0, in, tmp1fft); + rdft(nfft, 1, &tmp1fft[1], ip, w); + mp_mul_i2d(n, radix, nfft, shift, in, tmp2fft); + rdft(nfft, 1, &tmp2fft[1], ip, w); + mp_mul_cmul(nfft, tmp1fft, tmp2fft); + rdft(nfft, -1, &tmp2fft[1], ip, w); + mp_mul_d2i(n_h, radix, nfft, tmp2fft, tmp); + /* ---- out = 2 * tmp + ((upper) in)^2 ---- */ + mp_mul_csqu(nfft, tmp1fft); + rdft(nfft, -1, &tmp1fft[1], ip, w); + mp_mul_d2i(n, radix, nfft, tmp1fft, out); + if (tmp[0] != 0) { + mp_add(n_h, radix, tmp, tmp, tmp); + mp_add(n, radix, out, tmp, out); + } +} + + +void mp_mulh(int n, int radix, int in1[], int in2[], int out[], + int nfft, double in1fft[], double outfft[], int ip[], double w[]) +{ + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_i2d(int n, int radix, int nfft, int shift, + int in[], double dout[]); + void mp_mul_cmul(int nfft, double din[], double dinout[]); + void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]); + + mp_mul_i2d(n, radix, nfft, 0, in1, in1fft); + rdft(nfft, 1, &in1fft[1], ip, w); + mp_mul_i2d(n, radix, nfft, 0, in2, outfft); + rdft(nfft, 1, &outfft[1], ip, w); + mp_mul_cmul(nfft, in1fft, outfft); + rdft(nfft, -1, &outfft[1], ip, w); + mp_mul_d2i(n, radix, nfft, outfft, out); +} + + +void mp_mulh_use_in1fft(int n, int radix, double in1fft[], + int shift, int in2[], int out[], int nfft, double outfft[], + int ip[], double w[]) +{ + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_i2d(int n, int radix, int nfft, int shift, + int in[], double dout[]); + void mp_mul_cmul(int nfft, double din[], double dinout[]); + void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]); + int n_h; + + while (n > shift) { + if (in2[shift + 2] != 0) { + break; + } + shift++; + } + n_h = n / 2 + 1; + if (n_h < n - shift) { + n_h = n - shift; + } + mp_mul_i2d(n, radix, nfft, shift, in2, outfft); + rdft(nfft, 1, &outfft[1], ip, w); + mp_mul_cmul(nfft, in1fft, outfft); + rdft(nfft, -1, &outfft[1], ip, w); + mp_mul_d2i(n_h, radix, nfft, outfft, out); +} + + +void mp_squh(int n, int radix, int in[], int out[], + int nfft, double inoutfft[], int ip[], double w[]) +{ + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_i2d(int n, int radix, int nfft, int shift, + int in[], double dout[]); + void mp_mul_csqu(int nfft, double dinout[]); + void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]); + + mp_mul_i2d(n, radix, nfft, 0, in, inoutfft); + rdft(nfft, 1, &inoutfft[1], ip, w); + mp_mul_csqu(nfft, inoutfft); + rdft(nfft, -1, &inoutfft[1], ip, w); + mp_mul_d2i(n, radix, nfft, inoutfft, out); +} + + +void mp_squh_use_in1fft(int n, int radix, double inoutfft[], int out[], + int nfft, int ip[], double w[]) +{ + void rdft(int n, int isgn, double *a, int *ip, double *w); + void mp_mul_csqu(int nfft, double dinout[]); + void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]); + + mp_mul_csqu(nfft, inoutfft); + rdft(nfft, -1, &inoutfft[1], ip, w); + mp_mul_d2i(n, radix, nfft, inoutfft, out); +} + + +/* -------- mp_mul child routines -------- */ + + +void mp_mul_i2d(int n, int radix, int nfft, int shift, + int in[], double dout[]) +{ + int j, x, carry, ndata, radix_2, topdgt; + + ndata = 0; + topdgt = 0; + if (n > shift) { + topdgt = in[shift + 2]; + ndata = (nfft >> 1) + 1; + if (ndata > n - shift) { + ndata = n - shift; + } + } + dout[nfft + 1] = in[0] * topdgt; + for (j = nfft; j > ndata; j--) { + dout[j] = 0; + } + /* ---- abs(dout[j]) <= radix/2 (to keep FFT precision) ---- */ + if (ndata > 1) { + radix_2 = radix / 2; + carry = 0; + for (j = ndata + 1; j > 3; j--) { + x = in[j + shift] - carry; + carry = x >= radix_2 ? -1 : 0; + dout[j - 1] = x - (radix & carry); + } + dout[2] = in[shift + 3] - carry; + } + dout[1] = topdgt; + dout[0] = in[1] - shift; +} + + +void mp_mul_cmul(int nfft, double din[], double dinout[]) +{ + int j; + double xr, xi, yr, yi; + + dinout[0] += din[0]; + dinout[1] *= din[1]; + dinout[2] *= din[2]; + for (j = 3; j < nfft; j += 2) { + xr = din[j]; + xi = din[j + 1]; + yr = dinout[j]; + yi = dinout[j + 1]; + dinout[j] = xr * yr - xi * yi; + dinout[j + 1] = xr * yi + xi * yr; + } + dinout[nfft + 1] *= din[nfft + 1]; +} + + +void mp_mul_cmuladd(int nfft, double din1[], double din2[], + double dinout[]) +{ + int j; + double xr, xi, yr, yi; + + dinout[1] += din1[1] * din2[1]; + dinout[2] += din1[2] * din2[2]; + for (j = 3; j < nfft; j += 2) { + xr = din1[j]; + xi = din1[j + 1]; + yr = din2[j]; + yi = din2[j + 1]; + dinout[j] += xr * yr - xi * yi; + dinout[j + 1] += xr * yi + xi * yr; + } + dinout[nfft + 1] += din1[nfft + 1] * din2[nfft + 1]; +} + + +void mp_mul_csqu(int nfft, double dinout[]) +{ + int j; + double xr, xi; + + dinout[0] *= 2; + dinout[1] *= dinout[1]; + dinout[2] *= dinout[2]; + for (j = 3; j < nfft; j += 2) { + xr = dinout[j]; + xi = dinout[j + 1]; + dinout[j] = xr * xr - xi * xi; + dinout[j + 1] = 2 * xr * xi; + } + dinout[nfft + 1] *= dinout[nfft + 1]; +} + + +void mp_mul_d2i(int n, int radix, int nfft, double din[], int out[]) +{ + int j, carry, carry1, carry2, shift, ndata; + double x, scale, d1_radix, d1_radix2, pow_radix, topdgt; + + scale = 2.0 / nfft; + d1_radix = 1.0 / radix; + d1_radix2 = d1_radix * d1_radix; + topdgt = din[nfft + 1]; + x = topdgt < 0 ? -topdgt : topdgt; + shift = x + 0.5 >= radix ? 1 : 0; + /* ---- correction of cyclic convolution of din[1] ---- */ + x *= nfft * 0.5; + din[nfft + 1] = din[1] - x; + din[1] = x; + /* ---- output of digits ---- */ + ndata = n; + if (n > nfft + 1 + shift) { + ndata = nfft + 1 + shift; + for (j = n + 1; j > ndata + 1; j--) { + out[j] = 0; + } + } + x = 0; + pow_radix = 1; + for (j = ndata + 1 - shift; j <= nfft + 1; j++) { + x += pow_radix * din[j]; + pow_radix *= d1_radix; + if (pow_radix < DBL_EPSILON) { + break; + } + } + x = d1_radix2 * (scale * x + 0.5); + carry2 = ((int) x) - 1; + carry = (int) (radix * (x - carry2) + 0.5); + for (j = ndata; j > 1; j--) { + x = d1_radix2 * (scale * din[j - shift] + carry + 0.5); + carry = carry2; + carry2 = ((int) x) - 1; + x = radix * (x - carry2); + carry1 = (int) x; + out[j + 1] = (int) (radix * (x - carry1)); + carry += carry1; + } + x = carry + ((double) radix) * carry2 + 0.5; + if (shift == 0) { + x += scale * din[1]; + } + carry = (int) (d1_radix * x); + out[2] = (int) (x - ((double) radix) * carry); + if (carry > 0) { + for (j = n + 1; j > 2; j--) { + out[j] = out[j - 1]; + } + out[2] = carry; + shift++; + } + /* ---- output of exp, sgn ---- */ + x = din[0] + shift + 0.5; + shift = ((int) x) - 1; + out[1] = shift + ((int) (x - shift)); + out[0] = topdgt > 0.5 ? 1 : -1; + if (out[2] == 0) { + out[0] = 0; + out[1] = 0; + } +} + + +double mp_mul_d2i_test(int radix, int nfft, double din[]) +{ + int j, carry, carry1, carry2; + double x, scale, d1_radix, d1_radix2, err; + + scale = 2.0 / nfft; + d1_radix = 1.0 / radix; + d1_radix2 = d1_radix * d1_radix; + /* ---- correction of cyclic convolution of din[1] ---- */ + x = din[nfft + 1] * nfft * 0.5; + if (x < 0) { + x = -x; + } + din[nfft + 1] = din[1] - x; + /* ---- check of digits ---- */ + err = 0; + carry = 0; + carry2 = 0; + for (j = nfft + 1; j > 1; j--) { + x = d1_radix2 * (scale * din[j] + carry + 0.5); + carry = carry2; + carry2 = ((int) x) - 1; + x = radix * (x - carry2); + carry1 = (int) x; + x = radix * (x - carry1); + carry += carry1; + x = x - 0.5 - ((int) x); + if (x > err) { + err = x; + } else if (-x > err) { + err = -x; + } + } + return err; +} + + +/* -------- mp_inv routines -------- */ + + +int mp_inv(int n, int radix, int in[], int out[], + int tmp1[], int tmp2[], int nfft, + double tmp1fft[], double tmp2fft[], int ip[], double w[]) +{ + int mp_get_nfft_init(int radix, int nfft_max); + void mp_inv_init(int n, int radix, int in[], int out[]); + int mp_inv_newton(int n, int radix, int in[], int inout[], + int tmp1[], int tmp2[], int nfft, double tmp1fft[], + double tmp2fft[], int ip[], double w[]); + int n_nwt, nfft_nwt, thr, prc; + + if (in[0] == 0) { + return -1; + } + nfft_nwt = mp_get_nfft_init(radix, nfft); + n_nwt = nfft_nwt + 2; + if (n_nwt > n) { + n_nwt = n; + } + mp_inv_init(n_nwt, radix, in, out); + thr = 8; + do { + n_nwt = nfft_nwt + 2; + if (n_nwt > n) { + n_nwt = n; + } + prc = mp_inv_newton(n_nwt, radix, in, out, + tmp1, tmp2, nfft_nwt, tmp1fft, tmp2fft, ip, w); + if (thr * nfft_nwt >= nfft) { + thr = 0; + if (2 * prc <= n_nwt - 2) { + nfft_nwt >>= 1; + } + } else { + if (3 * prc < n_nwt - 2) { + nfft_nwt >>= 1; + } + } + nfft_nwt <<= 1; + } while (nfft_nwt <= nfft); + return 0; +} + + +int mp_sqrt(int n, int radix, int in[], int out[], + int tmp1[], int tmp2[], int nfft, + double tmp1fft[], double tmp2fft[], int ip[], double w[]) +{ + void mp_load_0(int n, int radix, int out[]); + int mp_get_nfft_init(int radix, int nfft_max); + void mp_sqrt_init(int n, int radix, int in[], int out[], int out_rev[]); + int mp_sqrt_newton(int n, int radix, int in[], int inout[], + int inout_rev[], int tmp[], int nfft, double tmp1fft[], + double tmp2fft[], int ip[], double w[], int *n_tmp1fft); + int n_nwt, nfft_nwt, thr, prc, n_tmp1fft; + + if (in[0] < 0) { + return -1; + } else if (in[0] == 0) { + mp_load_0(n, radix, out); + return 0; + } + nfft_nwt = mp_get_nfft_init(radix, nfft); + n_nwt = nfft_nwt + 2; + if (n_nwt > n) { + n_nwt = n; + } + mp_sqrt_init(n_nwt, radix, in, out, tmp1); + n_tmp1fft = 0; + thr = 8; + do { + n_nwt = nfft_nwt + 2; + if (n_nwt > n) { + n_nwt = n; + } + prc = mp_sqrt_newton(n_nwt, radix, in, out, + tmp1, tmp2, nfft_nwt, tmp1fft, tmp2fft, + ip, w, &n_tmp1fft); + if (thr * nfft_nwt >= nfft) { + thr = 0; + if (2 * prc <= n_nwt - 2) { + nfft_nwt >>= 1; + } + } else { + if (3 * prc < n_nwt - 2) { + nfft_nwt >>= 1; + } + } + nfft_nwt <<= 1; + } while (nfft_nwt <= nfft); + return 0; +} + + +/* -------- mp_inv child routines -------- */ + + +int mp_get_nfft_init(int radix, int nfft_max) +{ + int nfft_init; + double r; + + r = radix; + nfft_init = 1; + do { + r *= r; + nfft_init <<= 1; + } while (DBL_EPSILON * r < 1 && nfft_init < nfft_max); + return nfft_init; +} + + +void mp_inv_init(int n, int radix, int in[], int out[]) +{ + void mp_unexp_d2mp(int n, int radix, double din, int out[]); + double mp_unexp_mp2d(int n, int radix, int in[]); + int outexp; + double din; + + out[0] = in[0]; + outexp = -in[1]; + din = 1.0 / mp_unexp_mp2d(n, radix, &in[2]); + while (din < 1) { + din *= radix; + outexp--; + } + out[1] = outexp; + mp_unexp_d2mp(n, radix, din, &out[2]); +} + + +void mp_sqrt_init(int n, int radix, int in[], int out[], int out_rev[]) +{ + void mp_unexp_d2mp(int n, int radix, double din, int out[]); + double mp_unexp_mp2d(int n, int radix, int in[]); + int outexp; + double din; + + out[0] = 1; + out_rev[0] = 1; + outexp = in[1]; + din = mp_unexp_mp2d(n, radix, &in[2]); + if (outexp % 2 != 0) { + din *= radix; + outexp--; + } + outexp /= 2; + din = sqrt(din); + if (din < 1) { + din *= radix; + outexp--; + } + out[1] = outexp; + mp_unexp_d2mp(n, radix, din, &out[2]); + outexp = -outexp; + din = 1.0 / din; + while (din < 1) { + din *= radix; + outexp--; + } + out_rev[1] = outexp; + mp_unexp_d2mp(n, radix, din, &out_rev[2]); +} + + +void mp_unexp_d2mp(int n, int radix, double din, int out[]) +{ + int j, x; + + for (j = 0; j < n; j++) { + x = (int) din; + if (x >= radix) { + x = radix - 1; + din = radix; + } + din = radix * (din - x); + out[j] = x; + } +} + + +double mp_unexp_mp2d(int n, int radix, int in[]) +{ + int j; + double d1_radix, dout; + + d1_radix = 1.0 / radix; + dout = 0; + for (j = n - 1; j >= 0; j--) { + dout = d1_radix * dout + in[j]; + } + return dout; +} + + +int mp_inv_newton(int n, int radix, int in[], int inout[], + int tmp1[], int tmp2[], int nfft, double tmp1fft[], + double tmp2fft[], int ip[], double w[]) +{ + void mp_load_1(int n, int radix, int out[]); + void mp_round(int n, int radix, int m, int inout[]); + void mp_add(int n, int radix, int in1[], int in2[], int out[]); + void mp_sub(int n, int radix, int in1[], int in2[], int out[]); + void mp_mulh(int n, int radix, int in1[], int in2[], int out[], + int nfft, double in1fft[], double outfft[], + int ip[], double w[]); + void mp_mulh_use_in1fft(int n, int radix, double in1fft[], + int shift, int in2[], int out[], int nfft, double outfft[], + int ip[], double w[]); + int n_h, shift, prc; + + shift = (nfft >> 1) + 1; + n_h = n / 2 + 1; + if (n_h < n - shift) { + n_h = n - shift; + } + /* ---- tmp1 = inout * (upper) in (half to normal precision) ---- */ + mp_round(n, radix, shift, inout); + mp_mulh(n, radix, inout, in, tmp1, + nfft, tmp1fft, tmp2fft, ip, w); + /* ---- tmp2 = 1 - tmp1 ---- */ + mp_load_1(n, radix, tmp2); + mp_sub(n, radix, tmp2, tmp1, tmp2); + /* ---- tmp2 -= inout * (lower) in (half precision) ---- */ + mp_mulh_use_in1fft(n, radix, tmp1fft, shift, in, tmp1, + nfft, tmp2fft, ip, w); + mp_sub(n_h, radix, tmp2, tmp1, tmp2); + /* ---- get precision ---- */ + prc = -tmp2[1]; + if (tmp2[0] == 0) { + prc = nfft + 1; + } + /* ---- tmp2 *= inout (half precision) ---- */ + mp_mulh_use_in1fft(n_h, radix, tmp1fft, 0, tmp2, tmp2, + nfft, tmp2fft, ip, w); + /* ---- inout += tmp2 ---- */ + if (tmp2[0] != 0) { + mp_add(n, radix, inout, tmp2, inout); + } + return prc; +} + + +int mp_sqrt_newton(int n, int radix, int in[], int inout[], + int inout_rev[], int tmp[], int nfft, double tmp1fft[], + double tmp2fft[], int ip[], double w[], int *n_tmp1fft) +{ + void mp_round(int n, int radix, int m, int inout[]); + void mp_add(int n, int radix, int in1[], int in2[], int out[]); + void mp_sub(int n, int radix, int in1[], int in2[], int out[]); + void mp_idiv_2(int n, int radix, int in[], int out[]); + void mp_mulh(int n, int radix, int in1[], int in2[], int out[], + int nfft, double in1fft[], double outfft[], + int ip[], double w[]); + void mp_squh(int n, int radix, int in[], int out[], + int nfft, double inoutfft[], int ip[], double w[]); + void mp_squh_use_in1fft(int n, int radix, double inoutfft[], int out[], + int nfft, int ip[], double w[]); + int n_h, nfft_h, shift, prc; + + nfft_h = nfft >> 1; + shift = nfft_h + 1; + if (nfft_h < 2) { + nfft_h = 2; + } + n_h = n / 2 + 1; + if (n_h < n - shift) { + n_h = n - shift; + } + /* ---- tmp = inout_rev^2 (1/4 to half precision) ---- */ + mp_round(n_h, radix, (nfft_h >> 1) + 1, inout_rev); + if (*n_tmp1fft != nfft_h) { + mp_squh(n_h, radix, inout_rev, tmp, + nfft_h, tmp1fft, ip, w); + } else { + mp_squh_use_in1fft(n_h, radix, tmp1fft, tmp, + nfft_h, ip, w); + } + /* ---- tmp = inout_rev - inout * tmp (half precision) ---- */ + mp_round(n, radix, shift, inout); + mp_mulh(n_h, radix, inout, tmp, tmp, + nfft, tmp1fft, tmp2fft, ip, w); + mp_sub(n_h, radix, inout_rev, tmp, tmp); + /* ---- inout_rev += tmp ---- */ + mp_add(n_h, radix, inout_rev, tmp, inout_rev); + /* ---- tmp = in - inout^2 (half to normal precision) ---- */ + mp_squh_use_in1fft(n, radix, tmp1fft, tmp, + nfft, ip, w); + mp_sub(n, radix, in, tmp, tmp); + /* ---- get precision ---- */ + prc = in[1] - tmp[1]; + if (in[2] > tmp[2]) { + prc++; + } + if (tmp[0] == 0) { + prc = nfft + 1; + } + /* ---- tmp = tmp * inout_rev / 2 (half precision) ---- */ + mp_round(n_h, radix, shift, inout_rev); + mp_mulh(n_h, radix, inout_rev, tmp, tmp, + nfft, tmp1fft, tmp2fft, ip, w); + *n_tmp1fft = nfft; + mp_idiv_2(n_h, radix, tmp, tmp); + /* ---- inout += tmp ---- */ + if (tmp[0] != 0) { + mp_add(n, radix, inout, tmp, inout); + } + return prc; +} + + +/* -------- mp_io routines -------- */ + + +void mp_sprintf(int n, int log10_radix, int in[], char out[]) +{ + int j, k, x, y, outexp, shift; + + if (in[0] < 0) { + *out++ = '-'; + } + x = in[2]; + shift = log10_radix; + for (k = log10_radix; k > 0; k--) { + y = x % 10; + x /= 10; + out[k] = '0' + y; + if (y != 0) { + shift = k; + } + } + out[0] = out[shift]; + out[1] = '.'; + for (k = 1; k <= log10_radix - shift; k++) { + out[k + 1] = out[k + shift]; + } + outexp = log10_radix - shift; + out += outexp + 2; + for (j = 3; j <= n + 1; j++) { + x = in[j]; + for (k = log10_radix - 1; k >= 0; k--) { + y = x % 10; + x /= 10; + out[k] = '0' + y; + } + out += log10_radix; + } + *out++ = 'e'; + outexp += log10_radix * in[1]; + sprintf(out, "%d", outexp); +} + + +void mp_sscanf(int n, int log10_radix, char in[], int out[]) +{ + char *s; + int j, x, outexp, outexp_mod; + + while (*in == ' ') { + in++; + } + out[0] = 1; + if (*in == '-') { + out[0] = -1; + in++; + } else if (*in == '+') { + in++; + } + while (*in == ' ' || *in == '0') { + in++; + } + outexp = 0; + for (s = in; *s != '\0'; s++) { + if (*s == 'e' || *s == 'E' || *s == 'd' || *s == 'D') { + if (sscanf(++s, "%d", &outexp) != 1) { + outexp = 0; + } + break; + } + } + if (*in == '.') { + do { + outexp--; + while (*++in == ' '); + } while (*in == '0' && *in != '\0'); + } else if (*in != '\0') { + s = in; + while (*++s == ' '); + while (*s >= '0' && *s <= '9' && *s != '\0') { + outexp++; + while (*++s == ' '); + } + } + x = outexp / log10_radix; + outexp_mod = outexp - log10_radix * x; + if (outexp_mod < 0) { + x--; + outexp_mod += log10_radix; + } + out[1] = x; + x = 0; + j = 2; + for (s = in; *s != '\0'; s++) { + if (*s == '.' || *s == ' ') { + continue; + } + if (*s < '0' || *s > '9') { + break; + } + x = 10 * x + (*s - '0'); + if (--outexp_mod < 0) { + if (j > n + 1) { + break; + } + out[j++] = x; + x = 0; + outexp_mod = log10_radix - 1; + } + } + while (outexp_mod-- >= 0) { + x *= 10; + } + while (j <= n + 1) { + out[j++] = x; + x = 0; + } + if (out[2] == 0) { + out[0] = 0; + out[1] = 0; + } +} + + +void mp_fprintf(int n, int log10_radix, int in[], FILE *fout) +{ + int j, k, x, y, outexp, shift; + char out[256]; + + if (in[0] < 0) { + putc('-', fout); + } + x = in[2]; + shift = log10_radix; + for (k = log10_radix; k > 0; k--) { + y = x % 10; + x /= 10; + out[k] = '0' + y; + if (y != 0) { + shift = k; + } + } + putc(out[shift], fout); + putc('.', fout); + for (k = 1; k <= log10_radix - shift; k++) { + putc(out[k + shift], fout); + } + outexp = log10_radix - shift; + for (j = 3; j <= n + 1; j++) { + x = in[j]; + for (k = log10_radix - 1; k >= 0; k--) { + y = x % 10; + x /= 10; + out[k] = '0' + y; + } + for (k = 0; k < log10_radix; k++) { + putc(out[k], fout); + } + } + putc('e', fout); + outexp += log10_radix * in[1]; + sprintf(out, "%d", outexp); + for (k = 0; out[k] != '\0'; k++) { + putc(out[k], fout); + } +} + + |