来源
类欧几里德算法由洪华敦大佬在 2016 年 NOI 冬令营营员交流中提出,其本质可以理解为,一个类似辗转相除法的函数求和的过程。截止 2020 年 4 月 29 日,百度百科和维基百科都尚未收录该词条 Orz。
前置知识
取整函数
⌊ba⌋ 为 ba 向下取整,⌈ba⌉ 为 ba 向上取整。
基与下面的类欧几里得算法的推导过程,需要掌握一些取整函数的技巧。
⌊ba⌋⩾c⇔a⩾bc ⌈ba⌉⩽c⇔a⩽bc ⌊ba⌋<c⇔a<bc ⌈ba⌉>c⇔a>bc
⌊ba⌋=⌈ba+b+1⌉ ⌈ba⌉=⌊ba+b−1⌋
证明要略
当 amodb=0 时,上述式子显然都成立;当 amodb=0 时,可以令 a=kb+m ,其中 k,m∈N∗,1⩽m<b,可以轻易得到证明。
引入
类欧几里得算法主要用来计算形如 i=0∑n⌊cai+b⌋,i=0∑ni⌊cai+b⌋,i=0∑n⌊cai+b⌋2的表达式。
记 f(a,b,c,n)=i=0∑n⌊cai+b⌋,g(a,b,c,n)=i=0∑ni⌊cai+b⌋,h(a,b,c,n)=i=0∑n⌊cai+b⌋2。
其结果由递归函数得到。
由于推导过程类似欧几里得算法,时间复杂度同欧几里得算法一样也是 logn ,故称类欧几里得算法。
推导
推导 f(a,b,c,n)
f(a,b,c,n) 的形式最为简单,也最常出现;其几何意义为直角坐标系下,直线 y=cax+b 与 x,y 两轴包围区域内的整点个数。
当 a⩾c 或 b⩾c时,意味着可以将 a,b 对 c 取模简化问题。
当 a⩾c 或 b⩾c时
f(a,b,c,n)=i=0∑n⌊cai+b⌋=i=0∑n⌊c(⌊ca⌋c+amodc)i+(⌊cb⌋c+bmodc)⌋=i=0∑n[⌊ca⌋i+⌊cb⌋+⌊c(amodc)i+(bmodc)⌋]=i=0∑n(⌊ca⌋i+⌊cb⌋)+i=0∑n⌊c(amodc)i+(bmodc)⌋=2n(n+1)⌊ca⌋+(n+1)⌊cb⌋+f(amodc,bmodc,c,n)
那么问题就转化为 a<c 且 b<c 的情况。我们观察式子,发现只有 i 这一个变量,因此就只能从 i 下手。在推求和式时有一个常见的技巧,就是条件与贡献的放缩与转化。具体地说,在原式 f(a,b,c,n)=i=0∑n⌊cai+b⌋ 中, 0⩽i⩽n 是条件,而 ⌊cai+b⌋ 是对总和的贡献。要加快一个求和式的计算,所有的方法都可以都可以归约为 贡献合并计算。如果发现式子的贡献难以合并,就可以将贡献与条件作转化得到另一个形式的求和式。
请阁下看好接下来的神仙操作! Orz
f(a,b,c,n)=i=0∑n⌊cai+b⌋=i=0∑nj=1∑⌊cai+b⌋1
当 a<c 且 b<c 时
f(a,b,c,n)=i=0∑n⌊cai+b⌋=i=0∑nj=1∑⌊cai+b⌋1=i=0∑nj=1∑⌊can+b⌋[⌊cai+b⌋⩾j]=j=1∑⌊can+b⌋i=0∑n[⌊cai+b⌋⩾j]=j=1∑⌊can+b⌋i=0∑n[ai+b⩾cj]=j=1∑⌊can+b⌋i=0∑n[i⩾⌈acj−b⌉]=j=1∑⌊can+b⌋i=⌈acj−b⌉∑n1=j=1∑⌊can+b⌋(n−⌈acj−b⌉+1)=j=1∑⌊can+b⌋(n−⌊acj+a−b−1⌋+1)=n⌊can+b⌋−j=1∑⌊can+b⌋(⌊acj+a−b−1⌋−1)=n⌊can+b⌋−j=1∑⌊can+b⌋⌊acj−b−1⌋=n⌊can+b⌋−j=0∑⌊can+b⌋−1⌊ac(j+1)−b−1⌋=n⌊can+b⌋−j=0∑⌊can+b⌋−1⌊acj+c−b−1⌋=n⌊can+b⌋−f(c,c−b−1,a,⌊can+b⌋−1)
边界条件为,当 a=0 时,fa,b,c,n)=(n+1)⌊cb⌋。
推导 g(a,b,c,n)
当 a⩾c 或 b⩾c 时
g(a,b,c,n)=i=0∑ni⌊cai+b⌋=i=0∑ni[⌊c(amodc)i+(bmodc)⌋+⌊ca⌋i+⌊cb⌋]=i=0∑ni⌊c(amodc)i+(bmodc)⌋+i=0∑n⌊ca⌋i2+i=0∑n⌊cb⌋i=6n(n+1)(2n+1)⌊ca⌋+2n(n+1)⌊cb⌋+i=0∑ni⌊c(amodc)i+(bmodc)⌋=6n(n+1)(2n+1)⌊ca⌋+2n(n+1)⌊cb⌋+g(amodc,bmodc,c,n)
当 a<c 且 b<c 时
g(a,b,c,n)=i=0∑ni⌊cai+b⌋=i=0∑nj=1∑⌊cai+b⌋i=i=0∑nj=1∑⌊can+b⌋i[⌊cai+b⌋⩾j]=i=0∑nj=1∑⌊can+b⌋i[ai+b⩾cj]=i=0∑nj=1∑⌊can+b⌋i[i⩾⌈acj−b⌉]=i=0∑nj=1∑⌊can+b⌋i[i⩾⌊acj+a−b−1⌋]=j=1∑⌊can+b⌋i=0∑ni[i⩾⌊acj+a−b−1⌋]=j=1∑⌊can+b⌋i=⌊acj+a−b−1⌋∑ni=j=1∑⌊can+b⌋2(n+⌊acj+a−b−1⌋)(n−⌊acj+a−b−1⌋+1)=j=1∑⌊can+b⌋2(n+⌊acj−b−1⌋+1)(n−⌊acj−b−1⌋)=j=0∑⌊can+b⌋−12(n+⌊ac(j+1)−b−1⌋+1)(n−⌊ac(j+1)−b−1⌋)=j=0∑⌊can+b⌋−12(n+⌊acj+c−b−1⌋+1)(n−⌊acj+c−b−1⌋)=21j=0∑⌊can+b⌋−1[n(n+1)−⌊acj+c−b−1⌋−⌊acj+c−b−1⌋2]=21j=0∑⌊can+b⌋−1n(n+1)−j=0∑⌊can+b⌋−1⌊acj+c−b−1⌋−j=0∑⌊can+b⌋−1⌊acj+c−b−1⌋2=21[⌊can+b⌋n(n+1)−f(c,c−b−1,a,⌊can+b⌋−1)−h(c,c−b−1,a,⌊can+b⌋−1)]
边界条件为,当 a=0 时,g(a,b,c,n)=2n(n+1)⌊cb⌋ 。
推导 h(a,b,c,n)
当 a⩾c 或 b⩾c 时
h(a,b,c,n)=i=0∑n⌊cai+b⌋2=i=0∑n[⌊c(amodc)i+(bmodc)⌋+⌊ca⌋i+⌊cb⌋]2=i=0∑n⌊c(amodc)i+(bmodc)⌋2+i=0∑n⌊ca⌋2i2+i=0∑n⌊cb⌋2+2⌊ca⌋i=0∑n⌊c(amodc)i+(bmodc)⌋i+2⌊cb⌋i=0∑n⌊c(amodc)i+(bmodc)⌋+2⌊ca⌋⌊cb⌋i=0∑ni=h(amodc,bmodc,c,n)+2⌊ca⌋g(amodc,bmodc,c,n)+2⌊cb⌋f(amodc,bmodc,c,n)+6n(n+1)(2n+1)⌊ca⌋2+(n+1)⌊cb⌋2+n(n+1)⌊ca⌋⌊cb⌋
当 a<c 且 b<c 时
∀x∈N+,有如下骚操作 x2=2×2x(x+1)−x=2i=1∑xi−x
h(a,b,c,n)=i=0∑n⌊cai+b⌋2=i=0∑n2j=1∑⌊cai+b⌋j−⌊cai+b⌋=2i=0∑nj=1∑⌊cai+b⌋j−i=0∑n⌊cai+b⌋=2i=0∑nj=1∑⌊can+b⌋j[⌊cai+b⌋⩾j]−i=0∑n⌊cai+b⌋=2i=0∑nj=1∑⌊can+b⌋j[ai+b⩾cj]−i=0∑n⌊cai+b⌋=2i=0∑nj=1∑⌊can+b⌋j[i⩾⌈acj−b⌉]−i=0∑n⌊cai+b⌋=2j=1∑⌊can+b⌋ji=0∑n[i⩾⌈acj−b⌉]−i=0∑n⌊cai+b⌋=2j=1∑⌊can+b⌋ji=0∑n[i⩾⌊acj+a−b−1⌋]−i=0∑n⌊cai+b⌋=2j=1∑⌊can+b⌋[j(n−⌊acj+a−b−1⌋+1)]−i=0∑n⌊cai+b⌋=2j=1∑⌊can+b⌋[j(n−⌊acj−b−1⌋)]−i=0∑n⌊cai+b⌋=2j=0∑⌊can+b⌋−1[(j+1)(n−⌊ac(j+1)−b−1⌋)]−i=0∑n⌊cai+b⌋=2j=0∑⌊can+b⌋−1[(j+1)(n−⌊acj+c−b−1⌋)]−i=0∑n⌊cai+b⌋=2nj=0∑⌊can+b⌋−1(j+1)−2j=0∑⌊can+b⌋−1j⌊acj+c−b−1⌋−2j=0∑⌊can+b⌋−1⌊acj+c−b−1⌋−i=0∑n⌊cai+b⌋=⌊can+b⌋(⌊can+b⌋+1)n−2g(c,c−b−1,a,⌊can+b⌋−1)−2f(c,c−b−1,a,⌊can+b⌋−1)−f(a,b,c,n)
边界条件为,当 a=0 时,h(a,b,c,n)=(n+1)⌊cb⌋2
题目描述
给定 n,,a,b,c,分别求 i=0∑n⌊cai+b⌋,i=0∑n⌊cai+b⌋2,i=0∑ni⌊cai+b⌋ ,答案对 998244353 取模。多组数据。
输入格式
第一行给出数据组数 t 。
接下来 t 行,每行有四个整数,分别为每组数据的 n,a,b,c 。
输出格式
对于每组数据,输出一行三个整数,为三个答案对 998244353 取模的结果。
输入输出样例
输入 #1
2
2 1 0 2
4 3 9 6
输出 #1
1 1 2
11 27 27
说明/提示
本题采用 Special Judge 。
答对所有第一问可以获得测试点 40% 的分数,答对所有第二问、第三问可以分别获得另外 30% 的分数。
测试点编号 |
特殊性质 |
1 |
n,a,b,c⩽10 |
2∼3 |
n,a,b,c⩽100 |
4 |
n,a,b,c⩽106 |
5 |
t=1 |
6∼10 |
无 |
对于所有测试点,有 1⩽t⩽105, 0⩽n,a,b,c⩽109,c=0 。
分析
在计算的时侯,因为三个函数各有交错递归,因此可以考虑三个函数打包成结构体一起整体递归,同步计算,否则有很多项会被多次计算。
类欧几里得算法时间复杂度为 O(logn) 。
代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
| #include<iostream> #include<cstdio> #define ll long long using namespace std; const ll mod=998244353;
const ll inv2=499122177; const ll inv6=166374059;
struct A { ll f; ll g; ll h; }; A cal(ll a,ll b,ll c,ll n) { A res,temp; if(!a) { res.f=(n+1)*(b/c)%mod; res.g=(b/c)*n%mod*(n+1)%mod*inv2%mod; res.h=(n+1)*(b/c)%mod*(b/c)%mod; } else if(a>=c||b>=c) { temp=cal(a%c,b%c,c,n); res.f=((a/c)*n%mod*(n+1)%mod*inv2%mod +(b/c)*(n+1)%mod+temp.f)%mod; res.g=((a/c)*n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod +(b/c)*n%mod*(n+1)%mod*inv2%mod+temp.g)%mod; res.h=((a/c)*(a/c)%mod*n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod +(b/c)*(b/c)%mod*(n+1)%mod+(a/c)*(b/c)%mod*n%mod*(n+1)%mod +temp.h+2*(a/c)%mod*temp.g+2*(b/c)%mod*temp.f)%mod; } else { ll m=(a*n+b)/c; temp=cal(c,c-b-1,a,m-1); res.f=(n*(m%mod)%mod-temp.f)%mod; res.g=(n*(n+1)%mod*(m%mod)%mod-temp.f-temp.h)%mod*inv2%mod; res.h=(n*(m%mod)%mod*((m+1)%mod)%mod-2*temp.g-2*temp.f-res.f)%mod; } return res; } int main() { int _; for(cin>>_;_;_--) { ll n,a,b,c; scanf("%lld%lld%lld%lld",&n,&a,&b,&c); A ans=cal(a,b,c,n); printf("%lld %lld %lld\n",(ans.f+mod)%mod,(ans.h+mod)%mod,(ans.g+mod)%mod); } return 0; }
|