【題目地址】

題意簡述

給你$n,d$，求下面式子在${\rm mod}\ 10^9+7$意義下的值。

$\sum_{i=1}^n[gcd(i,n)=1]i^d$

其中$n$特別大，所以我們給你一個$w$，然後給你兩個陣列$p_i,c_i$，其中：

$n=\prod_{i=1}^wp_i^{c_i}$

$d\leq 100,w\leq 1000,p_i,c_i\leq 10^9$

---

我們仍舊先推一波式子：

$\begin{aligned} ans=&\sum_{i=1}^n[gcd(i,n)=1]i^d \\ =& \sum_{i=1}^ni^d\sum_{j|i,j|n}\mu(j) \\ =& \sum_{j=1}^n\mu(j)\sum_{j|i}^ni^d \\ =& \sum_{j=1}^n\mu(j)\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}(ij)^d \\ =&\sum_{j=1}^n\mu(j)j^d\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}i^d \end{aligned}$

後面部分就是自然數冪的字首和，拉格朗日差值就可以了，在$O(dlogd)$時間內求出，而前面是可以杜教篩的：

$\begin{aligned} \sum_{j|n}\mu(j)j^d=&(\mu\cdot id^d)*\mathbf1 \end{aligned}$

我們令$f=(\mu\cdot id^d)*\mathbf1,g=id^d$，然後$f*g$就等於：

$\sum_{j|n}\mu(j)j^d(\frac{n}{j})^d=n^d\sum_{j|n}\mu(j)=n^d[n=1]=[n=1]=\epsilon$

所以$f=(\mu\cdot id^d)*\mathbf1,g=id^d,h=\epsilon$，然後套入杜教篩公式即可：

$S(n)=1-\sum_{i=2}^ni^dS(\lfloor\frac{n}{i}\rfloor)$

然後就可以在$O(n^{\frac{2}{3}}dlogd)$的時間內解決了。

---

等等，$n$好像有$10^{10000...}$左右，-_-||。

但是不是就這樣束手就擒，打暴力去了呢？

不，我們離正解還差一步了。

---

但是，我們根據拉格朗日差值的經驗可以知道，$\sum_{i=1}^ni^d$為一個$d+1$次多項式，那麼我們令$g(x)=\sum_{i=1}^xi^d$，同時我們也可以知道$g(x)=\sum_{i=1}^{d+1}a_ix^i$，假設我們求出了所有$a_i$，我們將其帶入式子，看看有沒有什麼奇妙的反應：

$\begin{aligned} ans=&\sum_{j=1}^n\mu(j)j^d\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}i^d \\ =&\sum_{j=1}^n\mu(j)j^d\sum_{i=1}^{d+1}a_i\left(\left\lfloor\frac{n}{j}\right\rfloor\right)^i \\ =&\sum_{i=1}^{d+1}a_i\sum_{j|n}\mu(j)j^d\left(\left\lfloor\frac{n}{j}\right\rfloor\right)^i \end{aligned}$

現在我們主要是要解決後半部分的快速求法，我們令$f(i)=\sum_{j|n}\mu(j)j^d\left(\left\lfloor\frac{n}{j}\right\rfloor\right)^i$，那麼我們知道$f=(\mu\cdot id^d)*id^i$，那麼根據$\mu,id^k$都為積性函式，而積性函式的狄利克雷卷積都為積性函式，所以我們可以知道$f$也是一個積性函式。

這裡，我們看，題目都將$n$的質因數幫我們分解好了凱森(ノ￣▽￣)，所以我們將質因數分開考慮，因為$f(n)=\prod f(p_i^{c_i})$，所以此時，我們只需知道$f(p^c)$如何計算：

我們看原式
$\sum_{j|p^c}\mu(j)j^d\left(\left\lfloor\frac{n}{j}\right\rfloor\right)^i=\sum_{k=0}^c\mu(p^k)(p^k)^d\left(\left\lfloor\frac{n}{p^k}\right\rfloor\right)^i$

而根據$\mu$的定義，我們知道當$k>1$時的時候$\mu(p^k)=0$，所以$k=0,1$，那麼$f(p^c)$就等於：

$(p^c)^i-p^d(p^{c-1})^i=p^{ci}-p^{ci-i+d}$

那麼列舉$a_i$，然後列舉$p_j$，然後算就用快速冪，複雜度就是$O(dwlog(10^9+6))$（因為根據尤拉定理和費馬小定理，我們可以將指數${\rm mod}\ 10^9+6$）。

好啦，就解決啦！ヾ(@^▽^@)ノ

---

等等辣QAQ，$a_i$好像還沒說怎麼求-_-||。

其實，我們知道它是一個$d+1$次多項式的係數，那麼高斯消元就好啦，我們列出$g(1)\sim g(d+1)$的方程即可。

前面回顧：$g(x)=\sum_{i=1}^{d+1}a_ix^i$

這個方程組的初始化，我們每次列舉$x$，然後用快速冪計算出當前的$g(x)$的值，然後列舉$a_i$，依次計算出$x^i$作為係數即可，然後就是高斯消元板子啦~~

---

醜陋的加調了半天的程式碼_(¦3」∠)_：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=1e3+10,M=1e2+10;
const ll Mod=1e9+7;
int d,w;
ll fpow(ll a,ll b){
	ll res=1;
	for(;b;b>>=1,a=(a*a)%Mod)
		if(b&1)res=(res*a)%Mod;
	return res;
}
ll A[M][M],X[M],P[N],C[N];
void Guass(){
	int D=d+1;
	for(ll i=0,sum=0,x;i<=D;i++){
		x=1;sum=(sum+fpow(i,d))%Mod;
		A[i][D+1]=sum;
		for(ll j=0;j<=D;j++){
			A[i][j]=x;x=(x*i)%Mod;	
		}
	}
	for(int now=0,pos;now<=D;now++){
		pos=now;
		if(!A[pos][now])
			for(int i=now+1;i<=D;i++)
				if(A[i][now]){pos=i;break;}
		if(pos!=now)swap(A[now],A[pos]);
		ll inv=fpow(A[now][now],Mod-2);
		for(int i=now+1;i<=D;i++){
			ll tt=(A[i][now]*inv)%Mod;
			for(int j=0;j<=D+1;j++){
				A[i][j]=(A[i][j]-A[now][j]*tt%Mod)%Mod;
				if(A[i][j]<0)A[i][j]+=Mod;
			}
		}
	}
	for(int i=D;i>=0;i--){
		for(int j=D;j>=i+1;j--){
			A[i][D+1]=(A[i][D+1]-A[i][j]*X[j]%Mod)%Mod;
			if(A[i][D+1]<0)A[i][D+1]+=Mod;
		}
		X[i]=A[i][D+1]*fpow(A[i][i],Mod-2)%Mod;
	}
}
ll GetH(ll a,ll b){
	ll p1=C[b]*a%(Mod-1);
	ll p2=(p1+d-a)%(Mod-1);if(p2<0)p2+=(Mod-1);
	ll ans=(fpow(P[b],p1)-fpow(P[b],p2))%Mod;
	if(ans<0)ans+=Mod;
	return ans;
}
ll ans;
int main(){
	scanf("%d%d",&d,&w);
	Guass();
	for(int i=1;i<=w;i++)scanf("%lld%lld",&P[i],&C[i]);
	for(int i=0,D=d+1;i<=D;i++){
		ll xx=X[i];
		for(int j=1;j<=w;j++)xx=(xx*GetH(i,j))%Mod;
		ans=(ans+xx)%Mod;
	}
	printf("%lld\n",ans);
	return 0;
}

---

End

元旦節只放一天，1月1日還要在學校度過，慘兮兮QWQ