Google Git
Sign in
fuchsia / third_party / github.com / petewarden / OouraFFT / c6fd2dd6d21397baa6653139d31d84540d5449a2 / . / sample2 / pi_fft.c
blob: c9a76bf88a1eee2f05de0987a117871ccb7fc3a7 [file] [log] [blame]
/*
---- 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);
}
}