wftavect.cpp

任意精度计算的实现

    x[8] = m1;
    x[1] = t8 + t16;
    x[15] = t8 - t16;
    x[2] = t0 + t2;
    x[14] = t0 - t2;
    x[13] = t11 + t19;
    x[3] = t11 - t19;
    x[4] = m2 + m10;
    x[12] = m2 - m10;
    x[5] = t10 + t18;
    x[11] = t10 - t18;
    x[6] = t1 + t3;
    x[10] = t1 - t3;
    x[9] = t9 + t17;
    x[7] = t9 - t17;
}

void (*wftas[]) (modint x[], modint w[]) =
{0, 0, wfta2, wfta3, wfta4, wfta5, 0, wfta7, wfta8, wfta9, 0, 0, 0, 0, 0, 0, wfta16};

void permute (modint x[], int p[], int n)
{
    int k;
    modint t[MAXN];

    for (k = 0; k < n; k++)
        t[k] = x[p[k]];

    for (k = 0; k < n; k++)
        x[k] = t[k];
}

void unpermute (modint x[], int p[], int n)
{
    int k;
    modint t[MAXN];

    for (k = 0; k < n; k++)
        t[p[k]] = x[k];

    for (k = 0; k < n; k++)
        x[k] = t[k];
}

extern void (*wftan[]) (modint x[], modint w[], int n, int s[], int m, int ms[]);

void wftan2 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, nn, nm;
    modint *p;

    p = x + n;

    vecaddsub (x, p, x, p, n);

    k = s[0];

    if (k == n)
    {
        wftas[k] (x, w);
        wftas[k] (p, w + m);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        wftan[k] (x, w, nn, s, nm, ms);
        wftan[k] (p, w + m, nn, s, nm, ms);
    }
}

void wftan3 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, nn, nm;
    modint *p[3];

    p[2] = (p[1] = (p[0] = x) + n) + n;

    vecaddsub (p[1], p[2], p[1], p[2], n);
    vecadd (p[0], p[0], p[1], n);

    k = s[0];

    if (k == n)
    {
        wftas[k] (x, w);
        wftas[k] (x + k, w + m);
        wftas[k] (x + k + k, w + m + m);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        wftan[k] (x, w, nn, s, nm, ms);
        wftan[k] (x + n, w + m, nn, s, nm, ms);
        wftan[k] (x + n + n, w + m + m, nn, s, nm, ms);
    }

    vecadd (p[1], p[1], p[0], n);
    vecaddsub (p[1], p[2], p[1], p[2], n);
}

void wftan4 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, na, ma;
    modint *p[4], *t[2] = {t0, t1};

    p[3] = (p[2] = (p[1] = (p[0] = x) + n) + n) + n;

    vecaddsub (t[0], p[2], p[0], p[2], n);
    vecaddsub (t[1], p[3], p[1], p[3], n);
    vecaddsub (p[0], p[1], t[0], t[1], n);

    k = s[0];

    if (k == n)
    {
        for (na = ma = j = 0; j < 4; j++, na += k, ma += m)
            wftas[k] (x + na, w + ma);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        for (na = ma = j = 0; j < 4; j++, na += n, ma += m)
            wftan[k] (x + na, w + ma, nn, s, nm, ms);
    }

    vecaddsub (t[0], p[3], p[2], p[3], n);
    veccopy (p[2], p[1], n);
    veccopy (p[1], t[0], n);
}

void wftan5 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint *p[6], *t[2] = {t0, t1};

    for (p[0] = x, k = 1; k < 5; k++)
        p[k] = p[k - 1] + n;

    p[5] = b;

    vecaddsub (t[0], p[5], p[1], p[4], n);
    vecaddsub (t[1], p[4], p[3], p[2], n);
    vecaddsub (p[1], p[2], t[0], t[1], n);
    vecadd (p[0], p[0], p[1], n);
    vecadd (p[3], p[4], p[5], n);

    k = s[0];

    if (k == n)
    {
        for (ma = j = 0; j < 6; j++, ma += m)
            wftas[k] (p[j], w + ma);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        for (ma = j = 0; j < 6; j++, ma += m)
            wftan[k] (p[j], w + ma, nn, s, nm, ms);
    }

    vecadd (p[1], p[1], p[0], n);
    vecaddsub (p[1], p[2], p[1], p[2], n);
    vecsub (p[4], p[3], p[4], n);
    vecadd (p[3], p[3], p[5], n);

    vecaddsub (p[1], p[4], p[1], p[4], n);
    vecaddsub (p[2], p[3], p[2], p[3], n);
}

void wftan7 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint *p[9], *t[7] = {t0, t1, t2, t3, t4, t5, t6};

    for (p[0] = x, k = 1; k < 7; k++)
        p[k] = p[k - 1] + n;

    p[8] = (p[7] = b) + n;

    vecaddsub (t[0], p[6], p[1], p[6], n);
    vecaddsub (p[2], p[5], p[2], p[5], n);
    vecaddsub (p[4], p[3], p[4], p[3], n);
    vecadd (p[1], t[0], p[2], n);
    vecadd (p[1], p[1], p[4], n);
    vecadd (p[0], p[0], p[1], n);

    vecsub (p[7], p[3], p[5], n);
    vecsub (p[8], p[5], p[6], n);
    vecaddsub (t[1], p[6], p[6], p[3], n);
    vecadd (p[5], p[5], t[1], n);

    vecsub (p[3], p[4], p[2], n);
    vecsub (t[1], p[2], t[0], n);
    vecsub (p[2], t[0], p[4], n);
    veccopy (p[4], t[1], n);

    k = s[0];

    if (k == n)
    {
        for (ma = j = 0; j < 9; j++, ma += m)
            wftas[k] (p[j], w + ma);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        for (ma = j = 0; j < 9; j++, ma += m)
            wftan[k] (p[j], w + ma, nn, s, nm, ms);
    }

    vecadd (t[0], p[0], p[1], n);
    vecaddsub (p[1], p[2], t[0], p[2], n);
    vecadd (t[1], p[1], p[3], n);
    vecsub (t[2], p[2], p[4], n);
    vecsub (t[3], t[0], p[3], n);
    vecadd (t[3], t[3], p[4], n);

    vecaddsub (t[4], t[5], p[5], p[6], n);
    vecadd (t[4], t[4], p[7], n);
    vecsub (t[5], t[5], p[8], n);
    vecsub (t[6], p[5], p[7], n);
    vecadd (t[6], t[6], p[8], n);

    vecaddsub (p[1], p[6], t[1], t[4], n);
    vecaddsub (p[2], p[5], t[2], t[5], n);
    vecaddsub (p[4], p[3], t[3], t[6], n);
}

void wftan8 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint *p[8], *t[8] = {t0, t1, t2, t3, t4, t5, t6, t7};

    for (p[0] = x, k = 1; k < 8; k++)
        p[k] = p[k - 1] + n;

    vecaddsub (t[1], p[6], p[2], p[6], n);
    vecaddsub (t[2], t[3], p[1], p[5], n);
    vecaddsub (t[4], t[5], p[3], p[7], n);
    vecaddsub (t[0], p[3], p[0], p[4], n);
    vecaddsub (t[6], p[2], t[0], t[1], n);
    vecaddsub (t[7], p[5], t[2], t[4], n);

    vecaddsub (p[0], p[1], t[6], t[7], n);
    vecaddsub (p[7], p[4], t[3], t[5], n);

    k = s[0];

    if (k == n)
    {
        for (ma = j = 0; j < 8; j++, ma += m)
            wftas[k] (p[j], w + ma);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        for (ma = j = 0; j < 8; j++, ma += m)
            wftan[k] (p[j], w + ma, nn, s, nm, ms);
    }

    vecaddsub (t[0], t[1], p[3], p[4], n);
    vecaddsub (t[2], t[3], p[6], p[7], n);

    veccopy (p[4], p[1], n);
    vecaddsub (p[1], p[7], t[0], t[2], n);
    vecaddsub (p[2], p[6], p[2], p[5], n);
    vecaddsub (p[5], p[3], t[1], t[3], n);
}

void wftan9 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint *p[11], *t[8] = {t0, t1, t2, t3, t4, t5, t6, t7};

    for (p[0] = x, k = 1; k < 9; k++)
        p[k] = p[k - 1] + n;

    p[10] = (p[9] = b) + n;

    vecaddsub (t[0], t[7], p[8], p[1], n);
    vecaddsub (t[1], t[6], p[7], p[2], n);
    vecaddsub (p[2], p[6], p[6], p[3], n);
    vecaddsub (t[3], t[4], p[5], p[4], n);

    vecaddsub (p[1], p[3], t[0], t[3], n);
    vecadd (p[1], p[1], t[1], n);
    vecsub (p[3], t[0], t[3], n);
    vecsub (p[4], t[1], t[3], n);
    vecsub (p[5], t[0], t[1], n);

    vecaddsub (p[7], p[8], t[7], t[4], n);
    vecsub (p[7], p[7], t[6], n);
    vecsub (p[8], t[7], t[4], n);
    vecadd (p[9], t[4], t[6], n);
    vecadd (p[10], t[7], t[6], n);

    k = s[0];

    if (k == n)
    {
        for (ma = j = 0; j < 11; j++, ma += m)
            wftas[k] (p[j], w + ma);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        for (ma = j = 0; j < 11; j++, ma += m)
            wftan[k] (p[j], w + ma, nn, s, nm, ms);
    }

    vecaddsub (t[0], t[2], p[0], p[2], n);
    vecadd (t[0], t[0], p[2], n);

    vecadd (t[1], p[3], p[4], n);
    vecsub (t[3], p[3], p[5], n);
    vecadd (t[6], p[4], p[5], n);
    vecsub (t[4], p[8], p[10], n);
    vecadd (t[5], p[9], p[8], n);
    vecadd (t[7], p[9], p[10], n);

    vecaddsub (p[0], t[0], t[0], p[1], n);
    vecadd (p[0], p[0], p[1], n);
    vecsub (p[2], t[2], t[1], n);
    vecadd (p[3], t[2], t[6], n);
    vecadd (p[4], t[2], t[3], n);
    vecadd (p[5], t[4], p[6], n);
    vecadd (t[6], t[7], p[6], n);
    vecsub (t[7], t[5], p[6], n);

    vecaddsub (p[8], p[1], p[3], t[6], n);
    vecaddsub (p[6], p[3], t[0], p[7], n);
    vecaddsub (p[7], p[2], p[2], t[7], n);
    vecaddsub (p[5], p[4], p[4], p[5], n);
}

void wftan16 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint *p[18], *t[22] = {t0, t1, t2, t3, t4, t5, t6, t7,
                             t8, t9, t10, t11, t12, t13, t14,
                             t15, t16, t17, t18, t19, t20, t21};

    for (p[0] = x, k = 1; k < 16; k++)
        p[k] = p[k - 1] + n;

    p[17] = (p[16] = b) + n;

    vecaddsub (t[1], p[12], p[4], p[12], n);
    vecaddsub (t[0], p[4], p[0], p[8], n);
    vecaddsub (t[2], t[3], p[2], p[10], n);
    vecaddsub (t[4], t[5], p[6], p[14], n);
    vecaddsub (t[6], t[7], p[1], p[9], n);
    vecaddsub (t[8], t[9], p[3], p[11], n);
    vecaddsub (t[10], t[11], p[5], p[13], n);
    vecaddsub (t[12], t[13], p[7], p[15], n);

    vecaddsub (t[14], p[3], t[0], t[1], n);
    vecaddsub (t[15], p[11], t[2], t[4], n);
    vecaddsub (t[17], t[18], t[6], t[10], n);
    vecaddsub (t[19], t[20], t[8], t[12], n);
    vecaddsub (p[16], p[8], t[7], t[13], n);
    vecaddsub (p[17], p[9], t[11], t[9], n);

    vecaddsub (t[16], p[2], t[14], t[15], n);
    vecaddsub (t[21], p[10], t[17], t[19], n);

    vecaddsub (p[0], p[1], t[16], t[21], n);
    vecaddsub (p[13], p[5], t[18], t[20], n);
    vecaddsub (p[14], p[6], t[3], t[5], n);

    vecadd (p[7], p[8], p[9], n);
    vecadd (p[15], p[16], p[17], n);

    k = s[0];

    if (k == n)
    {
        for (ma = j = 0; j < 18; j++, ma += m)
            wftas[k] (p[j], w + ma);
    }
    else
    {
        nn = n / k;
        nm = m / ms[0];
        s++;
        ms++;
        for (ma = j = 0; j < 18; j++, ma += m)
            wftan[k] (p[j], w + ma, nn, s, nm, ms);
    }

    vecaddsub (t[0], t[1], p[3], p[5], n);
    vecaddsub (t[2], t[3], p[13], p[11], n);
    vecaddsub (t[4], t[5], p[4], p[6], n);

    vecsub (t[6], p[8], p[7], n);
    vecsub (t[7], p[9], p[7], n);

    vecaddsub (t[8], t[9], t[4], t[6], n);
    vecaddsub (t[10], t[11], t[5], t[7], n);
    vecaddsub (t[12], t[13], p[12], p[14], n);

    vecadd (t[14], p[15], p[16], n);
    vecsub (t[15], p[15], p[17], n);

    vecaddsub (t[16], t[17], t[12], t[14], n);
    vecaddsub (t[18], t[19], t[13], t[15], n);

    veccopy (p[8], p[1], n);

    vecaddsub (p[1], p[15], t[8], t[16], n);
    vecaddsub (p[4], p[12], p[2], p[10], n);
    vecaddsub (p[2], p[14], t[0], t[2], n);
    vecaddsub (p[13], p[3], t[11], t[19], n);
    vecaddsub (p[5], p[11], t[10], t[18], n);
    vecaddsub (p[6], p[10], t[1], t[3], n);
    vecaddsub (p[9], p[7], t[9], t[17], n);
}

void (*wftan[]) (modint x[], modint w[], int n, int s[], int m, int ms[]) =
{0, 0, wftan2, wftan3, wftan4, wftan5, 0, wftan7, wftan8, wftan9, 0, 0, 0, 0, 0, 0, wftan16};

void wfta (modint x[], modint *w[], int nn)
{
    int t, k;

    for (t = 0; t < NALGS; t++)
        if (n[t] == nn) break;

    if (t == NALGS) return;

    k = s[t][0];

    if (k == nn)
        wftas[k] (x, w[t]);
    else
    {
        permute (x, p[t], nn);
        wftan[k] (x, w[t], nn / k, s[t] + 1, m[t] / ms[t][0], ms[t] + 1);
        unpermute (x, ip[t], nn);
    }
}