pi_fft.c

This is a package to calculate Discrete Fourier/Cosine/Sine Transforms of 1-dimensional sequences of

    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);
    }
}