gxx_slide之Evaluation

hdu4656:
链接：http://acm.hdu.edu.cn/showproblem.php?pid=4656

题意：$f(x)=\sum_{i=0}^{n-1} A_ix^i$给出$A_0,A_1,A_2,...A_{n-1}$,对于所有的$0<=k<n$ ,求解$f(B*C^{2k}+D)$,$n<=10^5$
分析：首先，尝试把$B*C^{2k}+D$直接带入化简
$f(x)=\sum_{i=0}^{n-1} A_i*(B*C^{2k}+d)^i$
用二项式定理展开
$f(x)=\sum_{i=0}^{n-1} A_i*\sum_{j=0}^{i}C_i^j(B*C^{2k})^j*d^{i-j}$
由于内层求和当中包含k,我们可以考虑交换求和顺序，即本来是一行一行求和，现在改成一列一列求和
$f(x)=\sum_{j=0}^{n-1} \sum_{i=j}^{n-1}A_i*C_i^j(B*C^{2k})^j*d^{i-j}$
先把组合数展开
$f(x)=\sum_{j=0}^{n-1} \sum_{i=j}^{n-1}A_i*i!*(j!)^{-1}*{(i-j)!}^{-1}(B*C^{2k})^j*d^{i-j}$
由于内层是只关于i的求和，我们可以把所有只包含j的项都提出去
$f(x)=\sum_{j=0}^{n-1}(j!)^{-1}*(B*C^{2k})^j \sum_{i=j}^{n-1}A_i*i!*{(i-j)!}^{-1}*d^{i-j}$
于是我们惊喜地发现内层$\sum$是一个只关于$i$和$i-j$的表达式，可以用傅立叶算出，记第j项的值是$P_j$
那么
$f(x)=\sum_{j=0}^{n-1}(j!)^{-1}*(B*C^{2k})^j *P_j$
按照刚才的思路，我们仍然想办法把这个式子变成一个卷积，由于要我们求的是第k项，而且表达式中存在$C^{2kj}$这样的项，我们可以想办法构造，使得右边的表达式是一个只关于j和k-j的函数于是
因为$2kj=k^2+j^2-(k-j)^2$
代入式子得
$f(x)=\sum_{j=0}^{n-1}(j!)^{-1}*B^j*C^{k^2+j^2-(k-j)^2} *P_j$
把$C^{k^2}$提出，就得到了我们想要的
$f(x)=C^{k^2}\sum_{j=0}^{n-1}(j!)^{-1}*B^j*C^{j^2}*P_j*C^{(k-j)^2}$
再用一次傅立叶就得到了我们要求的答案
本题的考验之处在于他的模数$M=10^6+3$并不是一个费马素数，翻看了kuangbin的代码，发现他是通过另选两个比较大的费马素数M1和M2，由于结果 <1017 <script type="math/tex" id="MathJax-Element-12495"><10^{17}</script>,在M1*M2的范围内只有惟一解，因此可以通过中国剩余定理来计算，由于当中有乘法操作，而M1*M2实在是太大，因此需要用快速加法替代乘法取模操作
这里我选用的两个素数是302776321,330301441,可以满足maxlen=262144,即可以满足数据量在$10^5$级别的题目，而且乘法不会溢出longlong，可以通过这道题目

总结：这道题目可以说是比较经典的题型，多项式展开后化简成卷积的形式来加速运算，同时，这道题模数的设置也让我复习了中国剩余定理（并不复杂，有时间可以写写关于这个定理的题目），按照这道题的做法，任何val<=$10^{18}$的fft的题目都可以用ntt来做（取两个素数），最后是我的代码，g++过了，c++tle了，不知道是什么原因

#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<vector>
using namespace std;
typedef long long Int;
const int M1=302776321,M2=330301441;
const int g1=17,g2=22;
const int M=1e6+3,Maxn=262145;
Int inv[100020];
Int fac[100020];
Int a[Maxn],b[Maxn],tp[Maxn],c[Maxn],d[Maxn];
int MM,gg;
int powmod(int x,int y,int mod)
{
    int ret=1,t=x;
    while(y){if(y&1)ret=ret*(Int)t%mod;t=t*(Int)t%mod;y>>=1;}
    return ret;
}
void rev(Int *a,int N)
{
    int i,j,k;
    for(i=1,j=N>>1;i<N-1;i++)
    {
        if(i<j)swap(a[i],a[j]);
        for(k=N>>1;j>=k;j-=k,k>>=1);j+=k;
    }
}
void dft(Int *a,int N,int flag=1)
{
    rev(a,N);
    for(int m=2;m<=N;m<<=1)
    {
        Int wm=powmod(gg,(MM-1)/m,MM);
        if(flag<0)wm=powmod(wm,MM-2,MM);
        for(int k=0;k<N;k+=m)
        {
            Int w=1;
            for(int j=k;j<k+(m>>1);j++,w=w*wm%MM)
            {
                Int u=a[j],v=a[j+(m>>1)]*w%MM;
                a[j]=(u+v)%MM;
                a[j+(m>>1)]=(u-v+MM)%MM;
            }
        }
    }
}
void mull(Int *a,Int *b,int N)
{
    dft(a,N);dft(b,N);
    for(int i=0;i<N;i++)a[i]=a[i]*b[i]%MM;
    dft(a,N,-1);
    int revn=powmod(N,MM-2,MM);
    for(int i=0;i<N;i++)a[i]=a[i]*revn%MM;
}
int op=0;
Int qmul(Int x,Int y,Int mod)
{
    Int ret=0,t=x;
    while(y){if(y&1)ret=(ret+t)%mod;t=(t+t)%mod;y>>=1;}
    return ret;
}
void mul(Int *a,Int *b,int N)
{
    MM=M1;
    gg=g1;
    for(int i=0;i<N;i++)c[i]=a[i],d[i]=b[i];
    mull(a,b,N);
    MM=M2;
    gg=g2;
    mull(c,d,N);
    int rev1=powmod(M2,M1-2,M1),rev2=powmod(M1,M2-2,M2);
    Int tt=M1*(Int)M2;
    for(int i=0;i<N;i++)a[i]=(qmul(a[i]*M2%tt,rev1,tt)+qmul(c[i]*M1%tt,rev2,tt))%tt%M;
}
int main()
{
    int n,B,C,D;
    inv[1]=inv[0]=fac[0]=fac[1]=1;
    for(int i=2;i<=100000;i++)inv[i]=inv[M%i]*(M-M/i)%M,fac[i]=fac[i-1]*i%M;
    for(int i=2;i<=100000;i++)inv[i]=inv[i-1]*inv[i]%M;
    while(scanf("%d%d%d%d",&n,&B,&C,&D)!=EOF)
    {
        int N=1;
        while(N<n)N<<=1;N<<=1;
        for(int i=0,curd=1;i<n;curd=curd*(Int)D%M,i++)
        {
            int x;scanf("%d",&x);
            a[i]=x*(Int)fac[i]%M;
            b[n-1-i]=curd*inv[i]%M;
        }
        for(int j=n;j<N;j++)a[j]=b[j]=0;
        mul(a,b,N);
        for(int i=0;i<n;i++)
        {
            tp[i]=powmod(C,(int)(i*(Int)i%(M-1)),M);
            a[i]=powmod(B,i,M)%M*tp[i]%M*inv[i]%M*a[i+n-1]%M;
        }
        for(int i=n;i<N;i++)a[i]=0;
        for(int i=0;i<n+n-1;i++)b[i]=powmod(tp[abs(i-(n-1))],M-2,M);
        for(int i=n+n-1;i<N;i++)b[i]=0;
        op=1;
        mul(a,b,N);
        for(int i=0;i<n;i++)printf("%d\n",(int)(a[i+n-1]*tp[i]%M));
    }
}