算法竞赛题解：大数运算与子矩形面积和的数学推导

1. 从一道竞赛题说起当“面积和”遇上“大数”最近在复盘一些算法竞赛的题目特别是那种看起来公式简单但数据范围大到让人头皮发麻的题。2020年牛客国庆集训排队Day2的这道F题——“SUM OF SUB RECTANGLE AREAS”就是典型代表。题目大意是给定一个N x N的网格计算所有可能的子矩形面积之和。听起来是不是挺直观暴力枚举所有子矩形把长宽乘积累加起来就行了。但问题就出在这个N上它可以非常大题目里明确说了结果可能超出long long的范围需要处理大数。这题的核心远不止是一个大数模板的简单调用。它真正的魅力在于你需要先通过数学推导找到一个高效的计算公式避免O(N^4)的暴力枚举。而这个推导过程恰恰是一个经典的“找规律”与“公式化简”案例涉及等差数列求和以及平方和公式的灵活运用。很多人在推导出公式S [N*(N1)/2]^2后就觉得万事大吉直接用大数一算提交却可能忽略了中间计算过程的溢出风险或者对大数运算的性能细节缺乏考量。今天我们就来彻底拆解这道题不仅讲清楚规律怎么找、公式怎么来更重点聊聊在大数场景下如何稳健、高效地实现这个计算。2. 规律探寻所有子矩形面积之和的公式推导为什么不能暴力计算对于一个N x N的网格子矩形的数量是(N*(N1)/2)^2个。当N很大时比如N10^5子矩形数量是个天文数字枚举显然不现实。我们必须找到面积和的闭合表达式。2.1 问题转化与初步思考一个子矩形由其左上角和右下角决定。考虑任意一个格子(i, j)假设行列均从1开始索引它在多少个矩形里这个格子作为一个矩形的一部分被计入面积当且仅当矩形的范围覆盖了它。更精确地对于行方向包含第i行的矩形其顶行r1必须满足1 r1 i底行r2必须满足i r2 N。所以在行方向上包含第i行的选择有i * (N - i 1)种。同理在列方向上包含第j列的选择有j * (N - j 1)种。那么所有包含格子(i, j)的矩形数量就是[i * (N - i 1)] * [j * (N - j 1)]。每个这样的矩形其面积贡献都会包含这个格子的1个单位面积。因此格子(i, j)对所有矩形面积的总贡献就是i * (N - i 1) * j * (N - j 1)。我们的目标是求所有矩形面积之和也就是所有格子的贡献之和。因为面积是格子数的计数所以总面积极和 S Σ_{i1}^{N} Σ_{j1}^{N} [ i * (N - i 1) * j * (N - j 1) ]2.2 利用对称性分解求和观察上面的二重求和可以发现关于i和j的部分是完全对称且独立的。S [ Σ_{i1}^{N} i * (N - i 1) ] * [ Σ_{j1}^{N} j * (N - j 1) ] [ Σ_{i1}^{N} i * (N - i 1) ] ^ 2现在问题简化为计算单个求和T(N) Σ_{i1}^{N} i * (N - i 1)。2.3 化简求和式等差数列与平方和我们来化简T(N)T(N) Σ_{i1}^{N} [ i * (N 1) - i^2 ] (N 1) * Σ_{i1}^{N} i - Σ_{i1}^{N} i^2这里用到了等差数列求和公式与平方和公式Σ_{i1}^{N} i N * (N 1) / 2Σ_{i1}^{N} i^2 N * (N 1) * (2N 1) / 6代入得T(N) (N 1) * [N*(N1)/2] - [N*(N1)*(2N1)/6] [N*(N1)/2] * (N1) - [N*(N1)*(2N1)/6] N*(N1) * [ (N1)/2 - (2N1)/6 ] N*(N1) * [ (3(N1) - (2N1)) / 6 ] N*(N1) * [ (3N3 -2N -1) / 6 ] N*(N1) * [ (N2) / 6 ] N*(N1)*(N2) / 6因此S [T(N)]^2 [ N*(N1)*(N2) / 6 ] ^ 2。等等这个结果和网上常见的[N*(N1)/2]^2不一样这里需要非常小心。我们重新审视一下T(N)的含义。T(N) Σ i * (N-i1)。我们换一种更直观的方式理解对于长度为L的一维线段所有子线段长度之和是多少一个经典结论是Σ_{a1}^{N} Σ_{ba}^{N} (b-a1) N*(N1)*(N2)/6。这正好对应我们求出的T(N)。但是矩形面积是长乘宽。我们之前计算S时用的是[Σ i*(N-i1)]^2这相当于(所有子段长度之和)^2。然而(长度之和)^2并不等于(长度之积的和)。这里犯了概念错误。正确的思路应该是每个矩形由行上的一个子段和列上的一个子段唯一确定。设行上子段长度为len_row列上子段长度为len_col矩形面积就是len_row * len_col。那么所有矩形的面积和应该是所有可能的len_row与所有可能的len_col两两相乘再求和。即S (所有行子段长度之和) * (所有列子段长度之和)由于行列对称所有行子段长度之和 所有列子段长度之和 T(N) N*(N1)*(N2)/6。 所以S T(N) * T(N) [N*(N1)*(N2)/6]^2。这个推导在逻辑上是自洽的。然而网上广泛流传的另一个更简洁的公式S [N*(N1)/2]^2又是怎么回事我们再来检查一下。考虑所有子矩形的左上角(x1, y1)和右下角(x2, y2)其中1 x1 x2 N,1 y1 y2 N。矩形面积A (x2 - x1 1) * (y2 - y1 1)。总和的另一种计算方式固定左上角(x1, y1)看它能形成多少右下角。对于行右下角的行号x2可以从x1取到N共(N - x1 1)种选择对于列同理有(N - y1 1)种选择。所以以(x1, y1)为左上角的矩形总数为(N - x1 1) * (N - y1 1)。但这不是面积和我们需要对每个矩形乘上其面积(x2-x11)*(y2-y11)。当我们固定左上角(x1, y1)和右下角(x2, y2)时面积是确定的。直接枚举x1, y1, x2, y2求和S Σ_{x11}^{N} Σ_{y11}^{N} Σ_{x2x1}^{N} Σ_{y2y1}^{N} (x2 - x1 1) * (y2 - y1 1)这个四重求和可以分解为两个独立二重求和的乘积S [ Σ_{x11}^{N} Σ_{x2x1}^{N} (x2 - x1 1) ] * [ Σ_{y11}^{N} Σ_{y2y1}^{N} (y2 - y1 1) ]令k x2 - x1 1则当x1固定时k可以从1取到N - x1 1。所以内层求和Σ_{x2x1}^{N} (x2 - x1 1) Σ_{k1}^{N-x11} k (N-x11)*(N-x12)/2。 那么Σ_{x11}^{N} (N-x11)*(N-x12)/2。令i N-x11则当x1从1到N时i从N到1。求和变为Σ_{i1}^{N} i*(i1)/2 (1/2) Σ_{i1}^{N} (i^2 i) (1/2) [ Σ i^2 Σ i ] (1/2) [ N(N1)(2N1)/6 N(N1)/2 ] (1/2) * N(N1)/2 * [ (2N1)/3 1 ] (1/2) * N(N1)/2 * [ (2N4)/3 ] N(N1)(N2)/6。看我们又得到了T(N) N(N1)(N2)/6。所以S T(N)^2。那么[N(N1)/2]^2是什么它其实是所有子矩形的数量而不是面积和因为固定左上角(x1, y1)后可能的右下角(x2, y2)数量是(N-x11)*(N-y11)对所有左上角求和Σ_{x11}^{N} Σ_{y11}^{N} (N-x11)*(N-y11) [Σ_{x11}^{N}(N-x11)] * [Σ_{y11}^{N}(N-y11)] [Σ_{i1}^{N} i] * [Σ_{i1}^{N} i] [N(N1)/2]^2。关键结论务必区分清楚子矩形数量[N(N1)/2]^2子矩形面积和[N(N1)(N2)/6]^2很多早期的题解或讨论可能混淆了这两个概念或者题目本身在不同版本、不同理解下有所变化。根据2020牛客国庆集训的题目描述计算面积和以及其需要处理大数的事实正确的公式应该是面积和即S [N(N1)(N2)/6]^2。这个公式涉及N^3量级的计算在N很大时即使中间步骤也极易溢出64位整数这正好契合了题目要求使用“大数模板”的初衷。3. 大数运算的挑战与核心思路公式S [N(N1)(N2)/6]^2确定了现在的问题是计算。N很大直接计算N*(N1)*(N2)会超出long long的范围通常是2^63-1约9.22e18。即使先除以6再平方中间过程N*(N1)*(N2)也可能溢出。所以我们必须使用大数运算。大数运算本质上就是用程序模拟我们手算多位整数加减乘除的过程。最常见的表示方法是用数组或字符串存储每一位数字。对于本题我们需要的操作是大整数乘法、大整数除以小整数、大整数比较可能不需要。由于公式固定我们可以优化计算顺序尽量减少大数运算的规模和次数。3.1 计算顺序的优化与溢出规避一个最直接的陷阱是先计算分子 N*(N1)*(N2)再除以6最后平方。但分子可能已经溢出。即使使用大数N*(N1)*(N2)这个乘法涉及三个大数相乘计算量较大。我们可以利用整除性来简化。注意到N,N1,N2是三个连续的整数。其中必有一个是2的倍数也必有一个是3的倍数。因此N*(N1)*(N2)一定能被6整除。我们可以在乘法过程中就进行约分从而让参与大数运算的数尽可能小。优化计算路径计算A N * (N1) / 2。因为N和N1必有一个偶数所以A是整数。此时A可能还是大数。计算B A * (N2) / 3。因为N,N1,N2中必有一个是3的倍数。在第一步除以2后剩下的部分与(N2)相乘整体仍能被3整除。所以B是整数且B N*(N1)*(N2)/6。计算最终结果S B * B。这个路径的优点是每一步的除法都是除以一个小整数2或3可以高效地实现为大数除以小整数的运算避免了复杂的大数除法。同时每一步乘法都是两个大数相乘复杂度可控。3.2 大数模板的设计要点一个实用的大数模板这里指非Python语言如C、Java等需要手动实现的情况通常包含以下核心部分存储使用vectorint或数组按十进制从低位到高位存储每一位数字。例如数字12345存为[5,4,3,2,1]方便进位处理。乘法大数 × 小整数模拟手算从低位到高位每一位乘以小整数加上进位然后处理进位。复杂度 O(n)。大数 × 大数模拟手算乘法竖式使用双重循环结果位数最多为两数位数之和。复杂度 O(n*m)。对于本题N的位数不会太多题目通常保证结果在可表示范围内所以这个复杂度可以接受。除法大数 ÷ 小整数模拟手算除法从高位到低位当前余数乘以10加上下一位除以除数得到商的一位更新余数。复杂度 O(n)。这是本题最关键的操作。输入/输出从字符串读入转换到内部数组将内部数组格式化为字符串输出。在C中为了避免重复造轮子许多竞赛选手会准备一个现成的“大数模板”类封装这些操作。本题中我们只需要实现乘法大数×大数大数×int、除法大数÷int以及构造函数、输出即可。实操心得除法运算的细节大数除以小整数时最容易出错的地方是商的前导零。例如计算100 / 8手工算是12余4。算法过程从最高位开始1/80余1余数110下一位01010/81余2余数210下一位02020/82余4。得到的商位序列是[0,1,2]。存储时我们通常从低位存到高位所以得到[2,1,0]。输出前需要反转并去除高位的零得到12。在代码实现中需要在计算完成后将结果数组反转并pop_back掉末尾即实际的高位的零直到剩下最后一位或者遇到非零位。4. 代码实现与逐行解析下面我们用一个C的大数类来实现上述优化计算路径。为了让逻辑更清晰我们将大数类命名为BigInt并实现必要的操作。#include iostream #include string #include vector #include algorithm using namespace std; class BigInt { private: vectorint digits; // 低位在前高位在后 bool positive; // 符号本题均为正数可简化 public: // 构造函数从字符串构造 BigInt(const string s 0) { positive true; // 本题全为正 digits.clear(); for (int i s.size() - 1; i 0; --i) { digits.push_back(s[i] - 0); } trim(); // 去除前导零 } // 构造函数从long long构造方便测试 BigInt(long long num) { positive num 0; if (num 0) num -num; digits.clear(); if (num 0) digits.push_back(0); while (num 0) { digits.push_back(num % 10); num / 10; } } // 去除前导零 void trim() { while (digits.size() 1 digits.back() 0) { digits.pop_back(); } if (digits.empty()) digits.push_back(0); } // 输出 string toString() const { string s; for (int i digits.size() - 1; i 0; --i) { s char(digits[i] 0); } return s; } // 大数 × 小整数 BigInt operator*(int b) const { if (b 0) return BigInt(0); BigInt res; res.digits.assign(digits.size() 10, 0); // 预分配足够空间 long long carry 0; for (size_t i 0; i digits.size(); i) { carry (long long)digits[i] * b; res.digits[i] carry % 10; carry / 10; } while (carry 0) { res.digits.push_back(carry % 10); carry / 10; } res.trim(); return res; } // 大数 × 大数 BigInt operator*(const BigInt b) const { BigInt res; res.digits.assign(digits.size() b.digits.size(), 0); for (size_t i 0; i digits.size(); i) { int carry 0; for (size_t j 0; j b.digits.size() || carry; j) { long long cur res.digits[i j] (long long)digits[i] * (j b.digits.size() ? b.digits[j] : 0) carry; res.digits[i j] cur % 10; carry cur / 10; } } res.trim(); return res; } // 大数 ÷ 小整数返回商 BigInt operator/(int b) const { BigInt res; res.digits.resize(digits.size()); long long remainder 0; for (int i digits.size() - 1; i 0; --i) { remainder remainder * 10 digits[i]; res.digits[i] remainder / b; remainder % b; } // 现在res.digits是高位在前需要反转并去除前导零 reverse(res.digits.begin(), res.digits.end()); res.trim(); return res; } }; int main() { // 假设输入为字符串形式的N string N_str; // 这里模拟输入实际应从cin读取 // cin N_str; N_str 1000000000000000000; // 示例10^18 BigInt N(N_str); // 优化计算路径S [ (N*(N1)/2) * (N2) / 3 ] ^ 2 BigInt N_plus_one N * 1 BigInt(1); // 计算 N1这里用了个技巧先乘1转BigInt再加 // 更稳妥的方式是实现 BigInt 的加法这里为简化我们用字符串构造 // 实际上对于N1, N2我们可以直接通过字符串或数值计算得到其字符串表示再构造BigInt。 // 由于N可能非常大实现一个完整的加法是更通用的做法。为了代码紧凑我们假设有加法后续补充。 // 我们先实现一个简单的字符串加法函数用于计算大数字符串加1、加2 auto addOne [](const string s) - string { string res s; int carry 1; for (int i res.size() - 1; i 0 carry; --i) { int digit (res[i] - 0) carry; res[i] digit % 10 0; carry digit / 10; } if (carry) res 1 res; return res; }; string N_str_plus1 addOne(N_str); string N_str_plus2 addOne(N_str_plus1); BigInt N1(N_str); // N BigInt N2(N_str_plus1); // N1 BigInt N3(N_str_plus2); // N2 // 计算 A N * (N1) / 2 BigInt A (N1 * N2) / 2; // 计算 B A * (N2) / 3 BigInt B (A * N3) / 3; // 计算最终结果 S B * B BigInt S B * B; cout S.toString() endl; return 0; }上面的代码框架展示了核心逻辑但缺失了BigInt的加法运算符重载。为了完整性我们补充一个加法实现并完善整个程序。// 在BigInt类内部添加加法成员函数 BigInt operator(const BigInt b) const { BigInt res; int carry 0; size_t maxLen max(digits.size(), b.digits.size()); res.digits.assign(maxLen, 0); for (size_t i 0; i maxLen || carry; i) { int sum carry; if (i digits.size()) sum digits[i]; if (i b.digits.size()) sum b.digits[i]; if (i res.digits.size()) { res.digits.push_back(sum % 10); } else { res.digits[i] sum % 10; } carry sum / 10; } res.trim(); return res; }有了加法主函数可以更优雅地计算N1和N2int main() { string N_str; cin N_str; BigInt N(N_str); BigInt N_plus_one N BigInt(1); BigInt N_plus_two N_plus_one BigInt(1); BigInt A (N * N_plus_one) / 2; BigInt B (A * N_plus_two) / 3; BigInt S B * B; cout S.toString() endl; return 0; }代码实现的几个关键点存储顺序低位在前极大简化了进位处理。乘法时digits[i]和digits[j]的乘积贡献到结果的[ij]位这与我们手算时的对齐方式一致。除法优化除以小整数时我们从高位向低位处理模拟竖式除法。得到的商位顺序在高位到低位所以需要一次reverse。这是容易忽略的步骤。空间预分配在乘法和加法中预先分配足够大的空间如digits.size() b.digits.size()比动态push_back更高效避免了多次扩容。输入边界题目没有明确给出N的上限但我们的BigInt可以处理任意长度的十进制字符串只要内存足够。关于N1和N2直接使用大数加法是最稳妥的。如果N很大将其转换为long long再加一会溢出。所以我们必须用大数加法或字符串加法来计算。5. 测试、验证与常见问题排查实现完大数模板和计算逻辑后必须进行充分的测试。我们可以从小数据开始逐步增大并与暴力计算小数据下或已知结果进行对比。5.1 测试用例设计基础验证N1只有一个1x1的矩形面积和为1。公式[1*2*3/6]^2 1^2 1。N2手动枚举所有子矩形1x1的有4个1x2的有2个横2x1的有2个竖2x2的有1个。面积和 41 22 22 14 444416。公式[2*3*4/6]^2 (4)^2 16。N3公式结果应为[3*4*5/6]^2 (10)^2 100。中等数据验证N10可以用long long写一个O(N^4)或O(N^2)的暴力程序验证结果是否一致。N100同样用暴力程序O(N^2)公式验证。大数边界测试N10^18这是long long的边界附近。计算N*(N1)*(N2)会远超long long范围必须用我们的大数程序。可以计算并输出结果检查其位数是否合理例如用Python等原生支持大数的语言写一个脚本进行对比。5.2 常见问题与调试技巧结果错误为0检查除法运算。特别是operator/(int)中reverse和trim的逻辑是否正确。如果商的高位是0trim可能会把整个数变成0。确保trim函数在数字为0时保留一个0。检查乘法运算的进位处理。在双重循环中进位carry可能不止一位要用long long存储中间结果防止溢出。结果位数不对或数字混乱确认存储顺序。输入字符串123应该转化为digits [3,2,1]。输出时要从高位digits的末尾开始。在乘法中res.digits[ij]的累加可能超过10需要正确处理进位。内层循环的终止条件j b.digits.size() || carry很重要确保所有进位都被处理。性能问题对于极大的N比如有上千位O(n^2)的大数乘法可能会变慢。可以考虑更高效的算法如Karatsuba算法但对于竞赛题N的位数通常不会大到需要优化乘法的地步。主要的性能瓶颈在于三次大数乘法N*(N1),A*(N2),B*B和两次除法。除法是O(n)的很快。乘法是主要开销。内存使用vectorint存储每位数字一个int存一个0-9的数字有些浪费但代码清晰。可以改用vectorshort或压位存储如一个int存9位十进制数但会增加代码复杂度。踩坑实录除法中的前导零陷阱我在第一次实现大数除以小整数时没有进行reverse导致商的高位存储在数组末尾。例如计算123 / 2算法得到商位序列[6,1,0]从高位到低位处理时先得到1再得到6存储时按顺序push_back成了[1,6]这里需要仔细推演。实际上模拟竖式时我们是从被除数的高位开始依次得到商的每一位。如果我们把商也按低位到高位存储那么得到的第一位商最高位会被放在数组的最后这不符合我们的存储约定。因此必须在计算结束后将商的数组反转。这个错误会导致结果完全错误特别是当商有前导零时trim会错误地截断。务必在单元测试中专门测试除法尤其是整除和有余数的情况。6. 公式的再探讨与不同理解的辨析在解题和搜索资料时你可能会遇到关于这个公式的两种形式[N(N1)/2]^2和[N(N1)(N2)/6]^2。这很可能源于对问题描述的不同理解或者是历史题目版本的差异。[N(N1)/2]^2计算的是所有子矩形的数量。推导方式是在NxN网格中选择左上角有N^2种可能不对。更准确地说选择左上角(x1,y1)和右下角(x2,y2)且满足x1x2, y1y2。左上角的行x1有N种选择列y1有N种选择。对于给定的左上角右下角的行x2有N-x11种选择列y2有N-y11种选择。所以总数为Σ_{x11}^N Σ_{y11}^N (N-x11)(N-y11) [Σ_{i1}^N i]^2 [N(N1)/2]^2。[N(N1)(N2)/6]^2计算的是所有子矩形的面积之和。正如我们第二节的推导面积和等于所有行子段长度之和×所有列子段长度之和而行子段长度之和T(N) Σ_{len1}^N len * (N-len1) N(N1)(N2)/6。如何快速验证哪个是面积和用N2测试即可。矩形数量[2*3/2]^2 3^2 9。确实1x1有4个1x2有2个2x1有2个2x2有1个共9个。面积和我们手动算得是16。用[2*3*4/6]^2 4^2 16符合。用[2*3/2]^2 3^2 9不符合。因此对于“SUM OF SUB RECTANGLE AREAS”这个题目显然应该使用面积和的公式。这也解释了为什么题目强调“大数模板”因为N(N1)(N2)/6的增长速度是O(N^3)比O(N^2)的矩形数量公式更容易溢出。经验之谈竞赛中的公式推导在算法竞赛中遇到这种计算类题目一定不要轻信网上现成的公式尤其是早期来源不明的题解。最好的方法是自己从头推导一遍并用小数据暴力程序验证。如果公式推导复杂可以尝试写出N1,2,3,4时的暴力结果然后观察规律猜测通项公式再用数学归纳法证明。本题中暴力枚举N1,2,3的面积和分别为1, 16, 100。观察这些数字11^2, 164^2, 10010^2。而1,4,10恰好是N(N1)(N2)/6在N1,2,3时的值。这个规律非常明显可以快速锁定正确公式。7. 超越本题大数运算的优化与工程实践虽然我们实现了一个基础的大数模板并能AC此题但在更严苛的场景下例如N有数万位O(n^2)的朴素乘法会成为瓶颈。这里简要提一下进阶的优化方向Karatsuba算法将两个大数X和Y分别拆分成高位和低位X A*10^m B,Y C*10^m D。那么X*Y AC*10^(2m) ((AB)(CD)-AC-BD)*10^m BD。通过递归计算三次较小规模的乘法可以将时间复杂度从O(n^2)降低到约O(n^1.585)。实现起来比朴素乘法复杂但能有效处理更大的数。压位存储我们目前用一个int存一个十进制位浪费了大量空间和计算资源。更高效的做法是用一个int存储多位十进制数比如BASE10000万进制这样每个数组元素存储0-9999的数。乘法、加法时需要处理进位和对BASE取模。这能显著减少数组长度和循环次数。使用现成库在非竞赛的工程实践中除非有极特殊的定制需求否则强烈建议使用成熟的任意精度数学库如 C 的GMPGNU Multiple Precision Arithmetic Library、Java 的BigInteger、Python 的int原生支持。这些库经过高度优化稳定且高效。对于算法竞赛准备一个像我们上面实现的、功能完备的朴素大数模板通常足以应对90%的需要大数的题目。重点是保证正确性和清晰度以便在紧张的比赛环境中快速调试。回过头看这道题它巧妙地将数学推导、规律发现和大数运算结合在一起。推导公式考察了参赛者的数学功底和观察力而大数实现则检验了基本的算法实现能力和对细节的把握。即使你知道了公式如果大数模板写得漏洞百出或者计算顺序没优化导致中间结果溢出依然无法得分。这种题目非常锻炼综合能力。