通过这道题学了伯努利数，写篇题解推一下
先推一下式子
[sum_{i=1}^ni^d[gcd(i,n)=1]
]
[sum_{i=1}^{n}i^dsum_{k|i}sum_{k|n}mu(k)
]
[sum_{k|n}mu(k)sum_{i=1}^{frac{n}{k}}(ik)^d
]
[sum_{k|n}mu(k)k^dsum_{i=1}^{frac{n}{k}}i^d
]
我们发现这个东西不好预处理，那么我们再简化一下。
设 (S_k(n)=sum_{i=1}^{n-1}i^k)
则原式等于
[sum_{k|n}mu(k)k^d(S_d(frac{n}{k})+(frac{n}{k})^d)
]
[sum_{k|n}mu(k)k^dS_d(frac{n}{k})+sum_{k|n}mu(k)n^d
]
(sum_{k|n}mu(k)) 等价与 ([n=1])，而本题中 (n) 不可能为 (1)，所以原式为
[sum_{k|n}mu(k)k^dS_d(frac{n}{k})
]
设 (f_i) 为在伯努利公式中 (i) 次幂的系数
[sum_{k|n}mu(k)k^dsum_{i=1}^{d+1}f_i(frac{n}{k})^i
]
[sum_{k|n}mu(k)sum_{i=1}^{d+1}f_in^ik^{d-i}
]
[sum_{i=1}^{d+1}f_in^isum_{k|n}mu(k)k^{d-i}
]
我们会发现式子后半部分就是 ((mu×id_{d-i})*I)，所以肯定是一个积性函数。
设 (F(n)=sum_{k|n}mu(k)k^{d-i},G(p)=mu(p)p^{d-i},F(n)=sum_{p|n}G(p))
考虑质数取值：
[G(p)=left{
egin{array}{lcl}
1kern 2.2em(p=1)\
-p^{d-i}kern 1.0em(pin prime)\
end{array}
ight.
]
由 (mu) 的性质，可以得出上式当 (nin prime) 时，(F(n)=1-n^{d-i})
所以总的复杂度即为 (O(d^2+dw))，(d^2) 为预处理 (f_i) 。
(AC kern 0.4emCODE:)
#include<bits/stdc++.h>
#define ri register int
#define p(i) ++i
using namespace std;
const int MOD=1e9+7,W=1010,D=107;
inline int read() {
int x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
return x*f;
}
inline int fpow(int x,int y) {
int res=1;
while(y) {
if (y&1) res=1ll*res*x%MOD;
x=1ll*x*x%MOD;y>>=1;
}
return res;
}
inline int inv(int x) {return fpow(x,MOD-2);}
int frac[D],invf[D];
void init(int k) {
frac[0]=1;
for (ri i(1);i<=k;p(i)) frac[i]=1ll*frac[i-1]*i%MOD;
invf[k]=inv(frac[k]);
for (ri i(k);i;--i) invf[i-1]=1ll*invf[i]*i%MOD;
}
inline int C(int n,int m) {return 1ll*frac[n]*invf[m]%MOD*invf[n-m]%MOD;}
struct Bernou{
int B[D],f[D];
inline void Bernoulli(int k) {
B[0]=1;
for (ri i(1);i<=k;p(i)) {
int sum=0;
for (ri j(0);j<i;p(j)) sum=(sum+1ll*C(i+1,j)*B[j]%MOD)%MOD;
B[i]=(0-1ll*sum*inv(i+1)%MOD+MOD)%MOD;
}
int INV=inv(k+1);
for (ri i(0);i<=k;p(i)) f[k+1-i]=1ll*C(k+1,i)*B[i]%MOD*INV%MOD;
}
}B;
int d,w,p[W],a[W];
inline int Sum(int cm) {
if (cm<0) cm=MOD-2;
int res=1;
for (ri i(1);i<=w;p(i)) res=1ll*res*(1-fpow(p[i],cm)+MOD)%MOD;
return res;
}
int main() {
d=read(),w=read();
init(d+3);B.Bernoulli(d);
for (ri i(1);i<=w;p(i)) p[i]=read(),a[i]=read();
int ans=0,n=1;
for (ri i(1);i<=w;p(i)) n=1ll*n*fpow(p[i],a[i])%MOD;
for (ri i(1),cm(1);i<=d+1;p(i)) {
cm=1ll*cm*n%MOD;
ans=(ans+1ll*B.f[i]*cm%MOD*Sum(d-i)%MOD)%MOD;
}
printf("%d
",ans);
return 0;
}