【BZOJ3601】一个人的数论（莫比乌斯反演+高斯消元）

莫比乌斯反演

$f_d(n)=\sum_{i=1}^ni^d[gcd(i,n)=1]$

$f_d(n)=\sum_{i=1}^ni^d\sum_{k|i,k|n}\mu(k)$

$f_d(n)=\sum_{k|n}\mu(k)\sum_{i=1}^{\frac nk}(ik)^d$

$f_d(n)=\sum_{k|n}\mu(k)k^d\sum_{i=1}^{\frac nk}i^d$

进一步推式子

$f_d(n)=\sum_{k|n}\mu(k)k^d\sum_{i=1}^{d+1}a_i(\frac nk)^i=\sum_{i=1}^{d+1}a_i\sum_{k|n}\mu(k)k^d(\frac nk)^i$

$f_d(n)=\sum_{i=1}^{d+1}a_ig_i(n)=\sum_{i=1}^{d+1}a_i\prod_{v=1}^wg_i(p_v^{\alpha_v})$

$g_i(p_v^{\alpha_v})=(p_v^{\alpha_v})^i-p_v^d(p_v^{\alpha_v-1})^i=p_v^{\alpha_v\times i}-p_v^{(\alpha_v-1)\times i+d}$

$f_d(n)=\sum_{i=1}^{d+1}a_i\prod_{v=1}^w(p_v^{\alpha_v\times i}-p_v^{(\alpha_v-1)\times i+d})$

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 1000
#define X 1000000007
#define swap(x,y) (x^=y^=x^=y)
#define Qinv(x) Qpow(x,X-2)
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
using namespace std;
int n,d,p[N+5],q[N+5];
I int Qpow(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
namespace Gauss//高斯消元
{
#define SZ 100
int a[SZ+5][SZ+5],v[SZ+5];
I void Find(CI x) {RI i=x;W(!a[i][x]) ++i;if(i^x) for(RI j=x;j<=n;++j) swap(a[x][j],a[i][j]);}
I void Solve(CI n)
{
RI i,j,k,t;for(i=1;i<=n;++i) for(Find(i),j=i+1;j<=n;++j)
{
t=X-1LL*a[j][i]*Qinv(a[i][i])%X,v[j]=(1LL*t*v[i]+v[j])%X;
for(k=i;k<=n;++k) a[j][k]=(1LL*t*a[i][k]+a[j][k])%X;
}
for(i=n;i;--i)
{
v[i]=1LL*v[i]*Qinv(a[i][i])%X;
for(j=i-1;j;--j) v[j]=(v[j]-1LL*v[i]*a[j][i]%X+X)%X;
}
}
}
int main()
{
using namespace Gauss;
RI i,j;for(scanf("%d%d",&d,&n),i=1;i<=n;++i) scanf("%d%d",p+i,q+i);
RI t=0;for(i=1;i<=d+1;++i) for(Inc(t,Qpow(i,d)),v[i]=t,j=1;j<=d+1;++j) a[i][j]=Qpow(i,j);//建立方程
RI ans=0,res;for(Solve(d+1),i=1;i<=d+1;++i)//解出系数，然后枚举每一个系数去求答案
{
for(res=j=1;j<=n;++j) res=1LL*res*//枚举每个质因子，统计函数值
(Qpow(p[j],1LL*q[j]*i%(X-1))-Qpow(p[j],(1LL*(q[j]-1)*i+d)%(X-1))+X)%X;
ans=(1LL*v[i]*res+ans)%X;//统计答案
}return printf("%d",ans),0;
}

posted @ 2020-02-03 13:14  TheLostWeak  阅读(...)  评论(...编辑  收藏