pi_fft.c

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

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