【BZOJ3601】一个人的数论 - 莫比乌斯反演+高斯消元
咕了inf天终于来补这题了……
题意：
题解：
常规操作，先套莫比乌斯反演大力推式子：
$$f_d(n)=\sum\limits_{i=1}^{n}i^d[gcd(i,n)=1]$$
$$=\sum\limits_{i=1}^{n}i^d\sum\limits_{k|gcd(i,n)}\mu(k)$$
$$=\sum\limits_{k|n}\sum\limits_{i=1}^{\frac{n}{k}}(ik)^d\mu(k)$$
$$=\sum\limits_{k|n}k^d\mu(k)\sum\limits_{i=1}^{\frac{n}{k}}i^d$$
设$S(n)=\sum\limits_{i=1}^{n}i^d$，则：
$$ans=\sum\limits_{k|n}k^d\mu(k)S(\frac{n}{k})$$
经过一些微小的数学分析我们可知$S(n)$是一个关于$n$的$d+1$次多项式；
证明详见【BZOJ3453】XLkxc那题……
于是可以暴力高斯消元或者拉格朗日插值来求出这个多项式的每一项系数；
设$S(n)=\sum\limits_{i=1}^{d+1}s_in^i$，则：
$$ans=\sum\limits_{k|n}k^d\mu(k)\sum\limits_{i=1}^{d+1}s_i(\frac{n}{k})^i$$
$$=\sum\limits_{i=1}^{d+1}s_i\sum\limits_{k|n}\mu(k)k^d(\frac{n}{k})^i$$
$$=\sum\limits_{i=1}^{d+1}s_i\sum\limits_{k|n}\mu(k)k^{d-i}n^i$$
设$g_i(n)=\sum\limits_{k|n}\mu(k)k^{d-i}n^i$，易得这是个积性函数；
题目中也很良心的给出的是$n$的质因数分解式，那么可以分开来考虑$n$的所有质数的整数次幂的因数（即$p_i^{\alpha_i}$）：
$$ans=\sum\limits_{i=1}^{d+1}s_i\prod\limits_{j=1}^{w}g_i(p_j^{\alpha_i})$$
只看其中的$g_i(p^\alpha)$一项：
$$g_i(p^\alpha)=\sum\limits_{j=0}^{\alpha}\mu(p^j)p^{j\times(d-i)}p^{\alpha i}$$
显然$j>1$时后面的都为0……因此只用考虑前两项即可：
$$g_i(p^\alpha)=p^{\alpha i}-p^{d-i}p^{\alpha i}$$
$$=p^{\alpha i}(1-p^{d-i})$$
于是就做完了。
代码：
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 #define mod 1000000007
10 using namespace std;
11 typedef long long ll;
12 typedef double db;
13 int d,w,ans=0,tmp,p[1001],a[1001],sq[105][105],s[101];
14 int fastpow(int x,int y){
15 int ret=1;
16 if(y<0)y+=mod-1;
17 for(;y;y>>=1,x=(ll)x*x%mod){
18 if(y&1)ret=(ll)ret*x%mod;
19 }
20 return ret;
21 }
22 void gauss(){
23 int sum=0,nw,tmp;
24 for(int i=0;i<=d+1;i++){
25 sum=(sum+fastpow(i+1,d))%mod;
26 sq[i][0]=1;
27 sq[i][d+2]=sum;
28 for(int j=1;j<=d+1;j++){
29 sq[i][j]=(ll)sq[i][j-1]*(i+1)%mod;
30 }
31 }
32 for(int i=0;i<=d+1;i++){
33 nw=i;
34 for(int j=i;j<=d+1;j++){
35 if(sq[i][j]){
36 nw=j;
37 break;
38 }
39 }
40 if(nw!=i){
41 for(int j=0;j<=d+2;j++){
42 swap(sq[i][j],sq[nw][j]);
43 }
44 }
45 for(int j=0;j<=d+1;j++){
46 if(i!=j&&sq[j][i]){
47 tmp=(ll)sq[j][i]*fastpow(sq[i][i],mod-2)%mod;
48 for(int k=0;k<=d+2;k++){
49 sq[j][k]=(sq[j][k]-(ll)tmp*sq[i][k]%mod+mod)%mod;
50 }
51 }
52 }
53 }
54 for(int i=0;i<=d+1;i++){
55 s[i]=(ll)sq[i][d+2]*fastpow(sq[i][i],mod-2)%mod;
56 }
57 }
58 int main(){
59 scanf("%d%d",&d,&w);
60 for(int i=1;i<=w;i++){
61 scanf("%d%d",&p[i],&a[i]);
62 }
63 gauss();
64 for(int i=1;i<=d+1;i++){
65 tmp=1;
66 for(int j=1;j<=w;j++){
67 tmp=(ll)tmp*fastpow(p[j],(ll)a[j]*i%(mod-1))%mod;
68 tmp=(ll)tmp*(mod-fastpow(p[j],d-i)+1)%mod;
69 }
70 ans=(ans+(ll)tmp*s[i]%mod)%mod;
71 }
72 printf("%d",ans);
73 return 0;
74 }
