本文代码已在 MATLAB R2025b 25.2.0.2998904、PCWIN64 环境原生运行，并完成独立原生复核。最终发布源码在 2026-11-06 01:06:40–01:07:13 UTC 再次完整回归。以下区分作者测试、独立测试和产物检查，不将重复运行叠加计数。
1. 大数旁边的 1，可能没有留下来
想象一个累计量先达到 2^54，期间又加了若干个单位增量，最后减去 2^54。数学上结果就是单位增量的个数；如果每个 1 都在与大数相加时被舍掉，最后得到的却可能是 0。
这类问题会出现在大基线上的微小修正、正负贡献的汇总和离线统计审计里。加长显示位数不会增加有效精度。MathWorks 的 浮点数说明介绍了 double 的有限表示与舍入；本文通过可精确核算的输入，区分输入已丢失的信息和计算过程中新增的误差。
2. 三种结果不要混在一起
第一种是理想实数数据的总和。第二种是已经存成 double 的这批数的精确总和。第三种是某个浮点求和算法的输出。
例如先在 double 中形成 2^54+1，这个 1 可能在调用求和函数之前就不见了。之后任何补偿循环都看不到它。本文的主要测试把大数和小数作为独立元素输入，测的是求和算法能否减少额外舍入误差，不声称恢复测量或转换时已损失的信息。
也不要假定 MATLAB 的内置 sum 就是自己写的从左到右循环。MathWorks 的开发者文章 A faster and more accurate sum讨论过内置实现的改进。因此本文专门写一个显式顺序循环作为负例，内置 sum 只记录实际返回值，不断言它必须失败，也不借本文样例做性能排名。
3. 补偿项记录了什么
普通循环只保留 s。Neumaier 形式还保留补偿 c：先计算 t=s+a，再根据 s 与 a 的绝对值大小，用较大的一项去构造这次加法的修正量。
- 若 abs(s)>=abs(a)，修正量为
(s-t)+a。 - 否则，修正量为
(a-t)+s。 - 将修正量累加进 c，把主和更新为 t，最后返回 s+c。
大小比较不是装饰。当新来的 a 比当前主和更大时，它决定了修正表达式的顺序。不要擅自把括号重排或把两个分支合并成“代数上等价”的式子，浮点舍入依赖实际运算过程。
算法有明确的历史来源，并非本文发明。Ogita、Rump、Oishi 的原始研究论文 Accurate Sum and Dot Product在相关工作部分讨论了大小排序的两数补偿和 Neumaier 变体。这里给出独立编写的教学实现，不将论文中的其他高精度算法或定理保证套到这个简短函数上。
主循环一次扫描、额外状态常数个。补偿项自身也用 double 累加，最后 s+c 仍会舍入。因此“所有输入都能精确求和”“必然正确舍入”“顺序完全无关”都不是本文承诺。
4. 用什么当作独立答案
测试不拿 MATLAB sum 当正确答案，而使用构造出来的精确整数总和：
[2^54,1,-2^54]的精确目标是 1。[2^54,1,-2^54,3]的 24 个排列，精确目标均为 4。[2^54,ones(1,N),-2^54]的精确目标为 N，本文只取 N=1、2、3、17、256、4096。
最后一组还用 int64 顺序累加核对。这里只能这么做，是因为输入全为可精确表示的整数，且每一步均远小于 int64 上限；对一般浮点数据强行转 int64 会改变问题，不能当作高精度参考。
负例要求朴素循环至少在部分排列中失败，用来证明测试能抓到目标问题。输入舍入负例要求已经丢失的 1 不会被“神奇找回”。溢出负例 [realmax,realmax,-realmax] 即使实数总和可表示，也会在中间求和溢出，本文函数明确报错。
5. 完整实现与测试入口
保存下面三个文件。实现只接受有限、实数、非稀疏 double 向量；空输入返回 0，不支持矩阵按维度求和、复数、single、NaN 忽略语义或 GPU。输入、补偿和最终结果的非有限情况会显式拒绝。
本次原生环境为 MATLAB R2025b，版本 25.2.0.2998904，PCWIN64。对最终 11 个作者 .m 文件执行 requiredFilesAndProducts，静态依赖结果仅识别 MATLAB 产品；本机安装有其他工具箱，本次未进行仅安装 MATLAB 的干净环境试验，也未做许可证调用跟踪。因此这里只报告已执行环境与静态依赖结果，不扩大成所有环境兼容保证。
neumaier_sum.m
function total=neumaier_sum(x)
%NEUMAIER_SUM Compensated summation for finite real full double vectors.
% Returns zero for empty input. Rejects intermediate nonfinite results.
% This is not an exact accumulator or a correctly-rounded-sum guarantee.
if ~(isa(x,'double') && isreal(x) && ~issparse(x) && ...
(isvector(x) || isempty(x)) && all(isfinite(x(:))))
error('neumaier:Input','x must be a finite real full double vector, or empty.');
end
s=0; c=0;
for k=1:numel(x)
a=x(k); t=s+a;
if ~isfinite(t), error('neumaier:Overflow','Intermediate sum overflowed.'); end
if abs(s)>=abs(a)
correction=(s-t)+a;
else
correction=(a-t)+s;
end
c=c+correction;
if ~isfinite(c), error('neumaier:Overflow','Compensation overflowed.'); end
s=t;
end
total=s+c;
if ~isfinite(total), error('neumaier:Overflow','Final sum overflowed.'); end
end
test_neumaier_sum.m
function result=test_neumaier_sum()
% Tests use exact integer fixtures rather than sum as the oracle.
n=0; big=2^54;
check(neumaier_sum([big 1 -big])==1,'lost unit is recovered');
check(neumaier_sum([1 big -big])==1,'magnitude swap branch');
check(neumaier_sum([big -big 1])==1,'cancellation before small term');
check(neumaier_sum([])==0 && neumaier_sum(7)==7,'empty and singleton');
check(neumaier_sum([big;1;-big])==1,'column vector');
check(neumaier_sum([1 -1])==0,'exact zero sum');
P=perms([big 1 -big 3]); naiveFailures=0;
for k=1:size(P,1)
check(neumaier_sum(P(k,:))==4,'all 24 four-term permutations');
naiveFailures=naiveFailures+(sequential_sum(P(k,:))~=4);
end
check(naiveFailures>0,'plain left-to-right mutant is detected');
for count=[1 2 3 17 256 4096]
x=[big ones(1,count) -big];
exact=int64(0);
for k=1:numel(x), exact=exact+int64(x(k)); end
check(neumaier_sum(x)==double(exact),'bounded int64 reference');
end
% The unit is already lost before the function receives its input.
x=[big+1 -big];
check(x(1)==big && neumaier_sum(x)==0,'cannot recover input rounding');
expect_error(@()neumaier_sum([NaN 1]),'neumaier:Input');
expect_error(@()neumaier_sum([Inf 1]),'neumaier:Input');
expect_error(@()neumaier_sum([1 1i]),'neumaier:Input');
expect_error(@()neumaier_sum(single([1 2])),'neumaier:Input');
expect_error(@()neumaier_sum(ones(2)),'neumaier:Input');
expect_error(@()neumaier_sum(sparse([1;2])),'neumaier:Input');
expect_error(@()neumaier_sum([realmax realmax -realmax]),'neumaier:Overflow');
result=struct('status','PASS','checks',n,'naive_failed_permutations',naiveFailures, ...
'builtin_sum_observation_only',sum([big 1 -big]));
fprintf('neumaier_sum PASS checks=%d; builtin sum observation=%.17g\n',n, ...
result.builtin_sum_observation_only);
function check(ok,label)
assert(ok,'test:Failed','%s',label); n=n+1;
end
function expect_error(fn,id)
found=false;
try, fn(); catch ME
if ~strcmp(ME.identifier,id), rethrow(ME); end
found=true;
end
check(found,['expected error ' id]);
end
end
function s=sequential_sum(x)
s=0;
for k=1:numel(x), s=s+x(k); end
end
demo_neumaier_sum.m
function demo_neumaier_sum(outdir)
counts=[1 2 3 17 256 4096]; naive=zeros(size(counts)); compensated=naive;
builtin=naive; big=2^54;
for j=1:numel(counts)
x=[big ones(1,counts(j)) -big]; s=0;
for k=1:numel(x), s=s+x(k); end
naive(j)=s; compensated(j)=neumaier_sum(x); builtin(j)=sum(x);
end
P=perms([big 1 -big 3]); perNaive=zeros(size(P,1),1); perComp=perNaive;
for j=1:size(P,1)
for k=1:4, perNaive(j)=perNaive(j)+P(j,k); end
perComp(j)=neumaier_sum(P(j,:));
end
fig=figure('Visible','off','Color','w','Position',[100 100 1000 720]);
clean=onCleanup(@()close(fig));
subplot(2,1,1); semilogx(counts,counts,'k-',counts,naive,'rs--', ...
counts,compensated,'bo-',counts,builtin,'md:','LineWidth',1.3); grid on;
xlabel('Number of unit increments'); ylabel('Final sum after cancelling 2^{54}');
legend('Exact integer target','Sequential','Compensated','Built-in: observed only', ...
'Location','northwest'); title('Small increments beside a large running total');
subplot(2,1,2); plot(1:24,perNaive-4,'rs-',1:24,perComp-4,'bo-','LineWidth',1.2);
grid on; xlabel('Permutation index'); ylabel('Error relative to exact sum 4');
legend('Sequential','Compensated','Location','best');
title('Order sensitivity in a deliberately ill-conditioned fixture');
print(fig,fullfile(outdir,'03_compensated_sum.png'),'-dpng','-r160');
save(fullfile(outdir,'03_compensated_sum_data.mat'),'counts','naive','compensated', ...
'builtin','P','perNaive','perComp');
end
运行命令
在以上文件所在目录运行以下命令。使用新目录保留每次证据；出错时保留失败日志，不把失败记录改成通过。
outdir=fullfile(pwd,['native_run_' datestr(now,'yyyymmddTHHMMSS')]);
assert(exist(outdir,'dir')~=7,'Please use a new evidence directory.');
mkdir(outdir);
diary(fullfile(outdir,'native.log'));
fprintf('MATLAB: %s\nPlatform: %s\n',version,computer);
ver
result=test_neumaier_sum();
demo_neumaier_sum(outdir);
save(fullfile(outdir,'native_result.mat'),'result');
diary off
6. 看绝对误差，也看问题本身的敏感性
图一比较不同单位增量个数下的结果，图二显示 24 个排列相对于精确目标 4 的误差。本次内置 sum 数据来自上述 MATLAB 环境的实际运行，没有用 Python 或其他版本结果代填。图像已完成原生回读和独立视觉复核，MAT 数据完成独立重算；本文保留生成图代码，不嵌入实验图片。
如果精确总和接近 0，直接计算相对误差可能非常大或没有定义。工程报告应同时列出绝对误差、数据尺度及允许误差。常见的求和条件量是 sum(abs(x))/abs(exact_sum)，它解释正负巨大抵消为什么敏感；分母为 0 时不能给出一个普通有限值，更不能因为相对误差不好看就删掉该用例。
- 作者原生检查：45 项 API 断言通过，另有 3 项导出检查通过。24 个指定排列均满足精确整数目标，6 组单位累加用例也通过。
- 独立原生检查：1038 项 API 断言、9 项产物断言通过。
- 独立取消补偿的变异，在指定的三元素求和用例中得到 0，而精确目标为 1，运行时断言按预期失败。
- 本机实际观察：sum([254,1,-254]) 返回 0；这只是 MATLAB R2025b 本次环境的记录，不承诺其他实现或版本返回相同值，也不代表内置 sum 在所有数据上较差。
- 最终发布副本的整批回归退出码为 0，stderr 为空。
7. 什么时候该换更强的方法
补偿求和适合希望在普通浮点精度内减少累加误差、并且能接受额外运算的场景。若需求是跨硬件逐比特一致、严格正确舍入、任意精度、跨越溢出范围，或要给出可靠的误差界，就应选择满足这些要求的专门累加器或更高精度方法，并单独验证。
本文也没有把误差更小等同于测量更准确。传感器偏置、量化和模型误差仍然存在。先明确要对“哪一批已经存储的数据”求和，再决定是否需要补偿，才不会让更漂亮的数字掩盖数据本身的问题。
182
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
