🎢 锯齿过山车·续：数量飙十亿，矩阵快速幂！
上回说到，我用前缀和加原地滚动，把计算锯齿过山车轨道方案数的复杂度从 O(nk²) 压到了 O(nk)，村长乐颠颠地去游乐园领赏金了。老勇者临走时扔下一句：“如果n飙到10^6，前缀和也不够用，那就得上矩阵快速幂了。”
我本以为这只是老勇者随口一提的“远期预告”，毕竟谁能闲着没事建一条十亿个节点的过山车？
没想到今天一大早，村长灰头土脸地就撞了进来：“小白！！” 他双手撑在桌子上，眼神里写满了绝望：“游乐园没给奖金！他们说上次的方案只适用于小规模轨道，现在要建一条真正的‘超级锯齿过山车’！轨道节点数量从 2000 暴涨到……到……”
他颤抖着伸出手指，声音都在劈叉：“十亿！十亿个节点！从村头一直铺到海岸线上的那种！你上次那套 O(n×k) 算法直接跪在地上摩擦，连渣都不剩啊！”
我接过他手里皱巴巴的需求文档，还是那座熟悉的锯齿过山车，但设计图上的参数变了——
锯齿过山车·超长版
- 轨道点数 n：3 ≤ n ≤ 10^9（昨天才 2000，今天直接通货膨胀？）
- 高度范围 [l, r]：1 ≤ l < r ≤ 75（昨天是 2000，今天反而缩水了？）
- 其他规则不变：相邻不等，连续三个不能严格递增或递减
- 目标：求方案总数，对 10^9+7 取模
10^9，**十亿次循环。**就算每次只做一次加法，O(nk) 在 n=10^9、k=75 时也是 750 亿次运算！就算评测机是超算，每秒处理 10^9 次，也得跑 75 秒！况且它的耐心只有 2 秒——暴力是死路，前缀和也是死路——前缀和再快也甩不掉外层那十亿次迭代啊！
“矩阵快速幂……“我喃喃自语，脑子里回放着老勇者昨天的话，“能把 O(n) 的线性递推优化成 O(log n)……”
log₂(10^9) ≈ 30。也就是说，如果能用上矩阵快速幂，原本需要十亿次的迭代，现在只需要 30 次！这压缩比，简直是用了量子传送门直接从起点跳到终点！
矩阵快速幂？听起来像是高维魔法
我深吸一口气，强迫自己冷静下来，开始拆解这个听起来很唬人的概念。
矩阵快速幂，名字听着玄乎，像是什么 SSR 级角色的终极技能。但剥开外壳看本质，它其实就是把快速幂的思想从标量扩展到矩阵，从而把 n 次矩阵乘法压缩到 log n 次。
先来个快速幂的极简回顾，免得待会儿在矩阵的迷宫里迷路：
快速幂的核心思想：计算 a^n 不需要傻乎乎地乘 n 次，而可以利用二进制分解。
比如要计算 a^13，把 13 写成二进制形式：13 = 1101(₂)。这意味着：a^13 = a^(8+4+1) = a^8 × a^4 × a^1
具体做法就像是在玩"翻倍游戏”：
- 初始化结果
res = 1，底数base = a- 从低位到高位遍历 n 的二进制位
- 如果当前位是 1，就把当前的
base累乘到res中- 每步让
base自乘（base = base × base），n 右移一位通过二进制分解，只需要计算 log n 次乘法，而不是 n 次。 时间复杂度从 O(n) 降到 O(log n)，压缩比，简直像是用了空间折叠技术！
def pow_mul(a, n):
res = 1
while n:
if n & 1: # 如果当前二进制位是 1
res *= a # 累乘到结果
a *= a # 底数自乘
n >>= 1 # n 右移一位
return res
温习完毕，再回到这次的任务上来。
我的 DP 转移是线性的——每一步都只做加法和数乘，没有非线性操作（比如取 max、min 或者条件判断）。这意味着，整个转移过程完全可以表示成一个矩阵乘法！
如果我能构造一个转移矩阵 M 和一个初始状态向量 F₁，使得：
F₂ = M × F₁ F₃ = M × F₂ = M² × F₁ … Fₙ = M^(n-1) × F₁
那么我只需要计算 M^(n-1)，然后乘以 F₁，就能直接得到第 n 层的状态！而 M^(n-1) 可以用快速幂在 O(log n) 次矩阵乘法内完成！
我在纸上又画了个对比图：
普通 DP: F₁ → F₂ → F₃ → … → Fₙ (n-1 步，O(n) 时间) 矩阵快速幂: F₁ ⇒ M^(n-1)×F₁ = Fₙ (log n 次矩阵乘法，O(log n) 时间)
这简直像是给你的 DP 装上了火箭推进器——原本要一步步走的 n 步，现在只要 log n 次"瞬移"就能到达！从村头到海岸线，以前要走十亿步，现在只要 30 次传送门跳跃！
思路通了！我兴奋地掏出键盘，准备大干一场。
第一步：把 DP 状态转化成向量形式
“首先，我需要把 DP 状态转化成向量形式。“我自言自语道，手指在桌面上敲着节拍。
“昨天的状态是 f0[j] 和 f1[j]，分别表示最后一个元素是 j 且最后两个元素递增/递减的方案数。那我可以把这些状态拼成一个长度为 2k 的列向量……”
我在草稿纸上写下状态向量的结构：
状态向量 F = [f0[0], f0[1], …, f0[k-1], f1[0], f1[1], …, f1[k-1]]ᵀ
长度为 2k，当 k=75 时就是 150 维。这个维度完全可控，既不会大到内存爆炸，也不会小到失去表达能力。完美！
第二步：构造转移矩阵 M
“然后，我需要构造一个 2k × 2k 的转移矩阵 M，使得 F_{i+1} = M × F_i。”
我盯着纸上那个 150×150 的空矩阵格子，感觉像是在面对一个巨大的数独谜题。怎么填？
先回顾一下昨天的状态转移公式：
对于 f0_new[j]（以 j 结尾且最后两步递增）： f0_new[j] = sum(f1_old[p]) where p < j 对于 f1_new[j]（以 j 结尾且最后两步递减）： f1_new[j] = sum(f0_old[p]) where p > j
这意味着，F{i+1} 中的每一个分量，都依赖于 F_i 中多个分量的累加和。
我皱起眉头：“每个新状态都依赖于多个旧状态的累加，这怎么用矩阵表示？区间求和的关系怎么转化成矩阵的系数？这难道要我给每个位置都写一个求和公式？”
就在我抓耳挠腮之际，脑海中突然闪过一道灵光——矩阵乘法的本质就是’加权求和’！
当我计算 M × F 时，结果向量的第 i 个元素等于：
(M × F)[i] = Σ M[i] [j] × F[j] （对所有 j 求和）
这不就是矩阵 M 的第 i 行与向量 F 的点积吗？
也就是说，如果我想让：
f0_new[j] = sum(f1_old[k]) where k < j
我只需要在矩阵 M 的第 j 行，把所有对应 f1_old[k]（k < j）的位置填成 1，其他位置填 0！这样，当矩阵乘法执行时，这些 1 就会把对应的 f1_old[k] “挑选"出来并累加，而 0 会把不相关的元素过滤掉！
这就像是给矩阵装上了过滤器——只让需要的元素通过，其他的统统拦截！
我激动地在草稿纸上画出了转移矩阵的结构：
对于f0[j] = f1[0] + f1[1] + … + f1[j-1] ， 转移矩阵里对应的那部分是一个下三角全1矩阵。 对于f1[j] = f0[j+1] + f0[j+2] + … + f0[k-1]，对应的是一个上三角全1矩阵。
转移矩阵 M 的结构（2k × 2k）：
| 0 B | ← 上半部分：f0_new 依赖 f1_old | C 0 | ← 下半部分：f1_new 依赖 f0_old
其中： B 是 k×k 的下三角全 1 矩阵（B[i] [j]=1 当且仅当 j < i） C 是 k×k 的上三角全 1 矩阵（C[i] [j]=1 当且仅当 j > i） 0 是 k×k 的全零矩阵
我盯着这个分块反对角矩阵，感觉自己像是破解了一个古老的密码。原来，看似复杂的区间求和，在矩阵的世界里，不过是填 0 和 1 的游戏罢了！
# 转移矩阵 M，大小2k × 2k
m = [[0] * (k*2) for _ in range(k*2)]
# 上半部分：f0[i] 依赖所有 f1[j] where j < i
for i in range(k):
for j in range(i):
m[i][k+j] = 1 # 第 i 行，第 k+j 列为 1
# 下半部分：f1[i] 依赖所有 f0[j] where j > i
for i in range(k):
for j in range(i+1, k):
m[k+i][j] = 1 # 第 k+i 行，第 j 列为 1
第三步：组装矩阵快速幂——把零件拼成火箭
想清楚矩阵的构造，剩下的就是体力活了。我需要实现三个组件：
- 矩阵乘法
mul：计算两个 2k × 2k 矩阵的乘积 - 矩阵快速幂
pow_mul：计算 M^(n-1) × F₁ - 主函数
zigZagArrays：构造矩阵、调用快速幂、统计结果
我开始写完整的代码，就像是在组装一台精密的仪器：
MOD = 10**9 + 7
def mul(a: List[List[int]], b: List[List[int]]) -> List[List[int]]:
"""矩阵乘法：C = A × B"""
return [[sum(x*y for x, y in zip(row, col))%MOD for col in zip(*b)]
for row in a]
def pow_mul(a: List[List[int]], n: int, f1: List[List[int]]) -> List[List[int]]:
"""矩阵快速幂：计算 A^n × F1"""
res = f1
while n:
if n & 1:
res = mul(a, res)
a = mul(a, a)
n >>= 1
return res
class Solution:
def zigZagArrays(self, n: int, l: int, r: int) -> int:
k = r - l + 1 # 高度范围的长度
# 构造转移矩阵 M，大小 2k × 2k
m = [[0] * (k*2) for _ in range(k*2)]
for i in range(k):
for j in range(i):
m[i][k+j] = 1
for i in range(k):
for j in range(i+1, k):
m[k+i][j] = 1
# 初始状态向量 F1，所有值为 1（第一个元素可以是任意值）
f1 = [[1] for _ in range(k*2)]
# 计算 M^(n-1) × F1
fn = pow_mul(m, n-1, f1)
# 统计所有状态的方案数之和
return sum(row[0] for row in fn) % MOD
小样例先跑下，完美通过！再算算复杂度：矩阵乘法 O((2k)³) = O(8k³) ≈ O(k³)。快速幂需要 log n 次矩阵乘法，所以总复杂度是 O(k³ log n)。代入数值：k=75, n=10^9，log₂(10^9) ≈ 30。75³ × 30 ≈ 1.2 × 10^7 次运算。也就是说，总操作量大约一千两百万次。怎么说呢，这个量级就像是早高峰的地铁——挤是挤了点，但只要不被卡住，还是能在规定时间内到达终点。
我深吸一口气，点击提交。屏幕闪烁了两秒——绿色的Accepted弹了出来！运行时间1.8秒，差点就超时了。
“过了！”我差点从椅子上跳起来，“虽然跑得慢了点儿，但能在2秒内ac，就是胜利！” 十亿次迭代被压缩到三千万次矩阵运算，这压缩比简直像是把一整条过山车轨道卷成了一个弹簧！
就在我得意洋洋，准备发朋友圈炫耀这波“极限压线”操作时，后脑勺又挨了精准的一下——力道刚好能把我从“我真牛逼”的幻觉里敲回现实。
老勇者登场——优化的边界
“不错嘛，”老勇者站在我身后，嘴角带着一丝笑意，“学会用朴素矩阵快速幂 AC 了。1.8秒，压线过——看来评测机今天心情不错，没卡你的常数因子。”
我揉着后脑勺，有点不好意思地嘿嘿一笑：“运气好，运气好。不过您看这1.8秒，心悬得很。有没有什么优化方法，能让它过得更稳一点？”
老勇者点了点屏幕：“你的解法很标准了——构造转移矩阵，实现矩阵快速幂，计算 M^(n-1) × F₁。这套流程你掌握得很扎实。矩阵快速幂的适用条件、转移矩阵的构造技巧，这两个核心点你都踩准了。”
他顿了顿，扫了一眼我屏幕上那个矩阵乘法的代码：“不过，既然你问到了优化，那我就提供一个思路。”
“你的转移矩阵是分块反对角矩阵，一半下三角全1，一半上三角全1。这种特殊结构，乘向量的时候其实不需要 O(k²) 的三重循环——下三角全1乘向量就是做前缀和，上三角全1乘向量就是做后缀和。跟昨天的原地滚动一个道理。一趟 O(k) 搞定，不用构建完整的转移矩阵。
我脑子“嗡”了一下，有点懵：“您的意思是……矩阵快速幂里每一轮乘向量都能压到 O(k)，总复杂度压到 O(k² log n)？”
“对。”老勇者点点头，话锋却一转，“但我只是提供这个思路，不要求你现在就去实现。”
我愣了一下，满脸写着“为什么”：“为什么？这不是能让代码跑得更快吗？”
“因为你已经 AC 了。”老勇者瞥了我一眼，“1.8秒和0.5秒，在判题机眼里都是绿色的‘通过’。但如果现在去实现这个优化，你得把矩阵乘法拆成两种情况——矩阵乘矩阵和矩阵乘向量。代码量翻倍，边界条件翻倍，调试时间也翻倍。为了省这1.3秒，把原本清晰的代码搅成一锅粥，不值得。”
他敲了敲我的屏幕：“刷题不是打比赛，AC 之后多花一小时抠常数，不如把这一小时用来消化今天学到的矩阵快速幂骨架。优化是锦上添花，骨架才是雪中送炭。k=75 的时候 O(k³ log n) 够用，下次 k 飙到 500，再掏出前缀和优化不迟——到时候你自然会感谢今天先把骨架练熟了。”
尾声：
我靠在椅背上，盯着屏幕上那个绿色的 Accepted，紧绷的肩膀终于放松下来。
为了凑出那个 1.8 秒的极限压线，我敲代码时手心全是汗，生怕评测机稍微打个盹就给我个 TLE。现在想想，从十亿次循环压到三十次矩阵乘法，已经是一波酣畅淋漓的战斗了。至于前缀和优化矩阵乘法——遇到k飙到500的时候再说，毕竟骨架学透才是正经事——
十亿节点莫惊慌，线性递推矩阵扛。
快速幂里藏乾坤，有缘再战五百强。
本作品采用
知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议
进行许可。