万能欧几里得学习笔记
万欧很好的点出了类欧的本质,告诉我们推那一串式子和推两个操作并没有任何区别。
目录:
- 万欧的概念、套路与代码实现
- 万欧解决的问题形式及例题
- 参考资料
万欧的概念、套路与代码实现
有一些概念性的东西:
定义 操作 为修改某个变量状物,它满足 结合律,且两个操作可压缩为一个操作,下文提到的 U 和 R 还有 solve 的返回值都是操作,定义操作之间的运算 a×b 或 ab 得到的操作就是先执行 a 再执行 b 这个操作,类似的我们可以定义操作的自然数次幂为重复执行该操作。
万能欧几里得在笛卡尔坐标系上考虑 y=qpx+b 这条直线,我们规定向右走的过程中每与横坐标相交一次就执行 U 操作,每与纵坐标相交一次就执行 R 操作,同时相交先执行 U 再执行 R。
只要我们把原问题转化为在平面上执行操作这件事那就可以用万欧秒了。
我们需要在这样一个框架下思考:solve(p,q,b,n,U,R) 为处理 x∈(0,n] 的 y=qpx+b 这条线段,U,R 表示遇到横线和遇到竖线之后的操作,该函数返回的也是一个(压缩后的)操作。
仔细想欧几里得算法的过程:a=0 时结束,a≥b 时变为 (amodb,b),a<b 时变为 (b,a)。
那我们也就仿照这个过程来。
欸你发现这个时候不会执行 U(即使与横线重合也不认为“遇到”横线),于是直接返回 Rn。
不难发现此时斜率 ≥1,每往右走一步就会遇到至少 ⌊qp⌋ 个横线,于是我们可以把这 ⌊qp⌋ 个 U 和一个 R 合并起来,递归 solve(pmodq,q,b,n,U,U⌊qp⌋×R)。
为什么“至少”可以等价为“必定”,为什么 p→pmodq?
其实这俩是一个问题啦,就是因为有可能有超过 ⌊qp⌋ 个 U,所以我们把整点部分取出来,剩下的 pmodq 多走一步还是会多一个 U 的。
此时想要继续分析稍有困难,但根据欧几里得算法的过程,不难猜到我们需要用某种方式交换 p,q,我们需要将其按照 y=x 对称一下,即直线变为 y=pqx−b,同时需要 (U,R)′←(R,U)。
但朴素的翻转会带来一些问题:
- 此时经过整点时变成了先加 R 再加 U,这与我们规定好的顺序不符。
- 翻折前的问题从 y 轴交点开始计算,那么翻折后的问题理应从 x 轴交点开始计算,但 solve 的定义为处理 (0,n] 的问题。
- 原问题 n 处可能并不是一个整点,也就是说,翻转后 n 可能不是一个整数。
我们一个一个来解决,首先是先加 R 再加 U 的问题,我们来分析一些事实:
已知第 a 个 R 前面有 ⌊qap+b⌋ 个 U,那么我们设第 a 个 R 在第 c 个 U 前面,看能否得到第 c 个 U 前面有多少个 R。
(由于此时 p<q,所以第 c 个 U 就是 y=c。)
推导过程
ccapaaa>⌊qap+b⌋>qap+b<qc−b<pqc−b<⌈pqc−b⌉=⌊pqc−b+p−1⌋=⌊pqc−b−1⌋−1≤⌊pqc−b−1⌋
于是我们知道了第 c 个 U 前面有 ⌊pqc−b−1⌋ 个 R,那我们如果把坐标系按照 y=x 对称的话,第 c 个 U 其实就是 x=c,所以实际的直线应该是 y=pqx−b−1。
第二个问题:我们的问题对称过来理应从与 x 轴交点开始执行,但直接这样写变成了从纵截距开始。
首先原问题的 y=qpx+b 不难看出 b 其实是在 mod q 意义下的,因为把 b 增加 q 操作序列不会有任何变化,我们不妨认为 b∈[0,q)。
那你发现 f(x)=pqx−b−1,f(1)≥0,换句话说我们只要把第一个 U 及其前面的 R 全部处理掉,再从第一位开始递归即可。
第三个问题:新的 n 不是整数。
这个很好办,仿照第二个问题,把最后一个 U 后面的 R 单独挑出来执行即可。
接下来我们明确一下细节:由前面的推导第一个 U 前面有 ⌊pq−b−1⌋ 个 R,新的 f(x)=pqx+q−b−1,共有 cntU=⌊qpn+b⌋ 个 U,最后一个 U 后面有 cntR=n−⌊pqcntU−b−1⌋ 个 R。
故递归到:
R⌊pq−b−1⌋×U×solve(q,p,(q−b−1)modp,cntU−1,R,U)×RcntR
代码实现:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
| template <typename T> T ksm(T a,ll n) { T res{}; for(;n;n>>=1,a=a*a) if(n&1) res=res*a; return res; } template <typename T> T sol(ll p,ll q,ll b,ll n,T U,T R) { if(!n) return T{}; if(p>=q) return sol(p%q,q,b,n,U,ksm(U,p/q)*R); ll cntU=(p*n+b)/q; if(!cntU) return ksm(R,n); return ksm(R,(q-b-1)/p)*U* sol(q,p,(q-b-1)%p,cntU-1,R,U)* ksm(R,n-((__int128)q*cntU-b-1)/p); }
|
注:由于推导过程使用了字母 p,所以代码里喜欢用 mod p 的同学可能要注意一下不要在 solve 中使用对答案的取模操作,可以封装在结构体内。
万欧解决的问题形式及例题
一般如果题目中出现了 ⌊qpx+b⌋ 那就要警觉了,剩下的与万欧关系不大。
对于常数 p,q,b,n,快速求 i=0∑n⌊qpi+b⌋,i=0∑n⌊qpi+b⌋2,i=0∑ni⌊qpi+b⌋。
T≤105,n,a,b,c,d≤109。
第一个最菜,先向第一个咋做。
来套用一下我们的概念,我们想要在平面上定义两个可合并的操作,仅此而已。
欸不难发现,遇到一个横线我们把 c 自增,遇到一个竖线我们把答案加上 c 就好了,很完美。
于是就可以把操作写成对向量 [c1ans] 的线性变换(矩阵)。
U=110010001,R=100010101
通常矩阵的常数是大到离谱的,所以我们也可以利用线段树那套思考标记信息的方式(等价于把矩阵的所有转移拆出来)来思考:需要维护什么?
维护 (k0,k,a) 表示给答案加了 k0 的常数,又给答案加了 k 倍的 c,最后把 c 加上了 a:
(k0,k,a)×(K0,K,A)→(k0+K0+Ka,k+K,a+A)
你看后面这个形式看起来就常数很小。
然后是后面两个问题,平方也是好做的,思考区间平方和怎么做。
维护 (k0,k1,k2,a) 表示给答案加了 k0 的常数,加了 k1 倍的 c,加了 k2 倍的 c2,最后把 c 加上了 a:
(k0,k1,k2,a)×(K0,K1,K2,A)→(k0+K0+K1a+K2a2,k1+K1+2K2a,k2+K2,a+A)
然后是 i=0∑ni⌊qpi+b⌋。
(其实矩阵想这些东西都超级好想,场上肯定要先写矩阵做法再慢慢优化。)
维护 (k0,ki,kc,kic,ai,ac) 表示给答案加了 k0 的常数,加了 ki 倍的 i,加了 kc 倍的 c,加了 kic 倍的 ic,最后把 i 加了 ai,把 c 加了 ac。
→(k0,ki,kc,kic,ai,ac)×(K0,Ki,Kc,Kic,Ai,Ac)k0+K0+Kiai+Kcac+Kicaiac,ki+Ki+Kicac,kc+Kc+Kicai,kic+Kic,ai+Ai,ac+Ac
这玩意要是写成矩阵后果不堪设想。
本题其实可以三个一块算(ix35 的题解),但是通常用到的时候都是求单个,且那种方式不如直接推直观(虽然能少几个变量),所以这里就不涉及了。
注意本题是从 0 开始求和,最后处理一下即可,还有求解 i=0∑ni⌊qpi+b⌋ 的时候其实应该先加 i 再加 ans,很遗憾我们的标记系统先加了 ans,所以初始应该把 i 加一。
实现没什么亮点,就不放了。
类似的东西其实都比较能做,只需要思考矩阵设计的状态即可。
求:
x=1∑LAxB⌊qpx+b⌋
其中 A,B 是 n×n 的矩阵。
L≤1018,n≤20。
你发现其实没什么新意,把矩阵当作一个数,其实就是动态维护 [AxB⌊qpx+b⌋,ans],你需要实现对前者的 ×A,×B,对后者加前者,这都弱爆了,矩阵可以直接做。
那来一道比较有新意的。
当时看懂了,但是没记下来,再次看到已经退役了/kk。
参考资料
万能欧几里得 学习笔记 - ycx060617
万能欧几里得还是很万能的! - Querainy