【CF1103E】Radix sum

foreverlasting

2019-02-21 08:22:40

题解

同步发表在我的博客里哦

多项式。

神题。

首先拿到题目，看懂题目后，你就会知道这题要怎么做。显然的卷积。但问题在于它给出了一个神仙模数2^{58}。然后我就不会了，因此当时就咕了这道题。不过后来这位神仙给出了一个中文题解，研究了一下，突然找到了做法。不过还是好神仙啊！

我们考虑DFT和IDFT我们要做什么。

$2$、找到$10^5$的逆元，这里是在$2^{58}$下的逆元。首先看向第二个问题，这里看起来比较好解决。首先$5$在$2^{58}$下是有逆元的，所以只用考虑$2^5$的逆元。这里就有一步神仙转换，即发现$x\ast 2^5$在$2^{64}$下的结果可以直接计算，假设结果为$y$。那我们就有$\frac {y}{2}=x\ast 2^4$在模$2^{63}$下。$\frac {y}{2^5}=x$在模$2^{59}$下。而模数是$2^{58}$，这没关系，再取个模就可以了。现在再去处理第一个问题。我们可以发现的是$\omega_{10}^5$在模$2^{64}$下是有意义的，且答案为$-1$。并且有$\omega_{10}^2=\omega_{5}$。于是得到了$\omega_{5}^5-1=0$。这有什么用呢？于是我们把模$x^5-1$的多项式环作为一种数据类型，以此来做乘法等操作。（这里我和$Mr-Spade$的思路有点不同）这里先讲一下怎么做乘法，再证明一下正确性。做乘法就是把两个最高次为$x^4$的多项式相乘，然后模$x^5-1$即可。于是$x^8$就变成$+x^3$，$x^7$就变成$+x^2$以此类推，用个数组记录一下，然后像$FFT$原理那样处理一下单位根就好了。最后考虑正确性。其实在模$x^5-1$下直接输出常数项是错误的。我们按照$DFT$的正确性，$a_1+a_2\ast x+a_3\ast x^2+a_4\ast x^3+a_5\ast x^4$应该是个整数，但事实上原方程$x^5-1=0$的根中$\omega_5\omega_5^2\omega_5^3\omega_5^4$和$1$都可能令这个多项式为整数。那怎么处理呢？试根显然太麻烦，我们考虑证明$a_2=a_3=a_4=a_5$。首先说明为什么要证明这个。因为如果这个成立，我们只要再模$x-1$就可以消除其他系数，只留下常数项。这个怎么证明呢。首先我们要知道$1$不是解。因为如果$1$是解，$\omega_{10}=\pm 1$。这个显然错误的。我们令$a_1+a_2\ast x+a_3\ast x^2+a_4\ast x^3+a_5\ast x^4=y$，则$y$为整数。我们假设有四个根分别为$x_1,x_2,x_3,x_4$，移项因式分解得到 $k\ast \prod_{i=1}^4 (x-x_i)=0

拆掉，就会得到

a_5=k

a_4=-k\ast \sum_{i=1}^4 x_i

a_3=k\ast \sum_{i=1}^4 \sum_{j=i+1}^4 x_i\ast x_j

a_2=-k\ast \sum_{i=1}^4 \sum_{j=i+1}^4 \sum_{l=j+1}^4 x_i\ast x_j\ast x_l

当你去将上述四个根带入检验的时候，根据x^5-1=0这条方程，你会得到a_5=a_4=a_3=a_2=k。因此就得知了正确性。

这道神仙题最终还是做完了。

code:

//2019.2.20 by ljz
#include<bits/stdc++.h>
using namespace std;
#define res register int
#define LL long long
#define inf 0x3f3f3f3f
#define eps 1e-10
inline int read() {
    res s=0;
    bool w=0;
    char ch=getchar();
    while(ch<'0'||ch>'9') {
        if(ch=='-')w=1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9')s=s*10+ch-'0',ch=getchar();
    return w?-s:s;
}
inline void _swap(res &x,res &y) {
    x^=y^=x^=y;
}
inline int _abs(const res &x) {
    return x>0?x:-x;
}
inline int _max(const res &x,const res &y) {
    return x>y?x:y;
}
inline int _min(const res &x,const res &y) {
    return x<y?x:y;
}
#define ull unsigned LL
const int pw[]={1,10,100,1000,10000,100000,1000000};
const ull INV5=14757395258967641293ull;
const ull kcz=1ull<<58;
namespace MAIN {
    int n;
    struct cp{
        ull a[5];
        cp() {a[0]=a[1]=a[2]=a[3]=a[4]=0;}
    }w10[10],a[100000+10];
    inline cp operator + (const cp &x,const cp &y){
        cp z;
        for(res i=0;i<=4;i++)z.a[i]=x.a[i]+y.a[i];
        return z;
    }
    inline cp operator - (const cp &x,const cp &y){
        cp z;
        for(res i=0;i<=4;i++)z.a[i]=x.a[i]-y.a[i];
        return z;
    }
    inline cp operator * (const cp &x,const ull &y){
        cp z;
        for(res i=0;i<=4;i++)z.a[i]=x.a[i]*y;
        return z;
    }
    inline cp operator * (const cp &x,const cp &y){
        ull b[9];
        for(res i=0;i<=8;i++)b[i]=0;
        for(res i=0;i<=4;i++)
            for(res j=0;j<=4;j++)b[i+j]+=x.a[i]*y.a[j];
        for(res i=3;~i;i--)b[i]+=b[i+5],b[i+5]=0;
        cp z;
        for(res i=0;i<=4;i++)z.a[i]=b[i];
        return z;
    }
    inline void pre(){
        w10[0].a[0]=1,w10[1].a[3]=-1;
        for(res i=2;i<10;i++)w10[i]=w10[i-1]*w10[1];
    }
    inline void FFT(cp *A,const res &type,const res &lim){
        for(res mid=0;mid<5;mid++){
            for(res i=0;i<lim;i++)
                if((i/pw[mid])%10==0){
                    cp b[10];
                    for(res j=0;j<10;j++)
                        for(res k=0;k<10;k++)
                            b[j]=b[j]+A[i+k*pw[mid]]*w10[(type*j*k%10+10)%10];
                    for(res j=0;j<10;j++)
                        A[i+j*pw[mid]]=b[j];
                }
        }
        if(type==-1){
            ull INV=INV5*INV5*INV5*INV5*INV5;
            for(res i=0;i<lim;i++)A[i]=A[i]*INV;
        }
    }
    inline cp qpow(cp x,res y){
        cp ret=w10[0];
        while(y){
            if(y&1)ret=ret*x;
            x=x*x,y>>=1;
        }
        return ret;
    }
    inline void MAIN(){
        pre();
        n=read();
        for(res i=1;i<=n;i++)a[read()].a[0]++;
        FFT(a,1,100000);
        for(res i=0;i<100000;i++)a[i]=qpow(a[i],n);
        FFT(a,-1,100000);
        for(res i=0;i<n;i++)printf("%llu\n",((a[i].a[0]-a[i].a[4])/32)%kcz);
    }
}
int main() {
    MAIN::MAIN();
    return 0;
}