题解 I. gcd -“山大地纬杯”第十二届山东省ICPC大学生程序设计竞赛(正式赛）

传送门

VP 的时候失误推错太多次了，写个博客总结一下

---

【大意】

求所有长度为 $m$ 且和为 $n$ 的正整数序列 $a$ 的贡献和。其中，每个数列的贡献为 $\displaystyle \sum_{i>1} \gcd(a_{i-1}, a_i)\cdot w(a_i)$

---

【分析】

考虑从 $m$ 个位置中，选择两个相邻的位置分别填写 $i, j$ ，则方案数为 $(m-1)$

剩余的 $m-2$ 个正整数的和为 $n-i-j$ ，则等价于模型“$n-i-j$ 个小球放入 $m-2$ 个盒子里，要求盒子非空”。由隔板法，方案数为 $\dbinom {n-i-j-1} {m-3}$

因此，选中的两个位置分别填写 $i, j$ 后，对答案的贡献为 $\displaystyle (m-1)\cdot \dbinom {n-i-j-1} {m-3} \cdot \gcd(i, j)\cdot w(j)$

我们只需要枚举正整数 $i, j$ 即可：

$\displaystyle res=\sum_{i, j\geq 1}(m-1)\cdot \dbinom {n-i-j-1} {m-3} \cdot \gcd(i, j)\cdot w(j)$

注意到公式中重复出线的 $i+j$ 项，我们可以直接枚举 $k=i+j$ ：

$\displaystyle res=(m-1)\sum_{k=2}^n \dbinom {n-k-1} {m-3} \cdot \sum_{i+j=k}\gcd(i, j)\cdot w(j)$

考虑 $\gcd(i, j)=\gcd(k-j, j)=\gcd(k, j)$ ，因此有：

$\displaystyle \sum_{i+j=k}\gcd(i, j)\cdot w(j)=\sum_{j=1}^{k-1} \gcd(k, j)\cdot w(j)=\sum_{j=1}^k \gcd(k, j)\cdot w(j)-k\cdot w(k)$

因此，现在仅需要求解 $\displaystyle \sum_{j=1}^k\gcd(k, j)\cdot w(j)$ 即可求得答案

而 $\displaystyle w(j)=\sum_t c_tj^t$ 有 $\displaystyle \sum_{j=1}^k\gcd(k, j)\cdot w(j)=\sum_t c_t\sum_{j=1}^k\gcd(k, j)\cdot j^t$

于是问题简化为如何求解 $\displaystyle \sum_{j=1}^k\gcd(k, j)\cdot j^t$

---

对于代求式子，我们枚举其 gcd ，并莫比乌斯反演，得到：

$\begin{aligned} &\sum_{j=1}^k\gcd(k, j)\cdot j^t \\=&\sum_{g\mid k} g \sum_{j=1}^k[\gcd(k, j)=g]\cdot j^t \\=&\sum_{g\mid k} g \sum_{j=1}^{k\over g}[\gcd({k\over g}, j)=1]\cdot (jg)^t \\=&\sum_{g\mid k} g^{t+1} \sum_{j=1}^{k\over g}\sum_{d\mid {k\over g}\wedge d\mid j}\boldsymbol \mu(d)\cdot j^t \\=&\sum_{gd\mid k} g^{t+1}\boldsymbol \mu(d)\sum_{j=1}^{k\over gd} (jd)^t \\=&\sum_{gd\mid k} g^{t+1}\cdot (\boldsymbol \mu(d)\cdot d^t)\cdot \sum_{j=1}^{k\over gd} j^t \\=&\sum_{gdq=k}g^{t+1}\cdot (\boldsymbol \mu(d)\cdot d^t)\cdot S_t(q) \end{aligned}$

其中 $\displaystyle S_t(n)=\sum_{i=1}^n i^t=\sum_{i=0}^{t+1}A_{t,i}n^i$ ， $A_{t,i}$ 是自然数 $t$ 次幂和的第 $i$ 项系数

代入得到：

$\begin{aligned} &\sum_{j=1}^k\gcd(k, j)\cdot j^t \\=&\sum_{gdq=k}g^{t+1}\cdot (\boldsymbol \mu(d)\cdot d^t)\cdot S_t(q) \\=&\sum_{gdq=k}g^{t+1}\cdot (\boldsymbol \mu(d)\cdot d^t)\cdot \sum_{i=0}^{t+1}A_{t, i} q^i \\=&\sum_{i=0}^{t+1}A_{t, i} \sum_{gdq=k} g^{t+1}\cdot (\boldsymbol \mu(d)\cdot d^t)\cdot q^i \end{aligned}$

最后这个求和项 $\displaystyle \sum_{gdq=k} g^{t+1}\cdot (\boldsymbol \mu(d)\cdot d^t)\cdot q^i$ 是积性函数 $\boldsymbol {id}^{t+1}, \boldsymbol \mu\cdot \boldsymbol {id}^t, \boldsymbol {id}^i$ 的迪利克雷卷积，显然也是积性函数

不妨记上述积性函数为 $\boldsymbol f_{t,i}$ ，则有：

$\begin{aligned} &\boldsymbol f_{t,i}(p^e) \\=&\sum_{x+y+z=e}(p^x)^{t+1}\cdot (\boldsymbol \mu(p^y)\cdot (p^y)^t)\cdot (p^z)^i \\=&\sum_{x+z=e}(p^x)^{t+1}\cdot (p^z)^i-\sum_{x+z=e-1} (p^x)^{t+1}\cdot (p^z)^i\cdot p^t \\=&p^{ei}\sum_{x=0}^e (p^x)^{t+1-i}-p^{(e-1)i+t}\sum_{x=0}^{e-1}(p^x)^{t+1-i}&(i\leq t+1) \end{aligned}$

故 $t+1-i>0$ 时，有：

$\begin{aligned} &\boldsymbol f_{t,i}(p^e) \\=&p^{ei}\cdot {1-p^{(e+1)(t+1-i)}\over 1-p^{t+1-i}}-p^{(e-1)i+t}\cdot {1-p^{e(t+1-i)}\over 1-p^{t+1-i}} \\=&{p^{ei}\cdot (1-p^{t-i}) - p^{e(t+1)+t-i}\cdot (p-1)\over 1-p^{t-i+1}} \end{aligned}$

$t+1-i=0$ 时，有：

$\begin{aligned} &\boldsymbol f_{t,i}(p^e) \\=&p^{ei}\cdot (e+1)-p^{ei+t-i}e \end{aligned}$

这种式子显然可以欧拉筛：

对于某个质数 $p$ ，直接算出它所有幂次的值，然后欧拉筛过程中额外维护最小质因数出线的次数。转移的时候直接查表就行。

复杂度为原欧拉筛复杂度 $O(n)$ ，加上处理质数幂次位置的结果。经一些不严谨的半猜半分析，后半部分的复杂度应该也是 $O(n)$ 。

由于题目所给的式子仅有 $1$~$3$ 次项，因此我们只要预处理 $t=1,2,3;i\leq t+1$ 的 $\boldsymbol f_{t,i}$ 即可筛出所有答案

考虑到 $\begin{cases} \begin{aligned} S_1(n)&={n(n+1)\over 2}&={1\over 2}n^2+{1\over 2}n \\S_2(n)&={n(n+1)(2n+1)\over 6}&={1\over 3}n^3+{1\over 2}n^2+{1\over 6}n \\S_3(n)&=[{n(n+1)\over 2}]^2&={1\over 4}n^4+{1\over 2}n^3+{1\over 4}n^2 \end{aligned} \end{cases}$

因此，仅需要处理 $(t,i)\in\{(1, 1), (1,2), (2,1), (2,2), (2,3), (3,2), (3,3), (3,4)\}$ 共 $8$ 项即可

---

【代码】

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int Lim=1e7, MAXN=Lim+10, P=1e9+7;

inline int kpow(int a, int x) { int ans=1; for(; x; x>>=1, a=(ll)a*a%P) if(x&1) ans=(ll)ans*a%P; return ans; }
inline int inv(int a) { return kpow(a, P-2); }
inline int norm(int v) { return v>=P?v-P:v; }
struct Z {
	int v;
	inline Z(int v_=0):v(norm(v_)) {}
	inline Z& operator += (const Z& x) { v=norm(v+x.v); return *this; }
	inline Z& operator -= (const Z& x) { v=norm(v+P-x.v); return *this; }
	inline Z& operator *= (const Z& x) { v=(ll)v*x.v%P; return *this; }
	inline Z operator + (const Z& x) const { Z y=*this; return y+=x; }
	inline Z operator - (const Z& x) const { Z y=*this; return y-=x; }
	inline Z operator * (const Z& x) const { Z y=*this; return y*=x; }
	inline Z operator ! () const { return inv(v); }
	inline operator int() const { return v; }
};
const int SIZ=8;
struct multi_F : public array<Z, SIZ> {
	inline multi_F() { fill(0); }
	inline multi_F& operator += (const multi_F &x) {
		for(int i=0; i<SIZ; ++i) at(i)+=x[i];
		return *this;
	}
	inline multi_F& operator -= (const multi_F &x) {
		for(int i=0; i<SIZ; ++i) at(i)-=x[i];
		return *this;
	}
	inline multi_F& operator *= (const multi_F &x) {
		for(int i=0; i<SIZ; ++i) at(i)*=x[i];
		return *this;
	}
	inline multi_F operator + (const multi_F &x) const { multi_F y=*this; return y+=x; }
	inline multi_F operator - (const multi_F &x) const { multi_F y=*this; return y-=x; }
	inline multi_F operator * (const multi_F &x) const { multi_F y=*this; return y*=x; }
};
inline multi_F f(int p, int k, ll pk, multi_F lst) {
	static int t[] = {1, 1, 2, 2, 2, 3, 3, 3};
	static int i[] = {1, 2, 1, 2, 3, 2, 3, 4};
	for(int j=0; j<SIZ; ++j) {
		int dif = t[j] + 1 - i[j];
		int a = kpow(p, dif+P-2);
		if(dif) {
			int b = (kpow(pk, i[j]) * (1ll-a) - kpow(pk, t[j]+1) * (p-1ll)%P* a)%P;
			lst[j] = (ll)b * inv(1-(ll)a*p%P)%P + P;
		}
		else
			lst[j] = kpow(pk, i[j]) * (k + 1ll - (ll)a*k%P) %P+P;
	}
	return lst;
}
int minp[MAXN], fr[MAXN], prime[MAXN/10], cntprime;
multi_F fi[MAXN];
inline void sieve() {
	fi[1].fill(1);
	for(int i=2; i<=Lim; ++i) {
		if(!minp[i]) {
			prime[++cntprime]=i;
			multi_F val=fi[1];
			for(ll pk=i, k=1; pk<=Lim; pk*=i, ++k) {
				fi[pk]=val=f(i, k, pk, val);
				minp[pk]=i;
				fr[pk]=1;
			}
		}
		for(int j=1; j<=cntprime; ++j) {
			int x=i*prime[j];
			if(prime[j] > minp[i] || x>Lim) break;
			minp[x]=prime[j];
			if(prime[j]==minp[i]) fr[x]=fr[i];
			else fr[x]=i;
			fi[x]=fi[fr[x]]*fi[x/fr[x]];
		}
	}
}

int n, m, p, q, r;
Z fact[MAXN], inft[MAXN];
inline Z C(int n, int m) { return m<0||m>n?Z(0):fact[n]*inft[m]*inft[n-m]; }
Z h[MAXN];
inline void init() {
	sieve();
	cin>>n>>m;
	cin>>p>>q>>r;
	Z rr=Z(r)*!Z(2);
	Z qq=Z(q)*!Z(6);
	Z pp=Z(p)*!Z(4);
	for(int i=1; i<=Lim; ++i) {
		Z g1 = (fi[i][0]*Z(1) + fi[i][1]*Z(1)) * rr;
		Z g2 = (fi[i][2]*Z(1) + fi[i][3]*Z(3) + fi[i][4]*Z(2)) * qq;
		Z g3 = (fi[i][5]*Z(1) + fi[i][6]*Z(2) + fi[i][7]*Z(1)) * pp;
		h[i]=g1+g2+g3;
	}
	
	fact[0]=1;
	for(int i=1; i<=Lim; ++i) fact[i]=fact[i-1]*Z(i);
	inft[Lim]=!fact[Lim];
	for(int i=Lim; i; --i) inft[i-1]=inft[i]*Z(i);
}
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	init();
	Z res=0;
	for(int k=2; k<n; ++k) {
		Z tmp=((Z(p)*Z(k)+Z(q))*Z(k)+Z(r))*Z(k)*Z(k);
		res+=C(n-k-1, m-3)*(h[k]-tmp);
	}
	res*=Z(m-1);
	cout<<res;
	return 0;
}