【题解】【原创题目】第〇类循环数•行
【题解】【原创题目】第〇类循环数·行
传送门：第〇类循环数·行
【题目描述】
给出 \(n,m,l,k,x\)，求 \(\frac{1}{nml}\sum\limits_{a=1}^{n}\sum\limits_{b=1}^{m}\sum\limits_{c=1}^{l}f(k,(ax^2+bx+c)\bmod k) \pmod{1004535809}\)，其中 \(f(i,j)\) 定义如下：
【分析】
【Subtask #1】
抄题目描述总会吧。
本测试包轻微卡常，需要提前把 \(f\) 全部预处理出来。
复杂度 \(O(k^2+n^3T)\)，期望得分：\(1pts\)。
【Code #1】
in(T),inv[1]=1;
for(Re i=2;i<=M-3;++i)inv[i]=(LL)inv[P%i]*(P-P/i)%P;
for(Re i=1;i<=M-3;++i){
f[i][0]=i*i+i;
for(Re j=1;j<=i-1;++j)
f[i][j]=((j?(
(LL)inv[i]*inv[j]%P*f[i-1][j-1]%P+
(LL)inv[i]*f[i-1][j-1]%P+
(LL)inv[j]*f[i-1][j-1]%P+
f[i-1][j-1]
)%P:0)+
(j<i-1?((LL)inv[i]*f[i-1][j]%P+f[i-1][j])%P:0))%P;
}
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
for(Re i=1;i<=n;++i)
for(Re j=1;j<=m;++j)
for(Re k=1;k<=l;++k)
(ans+=f[K][((LL)i*x2%K+(LL)j*x%K+k)%K])%=P;
printf("%lld\n",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
【Subtask #2】
存在特殊性质 \(k|x\)，说明 \(q\) 与 \(a,b\) 无关，即求 \(\frac{1}{l}\sum_{c=1}^{l}f(k,c\bmod k)\)，发现暴力枚举有很多重复计算，让 \(l\) 对 \(k\) 取个模即可降到 \(O(k)\) 。
本测试包轻微卡常。
复杂度 \(O(k^2+kT)\)，期望得分：\(1+9=10pts\)。
【Code #2】
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
Re tmp=0;
for(Re k=0;k<K;++k)(tmp+=f[K][k])%=P;
(ans+=(LL)l/K%P*tmp%P)%=P;
for(Re k=1;k<=(l%K);++k)(ans+=f[K][k])%=P;
printf("%lld\n",(LL)ans*Inv(l)%P);
}
【Subtask #3】
用类似 \(\text{Subtask 2}\) 的降级做法，依次处理两个东西：\(g_1(a)=\sum_{i=1}^{l}f_{k}(a+i),\) \(g_2(a)=\sum_{i=1}^{m}g_1(a+ix)\)，最后答案为 \(\sum_{i=1}^{n}g_{2}(ix^2)\) 。
本做法较卡常，需要把加法取模优化掉。
复杂度 \(O(k^2+k^2T)\)，期望得分：\(1+9+11=21pts\)。
【Code #3】
inline void mo(Re &x,Re mod=P){x=x>=mod?x-mod:x;}
inline int Mo(Re x,Re mod=P){return x>=mod?x-mod:x;}
int main(){
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
for(Re a=0;a<K;++a){
Re tmp=0;
if(l>=K)for(Re k=0,Tmp=a;k<K;++k)mo(tmp+=f[K][Tmp]),mo(++Tmp,K);
g1[a]=(LL)l/K%P*tmp%P;
for(Re k=1,Tmp=Mo(a+1,K),ed=l%K;k<=ed;++k)mo(g1[a]+=f[K][Tmp]),mo(++Tmp,K);
}
for(Re a=0;a<K;++a){
Re tmp=0;
if(m>=K)for(Re k=0,Tmp=a;k<K;++k)mo(tmp+=g1[Tmp]),mo(Tmp+=x,K);
g2[a]=(LL)m/K%P*tmp%P;
for(Re k=1,Tmp=Mo(a+x,K),ed=m%K;k<=ed;++k)mo(g2[a]+=g1[Tmp]),mo(Tmp+=x,K);
}
Re tmp=0;
if(n>=K)for(Re k=0,Tmp=0;k<K;++k)mo(tmp+=g2[Tmp]),mo(Tmp+=x2,K);
ans=(LL)n/K%P*tmp%P;
for(Re k=1,Tmp=x2,ed=n%K;k<=ed;++k)mo(ans+=g2[Tmp]),mo(Tmp+=x2,K);
printf("%lld\n",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
}
【Subtask #4】
观察 \(f\) 递推式，首先通分合并同类项：
两边同时除以 \((i+1)(j+1)\)：
发现是组合数的杨辉三角递推式。设 \(g(i,j)=\frac{f(i,j)}{(i+1)(j+1)}\)，则有 \(g(i,j)=g(i-1,j-1)+[j\neq i-1]g(i-1,j)\) 。
注意此处的细节：通常组合数递推为 \(\dbinom{n}{m}=\dbinom{n-1}{m-1}+\dbinom{n-1}{m}\)，当 \(m=n\) 时 \(\dbinom{n-1}{m}\) 才为 \(0\)，但 \(g(i-1,j)\) 前面的系数为 \([j\neq i-1]\)，因此 \(g(i,j)\) 应该等于\(\dbinom{i}{j+1}\) 而不是 \(\dbinom{i}{j}\) ，否则就与递推式不相符。
得出 \(f(i,j)=(i+1)(j+1)\dbinom{i}{j+1}\)，这个东西算出来的 \(f(i,0)\) 和前面定义是符合的。
注意到 \(n\) 比较小，预处理出组合数后直接沿用 \(\text{Subtask 3}\) 的做法即可。
复杂度 \(O(\min(n,k)kT)\)，期望得分：\(1+9+11+19=40pts\)。
【Code #4】
代码略。
【Subtask #5】
推出 \(f\) 通项式后直接沿用 \(\text{Subtask 2}\) 的做法。
此做法良心出题人没有卡常。
复杂度 \(O(kT)\)，期望得分：\(9+18=27pts\)，结合前面算法可得 \(58pts\)。
【Code #5】
代码略。
【Subtask #6】
前置知识：单位根反演。【常见公式汇总】
按照套路，枚举 \(f\) 的参数并交换和式：
\([(ax^2+bx+c)\bmod k=d]\) 可化为 \([k|(ax^2+bx+c-d)]\)，套用单位根反演的柿子，得到：
提一个 \(w_{k}^{-id}\) 出来，并交换和式：
右边那一坨显然可以简化为三个等比数列和的乘积，记一个 \(calc\) 函数：
用组合数吸收公式去掉 \(f\) 柿子里的 \((j+1)\)，即 \(f(i,j)=(i+1)(j+1)\dbinom{i}{j+1}=i(i+1)\dbinom{i-1}{j}\)，然后代进去用二项式定理合并：
本测试包不卡常，可以直接暴算。
复杂度 \(O(Tk\log n))\)，期望得分：\(91pts\) 。
【Code #6】
inline int calc(Re x,Re n){//注意等比数列求和要特判公比等于1的情况
return x==1?n:(LL)x*(1-mi(x,n)+P)%P*Inv((1-x+P)%P)%P;
}
inline int calc(Re w){
return (LL)calc(mi(w,x2),n)*calc(mi(w,x),m)%P*calc(w,l)%P;
}
int main(){
in(T);
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
omega[0]=1,omega[1]=mi(G,(P-1)/K);
for(Re i=2;i<=K;++i)omega[i]=(LL)omega[i-1]*omega[1]%P;
for(Re i=0;i<=K-1;++i)
(ans+=(LL)mi(omega[K-i]+1,K-1)*calc(omega[i])%P)%=P;
ans=(LL)ans*(K+1)%P;
printf("%lld\n",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
}
【Subtask #7】
最后两个测试点需要大力卡常。
大常数主要来源于快速幂，枚举一个 \(i\) 要调用 \(9\) 次，极限数据跑了约 \(5.8s\)（洛谷上要快 \(0.6s\)）。
实际上其中有 \(5\) 次都在算单位根 \(w_{k}^{1}\) 的若干次幂，由于 \(w_{k}^{i}=w_{k}^{i\bmod k}\)，只要预处理出 \(0\sim k-1\) 次幂就能 \(O(1)\) 查询了。
然后是 \(3\) 次计算逆元，同样可以预处理，具体做法见 乘法逆元 \(2\) \(\text{[P5431]}\) 。
其实单位根还可以放到询问前面预处理，对于不同的 \(k\) 根据 \(w_{k}^{a}=w_{kn}^{an}\) 这一性质直接得到。
但瓶颈不在于此，常数影响不大。
复杂度 \(O(Tk\log k))\)，期望得分：\(100pts\) 。
【Code #7】
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#define LD double
#define LL long long
#define Re register int
using namespace std;
const int N=1e9+3,M=1048576+3,P=1004535809,G=3;
int n,m,l,K,T,x2,ans,Sm[M],invo[M],omega[M];LL x;
template<typename T>inline void in(T &x){
int f=0;x=0;char ch=getchar();
while(ch<'0'||ch>'9')f|=ch=='-',ch=getchar();
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
x=f?-x:x;
}
inline int mi(Re x,Re k){
Re s=1;
while(k){
if(k&1)s=(LL)s*x%P;
x=(LL)x*x%P,k>>=1;
}
return s;
}
inline int Inv(Re x){return mi(x,P-2);}
inline void mo(Re &x,Re mod=P){x=x>=mod?x-mod:x;}
inline int Mo(Re x,Re mod=P){return x>=mod?x-mod:x;}
inline int calc(Re i,Re n){//注意等比数列求和要特判公比等于1的情况
return (omega[i]==1?n:(LL)omega[i]*(omega[(LL)i*n%K]-1+P)%P*invo[i]%P);
}
int main(){
in(T);
while(T--){
in(n),in(m),in(l),in(K),in(x),x%=K,x2=(LL)x*x%K,ans=0;
omega[0]=Sm[0]=1,Sm[1]=(omega[1]=mi(G,(P-1)/K))-1;
for(Re i=2;i<=K;++i)omega[i]=(LL)omega[i-1]*omega[1]%P,Sm[i]=(LL)Sm[i-1]*(omega[i]-1)%P;
Re tmp=Inv(Sm[K-1]);
for(Re i=K-1;i>=1;--i)invo[i]=(LL)tmp*Sm[i-1]%P,tmp=(LL)tmp*(omega[i]-1)%P;
for(Re i=0,mp1=0,mp2=0;i<K;++i)
mo(ans+=(LL)mi(omega[K-i]+1,K-1)*calc(mp1,n)%P*calc(mp2,m)%P*calc(i,l)%P),mo(mp1+=x2,K),mo(mp2+=x,K);
ans=(LL)ans*(K+1)%P;
printf("%lld\n",(LL)ans*Inv(n)%P*Inv(m)%P*Inv(l)%P);
}
}
【题外话—命题报告】
出题人：辰星凌
验题人：鏡音リン
（感谢神仙 \(\text{karry}\) 解惑，\(\text{ezoixx130}\) 划水验题）
\(7\) 月 \(28\) 晚，翻看自己整理的数学知识点汇总时，发现单位根反演这个东西题不多，就想搞一个出来玩玩儿。但数学题并不好出，直接放柿子的话会被喷惨。但我又比较菜，斗不出组合意义，顿时感觉希望渺茫。
抱着试一试的心态，随便选了个方向：求 \(\sum_{i,j}f(i\bmod j)\)，其中 \(f\) 先待定成普通多项式和下降幂多项式两种情况。推推推，最后弄出来一长串柿子，用了各种各样的技巧，结果最后还是 \(O(n^2)\) ！甚至常数还更大了！
人傻了....
后来调整数据范围，另外搞了个柿子，并斗出一个简单的 \(f\) 函数，消去了很多东西，似乎能出成一道题了。
但这还不够，因为只是给出了一个柿子让做题者去推，这样是没有意义的。而且单反就那两种变化，由于白兔之舞的存在，似乎无论怎么出都会感觉很套路，多半会被狂 \(\text{D}\) 。
那就....加上构造！
把 \(f\) 函数里的组合数拆开，弄出一个递推式，在题面上加点迷惑元素，让出题者去尝试构造 \(f\)，并在边界处设置一个坑人的细节问题。
这样看起来似乎就不那么板了(≥▽≤)/
（其实还是很套路，只要对组合数够熟悉就能一眼秒）。
另外，关于用 \(f\) 递推式求通项表示式，一开始是想用生成函数暴力解的，但貌似要求导积分，小蒟蒻对这一类的东西不太熟悉，如果有直接硬推成功解出来的巨佬，请接受蒟蒻真诚的膜拜stO
emm....从巨佬那里得知 \(f\) 可以在 \(\text{oeis}\) 上找到，自闭....
还是不加入到公开赛了，免得被骂搬原式。
\(\text{updata:}\) 经 \(\text{jklover}\) 巨佬提醒，\(\text{Subtask 3,4}\) 里的做法是可以直接用 \(\text{DFT}\) 优化的，总复杂度也为 \(O(n\log n)\)，虽然常数大一些，但实际运行效果只比前面的无常数单 \(log\) 做法慢 \(0.3s\)，过于柯啪.....
懒得写新 std 和重新造 Subtask 卡常了，就酱吧。
