OCaml浮点拆分求和与全量求和结果不一致问题咨询

问题描述 #

我目前正在进行OCaml代码开发，首先编写了用于计算正整数倒数和的函数inverseSum：

let inverseSum n =
  let rec sI n acc =
    match n with
    | 1 -> acc
    | _ -> sI (n - 1) ((1.0 /. float n) +. acc)
  in sI n 1.0;;

例如inverseSum 2的计算逻辑为1/2 + 1 = 3/2 = 1.5，使用参数2、5测试该函数，运行结果符合预期：

inverseSum 2;;
inverseSum 5;;

(* 运行返回 *)
inverseSum 2;;
- : float = 1.5
inverseSum 5;;
- : float = 2.28333333333333321

该阶段运行无异常。随后我初始化了包含1到10000所有整数的列表initList：

let initList = List.init 10000 (fun n -> n + 1);;

该步骤运行正常。之后编写invSumLst函数，将列表中每个元素替换为对应值的倒数和计算结果（例：[1;2;3]转换为[inverseSum 1; inverseSum 2; inverseSum 3]）：

let rec invSumLst lst =
  match lst with
  | [] -> []
  | h::t -> (inverseSum h) :: invSumLst t;;

将该函数应用于initList得到invInit列表，此阶段运行无异常。后续我执行了如下操作：

• 筛选invInit中严格小于5.0的元素得到listLess5，通过List.fold_left对该列表求和得到foldLess5
• 筛选invInit中大于等于5.0的元素得到moreEg5，通过List.fold_left对该列表求和得到foldMore5
• 直接通过List.fold_left对全量invInit列表求和得到foldInvInit

对应代码如下：

let listLess5 = List.filter (fun n -> n < 5.0) invInit;;
let foldLess5 = List.fold_left (+.) 0.0 listLess5;;
let moreEg5 = List.filter (fun n -> n >= 5.0) invInit;;
let foldMore5 = List.fold_left (+.) 0.0 moreEg5;;
let foldInvInit = List.fold_left (+.) 0.0 invInit;;

最后计算两个拆分和相加与全量和的绝对误差时，得到了意料之外的结果：

Float.abs ((foldLess5 +. foldMore5) -. foldInvInit);;
Printf.printf "%f\n" (Float.abs ((foldLess5 +. foldMore5) -. foldInvInit));;
Printf.printf "%b\n" ((foldLess5+.foldMore5) = foldInvInit);;

运行返回结果如下：

let foldMore5 = List.fold_left (+.) 0.0 moreEg5;;
val foldMore5 : float = 87553.6762998474733
let foldInvInit = List.fold_left (+.) 0.0 invInit;;
val foldInvInit : float = 87885.8479664799379
Float.abs ((foldLess5 +. foldMore5) -. foldInvInit);;
- : float = 1.45519152283668518e-11
Printf.printf "%f\n" (Float.abs ((foldLess5 +. foldMore5) -. foldInvInit));;
0.000000
- : unit = ()
Printf.printf "%b\n" ((foldLess5+.foldMore5) = foldInvInit);;
false
- : unit = ()

我推测该现象可能和舍入问题相关，但想明确误差的具体来源：在解释器环境中可观测到误差值为1.45519152283668518e-11，但如果使用ocamlpro这类编译器运行，终端只会输出0.000000和相等判断为false的结果，难以定位问题。我想确认该问题是源于代码中某个函数的逻辑缺陷，还是浮点数本身的运算舍入问题，亦或是Printf.printf使用非科学计数法输出时产生的格式舍入导致。

问题结论 #

代码不存在逻辑缺陷，观测到的现象来自两个浮点数相关的常规机制：

• 浮点数加法不满足结合律，求和顺序不同会产生不同的舍入累积误差
  OCaml的float类型是标准64位双精度浮点数，仅能保证约15-17位十进制有效精度，每次加法运算后都会对结果做舍入处理。全量求和是按列表原有顺序依次累加所有元素，拆分求和是先累加所有小于5.0的元素、再累加大于等于5.0的元素，两种路径的加法顺序完全不同，累积的舍入误差自然存在差异。观测到的1.45e-11的误差相对于8万多的求和总值，相对误差不到1e-16，已经接近双精度浮点数的理论精度极限，属于完全正常的现象。
• %f格式符默认输出精度不足，会隐藏微小误差
  Printf.printf的%f格式符默认仅保留小数点后6位，1.45e-11远小于0.0000005，四舍五入后自然会显示为0.000000。但这只是输出显示层面的格式舍入，不会修改内存中存储的实际浮点值，所以相等判断依然会返回false。如果需要观测到真实的误差值，可以换用%e科学计数法格式符，或者指定更高的输出精度比如%.17f，就能看到完整的误差数值。

内容的提问来源于stack exchange，提问作者mxbr236