[cf1103E]Radix sum
类似于uoj272,即$B=10$的情况,然后有以下几个细节问题:
1.答案对$2^{58}$取模可以先使用自然溢出模$2^{64}$,最后对$2^{58}$取模即可
2.为了避免实数,令$\omega=\cos\frac{2\pi}{10}+\sin\frac{2\pi}{10}i$,初始每一个数必然是$\omega^{i}$,相乘也就是多项式乘法
根据$\omega^{10}=1$,可以将其幂次对10取模,即是一个9次多项式
又因为$\omega^{5}=-1$,因此$\omega^{i+5}=-\omega^{i}$,即可以将5-9次项都降为0-4次项
再根据$\omega^{4}-\omega^{3}+\omega^{2}-\omega^{1}+1=\frac{1+\omega^{5}}{1+\omega}=0$,可以将4次项降为0-3次项
接下来,考虑若最后1-3次项存在非0,由于最终结果是实数,因此这些次项的虚部带权和为0
其实部带权和一定可以用$\omega^{0}$来表示,也就是可以继续降幂,但事实上我们无法再找到可以降$\omega^{3}$的式子,因此最终必然只有$\omega^{0}$系数非0,即答案(严格地证明并不会证)
3.关于10在模$2^{64}$下不一定没有逆元,可以将10放在最外面除以$10^{5}$,之后由于2一定是可以直接除掉的(因为最终结果是整数),再求出5的逆元$inv$,乘上$inv^{5}$即可
时间复杂度为$o(4^{2}Bn+n\log_{2}n)$(这里$n=10^{5}$,不与数字数量区分)
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 100005 4 #define M 5 5 #define B 10 6 #define ll long long 7 int n,m,x,base[N][M]; 8 struct Complex{ 9 unsigned ll a[4]; 10 Complex(){ 11 memset(a,0,sizeof(a)); 12 } 13 Complex(int x){ 14 memset(a,0,sizeof(a)); 15 a[0]=x; 16 } 17 Complex(int x,int y){ 18 memset(a,0,sizeof(a)); 19 a[0]=x,a[1]=y; 20 } 21 Complex operator + (const Complex &k)const{ 22 Complex o; 23 for(int i=0;i<4;i++)o.a[i]=a[i]+k.a[i]; 24 return o; 25 } 26 Complex operator * (const Complex &k)const{ 27 Complex o; 28 for(int i=0;i<4;i++) 29 for(int j=0;j<4;j++){ 30 if (i+j<4)o.a[i+j]=o.a[i+j]+a[i]*k.a[j]; 31 if (i+j==4){ 32 unsigned ll s=a[i]*k.a[j]; 33 o.a[0]-=s,o.a[1]+=s,o.a[2]-=s,o.a[3]+=s; 34 } 35 if (i+j>4)o.a[i+j-5]=o.a[i+j-5]-a[i]*k.a[j]; 36 } 37 return o; 38 } 39 }inv,A[B][B],invA[B][B],a[N]; 40 Complex pow(Complex n,ll m){ 41 Complex s=n,ans=Complex(1); 42 while (m){ 43 if (m&1)ans=ans*s; 44 s=s*s; 45 m>>=1; 46 } 47 return ans; 48 } 49 void DFT(Complex *a){ 50 Complex aa[B]; 51 for(int i=0,s=1;i<M;i++,s*=B) 52 for(int j=0;j<n;j++) 53 if (!base[j][i]){ 54 for(int k=0;k<B;k++)aa[k]=Complex(); 55 for(int k=0;k<B;k++) 56 for(int l=0;l<B;l++)aa[k]=aa[k]+a[j+l*s]*A[l][k]; 57 for(int k=0;k<B;k++)a[j+k*s]=aa[k]; 58 } 59 } 60 void IDFT(Complex *a){ 61 Complex aa[B]; 62 for(int i=0,s=1;i<M;i++,s*=B) 63 for(int j=0;j<n;j++) 64 if (!base[j][i]){ 65 for(int k=0;k<B;k++)aa[k]=Complex(); 66 for(int k=0;k<B;k++) 67 for(int l=0;l<B;l++)aa[k]=aa[k]+a[j+l*s]*invA[l][k]; 68 for(int k=0;k<B;k++)a[j+k*s]=aa[k]; 69 } 70 } 71 int main(){ 72 n=1e5; 73 inv=pow(Complex(5),(1LL<<57)-1); 74 for(int i=0;i<n;i++){ 75 base[i][0]=i%B; 76 for(int j=1;j<M;j++)base[i][j]=base[i/B][j-1]; 77 } 78 for(int i=0;i<B;i++) 79 for(int j=0;j<B;j++){ 80 A[i][j]=pow(Complex(0,1),i*j); 81 invA[i][j]=pow(Complex(0,1),B*B-i*j)*inv; 82 } 83 scanf("%d",&m); 84 for(int i=1;i<=m;i++){ 85 scanf("%d",&x); 86 a[x]=a[x]+Complex(1); 87 } 88 DFT(a); 89 for(int i=0;i<n;i++)a[i]=pow(a[i],m); 90 IDFT(a); 91 for(int i=0;i<m;i++){ 92 a[i].a[0]>>=M; 93 a[i].a[0]&=((1LL<<58)-1); 94 printf("%llu\n",a[i].a[0]); 95 } 96 }