题解 P6271 【[湖北省队互测2014]一个人的数论】
题意
\(\sum_{i=1}^n[\gcd(i,n)=1]i^d\)
由于 \(n\) 可能很大，给出\(n\)的质因数分解式。
\(n=\prod_{i=1}^{w} p_{i}^{\alpha_{i}}\)
\(1≤w≤10^3,2\leq p_i \leq 10^9,1 \leq a_i \leq 10^9\)
题解
定位：莫反板子题莫反神仙题
\(\sum_{i=1}^n[\gcd(i,n)=1]i^d\)
\(\sum_{i=1}i^d\sum_{x|\gcd(i,n)}\mu(x)\)
\(\sum_{x|n}\mu(x)\sum_{i=1}^{\frac n x}(ix)^d\)
\(\sum_{x|n}\mu(x)x^d\sum_{i=1}^{\frac n x}i^d\)
根据经验我们得到，\(f(x)=\sum_{i=1}^xi^d\)是一个关于\(x\)的\(d+1\)次的多项式。
因此，设，\(f(x)=a_0+a_1x+a_2x^2+a_3x^3+\ldots+a_{d+1}x^{d+1}\)
\(\sum_{x|n}\mu(x)x^d f(\frac n x)\)
\(\sum_{x|n}\mu(x)x^d \sum_{i=0}^{d+1}a_i(\frac n x)^i\)
\(\sum_{i=0}^{d+1}a_i\sum_{x|n}\mu(x)x^d (\frac n x)^i\)
设\(g(n)=\sum_{x|n}\mu(x)x^d (\frac n x)^i\)，为\(\mu\times Id_d\)与\(Id_i\)卷起来得到的，前者与后者显然都为积性函数，那么\(g\)也必为积性函数。
\(\sum_{i=0}^{m+1}a_ig(n)\)
那么根据积性函数的性质，我们知道：
\(g(n)=\prod_{j} g(p_j^{\alpha_j})\)
再考虑当\(n=p^{\alpha}\)时，\(\mu\)当且仅当\(x=1\ \text{or}\ p\)时值不为零，因此：
\(g(p^{\alpha})=\mu(1)1^d (\frac {p^\alpha} 1)^i+\mu(p)p^d (\frac {p^\alpha} p)^i=p^{i\alpha}-p^{d+(\alpha-1)\times i}\)
再带回去，\(g(n)=\prod_{j} (p_j^{i\alpha_j}-p_j^{d+(\alpha_j-1)\times i})\)
再带回去，\(\sum_{i=0}^{d+1}a_i\prod_{j} (p_j^{i\alpha_j}-p_j^{d+(\alpha_j-1)\times i})\)
于是我们通过莫比乌斯函数，得出了一个与它无关的柿子。
\(a_i\)可以用高斯消元求出。
复杂度：\(\mathcal{O}(d^3+d\times w)\)可能还有快速幂的\(\log mod\)
代码
#include<bits/stdc++.h>
using namespace std;
template<class T,const int M>
struct Gauss{
//解n元方程组
int n;
T B[M][M+1];
T ans[M];
bool run(){
for(int i=0;i<n;i++){
int pivot=i;
for(int j=i;j<n;j++)
if(B[j][i]==B[pivot][i])
pivot=j;
for(int j=0;j<=n;j++)swap(B[i][j],B[pivot][j]);
if(B[i][i]==0)return false;
for(int j=i+1;j<=n;j++)B[i][j]/=B[i][i];
for(int j=0;j<n;j++)
if(i!=j)
for(int k=i+1;k<=n;k++)
B[j][k]-=B[j][i]*B[i][k];
}
for(int i=0;i<n;i++)ans[i]=B[i][n];
return true;
}
};
typedef long long ll;
const ll mod=1e9+7;
ll ksm(ll a,ll b){
ll res=1;
while(b){
if(b&1)res=res*a%mod;
a=a*a%mod;b>>=1;
}return res;
}
struct node{
ll x;
node(ll _=0){x=_;}
bool operator==(node b){return x==b.x;}
bool operator==(ll b){return x==b;}
node operator*(node b){return node(x*b.x%mod);}
node *operator/=(node b){x*=ksm(b.x,mod-2),x%=mod;return this;}
node *operator-=(node b){x-=b.x,x%=mod,x+=mod,x%=mod;;return this;}
node *operator+=(node b){x+=b.x,x%=mod,x+=mod,x%=mod;;return this;}
};
Gauss<node,200>G;
int d,w,p[1200],a[1200];
ll n=1;
signed main(){
cin>>d>>w;
G.n=d+2;
for(int i=0;i<G.n;i++){
for(int j=0;j<G.n;j++)
G.B[i][j].x=ksm(i+1,j);
for(int j=1;j<=i+1;j++)
(G.B[i][G.n].x+=ksm(j,d))%=mod;
}
G.run();
for(int i=1;i<=w;i++)
cin>>p[i]>>a[i],(n*=ksm(p[i],a[i]))%=mod;
ll ans=0,res;
for(int i=0;i<G.n;i++){
res=G.ans[i].x;
for(int j=1;j<=w;j++)
res*=ksm(p[j],1ll*i*a[j])-ksm(p[j],d+1ll*(a[j]-1)*i),res%=mod;
ans+=res;ans%=mod;
}
cout<<(ans+mod)%mod;
}
