双精度浮点数幂运算舍入误差问题及Common Lisp解决方案咨询

问题背景 #

我正在为LeetCode第50题“Pow(x, n)”（实现pow(x, n)，计算x的n次幂xⁿ）开发Common Lisp解决方案，当前的直接实现代码如下：

(defun pow (x n)
  (declare (ftype (function (double-float fixnum) double-float) pow))
  (if (zerop n)
      1d0
      (let ((r 1d0))
        (dotimes (i (abs n)) (setf r (* r x)))
        (if (plusp n)
            r
            (/ 1d0 r)))))

(pow 2.1d0 3)
=> 9.261000000000001d0

为何结果小数点后第15位会出现1？如何消除该舍入误差？运行环境为Linux下的SBCL 2.3.4。

误差原因 #

这个问题的核心在于双精度浮点数（double-float）的二进制存储特性：

• 十进制的2.1无法被二进制浮点数精确表示，它在二进制中是无限循环小数，存储时会被截断为最接近的可表示近似值。
• 连续三次乘以这个不精确的近似值时，每次乘法的微小误差会累积，最终导致结果末尾出现额外的1。

解决方法 #

方法1：使用精确十进制运算（适合高精度场景） #

Common Lisp支持十进制浮点数类型（decimal-float），可以避免二进制浮点数的精度问题。修改代码使用十进制类型：

(defun pow-decimal (x n)
  (declare (ftype (function (decimal-float fixnum) decimal-float) pow-decimal))
  (if (zerop n)
      1.0d0M  ; 十进制浮点数字面量
      (let ((r 1.0d0M))
        (dotimes (i (abs n)) (setf r (* r x)))
        (if (plusp n)
            r
            (/ 1.0d0M r)))))

(pow-decimal 2.1d0M 3)
=> 9.261d0M

注意：十进制浮点数运算速度通常比二进制双精度慢，仅适合对精度要求极高的场景。

方法2：格式化输出时截断或四舍五入 #

如果只需要在展示时消除误差，可以对结果进行格式化，保留合理的小数位数：

; 保留12位小数输出
(format t "~,12F~%" (pow 2.1d0 3))
=> 9.261000000000

; 四舍五入到指定精度
(round (* (pow 2.1d0 3) 1d12)) 1d12
=> 9.261d0

方法3：用快速幂减少误差累积（同时优化效率） #

当前线性遍历的实现不仅效率低（时间复杂度O(n)），还会增加误差累积次数。使用快速幂算法（时间复杂度O(log n)）可以减少乘法次数，从而降低误差累积：

(defun fast-pow (x n)
  (declare (ftype (function (double-float fixnum) double-float) fast-pow))
  (cond
    ((zerop n) 1d0)
    ((minusp n) (/ 1d0 (fast-pow x (- n))))
    ((oddp n) (* x (fast-pow x (- n 1))))
    (t (let ((half (fast-pow x (/ n 2))))
         (* half half)))))

(fast-pow 2.1d0 3)
=> 9.261000000000001d0

注：快速幂虽减少了乘法次数，但由于2.1本身的二进制表示不精确，结果仍会有微小误差，结合格式化输出即可解决展示问题。

内容的提问来源于stack exchange，提问作者Jérôme Radix