二分查找与NTC查表:边界、插值与C工程避坑
发布时间:2026/10/1 9:24:11来源:尧图网络
凌晨两点温控板上那支 NTC 又把 62 度的水温读成了 76 度而同一块板子在常温下准得让人放心。那一刻我盯着屏幕上的二分查找代码明明教科书上三行就能写完、二分法查找的复杂度还漂亮得不行问题却出在最不该出问题的地方。后来把链路一层层剥开才发现二分法本身没算错任何一个下标真正捅刀子的是折半定位之后拿什么当结果这件事——用最近邻查表去逼近一条指数曲线误差被放大到十几度完全不奇怪。这篇就把二分法从算法和工程两个角度一起捋一遍前半段讲清楚折半的边界怎么写、为什么那么写后半段专治 NTC 查表里那些二分法不准的疑难杂症顺带把数学课本里那个求方程近似根的二分法也一并收了。代码给 C因为真正在板子上跑的、真正会踩坑的基本都是 C。1. 二分查找的两条铁律与一次手算推演1.1 折半的本质用有序性换对数时间二分法最原始的形态就是猜数字我心里想一个 1 到 100 之间的数你猜我告诉你大了还是小了。笨办法是从 1 开始一个一个试运气差要试 100 次聪明办法是先猜 50一刀砍掉一半剩下 25 次的机会都没有最多 7 次就能锁定。这就是二分法的全部魅力——它不是更快一点的查找而是把一个 O(n) 的问题硬生生压成了 O(log n)。这个压缩能成立靠的是一个非常强的不变式每次比较之后我们能确定答案一定不在刚刚排除掉的那一半里。换句话说有序性把一个全局未知的问题变成了局部可判定的问题。我只要知道中点的值和目标的大小关系就能给整个左半边或者右半边判死刑不需要看它们具体是什么。这个思维方式比代码本身重要得多后面你会看到 NTC 查表之所以能用二分也是因为它满足同样的单调结构。具体到次数n 个元素最多需要 ⌈log2(n1)⌉ 次比较。1000 个元素 10 次100 万个元素 20 次40 亿个元素 32 次。我做过一个粗略对比在一个 4096 项的整型表里做线性查找平均值大约 2048 次比较每次几个周期算下来几万个时钟周期换成二分只要 12 次比较几百个周期就结束了。对那种每秒要采几千次温度、每次都要查表的嵌入式场景这个差距是实打实的。1.2 单调性、可随机访问、边界可比缺一不可很多人写二分只在数组上练过于是下意识以为二分是数组专属技巧。其实二分的前提条件只有三条数组只是恰好最方便满足它们而已理解这三条能帮你在遇到二分结果不对时快速定位问题出在哪一层。前提条件含义不满足时会怎样典型反例序列单调整体升序或整体降序且比较方向与代码假设一致排除掉的一侧可能真的藏着答案结果直接错NTC 阻值表是降序套用升序模板O(1) 随机访问能在常数时间内取到任意下标的元素每次跳一半都要走一遍链表复杂度退化成 O(n log n)单链表、只读流式数据比较关系良序任意两元素可比且无歧义收缩方向不确定可能死循环或漏掉元素含 NaN 的浮点数组、结构体只比了部分字段还有一个容易被忽略的隐含条件每一步都必须让搜索区间严格变小。也就是说算出来的 mid 必须真的落在 (left, right) 内部否则 left 和 right 不动while 条件永远为真程序就挂在那儿了。这个坑在找下界这类变体里特别容易踩因为那种写法里有一侧是保留 mid而不是跳过 mid的一旦边界公式写错就死循环。提示判断一段数据能不能二分最快的自检方式是问自己如果我排除掉左半边我能保证答案不在这里面吗。答不上来就别用二分。1.3 十个元素的手算全过程光看代码很难建立手感我习惯让新人先手算一遍。取数组[3, 7, 11, 18, 22, 30, 41, 55, 63, 90]下标 0 到 9找目标 41用最朴素的闭区间写法轮次leftrightmida[mid]与目标比较下一轮区间109422小于 41答案在右侧[5, 9]259755大于 41答案在左侧[5, 6]356530小于 41答案在右侧[6, 6]466641命中结束四轮就把 10 个元素锁死了。注意第 3 轮里 mid 算出来是 5而不是 5.5 或者 6这个取整方向决定了区间能不能严格收缩。再看不存在的目标比如找 40最后一轮会走到 left6、right5循环条件left right不成立退出返回未找到。此时 left 的语义是第一个大于 40 的位置也就是 6值是 41而 right 是最后一个小于 40 的位置也就是 5值是 30。这两个位置的语义在后面插值查表时会变成主角所以现在就把它们记牢。2. 三种二分写法的边界差异与选型依据2.1 闭区间模板left right这是最常见的写法区间语义是答案可能在 [left, right] 里两端都还没被排除。int bsearch_closed(const int *a, int n, int target) { int left 0, right n - 1; while (left right) { int mid left (right - left) / 2; if (a[mid] target) return mid; if (a[mid] target) left mid 1; /* mid 已排除 */ else right mid - 1; /* mid 已排除 */ } return -1; /* 没找到 */ }三个细节值得说清楚。第一mid写成left (right - left) / 2而不是(left right) / 2不是为了好看是因为后者在 left 和 right 都很大时会整数溢出。32 位平台上 left 和 right 都接近 2^31 时相加直接翻负mid 变成负数然后你就眼睁睁看着数组越界访问。第二两个分支都写了mid ± 1因为闭区间写法里 mid 自己已经被比较过了能排除就一定要排除否则区间不收缩。第三left right里的等号不能丢因为它对应区间里还剩一个元素的情况丢了等号就会漏掉最后一个候选。2.2 左闭右开模板left right把右端改成开区间语义变成[left, right) 是候选范围right 本身不在范围里。int bsearch_halfopen(const int *a, int n, int target) { int left 0, right n; /* 注意这里是 n不是 n-1 */ while (left right) { int mid left (right - left) / 2; if (a[mid] target) left mid 1; else right mid; /* mid 保留 */ } return (left n a[left] target) ? left : -1; }这个写法我一开始很不习惯觉得右边都不在范围里了还怎么查。真正上手之后才发现它有两个好处一是区间收缩的公式只有一条right mid不容易写错方向二是循环结束时 left 和 right 相等而且这个位置天然就是第一个不小于 target 的下标语义非常干净拿来插值正合适。代价是循环结束后必须再判断一次a[left]到底等不等于目标不能直接返回。2.3 下界模板找第一个不小于目标的位置把 2.2 里的相等判断去掉就是标准的下界函数。它不告诉你有没有这个元素只告诉你从哪个位置开始元素就不小于目标了。/* 返回第一个满足 a[i] target 的下标可能等于 n */ int lower_bound(const int *a, int n, int target) { int left 0, right n; while (left right) { int mid left (right - left) / 2; if (a[mid] target) left mid 1; else right mid; } return left; }不变式是[0, left)里的元素全部小于 target[left, n)里的元素全部不小于 target。返回值等于 n 表示整张表都小于 target。这个语义为什么重要因为插值查表要的恰恰就是目标值夹在哪两个表项之间而下界给出的正是这个分界点的下标。找到下标 i 之后i-1 和 i 就是夹住目标的两个表项可以直接代入插值公式。这一段先记住第 4 节会直接用到。2.4 循环结束后 left 与 right 到底落在哪初学者最容易蒙的就是退出循环之后这两个变量还有意义吗。有而且意义很明确。以左闭右开模板为例循环结束的那一刻 left 一定等于 right因为区间长度为 0 才会退出并且这个位置满足左边所有元素都小于 target右边所有元素都不小于 target。所以 left 就是插入点right 也是插入点它俩指的是同一个位置。用这个视角回看闭区间模板退出时 left right 1left 是第一个大于 target 的位置right 是最后一个小于 target 的位置两者刚好把 target 应该待的缝夹在中间。搞清楚了这一点你就会发现很多找不到就返回 -1的写法其实是把有用信息丢掉了。我用一个表把三种模板的差异压缩一下方便贴在手边模板初始区间循环条件收缩方式退出后 left 的含义常用场景闭区间[0, n-1]left rightmid ± 1 双排除第一个大于目标的位置单纯判断存在性左闭右开[0, n]left right单侧保留插入点需要插入位置时下界[0, n]left right单侧保留第一个不小于目标的位置查表插值、分段定位3. NTC查表二分法不准的根因逐条拆解3.1 从ADC值到温度查表这一步到底在做什么先说清楚整条链路不然很容易把锅甩给二分。典型电路是 NTC 和一只上拉电阻串联分压NTC 在下、上拉在上、分压点接 ADC。上拉电阻取 Rp、ADC 位数为 N、满量程码值是FULL (1N) - 1那么从采样码值反算 NTC 阻值的式子是/* code 是 ADC 采样值Rp 是上拉电阻FULL 是满量程码值 */ int32_t r_ntc (int32_t)((int64_t)Rp * code / (FULL - code));拿到阻值之后第五步才是查表在温度-阻值表里找到这个阻值对应的温度。二分法只负责第五步中的定位动作它不产生任何数值误差。你说二分法不准其实是定位之后拿什么当结果不准。把这句话想明白排查方向就清晰了——问题一定出在定位的方向、定位的精度、或者数值类型这三处中的某一处。3.2 方向反了阻值-温度表是降序的NTC 的名字就叫负温度系数温度升高阻值下降。所以按温度升序排列的表阻值一列是严格递减的。而所有教科书模板、几乎所有标准库的bsearch()都默认数据是升序的。你要是把这只表直接丢给bsearch()比较函数里还写着a b返回负数那结果就是随机的——可能返回一个莫名其妙的下标也可能干脆返回 NULL。更隐蔽的是看起来能用。因为二分的排除逻辑在降序表上并不是完全随机崩溃有些查询会碰巧返回正确位置让你误以为代码没问题直到某一段温度区间开始偏才意识到方向错了。所以降序表必须用降序模板我一般写成这样/* r[] 按温度升序排列同时阻值严格递减返回第一个满足 r[i] rx 的下标 */ static int ntc_lower_bound_desc(const int32_t *r, int n, int32_t rx) { int left 0, right n; /* 答案范围 [0, n] */ while (left right) { int mid left (right - left) / 2; if (r[mid] rx) left mid 1; /* 中点在冷侧往热端走 */ else right mid; /* 中点在热侧保留 */ } return left; /* 可能是 n表示比所有表项都热 */ }判断方向的诀窍是不看符号看物理把r[mid]和rx的关系翻译成谁更冷比在脑子里推大小号靠谱得多。r[mid] rx表示表项阻值更大也就是表项更冷而我们要找的温度比它热所以往热端下标增大方向走。3.3 返回最近邻把连续的物理量当成量化的档位这是不准最主要的来源。很多老代码的逻辑是二分找到某个下标 i直接return temp_table[i]。这在温度上看就是档位制——表的步长是 25 度那读数只能取到 25 的倍数步长是 1 度读数就只能取整。表面上看 1 度精度也够用但下面的算例会告诉你即使在阻值域里找最近误差也远比你想的夸张。用常见的 10K、B3950 的 NTC 举例子理想模型下的阻值表25 度步长是这样的温度0 °C25 °C50 °C75 °C100 °C阻值33620 Ω10000 Ω3587 Ω1492 Ω698 Ω假设真实水温 60 °C模型阻值是 2487 Ω。二分在阻值域里找最近邻到 3587 Ω50 °C的距离是 1100 Ω到 1492 Ω75 °C的距离是 995 Ω最近的居然是 75 °C 那一项。于是程序老老实实输出 75 °C而实际是 60 °C偏差 15 度。这不是编程错误是数学上的必然——阻值随温度是近似指数下降的中间段极其陡峭在阻值域里等距地看温度域上是严重偏斜的。曲线越弯最近邻越容易往热端倒。第二层原因是灵敏度不均匀。高温段阻值变化慢同样的阻值差对应更大的温度差所以最近邻误差在热端被进一步放大。这解释了一个很典型的现象常温准得不行、一到 70 度以上就开始飘。不是你运气不好是物理特性决定的。3.4 整数与截断误差被放大到肉眼可见还有一类不准纯粹是数值处理惹的祸而且症状非常有辨识度——读数只在几个整数值之间跳永远看不到小数。第一个元凶是先除后乘。有人写插值的时候顺手写成temp t0 (r0 - rx) / (r0 - r1) * (t1 - t0)整数运算下(r0 - rx) / (r0 - r1)先被截断成 0结果恒等于t0。整段温度全变成表项的整数度而且永远取冷端那个值。修法很简单先乘后除中间结果用 64 位承住。第二个元凶是阻值类型选得太小。NTC 在低温段阻值可以是几十千欧加上上拉电阻后分压算出来的中间量很容易超过 16 位。我用int16_t存过一次阻值-10 °C 附近直接翻负查表拿到一个负数去二分输出一个疯狂的高温值排查了半天才反应过来是溢出。第三个元凶是表里出现重复或近重复的阻值。把厂家表按四舍五入取整到欧姆时低温段相邻两度可能取到同一个数于是插值公式的分母变成 0除法直接崩。这类问题在低温段尤其常见因为那里的曲线太平了。我的做法是生成表的时候加一道校验相邻两项差值为 0 就报警宁可手工改一个欧姆也不能留个地雷。4. 插值查表的实现把误差从十几度压到零点几度4.1 线性插值的公式推导与定点化思路很朴素二分只负责找到夹住目标的两个表项剩下的交给线性插值。设找到的下标关系是r[i] rx r[i1]降序表对应温度是t[i] t[i1]那么temp t[i] (r[i] - rx) * (t[i1] - t[i]) / (r[i] - r[i1])分子(r[i] - rx)非负分母(r[i] - r[i1])恒正公式在数值上是安全的。物理意义也直观把 rx 在[r[i1], r[i]]这段阻值区间里的相对位置映射到温度区间[t[i], t[i1]]上。因为阻值和温度在局部区间内是一一对应的单调关系这种映射不会出现多值。定点化有两个要点。第一温度统一用 0.01 °C 为单位也就是把 25.00 °C 存成 2500这样输出天然带两位小数不用浮点。第二中间乘积用int64_t。取个最坏情况阻值差最大约 3 万低温段温度差最大约 2500百摄氏度步长乘起来 7500 万32 位放得下但为了不给自己留隐患我直接上 64 位代价可以忽略。4.2 完整查表代码下面是能直接编译运行的版本包括降序二分、插值、边界钳位和异常兜底。表数据按 10 °C 步长给出用的是 B3950 的理想模型值实际项目请换成厂家给的 R-T 表。#include stdint.h #include stddef.h /* 温度单位 0.01 摄氏度阻值单位欧姆均按温度升序排列 */ static const int16_t s_temp_c100[] { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; static const int32_t s_res_ohm[] { 33620, 20175, 12536, 8037, 5302, 3587, 2487, 1760, 1270, 934, 698 }; #define NTC_TABLE_LEN (sizeof(s_res_ohm) / sizeof(s_res_ohm[0])) /* 返回 0 表示成功-1 表示阻值超出表的覆盖范围结果已钳位 */ int ntc_res_to_temp(int32_t r_meas, int16_t *temp_out_c100) { int n (int)NTC_TABLE_LEN; if (r_meas s_res_ohm[0]) { /* 比最冷端还冷钳到表头 */ *temp_out_c100 s_temp_c100[0]; return -1; } if (r_meas s_res_ohm[n - 1]) { /* 比最热端还热钳到表尾 */ *temp_out_c100 s_temp_c100[n - 1]; return -1; } /* 降序下界找第一个满足 s_res_ohm[i] r_meas 的下标 */ int left 0, right n; while (left right) { int mid left (right - left) / 2; if (s_res_ohm[mid] r_meas) left mid 1; else right mid; } /* 前面已经排除两端越界所以 left 一定落在 [1, n-1] */ int i left - 1; int32_t num (s_res_ohm[i] - r_meas); /* 0 */ int32_t den (s_res_ohm[i] - s_res_ohm[i 1]); /* 0 */ int32_t dt (int32_t)s_temp_c100[i 1] - s_temp_c100[i]; int64_t delta (int64_t)num * dt / den; /* 先乘后除注意顺序 */ *temp_out_c100 (int16_t)(s_temp_c100[i] delta); return 0; }有两处容易看漏的地方。一是钳位分支放在二分之前这样后面的二分就不必处理 left 等于 0 或者等于 n 的边界代码短一截也不容易出错。二是delta的除法发生在乘法之后且被除数是 64 位不会出现 3.4 节里那种恒等于冷端温度的截断。4.3 表项疏密的取舍与实测误差表越密插值误差越小代价是 Flash 占用和生成表的工作量。我用理想 B 模型估算过不同步长下的最大插值误差在区间中段取点实际以厂家表为准表步长覆盖 -40~125 °C 的表项数仅最近邻的误差插值后的误差两个数组的 Flash 占用25 °C约 8 项可达 15 °C约 3 °C约 100 字节10 °C约 18 项约 5 °C约 0.5 °C约 220 字节5 °C约 35 项约 2.5 °C约 0.15 °C约 420 字节1 °C约 166 项约 0.5 °C小于 0.01 °C约 2 千字节验证方式和数据都能对得上25 °C 步长下 60 °C 插值出来是 63.1 °C误差 3.1 °C10 °C 步长下 65 °C 插值出来是 65.5 °C误差 0.5 °C5 °C 步长下 62.5 °C 插值出来是 62.64 °C误差 0.14 °C。趋势非常清楚步长减半误差大致降到四分之一因为插值误差和步长的平方成正比。这里有个很实际的取舍判断。当表步长降到 5 °C插值误差已经到 0.15 °C 量级而常见的 NTC 传感器本身公差是 R25 ±1%、B 值 ±1%换算到温度上是零点几度到一度。也就是说算法误差已经被传感器误差盖住了再加密表就是白费 Flash。我的经验线是插值后的算法误差做到传感器公差的三分之一以下就可以收手了5 °C 步长对绝大多数消费级和工业级场景刚好卡在这条线上。4.4 边界钳位与异常兜底查表函数最容易在边界上翻车。真实阻值可能因为以下原因跑到表外传感器断线阻值无穷大、短路阻值接近 0、接线氧化导致接触电阻偏大、低温环境超出表的下限。这些情况如果直接喂给插值公式分母可能为 0或者算出下标 -1 直接越界读数组程序的状态就完全不可预期了。我习惯在函数出口统一约定返回值语义0 表示正常-1 表示钳位。调用方看到 -1 就知道该走故障处理流程了——断线报故障、超范围报警、或者按场景降级使用。这样既不阻塞正常的温度读取又不会把异常数据悄悄混进控制回路里。还有一类兜底是排序稳定性。表是工具生成的万一手抖把某两行顺序写反单调性就破了二分会给出一个看起来正常的错误结果。我的做法是在初始化阶段跑一次全表校验检查阻值列是否严格递减、温度列是否严格递增、相邻阻值差是否都大于某个阈值比如 5 欧姆。这个校验只有几十次循环上电时跑一次几乎不花时间但能挡住绝大多数低级错误。5. 一次完整排障从常温准高温偏到定位修复5.1 现象记录与最小复现先把现象钉死不然排查会一直在猜。当时记录到的是室温 25 °C 附近读数偏差小于 0.5 °C55 °C 以上开始明显偏大且偏得没规律最诡异的线索是读数只在几个固定值之间跳比如 60、65、70 这种整五度中间值永远看不到。把这三条放一起指向性其实很强——只跳固定值说明根本没有插值只用表项温度当结果热端偏大说明最近邻在阻值域里往热端倒也就是 3.3 节那个机制。复现方式也很简单不用真的去搭恒温槽直接绕开 ADC把一组已知阻值的电阻或者高精度电阻箱接到分压点位置手工喂给ntc_res_to_temp()让函数把返回值打印出来。这是最快的定位手段因为它把 ADC 本身的误差、参考电压漂移、走线噪声全部剔除了剩下的只有查表逻辑。5.2 分层定位把链路切成三段分别验证链路是分压 → ADC → 反算阻值 → 查表 → 温度我把它切成三段独立验证每段只留一个变量。第一段是 ADC 到阻值。拿万用表直接量分压点的电压按分压公式算出理论阻值再和代码里反算出来的阻值对比。这一步的目的是确认 ADC 采样值和换算公式都没问题。当时量下来两者差不到 1%说明前两段是干净的。第二段是阻值到温度。把上一步得到的阻值丢到厂家提供的 R-T 表里用计算器手查得到的温度和标准温度计差 0.3 °C说明表数据本身没问题。第三段是查表代码。把ntc_res_to_temp()单独抠出来喂入几个已知阻值把left、right、mid、r[mid]全打印出来。马上就看到了问题函数返回的温度全是表项上的整度值delta恒等于 0。再翻代码delta那行的写法是(num / den) * dt——先除后乘整数除法把 num/den 截成了 0。与此同时二分定位到的下标也确实偏向了热端因为表步长是 25 °C 的粗表最近邻本来就会倒向热端。嵌入式环境里打印变量有时候不方便替代方案是用一个环形缓冲把left/right/mid记下来跑完之后一次性 dump或者用一根空闲 GPIO 翻转电平接逻辑分析仪看波形。都比盲猜快。5.3 修复与验证怎么确认不是看起来对了改动只有三处把升序模板换成 3.2 节的降序模板把先除后乘改成先乘后除并加 64 位中间量补上边界钳位和返回值语义。改完之后读数只跳整度和热端偏大两个症状同时消失这本身就说明方向找对了。但我不太信任改完看着对了这种验证方式。真正让我放心的是两组等价性测试。第一组是逐点比对写一个线性扫描版本从表头开始一项一项找配上同样的插值公式然后把二分版本和线性版本在 -50 到 150 °C 之间按 1 Ω 步长全部跑一遍要求两者输出完全一致一个 bit 都不能差。第二组是随机模糊测试随机生成十万个阻值包含越界值同样要求两个版本结果一致。二分和线性扫描在逻辑上是等价的任何不一致都说明二分模板有问题这个测试能覆盖掉绝大多数边界 bug。还有一组针对性的边界测试值得单独跑喂入比表头还冷的阻值比如 1 MΩ和比表尾还热的阻值比如 10 Ω检查返回的是钳位值加 -1而不是越界下标。这两种情况在传感器断线和短路时是真实会出现的测试里必须覆盖。5.4 同类问题清单排完这一趟我把这个项目里所有查表不准的坑都整理了一遍后来发现它们反复出现在几乎每一个用 NTC 的项目里症状最可能的原因快速验证方法读数只在固定值上跳没有插值直接返回表项看返回温度是否总是表步长的整数倍热端偏大、冷端还行表方向搞反或最近邻偏向热端用电阻箱喂已知阻值对照厂家表返回值恒等于冷端温度整数除法先除后乘看插值中间量是否为 0低温段直接算出巨大异常值阻值用 16 位溢出或下标越界查数据类型和数组下标范围插值算到一半崩溃表里有重复阻值分母为 0全表扫描相邻差值换了上拉电阻后整体偏移表按另一个上拉阻值生成核对生成表时的分压参数高低温都有零点几度固定偏移传感器公差不是算法问题用高精度电阻箱替掉传感器再测最后那一行值得多说一句。很多人把算法误差压到 0.01 °C 之后还在纠结为什么实测差 0.5 °C然后继续改代码越改越离谱。其实这时候误差已经转移到传感器公差和自发热上了——自发热尤其容易被忽略如果分压电流偏大比如上拉用了 1 kΩNTC 自身会被加热读数系统性偏低而且这个偏差和功耗、散热、封装都相关绝不是改代码能解决的。6. 换个场景二分法求方程近似根6.1 从查找问题到求根问题的转换数学教材里二分法这个名字更多指的是求方程近似根和查表定位是同一套思想的两个投影。给定连续函数 f(x)如果在区间 [a, b] 上 f(a) 和 f(b) 异号那么区间内至少有一个根。取中点 m看 f(m) 的符号它和 f(a) 同号就把 a 移到 m异号就把 b 移到 m区间长度每次减半若干轮之后中点就是根的一个近似值。本质上它和二分查找干的是同一件事靠一个可判定的性质符号每次砍掉一半的搜索空间。区别在于查找的判定是等于/大于/小于求根的判定是同号/异号但收缩区间的机制完全一致。6.2 迭代终止条件与精度换算精度和迭代次数是可以直接算出来的。初始区间长度是b - a每轮减半k 轮之后长度是(b-a)/2^k。要求最终精度不超过 eps解不等式得到k log2((b - a) / eps)举个数区间长度 2目标精度 1e-6那么 k log2(2×10^6) ≈ 20.93取 21 轮。二十来次迭代就能把一个两位数精度的答案直接压到六位小数这就是指数收敛的力量。反过来说如果你发现有人跑了一万次迭代那说明终止条件写错了不是收敛慢的问题。#include math.h double bisect_root(double (*f)(double), double a, double b, double eps, int max_iter) { double fa f(a), fb f(b); if (fa * fb 0.0) return NAN; /* 端点同号不能保证有根 */ for (int i 0; i max_iter (b - a) eps; i) { double m 0.5 * (a b); double fm f(m); if (fm 0.0) return m; /* 运气好正中根上 */ if (fa * fm 0.0) { /* 根在左半段 [a, m] */ b m; } else { /* 根在右半段 [m, b] */ a m; fa fm; /* 复用避免重复调 f */ } } return 0.5 * (a b); }代码里有两个设计细节值得展开。第一fa被缓存下来每次只需要算一次f(m)而不是两次f(a)和f(m)。函数调用次数减半在这个场景下是很划算的优化因为 f 往往是整个流程里最贵的一环。第二终止条件用的是区间长度(b - a) eps而不是fabs(fm) eps。这两者差别很大当函数在根附近非常平缓导数接近 0时x 偏离根很远但 f(x) 已经很小了用函数值做终止条件会提前退出给出的 x 完全不准。区间长度是自变量空间的度量用它才算真正控制了根的精度。用浮点做二分会遇到一个理论上的极限当区间长度小到相邻两个可表示浮点数之间的距离ulp时m算出来只会等于 a 或 b区间不再收缩就陷入死循环了。所以max_iter这个保险必须加不能只依赖 eps。工程上一般给 100 到 200 轮上限足够覆盖 double 的全部有效位。6.3 求根二分的局限与替代方案二分求根最大的优点是稳缺点是慢。它每轮只把误差减半属于线性收敛大约每 3.3 轮才多出一位有效数字。相比之下牛顿法在根附近是平方收敛——每轮有效位数翻倍通常四到六轮就能到机器精度。但牛顿法需要导数而且初值选得不好会发散、会跑到别的根上去甚至在迭代中跨过不连续点直接飞掉。我实际的做法是混合先用二分把区间收缩到足够小比如收缩到原区间的千分之一再切到牛顿法或者割线法快速收尾。这样既有二分的全局保证只要端点异号一定能收敛又有牛顿法的末段速度。交换机内燃机标定、传感器线性化系数反解、温控回路的稳态点求解这类场景这个策略都很好用。还有几个信号很明确的失效情形需要提前防住。如果 f(a) 和 f(b) 同号二分根本启动不了这时候得先扫描一遍找变号区间不能硬塞进去。如果区间内有偶数个根端点异号并不能保证你收敛到想要的那个根可能收敛到旁边的根上。如果函数在某处不连续二分的收缩可能直接把根跳过。这些都是数学前提的问题不是代码能补救的。最后再分享一个我在实际操作中的小体会无论用二分做哪种事我都强迫自己先把不变式写在注释里——循环开始时[0, left) 全部小于目标[left, n) 全部不小于目标或者循环开始时f(a) 与 f(b) 异号。写完注释再写代码边界条件的错误率能降一个数量级。二分法的代码短到可以背下来但正是因为它短一旦写错就特别难看出来而不变式是唯一能在不看运行结果的情况下判断对错的东西。这个习惯我从查表这件事上养起来后来写别的算法也一直在用。
网站建设高端定制企业官网