[cf1103E]Radix sum

类似于uoj272,即$B=10$的情况,然后有以下几个细节问题:

1.答案对$2^{58}$取模可以先使用自然溢出模$2^{64}$,最后对$2^{58}$取模即可

2.为了避免实数,令$omega=cosfrac{2pi}{10}+sinfrac{2pi}{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+nlog_{2}n)$(这里$n=10^{5}$,不与数字数量区分)

 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
",a[i].a[0]);
95     }
96 } 
View Code
原文地址:https://www.cnblogs.com/PYWBKTDA/p/14490882.html