题目
定义：
[egin{aligned}
f_d(n)=sum _{gcd(x,n)=1} x^d
end{aligned}
]
给定(d)和(n) ((n)以质因数分解的形式给出)，求(f_d(n) (mod 10^9+7))。
分析
本来是学莫比乌斯反演时候的题，拖到了高斯消元来做。
看到一个互质，即(gcd(x,n)=1)，马上就能想到莫比乌斯函数的那个性质：
[令 f(n)=sum _{d|n} mu(d) \
则有 f(n)= left{egin{aligned} 1，n=1\0，n>1end{aligned}
ight.
]
下面我们证明这个结论：
若(n=1)，结论显然。
若(n>1)，那么我们设
[egin{aligned}
s=prod _{i=1} ^k p_i^{a_i}，k|n
end{aligned}
]
当任意一个(a\_i>1)时，(mu(s)=0)
所以上面的(f(n))等价于(f(s))，其中(s=prod\_{i=1}^k p\_i)
根据莫比乌斯函数的定义，当(k)为奇数时，(f(s)=-1)，(k)为偶数时，(f(s)=1).
所以(f(n))可以写成：
[egin{aligned}
f(n)=C_k^0-C_k^1+C_k^2....=sum _{i=0}^k (-1)^i*C_k^i
end{aligned}
]
所以我们现在要证明的就是：
[egin{aligned}
sum _{i=0}^k (-1)^i*C_k^i=0
end{aligned}
]
注意到二项式定理：
[egin{aligned}
(x+y)^n=sum _{i=0}^k C_k^i*x^i*y^{k-i}
end{aligned}
]
与我们要证明的结论的相似性，我们令(x=-1)，(y=1)，即可得证。
这样，我们可以对原式进行变形：
[egin{aligned}
f_d(n)&=sum _{gcd(x,n)=1} x^d \
&=sum _{x=1}^n x^d[gcd(x,n)=1] \
&=sum _{x=1}^n x^dsum _{c|gcd(x,n)}mu(c) \
&=sum _{c|n} sum _{x=1}^{frac{n}{c}} (cx)^dmu(c) \
&=sum _{c|n} mu(c)c^d sum _{x=1} ^{frac{n}{c}} x^d \
end{aligned}
]
令：
[egin{aligned}
h(n)=sum _{i=1}^n i^d
end{aligned}
]
我们称(h)为自然数幂和。
则有：
[egin{aligned}
f_d(n)=sum _{c|n} mu(c)c^d h(frac{n}{c}) \
end{aligned}
]
这里有一个很奇妙的，并且正确的猜想：
[egin{aligned}
h(n)=sum _{i=1} ^{d+1}a_in^i
end{aligned}
]
也就是说，自然数幂和可以用(n)的多项式表示出来。
于是奇妙地用到了高斯消元，直接用(1)到(d+1)为例子，列出(d+1)个方程，强行把多项式的系数解出来。其实自然数幂和有通项公式，近期会有相关的文章。于是我们有：
[egin{aligned}
f_d(n)&=sum _{c|n} mu(c)c^d sum _{i=1} ^{d+1}a_i(frac{n}{c})^i \
&=sum _{i=1}^{d+1} a_i sum _{c|n} mu(c)c^d(frac{n}{c})^i
end{aligned}
]
令：
[egin{aligned}
g_i(n)=sum _{c|n} mu(c)c^d(frac{n}{c})^i
end{aligned}
]
那么有：
[egin{aligned}
f_d(n)=sum _{i=1}^{d+1} a_i *g_i(n)
end{aligned}
]
现在我们可以在(O(d^3))的时间内求出所有的(a_i)，只要对每个(i)求出(g_i(n))，我们就能(O(d))算出答案了。
若令(r(c)=mu(c)c^d)，那么有：
[egin{aligned}
g_i(n)=sum _{c|n} r(c)(frac{n}{c})^i
end{aligned}
]
这是一个狄利克雷卷积(即形如 (t(n)=sum _{d|n}f(d)g(frac{n}{d}))的函数)！
狄利克雷卷积函数有一个性质，如果原来两个函数都是积性函数，那么新函数也是积性函数。证明如下：
如果两个函数不完全积性，那么需要(gcd(x,y)=1)，否则不需要。
[egin{aligned}
t(x)t(y)&=sum _{d|x}f(d)g(frac{x}{d}) sum _{e|y}f(e)g(frac{y}{e}) \
&=sum _{de|xy}f(de)g(frac{xy}{de}) \
&=t(xy)
end{aligned}
]
所以我们可以通过(n)的质因数分解求出每个(g_i(n))：
[egin{aligned}
g_i(n)=prod _{j=1}^w g_i(p_j^{q_j})
end{aligned}
]
而对于只包含一个质数的(g)函数：
[egin{aligned}
g_i(p^q) &=sum _{j=0}^q mu(p^j)p^{jd}(frac{p^q}{p^j})i \
&= left{egin{aligned} p^{qi}，j=0\-p^{qi+d-i}，j=1\ 0, otherwiseend{aligned}
ight.
end{aligned}
]
所以有：
[egin{aligned}
g_i(p^q) =p^{qi}-p^{qi+d-i}
end{aligned}
]
直接计算即可。还是很简单的嘛！
代码
公式推错了四次～
#include<cstdio>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long giant;
giant read() {
giant x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const giant maxn=1e3+5;
const giant maxd=105;
const giant q=1e9+7;
giant p[maxn],s[maxn];
giant a[maxd][maxd],as[maxn];
giant mi(giant x,giant y) {
giant ret=1;
while (y) {
if (y&1) (ret*=x)%=q;
y>>=1,(x*=x)%=q;
}
return ret;
}
void elm(giant n) {
for (giant i=1;i<n;++i) {
if (!a[i][i]) {
for (giant j=i+1;j<=n;++j) if (a[j][i]) {
for (giant k=1;k<=n+1;++k) swap(a[j][k],a[i][k]);
break;
}
}
giant inv=mi(a[i][i],q-2);
for (giant j=i+1;j<=n;++j) if (a[j][i]) {
giant tmp=(a[j][i]*inv)%q;
for (giant k=1;k<=n+1;++k) a[j][k]=((a[j][k]-a[i][k]*tmp+q)%q+q)%q;
}
}
for (giant i=n;i>1;--i) {
if (!a[i][i]) {
for (giant j=i-1;j>0;--j) if (a[j][i]) {
for (giant k=1;k<=n+1;++k) swap(a[j][k],a[i][k]);
break;
}
}
giant inv=mi(a[i][i],q-2);
for (giant j=i-1;j>0;--j) if (a[j][i]) {
giant tmp=(a[j][i]*inv)%q;
for (giant k=1;k<=n+1;++k) a[j][k]=((a[j][k]-a[i][k]*tmp+q)%q+q)%q;
}
}
}
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
freopen("my.out","w",stdout);
#endif
giant d=read(),w=read();
for (giant i=1;i<=w;++i) p[i]=read(),s[i]=read();
giant sum=0;
for (giant i=1;i<=d+1;++i) {
(sum+=mi(i,d))%=q;
giant tmp=1;
for (giant j=1;j<=d+1;++j) (tmp*=i)%=q,a[i][j]=tmp;
a[i][d+2]=sum;
}
elm(d+1);
for (giant i=1;i<=d+1;++i) as[i]=(a[i][d+2]*mi(a[i][i],q-2))%q;
giant ans=0;
for (giant i=1;i<=d+1;++i) {
giant g=1;
for (giant j=1;j<=w;++j) {
giant tmp=(mi(p[j],s[j]*i)-mi(p[j],s[j]*i+d-i)+q)%q;
(g*=tmp)%=q;
}
(ans+=(as[i]*g))%=q;
}
printf("%lld
",ans);
}