题解 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;
}
