【CF1103E】Radix sum
foreverlasting
2019-02-21 08:22:40
[同步发表在我的博客里哦](https://foreverlasting1202.github.io/)
多项式。
神题。
首先拿到题目,看懂题目后,你就会知道这题要怎么做。显然的卷积。但问题在于它给出了一个神仙模数$2^{58}$。然后我就不会了,因此当时就咕了这道题。不过后来[这位神仙](https://www.cnblogs.com/Mr-Spade/p/10390667.html)给出了一个中文题解,研究了一下,突然找到了做法。不过还是好神仙啊!
我们考虑$DFT$和$IDFT$我们要做什么。
$1$、找到单位根,这里是$\omega_{10}$。
$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:
```cpp
//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;
}
```