从今天开始，蒟蒻_WA自动机也要学多项式辣！虽然我什么都不会，但是好在我们机房有~~HA A队队长~~LCH~~CTP_998244353~~巨佬啊！~~所以天天被巨佬爆踩~~

lg P3723 [AH/HNOI2017]礼物

题目大意

两个长度为n的环形整数数列，您可以进行两种操作： 1. 对一个数列的所有数加上一个值C 2. 将一个数列旋转一位
定义两个数列的差异值为
$\sum_{i = 1}^{n} (x_{i} - y_{i})^{2}$ 求差异值的最小值。

数据范围

30%的数据满足 $n \leq 500, m \leq 10$ ；
70%的数据满足 $n \leq 5000$ ；
100%的数据满足 $1 \leq n \leq 50000, 1 \leq m \leq 100, 1 \leq a_{i} \leq m$ 。 ## 题解：考虑如何使上式值最小
设给x序列加了一个数 $C \in Z$ ,则答案为 $(x_{i} - y_{i} + C)^{2}$
拆掉上面的柿子..
$(x_{i} - y_{i} + C)^{2} = x_{i}^{2} + y_{i}^{2} + C^{2} - 2 x_{i} y_{i} + 2 C (x_{i} - y_{i})$ 柿子中除了 $- 2 x_{i} y_{i}$ 都是定值，所以我们只需要最大化 $x_{i} y_{i}$ 这个非常像一个卷积..再看看数据范围，也是FFT的数据范围。。
对于这样的柿子我们通常的处理方法就是把一个序列翻转。这样就变成了卷积，可以 $O (n \log n)$ 算出值了。现在我们要最大化所有可能的这样的值，所以考虑断环为链，把A复制一遍接在末尾。算出来的 $a [n + 1. . n + n]$ 就是所有可能的答案。再考虑关于C的项，含有C的项总和为 $n C^{2} - 2 C \sum_{i = 1}^{n} x_{i} y_{i}$ 这是关于C的二次函数，用对称轴就好了。注意C是整数，所以算出对称轴后要四舍五入取接近的整数。 注意数学函数调用时一定要加std::，不然会出现奇怪的精度问题。。

我的代码（NTT实现

#include <cstdio>
#include <algorithm>
#include <cmath>

using std::reverse;

typedef long long ll;

const int P=1004535809;
const int g=3;
const int maxn=4e6+2000;

int A[maxn],B[maxn],limit=1,w[maxn],winv[maxn],rev[maxn];

inline void swap(int &a,int & b){a^=b^=a^=b;}

inline int qpow(int a,int b,int p)
{
    int ans=1%p;
    for (;b;b>>=1,a=(ll)a*a%p)
        if (b&1) ans=(ll)ans*a%p;
    return ans;
}

inline void prework(int n,int m)
{
    int l=0;
    while (limit<=(n+m+1)) limit<<=1,++l;
    w[0]=1;w[1]=qpow(g,(P-1)/limit,P),winv[0]=1,winv[1]=qpow(w[1],P-2,P);
    for (int i=2;i<limit;++i)
        w[i]=(ll)w[i-1]*w[1]%P,winv[i]=1ll*winv[i-1]*winv[1]%P;
    for (int i=1;i<limit;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
}

inline void NTT(int *A,int *w,int limit)
{
    for (int i=0;i<limit;++i)
        if (i<rev[i]) swap(A[i],A[rev[i]]);
    for (int mid=1;mid<limit;mid<<=1)
        for (int R=mid<<1,j=0;j<limit;j+=R)
            for (int k=0;k<mid;++k)
            {
                int x=A[j+k],y=(ll)A[j+k+mid]*w[limit/2/mid*k]%P;
                A[j+k]=(x+y)%P;A[j+mid+k]=(x-y+P)%P;
            }
}

inline void DFT(int *A){ NTT(A,w,limit); }

inline void IDFT(int *A)
{
    NTT(A,winv,limit);
    int inv=qpow(limit,P-2,P);
    for (int i=0;i<=limit;++i)
        A[i]=((ll)A[i]*inv)%P;
}

int main()
{
    // freopen("./4827/4827/10.in","r",stdin);
    int n,m;
    scanf("%d%d",&n,&m);
    for (int i=1;i<=n;++i)
        scanf("%d",A+i);
    for (int i=1;i<=n;++i)
        scanf("%d",B+i);
    ll k=0;
    for (int i=1;i<=n;++i)
        k+=A[i]-B[i];
    ll C=std::round((double)-k/n);
    ll tot=0;
    for (int i=1;i<=n;++i)
        tot+=(ll)A[i]*A[i]+(ll)B[i]*B[i]+2*C*(A[i]-B[i])+C*C;
    std::reverse(A+1,A+n+1);
    std::copy(A+1,A+n+1,A+n+1);
    prework(n,2*n);
    DFT(A);DFT(B);
    for (int i=0;i<limit;++i)
        A[i]=(ll)B[i]*A[i]%P;
    IDFT(A);
    ll ans=-0x3f3f3f3f;
    for (int i=1;i<=n;++i)
        if (A[i+n]>ans) ans=A[i+n];
    printf("%lld",tot-2*ans);
}

lg P3338 [ZJOI2014]力

题目大意：

给出 $n$ 个数 $q_{i}$ ，给出 $F_{j}$ 的定义如下：

$F_{j} = \sum_{i < j} \frac{q_{i} q_{j}}{(i - j)^{2}} - \sum_{i > j} \frac{q_{i} q_{j}}{(i - j)^{2}}$
令 $E_{i} = F_{i} / q_{i}$ ，求 $E_{i}$ .

数据范围：

对于 $30 %$ 的数据， $n \leq 1000$ 。对于 $50 %$ 的数据， $n \leq 60000$ 。对于 $100 %$ 的数据， $n \leq 100000$ ， $0 < q_{i} < 1000000000$ 。

题解：

首先，可以直接得到 $E_{i}$ 的表达式： $E_{i} = \sum_{j < i} \frac{q_{j}}{(i - j)^{2}} - \sum_{j > i} \frac{q_{j}}{(i - j)^{2}}$ 分开考虑每一项：两项的推导方式差不多，所以只用后面的那一项做栗子吧：我们观察到这很像一个卷积的形式，所以我们考虑化成卷积：
设 $G_{i} = \sum_{j > i} \frac{q_{j}}{(i - j)^{2}} = \sum_{j = i + 1}^{n} \frac{q_{j}}{(i - j)^{2}}$ 令 $A_{i} = q_{i}, B_{i} = \frac{1}{i^{2}}$
则上式珂变形为： $G_{i} = \sum_{j = i + 1}^{n} A_{j} B_{j - i}$ 看到A和B的下标和不是定值，考虑翻转系数： $G_{i} = \sum_{j = i + 1}^{n} A_{n - j}^{R} B_{j - i}$ 这时下标和是定值了，要化成卷积还要把sigma的上下标换一下：
平移j的取值范围：
$G_{i} = \sum_{j = 1}^{n - i} A_{n - j + i}^{R} B_{j}$ 令 $B_{0} = 0$ ，此时j可以取到0，并将 $G$ 翻转 $G_{i}^{R} = \sum_{j = 0}^{i} A_{i - j}^{R} B_{j}$ 然后做卷积就可以了！代码:

// luogu-judger-enable-o2
#pragma GCC optimize("-O3")
#pragma GCC optimize("-Ofast")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-ffast-math")
#include <cstdio>
#include <cmath>
#include <cctype>
#include <algorithm>

const double Pi=acos(-1.0);
const int maxn=4e5+100;

double q[maxn];
int limit=1,rev[maxn];

struct Complex
{
    double real,imag;
    Complex(double real,double imag):real(real),imag(imag){}
    Complex(){}
    Complex conj();
}w[maxn],winv[maxn],A[maxn],B[maxn],C[maxn];

inline Complex Complex::conj(){return Complex(real,-imag);}
inline Complex operator+(const Complex& a,const Complex& b){return Complex(a.real+b.real,a.imag+b.imag);}
inline Complex operator-(const Complex& a,const Complex& b){return Complex(a.real-b.real,a.imag-b.imag);}
inline Complex operator*(const Complex& a,const Complex& b){return Complex(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);}

template<typename T>
inline void swap(T& a,T& b){T t=a;a=b;b=t;}

inline void DFT(Complex* A,Complex* w,int limit)
{
    for (register int i=0;i<limit;++i)
        if (i<rev[i]) swap(A[i],A[rev[i]]);
    for (register int mid=1;mid<limit;mid<<=1)
        for (register int R=mid<<1,j=0;j<limit;j+=R)
            for (register int k=0;k<mid;++k)
            {
                Complex x=A[j+k],y=w[limit/mid/2*k]*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
}

inline void prework(int n)
{
    int l=0;
    while (limit<=(n<<1)+1) limit<<=1,++l;
    for (register int i=0;i<limit;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    for (register int i=0;i<limit;++i) 
        w[i]=Complex(cos(Pi*2/limit*i),sin(Pi*2/limit*i)),winv[i]=w[i].conj();
}

int main()
{
    int n;
    scanf("%d",&n);
    for (int i=1;i<=n;++i)
        scanf("%lf",q+i),A[i].real=q[i];
    prework(n);
    for (int i=1;i<=n;++i)
        B[i].real=1.0/i/i;
    DFT(A,w,limit);DFT(B,w,limit);
    static Complex G[maxn],H[maxn];
    for (int i=0;i<limit;++i)
        G[i]=A[i]*B[i];
    DFT(G,winv,limit);
    for (int i=0;i<limit;++i)
        G[i].real/=limit;
    std::reverse(q,q+n+1);
    for (int i=0;i<limit;++i)
        A[i]=Complex(q[i],0);
    // for (int i=0;i<=n;++i)
    //     B[i]=Complex(1.0/(i+1)/(i+1),0);
    DFT(A,w,limit);
    // DFT(B,w,limit);
    for (int i=0;i<limit;++i)
        H[i]=A[i]*B[i];
    DFT(H,winv,limit);
    for (int i=0;i<limit;++i)
        H[i].real/=limit;
    for (int i=1;i<=n;++i)
        printf("%.5lf\n",G[i].real-H[n-i].real);
    // while (~puts("\a"));
}

_WA自动机的blog

多项式算法做题记录