传送门

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

由于牛客限定了单份题解的字数上限，本文删去了一些公式的中间推导过程，仅留下了开头与结果。若无法推导出一致结果的，可以参考本人的 博客园博客

【大意】

求所有长度为 mm 且和为 nn 的正整数序列 aa 的贡献和。其中，每个数列的贡献为 ∑i>1gcd⁡(ai−1,ai)⋅w(ai)\displaystyle \sum_{i>1} \gcd(a_{i-1}, a_i)\cdot w(a_i)

【分析】

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

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

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

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

res=∑i,j≥1(m−1)⋅(n−i−j−1m−3)⋅gcd⁡(i,j)⋅w(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+ji+j 项，我们可以直接枚举 k=i+jk=i+j ：

res=(m−1)∑k=2n(n−k−1m−3)⋅∑i+j=kgcd⁡(i,j)⋅w(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)\gcd(i, j)=\gcd(k-j, j)=\gcd(k, j) ，因此有：

∑i+j=kgcd⁡(i,j)⋅w(j)=∑j=1k−1gcd⁡(k,j)⋅w(j)=∑j=1kgcd⁡(k,j)⋅w(j)−k⋅w(k)\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)

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

而 w(j)=∑tctjt\displaystyle w(j)=\sum_t c_tj^t 有 ∑j=1kgcd⁡(k,j)⋅w(j)=∑tct∑j=1kgcd⁡(k,j)⋅jt\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

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

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

∑j=1kgcd⁡(k,j)⋅jt=∑gdq=kgt+1⋅(μ(d)⋅dt)⋅St(q)\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) \end{aligned}

其中 St(n)=∑i=1nit=∑i=0t+1At,ini\displaystyle S_t(n)=\sum_{i=1}^n i^t=\sum_{i=0}^{t+1}A_{t,i}n^i ， At,iA_{t,i} 是自然数 tt 次幂和的第 ii 项系数

代入得到：

∑j=1kgcd⁡(k,j)⋅jt=∑i=0t+1At,i∑gdq=kgt+1⋅(μ(d)⋅dt)⋅qi\begin{aligned} &\sum_{j=1}^k\gcd(k, j)\cdot j^t \\=&\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}

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

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

ft,i(pe)=pei∑x=0e(px)t+1−i−p(e−1)i+t∑x=0e−1(px)t+1−i(i≤t+1)\begin{aligned} &\boldsymbol f_{t,i}(p^e) \\=&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>0t+1-i>0 时，有：

ft,i(pe)=pei⋅(1−pt−i)−pe(t+1)+t−i⋅(p−1)1−pt−i+1\begin{aligned} &\boldsymbol f_{t,i}(p^e) \\=&{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=0t+1-i=0 时，有：

ft,i(pe)=pei⋅(e+1)−pei+t−ie\begin{aligned} &\boldsymbol f_{t,i}(p^e) \\=&p^{ei}\cdot (e+1)-p^{ei+t-i}e \end{aligned}

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

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

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

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

考虑到 {S1(n)=n(n+1)2=12n2+12nS2(n)=n(n+1)(2n+1)6=13n3+12n2+16nS3(n)=[n(n+1)2]2=14n4+12n3+14n2\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)∈{(1,1),(1,2),(2,1),(2,2),(2,3),(3,2),(3,3),(3,4)}(t,i)\in\{(1, 1), (1,2), (2,1), (2,2), (2,3), (3,2), (3,3), (3,4)\} 共 88 项即可

【代码】

#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;
}