【题解】【原创题目】第〇类循环数•行

【题解】【原创题目】第〇类循环数·行

传送门:第〇类循环数·行

【题目描述】

给出 \(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)\) 定义如下:

\[\begin{cases}i^2+i&i\in[1,\infty],j=0\\\frac{f(i-1,j-1)}{ij}+\frac{f(i-1,j-1)+[j<i-1]f(i-1,j)}{i}+\frac{f(i-1,j-1)}{j}+f(i\!-\!1,j\!-\!1)+[j\!\neq\!i\!-\!1]f(i\!-\!1,j)&i\in[2,\infty],j\in[1,i\!-\!1]\end{cases} \]

【分析】

【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\) 递推式,首先通分合并同类项:

\[f(i,j)=\frac{(i+1)(j+1)}{ij}f(i-1,j-1)+[j\neq i-1]\frac{(i+1)}{i}f(i-1,j) \]

两边同时除以 \((i+1)(j+1)\)

\[\frac{f(i,j)}{(i+1)(j+1)}=\frac{f(i-1,j-1)}{ij}+[j\neq i-1]\frac{f(i-1,j)}{i(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\) 的参数并交换和式:

\[\frac{1}{nml}\sum_{d=0}^{k-1}f(k,d)\sum\limits_{a=1}^{n}\sum\limits_{b=1}^{m}\sum\limits_{c=1}^{l}[(ax^2+bx+c)\bmod k=d] \]

\([(ax^2+bx+c)\bmod k=d]\) 可化为 \([k|(ax^2+bx+c-d)]\),套用单位根反演的柿子,得到:

\[\frac{1}{nml}\sum_{d=0}^{k-1}f(k,d)\sum\limits_{a=1}^{n}\sum\limits_{b=1}^{m}\sum\limits_{c=1}^{l}\frac{1}{k}\sum_{i=0}^{k-1}w_{k}^{i(ax^2+bx+c-d)} \]

提一个 \(w_{k}^{-id}\) 出来,并交换和式:

\[\frac{1}{nmlk}\sum_{i=0}^{k-1}\sum_{d=0}^{k-1}f(k,d)w_{k}^{-id}\sum\limits_{a=1}^{n}\sum\limits_{b=1}^{m}\sum\limits_{c=1}^{l}w_{k}^{i(ax^2+bx+c)} \]

右边那一坨显然可以简化为三个等比数列和的乘积,记一个 \(calc\) 函数:

\[calc(w)=\sum\limits_{a=1}^{n}w^{ax^2}\sum\limits_{b=1}^{m}w^{bx}\sum\limits_{c=1}^{l}w^{c}=\frac{w^{x^2}(1-(w^{x^2})^{n})}{1-w^{x^2}}\times\frac{w^{x}(1-(w^{x})^{m})}{1-w^{x}}\times\frac{w^{}(1-w^{l})}{1-w^{}} \]

用组合数吸收公式去掉 \(f\) 柿子里的 \((j+1)\),即 \(f(i,j)=(i+1)(j+1)\dbinom{i}{j+1}=i(i+1)\dbinom{i-1}{j}\),然后代进去用二项式定理合并:

\[\frac{1}{nmlk}\sum_{i=0}^{k-1}\sum_{d=0}^{k-1}k(k+1)\dbinom{k-1}{d}w_{k}^{-id}calc(w_{k}^{i})\\ =\frac{k+1}{nml}\sum_{i=0}^{k-1}calc(w_{k}^{i})\sum_{d=0}^{k-1}\dbinom{k-1}{d}w_{k}^{-id}\\ =\frac{k+1}{nml}\sum_{i=0}^{k-1}calc(w_{k}^{i})(w_{k}^{-i}+1)^{k-1} \]

本测试包不卡常,可以直接暴算。

复杂度 \(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 卡常了,就酱吧。


-----若要转载请私信作者获得许可并在文首标出转载来源-----

posted @ 2020-08-03 18:12  辰星凌  阅读(341)  评论(0编辑  收藏  举报