一、题目
二、解法
你会发现好像没有什么巧妙的算法，不如我们直接把答案形式化地写出来：
[sum_{i=1}^n[gcd(i,n)=1]space i^k
]
然后看到了熟悉的 (gcd) 结构，我们直接反演：
[sum_{i=1}^ni^ksum_{x|(i,n)}mu(x)
]
然后我们枚举 (x) ：
[sum_{x|n}mu(x)x^ksum_{i=1}^{n/x} i^k
]
你会发现还是很复杂，现在的契机在于后面的 (sum_{i=1}^{n/x}i^k) 是有很多方法的，但是我们的标准只有一个：把 (k) 融入到枚举中所以我们就能得到关于 (k) 的复杂度 ，我们可以试一试用第二类斯特林数展开它（但是我试过不行），还有一种方法，因为他是一个 (k+1) 次的关于 (n/x) 的多项式，我们假设已经计算了他的系数 (f_i) ，这样改写算式：
[sum_{x|n}mu(x)x^ksum_{i=1}^kf_i imes(frac{n}{x})^i
]
然后我们先枚举 (i) ，这也是我们希望看到的：
[sum_{i=1}^kf_isum_{x|n}mu(x)x^k(frac{n}{x})^i
]
[sum_{i=1}^kf_in^isum_{x|n}mu(x)x^{k-i}
]
现在要好好观察后面一部分了，他是我们的复杂度瓶颈，这道题很特殊的地方是题目给了我们一个乘积式来表示 (n) ，后面那一部分又一定是一个积性函数（因为 (mu) 和 (x^{k-i}) 都是可以拆分的），所以可以对每个质数分开考虑。
对于质数 (p_i^{a_i}) ，当 (x=1) 时答案是 (1) ，当 (x=p_i) 时答案是 (-p_i^{k-i}) ，对于其他情况 (mu=0) 所以不用考虑。
现在问题变成了算系数 (f_i) ，你可以用高斯消元爆操，但是有一种更优美的方法：伯努利数。伯努利数就是用来处理这种幂求和的问题，学会这种方法你只需要记住两个公式，(B_i) 表示伯努利的定义式是这样的：
[sum_{i=0}^n C(n+1,i)B_i=[n=0]
]
那么我们 (O(k^2)) 就可以求 (B_i) ，怎么求 (f_i) 呢？首先要知道伯努利公式：
[sum_{i=1}^ni^k=frac{1}{k+1}sum_{i=0}^kC(k+1,i)B_in^{k+1-i}
]
所以我们就知道了 (f_{k+1-i}=frac{C(k+1,i)B_i}{k+1}) ，证明我暂时不懂，以后再填坑吧。
#include <cstdio>
const int M = 1005;
const int MOD = 1e9+7;
#define int long long
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int k,w,ans,p[M],a[M],B[M],f[M],fac[M],inv[M];
int C(int n,int m)
{
if(n<m) return 0;
return fac[n]*inv[m]%MOD*inv[n-m]%MOD;
}
int qkpow(int a,int b)
{
int r=1;
while(b>0)
{
if(b&1) r=r*a%MOD;
a=a*a%MOD;
b>>=1;
}
return r;
}
void init()
{
fac[0]=inv[0]=inv[1]=1;
for(int i=1;i<=k+1;i++) fac[i]=fac[i-1]*i%MOD;
for(int i=2;i<=k+1;i++) inv[i]=inv[MOD%i]*(MOD-MOD/i)%MOD;
for(int i=2;i<=k+1;i++) inv[i]=inv[i-1]*inv[i]%MOD;
B[0]=1;
for(int i=1;i<=k;i++)
{
int sum=0;
for(int j=0;j<i;j++)
sum=(sum-C(i+1,j)*B[j])%MOD;
B[i]=sum*qkpow(i+1,MOD-2)%MOD;
}
for(int i=0;i<=k;i++)
f[k+1-i]=C(k+1,i)*B[i]%MOD*qkpow(k+1,MOD-2)%MOD;
}
signed main()
{
k=read();w=read();
for(int i=1;i<=w;i++)
{
p[i]=read();a[i]=read();
}
init();
for(int i=1;i<=k+1;i++)
{
int res=f[i];
for(int j=1;j<=w;j++)
{
res=res*qkpow(p[j],a[j]*i)%MOD;//这里我懒了，多了个log
res=res*(1-qkpow(p[j],k-i+MOD-1))%MOD;
}
ans=(ans+res)%MOD;
}
printf("%lld
",(ans+MOD)%MOD);
}