题解 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;
}
posted @ 2022-11-25 11:54  JustinRochester  阅读(51)  评论(0编辑  收藏  举报