wfta.cpp

任意精度计算的实现

        p[2][k] = t0 - t2;
        p[3][k] = t1 - t3;
    }

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

    for (k = 0; k < n; k++)
    {
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        p[1][k] = t2 + t3;
        p[3][k] = t2 - t3;
        p[2][k] = t1;
    }
}

void wftan5 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint t0, t1, t2, t3, t4, t5;
    modint r0, r1, r2, r3, r4;
    modint *p[6], b[B5];

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

    p[5] = b;

    for (k = 0; k < n; k++)
    {
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        r0 = t1 + t4;
        r2 = t1 - t4;
        r1 = t3 + t2;
        r3 = t3 - t2;
        r4 = r0 + r1;
        p[0][k] += r4;
        p[1][k] = r4;
        p[2][k] = r0 - r1;
        p[3][k] = r2 + r3;
        p[4][k] = r3;
        p[5][k] = r2;
    }

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

    for (k = 0; k < n; k++)
    {
        t0 = p[0][k];
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        r0 = t0 + t1;
        r1 = r0 + t2;
        r3 = r0 - t2;
        r2 = t3 - t4;
        r4 = t3 + t5;
        p[1][k] = r1 + r2;
        p[4][k] = r1 - r2;
        p[2][k] = r3 + r4;
        p[3][k] = r3 - r4;
    }
}

void wftan7 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint t0, t1, t2, t3, t4, t5, t6, t7, t8;
    modint r0, r1, r2, r3, r4, r5, r6;
    modint *p[9], b[2 * B7];

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

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

    for (k = 0; k < n; k++)
    {
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        r0 = t1 + t6;
        r4 = t1 - t6;
        r1 = t2 + t5;
        r5 = t2 - t5;
        r2 = t4 + t3;
        r6 = t4 - t3;
        r3 = r0 + r1 + r2;
        p[0][k] += r3;
        p[1][k] = r3;
        p[2][k] = r0 - r2;
        p[3][k] = r2 - r1;
        p[4][k] = r1 - r0;
        p[5][k] = r4 + r5 + r6;
        p[6][k] = r4 - r6;
        p[7][k] = r6 - r5;
        p[8][k] = r5 - r4;
    }

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

    for (k = 0; k < n; k++)
    {
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        t8 = p[8][k];
        r0 = p[0][k] + p[1][k];
        r1 = r0 + t2 + t3;
        r2 = r0 - t2 - t4;
        r3 = r0 - t3 + t4;
        r4 = t5 + t6 + t7;
        r5 = t5 - t6 - t8;
        r6 = t5 - t7 + t8;
        p[1][k] = r1 + r4;
        p[6][k] = r1 - r4;
        p[2][k] = r2 + r5;
        p[5][k] = r2 - r5;
        p[4][k] = r3 + r6;
        p[3][k] = r3 - r6;
    }
}

void wftan8 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint t0, t1, t2, t3, t4, t5, t6, t7;
    modint r0, r1, r2, r3, r4, r5, r6, r7;
    modint *p[8];

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

    for (k = 0; k < n; k++)
    {
        t0 = p[0][k];
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        r0 = t0 + t4;
        r1 = t2 + t6;
        r2 = t1 + t5;
        r3 = t1 - t5;
        r4 = t3 + t7;
        r5 = t3 - t7;
        r6 = r0 + r1;
        r7 = r2 + r4;
        p[0][k] = r6 + r7;
        p[1][k] = r6 - r7;
        p[2][k] = r0 - r1;
        p[3][k] = t0 - t4;
        p[5][k] = r2 - r4;
        p[6][k] = t2 - t6;
        p[7][k] = r3 + r5;
        p[4][k] = r3 - r5;
    }

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

    for (k = 0; k < n; k++)
    {
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        r0 = t3 + t4;
        r1 = t3 - t4;
        r2 = t6 + t7;
        r3 = t6 - t7;
        p[4][k] = t1;
        p[1][k] = r0 + r2;
        p[7][k] = r0 - r2;
        p[2][k] = t2 + t5;
        p[6][k] = t2 - t5;
        p[5][k] = r1 + r3;
        p[3][k] = r1 - r3;
    }
}

void wftan9 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    modint r0, r1, r2, r3, r4, r5, r6, r7;
    modint q0, q1, q2, q3, q4, q5, q6, q7;
    modint *p[11], b[2 * B9];

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

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

    for (k = 0; k < n; k++)
    {
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        t8 = p[8][k];
        r0 = t8 + t1;
        r7 = t8 - t1;
        r1 = t7 + t2;
        r6 = t7 - t2;
        r2 = t6 + t3;
        r5 = t6 - t3;
        r3 = t5 + t4;
        r4 = t5 - t4;
        p[2][k] = r2;
        p[1][k] = r0 + r1 + r3;
        p[3][k] = r0 - r3;
        p[4][k] = r1 - r3;
        p[5][k] = r0 - r1;
        p[6][k] = r5;
        p[7][k] = r7 - r6 + r4;
        p[8][k] = r7 - r4;
        p[9][k] = r4 + r6;
        p[10][k] = r7 + r6;
    }

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

    for (k = 0; k < n; k++)
    {
        t0 = p[0][k];
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        t8 = p[8][k];
        t9 = p[9][k];
        t10 = p[10][k];
        r0 = t2 + t2 + t0;
        r2 = t0 - t2;
        r1 = t3 + t4;
        r3 = t3 - t5;
        r6 = t4 + t5;
        r4 = t8 - t10;
        r5 = t9 + t8;
        r7 = t9 + t10;
        q0 = r0 + t1 + t1;
        q1 = r0 - t1;
        q2 = r2 - r1;
        q4 = r2 + r3;
        q3 = r6 + r2;
        q5 = r4 + t6;
        q6 = r7 + t6;
        q7 = r5 - t6;
        p[0][k] = q0;
        p[8][k] = q3 + q6;
        p[1][k] = q3 - q6;
        p[7][k] = q2 + q7;
        p[2][k] = q2 - q7;
        p[6][k] = q1 + t7;
        p[3][k] = q1 - t7;
        p[5][k] = q4 + q5;
        p[4][k] = q4 - q5;
    }
}

void wftan16 (modint x[], modint w[], int n, int s[], int m, int ms[])
{
    int k, j, nn, nm, ma;
    modint t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17;
    modint r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12,
           r13, r14, r15, r16, r17, r18, r19, r20, r21, r22, r23, r24, r25;
    modint *p[18], b[2 * B16];

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

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

    for (k = 0; k < n; k++)
    {
        t0 = p[0][k];
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        t8 = p[8][k];
        t9 = p[9][k];
        t10 = p[10][k];
        t11 = p[11][k];
        t12 = p[12][k];
        t13 = p[13][k];
        t14 = p[14][k];
        t15 = p[15][k];
        r0 = t0 + t8;
        r1 = t4 + t12;
        r2 = t2 + t10;
        r3 = t2 - t10;
        r4 = t6 + t14;
        r5 = t6 - t14;
        r6 = t1 + t9;
        r7 = t1 - t9;
        r8 = t3 + t11;
        r9 = t3 - t11;
        r10 = t5 + t13;
        r11 = t5 - t13;
        r12 = t7 + t15;
        r13 = t7 - t15;
        r14 = r0 + r1;
        r15 = r2 + r4;
        r16 = r14 + r15;
        r17 = r6 + r10;
        r18 = r6 - r10;
        r19 = r8 + r12;
        r20 = r8 - r12;
        r21 = r17 + r19;
        r22 = r7 + r13;
        r23 = r7 - r13;
        r24 = r11 + r9;
        r25 = r11 - r9;

        p[0][k] = r16 + r21;
        p[1][k] = r16 - r21;
        p[2][k] = r14 - r15;
        p[3][k] = r0 - r1;
        p[4][k] = t0 - t8;
        p[13][k] = r18 + r20;
        p[5][k] = r18 - r20;
        p[14][k] = r3 + r5;
        p[6][k] = r3 - r5;
        p[7][k] = r23 + r25;
        p[8][k] = r23;
        p[9][k] = r25;
        p[10][k] = r17 - r19;
        p[11][k] = r2 - r4;
        p[12][k] = t4 - t12;
        p[15][k] = r22 + r24;
        p[16][k] = r22;
        p[17][k] = r24;
    }

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

    for (k = 0; k < n; k++)
    {
        t1 = p[1][k];
        t2 = p[2][k];
        t3 = p[3][k];
        t4 = p[4][k];
        t5 = p[5][k];
        t6 = p[6][k];
        t7 = p[7][k];
        t8 = p[8][k];
        t9 = p[9][k];
        t10 = p[10][k];
        t11 = p[11][k];
        t12 = p[12][k];
        t13 = p[13][k];
        t14 = p[14][k];
        t15 = p[15][k];
        t16 = p[16][k];
        t17 = p[17][k];

        r0 = t3 + t5;
        r1 = t3 - t5;
        r2 = t13 + t11;
        r3 = t13 - t11;
        r4 = t4 + t6;
        r5 = t4 - t6;
        r6 = t8 - t7;
        r7 = t9 - t7;
        r8 = r4 + r6;
        r9 = r4 - r6;
        r10 = r5 + r7;
        r11 = r5 - r7;
        r12 = t12 + t14;
        r13 = t12 - t14;
        r14 = t15 + t16;
        r15 = t15 - t17;
        r16 = r12 + r14;
        r17 = r12 - r14;
        r18 = r13 + r15;
        r19 = r13 - r15;

        p[8][k] = t1;
        p[1][k] = r8 + r16;
        p[15][k] = r8 - r16;
        p[2][k] = r0 + r2;
        p[14][k] = r0 - r2;
        p[13][k] = r11 + r19;
        p[3][k] = r11 - r19;
        p[4][k] = t2 + t10;
        p[12][k] = t2 - t10;
        p[5][k] = r10 + r18;
        p[11][k] = r10 - r18;
        p[6][k] = r1 + r3;
        p[10][k] = r1 - r3;
        p[9][k] = r9 + r17;
        p[7][k] = r9 - r17;
    }
}

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