[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}$，不与数字数量区分）

 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 }

巴特西

[cf1103E]Radix sum

最新文章

热门文章