[luogu7092]计数题

由于$\mu(i)$，因此每一个素数最多存在1次，当$k=0$答案必然为0

根据莫比乌斯和欧拉函数的积性，答案与对素数的划分无关，仅与每一个素数是否出现有关，换言之枚举素数出现的集合$P'$，答案即为$\sum_{P'\subseteq P}(-1)^{|P'|}div(|P'|)\prod_{p\in P'}(p-1)$

（其中$div(n)$表示对$n$个数划分的方案数，当$k=1$时即$div(n)=1$）

令$f(x)=\sum_{i=0}^{\infty}(\sum_{|P'|=i}\prod_{p\in P'}(p-1))x^{i}$，答案即为$\sum_{i=0}^{\infty}(-1)^{i}div(i)f(x)[i]$

考虑求$f(x)$，当插入一个素数$p$，则$f(x)[i]=(p-1)f(x)[i-1]+f(x)[i]$，可以看作乘上$1+(p-1)x$这个多项式，重复此过程可得$f(x)=\prod_{p+1\in P}(1+px)$，分治fft即可

当$p=1$时$div(i)=1$，直接计算即可；当$p=2$时，考虑$div(i)$，即贝尔数，记作$B_{n}$

考虑其指数生成函数，即$f(x)=\sum_{i=0}^{\infty}\frac{B_{i}}{i!}x^{i}$

代入其递推式，即$f(x)=B_{0}+\sum_{i=1}^{\infty}\frac{\sum_{j=0}^{i-1}\frac{(i-1)!}{j!(i-j-1)!}B_{j}}{i!}x^{i}$

调换枚举顺序，即$f(x)=B_{0}+\sum_{j=0}^{\infty}\frac{B_{j}}{j!}\sum_{i=j+1}^{\infty}\frac{x^{i}}{i(i-j-1)!}$

对其求导，即$f'(x)=\sum_{j=0}^{\infty}\frac{B_{j}}{j!}\sum_{i=j}^{\infty}\frac{x^{i}}{(i-j)!}=\sum_{j=0}^{\infty}\frac{B_{j}x^{j}}{j!}\sum_{i=j}^{\infty}\frac{x^{i-j}}{(i-j)!}$

根据$e^{x}$泰勒展开的结果，后半部分即$e^{x}$，代入得$f'(x)=e^{x}\sum_{j=0}^{\infty}\frac{B_{j}x^{j}}{j!}=e^{x}f(x)$

考虑多项式$\frac{f'(x)}{f(x)}$，联系多项式ln的推导过程或对$\ln f(x)$求导，可得$\frac{f'(x)}{f(x)}=(\ln f(x))'=e^{x}$

对两边积分，即$\ln f(x)=e^{x}+C$，再同时exp，即$f(x)=e^{e^{x}+C}$，由于exp要求常数项为0，令$C=-1$即可，最终得到$f(x)=e^{e^{x}}-1$，多项式exp即可

由于素数个数为$\frac{n}{\ln n}$个，因此总复杂度为$o(n\log_{2}n)$，可以通过

  1 #include<bits/stdc++.h>

  2 using namespace std;

  3 #define N 1000005

  4 #define mod 998244353

  5 struct poly{

  6     vector<int>a;

  7 }a,b;

  8 int n,type,ans,p[N],vis[N];

  9 int ksm(int n,int m){

 10     if (!m)return 1;

 11     int s=ksm(n,m>>1);

 12     s=1LL*s*s%mod;

 13     if (m&1)s=1LL*s*n%mod;

 14     return s;

 15 }

 16 void ntt(poly &a,int n,int p){

 17     for(int i=0;i<(1<<n);i++){

 18         int rev=0;

 19         for(int j=0;j<n;j++)rev=rev*2+((i&(1<<j))>0);

 20         if (i<rev)swap(a.a[i],a.a[rev]);

 21     }

 22     for(int i=2;i<=(1<<n);i*=2){

 23         int s=ksm(3,(mod-1)/i);

 24         if (p)s=ksm(s,mod-2);

 25         for(int j=0;j<(1<<n);j+=i)

 26             for(int k=0,ss=1;k<(i>>1);k++,ss=1LL*ss*s%mod){

 27                 int x=a.a[j+k],y=1LL*ss*a.a[j+k+(i>>1)]%mod;

 28                 a.a[j+k]=(x+y)%mod;

 29                 a.a[j+k+(i>>1)]=(x+mod-y)%mod;

 30             }

 31     }

 32     if (p){

 33         int s=ksm((1<<n),mod-2);

 34         for(int i=0;i<(1<<n);i++)a.a[i]=1LL*a.a[i]*s%mod;

 35     }

 36 }

 37 poly dfs(int n,int l,int r){

 38     poly ans;

 39     if (n==1){

 40         ans.a.push_back(1);

 41         if ((!l)||(r>p[0]))ans.a.push_back(0);

 42         else ans.a.push_back(p[l]-1);

 43         ans.a.push_back(0);

 44         ans.a.push_back(0);

 45         return ans;

 46     }

 47     int mid=(l+r>>1);

 48     poly L=dfs(n-1,l,mid),R=dfs(n-1,mid+1,r);

 49     for(int i=0;i<(1<<n);i++)L.a.push_back(0);

 50     for(int i=0;i<(1<<n);i++)R.a.push_back(0);

 51     ntt(L,n+1,0);

 52     ntt(R,n+1,0);

 53     for(int i=0;i<(1<<n+1);i++)ans.a.push_back(1LL*L.a[i]*R.a[i]%mod);

 54     ntt(ans,n+1,1);

 55     return ans;

 56 }

 57 poly inv(int n,poly a){//返回为2^(n+1)次多项式

 58     poly ans;

 59     if (n==0){

 60         ans.a.push_back(ksm(a.a[0],mod-2));

 61         ans.a.push_back(0);

 62         return ans;

 63     }

 64     ans=inv(n-1,a);

 65     for(int i=0;i<(1<<n);i++)ans.a.push_back(0);

 66     poly b=a;

 67     for(int i=(1<<n);i<(1<<n+1);i++)b.a[i]=0;

 68     ntt(ans,n+1,0);

 69     ntt(b,n+1,0);

 70     for(int i=0;i<(1<<n+1);i++)ans.a[i]=1LL*ans.a[i]*(mod+2-1LL*ans.a[i]*b.a[i]%mod)%mod;

 71     ntt(ans,n+1,1);

 72     for(int i=(1<<n);i<(1<<n+1);i++)ans.a[i]=0;

 73     return ans;

 74 }

 75 poly ln(int n,poly a){//返回为2^(n+1)次多项式

 76     poly ans=inv(n,a);

 77     for(int i=0;i<(1<<n+1)-1;i++)a.a[i]=1LL*a.a[i+1]*(i+1)%mod;

 78     a.a[(1<<n+1)-1]=0;

 79     ntt(ans,n+1,0);

 80     ntt(a,n+1,0);

 81     for(int i=0;i<(1<<n+1);i++)ans.a[i]=1LL*ans.a[i]*a.a[i]%mod;

 82     ntt(ans,n+1,1);

 83     for(int i=(1<<n+1)-1;i;i--)ans.a[i]=1LL*ksm(i,mod-2)*ans.a[i-1]%mod;

 84     ans.a[0]=0;

 85     for(int i=(1<<n);i<(1<<n+1);i++)ans.a[i]=0;

 86     return ans;

 87 }

 88 poly exp(int n,poly a){//返回为2^(n+1)次多项式

 89     poly ans;

 90     if (!n){

 91         ans.a.push_back(1);

 92         ans.a.push_back(0);

 93         return ans;

 94     }

 95     ans=exp(n-1,a);

 96     for(int i=0;i<(1<<n);i++)ans.a.push_back(0);

 97     poly l=ln(n,ans);

 98     for(int i=0;i<(1<<n);i++)l.a[i]=(a.a[i]-l.a[i]+mod)%mod;

 99     for(int i=(1<<n);i<(1<<n+1);i++)l.a[i]=0;

100     l.a[0]++;

101     ntt(l,n+1,0);

102     ntt(ans,n+1,0);

103     for(int i=0;i<(1<<n+1);i++)ans.a[i]=1LL*ans.a[i]*l.a[i]%mod;

104     ntt(ans,n+1,1);

105     for(int i=(1<<n);i<(1<<n+1);i++)ans.a[i]=0;

106     return ans;

107 }

108 int main(){

109     scanf("%d%d",&n,&type);

110     if (!type){

111         printf("0");

112         return 0;

113     }

114     for(int i=2;i<=n;i++){

115         if (!vis[i]){

116             p[++p[0]]=i;

117             vis[i]=1;

118         }

119         for(int j=1;(j<=p[0])&&(i*p[j]<=n);j++){

120             vis[i*p[j]]=1;

121             if (i%p[j]==0)break;

122         }

123     }

124     a=dfs(18,0,(1<<17)-1);

125     if (type==1){

126         for(int i=1;i<=p[0];i++)

127             if (i&1)ans=(ans+mod-a.a[i])%mod;

128             else ans=(ans+a.a[i])%mod;

129         printf("%d",ans);

130         return 0;

131     }

132     b.a.push_back(0);

133     b.a.push_back(1);

134     for(int i=2;i<(1<<17);i++)b.a.push_back(1LL*b.a[i-1]*ksm(i,mod-2)%mod);

135     b=exp(17,b);

136     int fac=1;

137     for(int i=1;i<=p[0];i++){

138         fac=1LL*fac*i%mod;

139         int s=1LL*a.a[i]*b.a[i]%mod*fac%mod;

140         if (i&1)ans=(ans+mod-s)%mod;

141         else ans=(ans+s)%mod;

142     }

143     printf("%d",ans);

144 }

巴特西

[luogu7092]计数题

最新文章

热门文章