一、題目

點此看題

二、解法

你會發現好像沒有什麼巧妙的演算法，不如我們直接把答案形式化地寫出來：

\[\sum_{i=1}^n[\gcd(i,n)=1]\space i^k \]

然後看到了熟悉的 \(\gcd\) 結構，我們直接反演：

\[\sum_{i=1}^ni^k\sum_{x|(i,n)}\mu(x) \]

然後我們列舉 \(x\) ：

\[\sum_{x|n}\mu(x)x^k\sum_{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^k\sum_{i=1}^kf_i\times(\frac{n}{x})^i \]

然後我們先列舉 \(i\) ，這也是我們希望看到的：

\[\sum_{i=1}^kf_i\sum_{x|n}\mu(x)x^k(\frac{n}{x})^i \]

\[\sum_{i=1}^kf_in^i\sum_{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\n",(ans+MOD)%MOD);
}