有一类MATLAB题看起来毫无技术含量做完却让人记很久。“pizza”就是其中一道。我最早是在MathWorks的在线练习Cody上看到它的题目说明只有一句话大致意思是有一张圆形披萨你一刀一刀切下去每刀都是直的、必须横穿整张披萨问切n刀最多能切成多少块。很多人的第一反应是打开几何画板开始画直线甚至有人真的先画了一堆直线再逐个区域去数但把纸和笔拿出来推两分钟你会发现答案就是一个非常干净的公式。这篇文章就来讲清楚这道题从建模到MATLAB实现再到各种变体的完整过程适合刚学MATLAB的新手也适合想复习递推与算法分析的老手。1. 从一道披萨选择题说起题目本意与常见误解1.1 我是在哪遇到“Pizza”这道题的当时是随手刷题题名就叫 Pizza函数签名要求输入一个非负整数 n返回一个整数。评论区里很有意思前排清一色都是“这不就是切披萨吗”后半页则是各种觉得自己被小学数学骗了的回复。说实话这种题放在数学练习册里就是一道平平无奇的递推题但放在 MATLAB 里它会逼着你考虑很多平时写业务代码根本不会想的事向量化、整数类型、内存占用、递归深度。这也是我后来特别喜欢拿它做教学案例的原因。这道题不依赖任何工具箱不涉及图像处理、不涉及优化求解器甚至连 Simulink 都用不上。装好 MATLAB新建一个函数文件或者直接在命令行里写一行匿名函数就能跑通。对于刚下载完 MATLAB、还在纠结“我装这个软件到底能干嘛”的新手来说这几乎是完美的第一个练手函数。1.2 题目里那些容易忽略的“默认规则”这类题目的难点不在代码而在“把现实问题翻译成数学模型”时容易漏条件。我见过好几个同学卡在错误答案上都是因为对题目默认规则的理解有偏差刀必须是直线。不是锯齿刀不是曲线切法就是一条直线切下去。每一刀都要贯穿整张披萨。直直地从这一侧切到那一侧不能只切一个角不能切一半停在中间。不要求每块一样大。题目问的是“最多多少块”不是“能不能平均分”。切完不能移动、叠放。切完就放在原地下一刀还是在当前这块圆形披萨上切。这些约束缺一个数学模型都会变。比如“允许叠起来切”的话答案会变成完全不同的东西我后面会专门讲。先记住原题的标准场景就是“一张圆形披萨、直线贯穿、不移动、追求最大块数”。1.3 一个先入为主的错误答案遇到这种题人脑很容易蹦出几个看起来合理的答案。有人猜是 2^n理由是切一刀变两块、切两刀变四块、切三刀好像也能变七块但实在懒得画也有人猜是 2n理由是每切一刀多两块还有人会纠结“如果不经过圆心呢如果不平行呢”。真相是前几个小 n 勉强能骗人n1 是 2 块n2 是 4 块n3 是 7 块而不是 8 块n4 是 11 块而不是 16 块。所以 2^n 在 n3 就已经错了而 2n 只在每刀都过圆心的时候才成立。真正决定答案的是“新增一刀到底最多穿过了多少个旧区域”。2. 一张纸推完递推式为什么第 n 刀恰好增加 n 块2.1 从手工画刀开始n0 到 n4先别急着写代码拿笔画一画。设 P(n) 表示切 n 刀后最多得到的块数刀数最多块数直观解释01整张披萨12一刀切成两半24第二刀与第一刀相交37第三刀穿过前两刀形成的两个交点把路径上经过的 3 个区域各一分为二411第四刀最多与前三条各相交一次产生 3 个交点把新刀切成 4 段每段都切出一个新区块注意 n3 到 n4 的增幅是 4n2 到 n3 的增幅是 3n1 到 n2 的增幅是 2。看到规律了吗每多切一刀增幅都在递增。2.2 把“新刀增加多少块”变成计数问题核心观察其实只有一句话第 k 刀最多被前 k-1 刀交出 k-1 个交点这 k-1 个交点把新刀在披萨内部的那段线段分成了 k 小段每一小段都会把一个旧区域一分为二所以新增 k 块。这就是“递推增量分析”的典型套路。这里的每小段必须满足一个条件它完整地落在某个旧区域内部。如果新刀和之前的刀平行它可能少产生交点如果三刀交于同一点交点数会“合并”新刀被切出的段数也会变少如果一刀下去只穿过了披萨的一个角那它可能连一个贯穿的完整线段都没有增量自然更小。所以在追求“最多”时我们要求刀与刀之间都保持合理的一般位置。作为对比如果每刀都过圆心那么第 k 刀和前 k-1 刀都在圆心处相交k-1 个交点重合成了一个交点新刀在圆内被分成 2 段所以每一刀都只增加 2 块最终是 2n。这也是为什么很多人一开始会猜 2n——他们潜意识里把“不共点”这个条件忽略了。2.3 闭式解、等差数列与最终公式有了 P(k) P(k-1) k 之后剩下的就是求和了P(0) 1P(n) 1 1 2 3 ... nP(n) 1 n(n1)/2也写作 (n^2 n 2)/2。这个公式漂亮得不像话n10 的时候答案是 56 块。我每次给学生讲到这里都会补一句数学上你可以做到但现实里你大概率切不出 56 块因为人眼手配合根本保证不了“每条切线都避开所有三重交点”。这就是理论与实操的浪漫差距。2.4 一般位置条件的等价说法把“一般位置”翻译成更容易检查的条件其实就三条没有两条切线平行。没有三条切线交于同一点。每条切线都必须真正横穿圆形披萨内部。这里第三条很容易被忽略。举个例子n2 时如果两条直线都只擦过圆形边缘交点落在披萨外面那它们对披萨的分割效果甚至不如两条平行线。所以讨论最大块数时要保证所有交点都落在圆盘内部或者至少保证在圆内的切割线段没有被“浪费”。验证一条直线是否贯穿圆盘最简单的办法是计算圆心到直线的距离。设直线为 ax by c归一化后圆心到直线距离为 |c| / sqrt(a^2b^2)只有小于半径时直线才真正割到披萨。后面可视化验证那段我就是用这个条件来筛选割线的。3. 用 MATLAB 实现时别急着写循环3.1 四种主流写法循环、向量化、匿名函数、递归既然公式已经推出来了实现方式有很多种。下面是我见过的几种典型写法。最直白的循环版function p pizzaLoop(n) p 1; for k 1:n p p k; end end这个版本完全对应递推式 P(k) P(k-1) k可读性最好适合第一次接触这道题的人。适合批量计算多个 n 的向量化版本function p pizzaVec(n) p 1 sum(1:n); end利用 sum(1:n) 一次性完成等差数列求和。看着简洁但注意它会在内存里构造一个长度为 n 的向量。递归版展示数学递归定义function p pizzaRec(n) if n 0 p 1; else p pizzaRec(n - 1) n; end end这个版本最贴近递推式但性能最差n 到几千就会明显变慢而且有递归调用开销MATLAB 没有做尾递归优化不推荐在正式提交时用。最推荐的搞定版一行匿名函数pizza (n) 1 n .* (n 1) / 2;注意用点乘 .*这样输入向量时也能一次返回多个结果。想验证的话可以直接跑n 0:10; pizza(n)输出 1, 2, 4, 7, 11, 16, 22, 29, 37, 46, 56和递推结果完全一致。3.2 为可读性和性能我推荐哪种如果是在 Cody 这类平台上提交函数我会直接写公式版function p pizza(n) p 1 n * (n 1) / 2; end原因是它 O(1) 时间复杂度不依赖任何循环也不会在 n 很大时因为内存分配而崩掉。更重要的是公式版把“这道题到底在考什么”表达得特别清楚递推关系、求和、闭式解。如果你读别人的代码一行公式就能看懂如果写循环你还要顺着状态累积去脑内模拟。但如果是学习用途我反而建议先把循环版写对再用公式版提速。很多初学者跳过了“模拟递推”这一步直接背公式结果题目稍一变形就不知道从哪里入手。先看懂结构再追求简洁这是我一直坚持的顺序。3.3 一个会被 sum(1:n) 坑到的内存问题向量化看起来又短又“MATLAB”但这里有一个隐藏陷阱。sum(1:n) 是先构造 1:n 这个长度为 n 的 double 数组再求和。如果 n 1e7数组占 80MB勉强能跑n 1e8数组占 800MB很多机器直接内存不足或者开始疯狂交换内存。公式版永远没有这个烦恼。这提醒我们一件事向量化不一定等于高性能。MATLAB 的向量化优势在于“减少循环解释开销”但代价可能是额外的临时数组分配。对这种规模很简单的问题老老实实一行公式是最稳的做法。如果你碰到其它问题看到 sum 后面跟了一个很大的冒号表达式先停下来想一想有没有等效的闭式公式。4. 大数陷阱double、int64 与符号运算4.1 double 真的够用吗大多数人的第一版代码是这样写的p 1 n * (n 1) / 2;n 是 double 类型结果自然也是 double。对常规测试用例来说完全没问题。但如果把 n 拉到 1e8P(n) 大约等于 5e15n 拉到 1e9结果约 5e17。而 double 能够精确表示不超过 2^53 的整数2^53 约等于 9.007e15。所以 n 1e8 时已经在悬崖边n 1e9 时必然产生舍入误差。也就是说如果题目只保证 n 1e4double 完全够用但如果你在做数值实验或者想得到超大 n 的精确整数结果就必须换思路。4.2 用 int64/uint64 前先算一算溢出边界很多同学会想“那我用 int64 不就行了”。int64 的最大值是 9.223e18看起来很大但直接算 n^2 依然危险。int64 溢出时 MATLAB 会静默回绕不报错你得到的是一个看起来莫名其妙的负数。这类 bug 非常难排查因为大多数时候它不会触发一旦触发就是灾难式错误。同样地uint64 最大值约为 1.844e19能承受更大的平方但当 n 超过 4e9 时 n^2 也会爆掉。即使写成 n/2*(n1)虽然推迟了溢出也只是把边界从 n≈2^31 提升到 n≈4e9不可能一劳永逸。所以在整型方案里必须先算清楚你的 n 量级再决定类型。一般算法题的 n 都很小int64 足够但想写出“任何 n 都正确”的函数就得用下面这招。4.3 需要精确大整数时直接用 symMATLAB 的符号运算可以处理任意大整数只要把 n 转成 symp 1 sym(n) * (sym(n) 1) / 2;例如 n 12345678901234567890double 早就面目全非了sym 仍然能精确算出结果。代价是速度慢很多但在这种求公式值的场景下符号运算几乎瞬间完成。我个人处理这种边界问题时的习惯是先用 double 跑通逻辑再根据输入规模决定要不要加健壮性处理。不要一上来就上 sym因为符号运算对数组不友好写起来也啰嗦。等测试用例真的出现了大数再针对它加一层判断。5. 代码对了还不算完用图形动画和暴力计数验证5.1 画一张动态切披萨图我拿到任何算法题都喜欢在纸上画完推完公式后再用 MATLAB 画一张可视化图确认自己没有推导到另一个平行宇宙去。“pizza”这道题特别适合画单位圆表示披萨随机生成相互穿插的割线然后逐刀显示。下面这个函数我经常用来演示。它生成 n 条直线每条直线用“法向角度 theta 到圆心距离 d”表示即 xcos(theta) ysin(theta) d。归一化后d 的绝对值小于 1 表示直线穿过了圆。function showPizzaCuts(n) th linspace(0, 2*pi, 200); figure(Color, w); plot(cos(th), sin(th), k-, LineWidth, 1.5); axis equal; axis([-1.2 1.2 -1.2 1.2]); hold on; grid on; thetaList zeros(1, n); dList zeros(1, n); for k 1:n % 随机生成一条穿过圆盘的直线 while true theta rand * pi; d (rand * 2 - 1) * 0.95; % 半弦长 s sqrt(1 - d^2); % 两个端点 hx d * cos(theta); hy d * sin(theta); ux -sin(theta); uy cos(theta); x1 hx - s * ux; y1 hy - s * uy; x2 hx s * ux; y2 hy s * uy; % 简单的“一般位置”检测新交点不要离已有交点太近 ok true; for j 1:k-1 % 两直线交点 A [cos(theta), sin(theta); cos(thetaList(j)), sin(thetaList(j))]; b [d; dList(j)]; if abs(det(A)) 1e-12 ok false; break; end inter A \ b; if inter(1)^2 inter(2)^2 1 ... || min(abs(inter(1)-x1), abs(inter(2)-y1)) 1e-6 ok false; break; end end if ok break; end end thetaList(k) theta; dList(k) d; % 画这条割线 s sqrt(1 - dList(k)^2); hx dList(k) * cos(thetaList(k)); hy dList(k) * sin(thetaList(k)); ux -sin(thetaList(k)); uy cos(thetaList(k)); plot([hx - s*ux, hx s*ux], [hy - s*uy, hy s*uy], LineWidth, 1.2); % 标题显示理论块数 p 1 k * (k 1) / 2; title(sprintf(k %d, P(k) %d, k, p)); pause(0.8); end end运行 showPizzaCuts(6)你会看到前几刀还算规整后几刀开始交叉得很乱但标题上的 P(k) 完美对应 2、4、7、11、16、22。这个过程能非常直观地验证公式。5.2 为什么用“角度距离”而不是斜截式写直线时很多人习惯 y kx b但斜截式对竖直直线是灾难而且两条直线是否平行在数值上也不好判断。用 xcos(theta) ysin(theta) d 这个法线式表达天然规避了竖直直线问题裁剪到圆内也方便垂足是 (dcos(theta), dsin(theta))方向向量是 (-sin(theta), cos(theta))半弦长由 s sqrt(1 - d^2) 给出。这也是我在实际画图时比较喜欢的一种参数化方式。虽然多写了几行但代码的稳定性上了一个台阶。你以后如果在 MATLAB 里做任意直线与圆求交可以直接抄这个结构。5.3 按区域颜色编号暴力验证 n 较小时的答案光看图形有点“眼见为虚”尤其当 n 到 7、8 之后圈圈叉叉已经很难数了。这时可以用一个简单粗暴的办法像素法计数。把圆盘离散成网格每条割线附近的像素判定为边界其余像素判定为内部然后用连通域标记统计白色区域的个数。function cnt countPizzaRegions(thetaList, dList, gridN) % 生成网格 g linspace(-1.2, 1.2, gridN); [X, Y] meshgrid(g, g); % 距离所有割线的最短距离 Dmin inf(size(X)); for i 1:numel(thetaList) Di abs(X * cos(thetaList(i)) Y * sin(thetaList(i)) - dList(i)); Dmin min(Dmin, Di); end % 圆内且不是边界的像素视为披萨内部 inPizza (X.^2 Y.^2) 1; boundary Dmin * gridN / 2.4 1.2; % 大约1.2个像素宽度 mask inPizza ~boundary; % 统计连通域数量 L bwlabel(mask); cnt max(L(:)); end注意这段代码需要 Image Processing Toolbox 里的 bwlabel。如果没有这个工具箱可以自己写一个 BFS 连通域标记逻辑也不难只是代码更长一些。像素法虽然不如平面图法严谨但只要网格足够密、边界阈值足够小它对 n6 的验证是可靠的。我实测 n6 随机割线网格取 800x800大多数情况都能数出 22 块。这算是一个“工程验证”思路数学推导归数学推导代码实现归代码实现用一条完全独立的路径去交叉验证能躲掉绝大多数粗心错误。6. 这道披萨题的亲戚们切蛋糕、曲线切割与分形脑洞6.1 三维版本n 刀最多能切出多少块蛋糕把“pizza”从二维升级到三维就变成了经典的切蛋糕问题一个长方体蛋糕用 n 个平面切最多能切成多少块这里一刀是一个平面贯穿整个蛋糕。结论是S(n) (n^3 5n 6) / 6前几项是 1、2、4、8、15、26、42……注意 n3 时是 8n4 时是 15不是 16。这个“少一块”的结果让很多人意外但用它来理解递推关系非常舒服第 n 个平面与前 n-1 个平面最多相交出若干条直线。这些直线在该平面上划分了一些区域。每个区域会让这个新平面穿过的一个旧空间块分裂成两块。所以 S(n) S(n-1) L(n-1)其中 L(n-1) 是二维版本 n-1 条直线最多分割平面的数量。L(n-1) 1 (n-1)n/2累加后就得到上面那个公式。这个推导过程恰好是把披萨题又用了一遍只是每个维度往上升了一级。学递推的时候能把二维三维放一起看理解会深很多。6.2 如果必须过圆心、或者允许叠起来切这道题有三个常见变形答案完全不同我列个表方便对比切割规则最大块数例子 n4直线贯穿一般位置原题1 n(n1)/211每刀都必须过圆心2n8每次切完后叠起来再一刀全切2^n16“每次切完叠起来”这个变体特别有意思。它对应的是信息论里的二分逻辑每一刀都能把现有的每一块都切到所以块数翻倍。很多公司的面试题“一张纸对折 n 次后有多少层”或者“一根绳子折 n 次剪一刀得到多少段”本质都是这个模型。6.3 从平面分割到超平面排列把一维、二维、三维的结论放到一起你会发现一条清晰的规律一维n 个点最多把直线分成 1 n 段。二维n 条直线最多把平面分成 1 n C(n,2) 块。三维n 个平面最多把空间分成 1 n C(n,2) C(n,3) 块。也就是说d 维空间被 n 个超平面最多分割的区域数是前几项二项式系数之和R(n, d) sum_{i0}^{d} C(n, i)披萨题不过就是 d2 的特例。如果你感兴趣把这个公式写成 MATLAB 循环其实也就几行但那种“从一个披萨挖到高维几何”的延伸感才是这类小题目真正值钱的地方。最后说点个人体会。这道题给我最大的启发不是公式本身而是“先想清楚数学结构再动手写代码”的习惯。很多时候我们打开 MATLAB 第一件事就是画图、模拟、写循环但把递推关系摆在纸上之后代码量能压缩到一行。它也是我讲给学生的一道经典入门题因为一行公式背后有足够多的坑可以聊一般位置、整数类型、递归栈、内存占用、图形验证全都占齐了。如果你也碰到这种“看起来简单、推起来有趣”的题目建议按这个套路做一次完整复盘收获会比刷十道重复题都大。