从今天开始，蒟蒻_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≤500,m≤10\)；
70%的数据满足\(n≤5000\)；
100%的数据满足\(1≤n≤50000, 1≤m≤100, 1≤a_i≤m\)。 ## 题解：考虑如何使上式值最小
设给x序列加了一个数\(C\in \mathbb{Z}\),则答案为\((x_i-y_i+C)^2\)
拆掉上面的柿子..
\[(x_i-y_i+C)^2=x_i^2+y_i^2+C^2-2x_iy_i+2C(x_i-y_i)\] 柿子中除了\(-2x_iy_i\)都是定值，所以我们只需要最大化\(x_iy_i\) 这个非常像一个卷积..再看看数据范围，也是FFT的数据范围。。
对于这样的柿子我们通常的处理方法就是把一个序列翻转。这样就变成了卷积，可以\(O(n \log n)\)算出值了。现在我们要最大化所有可能的这样的值，所以考虑断环为链，把A复制一遍接在末尾。算出来的\(a[n+1..n+n]\)就是所有可能的答案。再考虑关于C的项，含有C的项总和为 \[nC^2-2C\sum_{i=1}^{n}x_iy_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≤1000\)。对于\(50\%\)的数据，\(n≤60000\)。对于\(100\%\)的数据，\(n≤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_jB_{j-i}\] 看到A和B的下标和不是定值，考虑翻转系数： \[G_i = \sum_{j=i+1}^nA^R_{n-j}B_{j-i}\] 这时下标和是定值了，要化成卷积还要把sigma的上下标换一下：
平移j的取值范围：
\[G_i=\sum_{j=1}^{n-i}A^R_{n-j+i}B_{j}\] 令\(B_0=0\)，此时j可以取到0，并将\(G\)翻转 \[G^R_i=\sum_{j=0}^{i}A^R_{i-j}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

多项式算法做题记录