DoorKickers的博客

博客

多项式入门推式子好题

2020-04-18 16:05:01 By DoorKickers

题目大意

给定n个数q1,q2qn 定义 Fj=j1i=1qi×qj(ij)2ni=j+1qi×qj(ij)2

同时 Ei=Fiqi

你的任务是对于所有的1in, 求出 Ei


Solution

推下公式: Ej=Fjqj=j1i=1qi(ij)2ni=j+1qi(ij)2 然后如果要用FFT加速的话, 要把这个推成卷积式

卷积式就是Ck=k0A[i]B[ki]

继续推 Ej=ji=1qi(ij)2ni=jqi(ij)2f[i]=qi, g[i]=1i2

原式变成 Ej=ji=1f[i]g[ji]ni=jf[i]g[ij] 然后设f[0]=g[0]=0 Ej=ji=0f[i]g[ji]ni=jf[i]g[ij] 发现左边已经是一个卷积式子了, 继续推右边的.

将他展开, 发现 ni=jf[i]g[ij]=f[j]g[0]+f[j+1]g[1]++f[j+(nj)]g[nj]=nji=0f[j+i]g[i] 然后发现f还不太行, 翻转一下.

f[i]=f[ni] nji=0f[j+i]g[i]=nji=0f[n(j+i)]g[i]=nji=0f[nji]g[i] 然后设t=nj

原式 ti=0f[ti]g[i] 答案就是 Ej=ji=0f[i]g[ji]ti=0f[ti]g[i](t=nj) 然后设A(x)=ni=0f[i],B(x)=ni=0g[i],C(x)=ni=0f[i]

D(x)=A(x)B(x),E(x)=B(x)C(x)

然后 Ej=D(j)E(nj)

就做完了

代码:

#include <bits/stdc++.h>

using namespace std;

#define int long long

const long double PI = 3.14159265358978323846264338;
const int N = 1e6 + 10;

complex<double> omega[N], omegaInv[N];
int rev[N];
double f[N], g[N], ff[N], res1[N], res2[N];

inline void init(const int n)
{
    int k = 0;
    while ((1 << k) < n) k++;
    for (int i = 0; i < n; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
        omega[i] = complex<double>(cos(2 * PI / n * i), sin(2 * PI / n * i));
        omegaInv[i] = conj(omega[i]);
    }
    rev[0] = 0;
}

void fft(complex<double> *a, const int n, const complex<double> *omega)
{
    for (int i = 0; i < n; i++)
    {
        if (i > rev[i]) continue;
        swap(a[i], a[rev[i]]);
    }
    for (int l = 2; l <= n; l *= 2)
    {
        int m = l / 2;
        for (complex<double> *p = a; p != a + n; p += l)
            for (int i = 0; i < m; i++)
            {
                complex<double> t = omega[n / l * i] * p[m + i];
                p[m + i] = p[i] - t;
                p[i] += t;
            }
    }
}

inline void dft(complex<double> *a, const int n)
{
    fft(a, n, omega);
}

inline void idft(complex<double> *a, const int n)
{
    fft(a, n, omegaInv);
    for (int i = 0; i < n; i++)
    {
        a[i] /= n;
    }
}

inline void squared(int *a1, const int n1)
{
    int n = 1;
    while (n < n1 + n1 + 2) n *= 2;
    static complex<double> c1[N];
    for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
    init(n);
    dft(c1, n);
    for (int i = 0; i < n; i++)
        c1[i] *= c1[i];
    idft(c1, n);
    for (int i = 0; i <= n1 * 2; i++)
        a1[i] = static_cast<int>(floor(c1[i].real() + 0.5));
}

inline void multiply(double *a1, const int n1, double *a2, const int n2, double *res)
{
    int n = 1;
    while (n < n1 + n2 + 2) n *= 2;
    static complex<double> c1[N], c2[N];
    memset(c1, 0, sizeof(c1));
    memset(c2, 0, sizeof(c2));
    for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
    for (int i = 0; i <= n2; i++) c2[i].real(a2[i]);
    init(n);
    dft(c1, n);
    dft(c2, n);
    for (int i = 0; i < n; i++)
        c1[i] *= c2[i];
    idft(c1, n);
    // idft(c2, n);
    for (int i = 0; i <= n1 + n2; i++)
        res[i] = c1[i].real();
}

signed main()
{
    int n;
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> f[i];
        ff[n - i] = f[i];
        g[i] = (double)(1.0 / i / i);
    }
    // for (int i = 1; i <= n; i++)
    //     cin >> g[i];
    multiply(f, n, g, n, res1);
    // for (int i = 0; i <= n * 2; i++) 
    //     cout << fixed << setprecision(10) << res1[i] << ' ';
    // cout << endl;
    // exit(0);
    multiply(ff, n, g, n, res2);
    for (int i = 1; i <= n; i++)
        cout << fixed << setprecision(3) << res1[i] - res2[n - i] << endl;
    return 0;
}

评论

  • 2020-04-18 16:09:32
  • Reply

发表评论

可以用@mike来提到mike这个用户,mike会被高亮显示。如果你真的想打“@”这个字符,请用“@@”。