一种简单的归约区间方法

文章背景与核心概要

计算大数字的正弦或余弦值需要进行“归约区间”（Range Reduction）——即将输入值 \(x\) 映射到一个更小、更易处理的区间（通常为 \([0, \pi/2]\)）中，同时不丢失精度。本文探讨了 Cody 和 Waite 方法，该方法利用高精度常数，在中等大小的自变量情况下实现准确计算，而这是朴素方法无法做到的。

文章深入浅出地介绍了三角函数计算中的浮点数精度问题。通过对比朴素的模运算与 Cody 和 Waite 提出的双常数拆分法，展示了如何通过精巧的操作顺序保留关键的有效数字。对于从事科学计算和数值优化的开发者而言，这是一篇理解底层浮点数精度处理机制的极佳入门指南。

模 \(\pi/2\) 的归约区间

为了高效计算 \(\sin(x)\) 或 \(\cos(x)\)，我们将 \(x\) 归约到 \([0, \pi/2]\) 区间。如果我们能找到一个整数 \(k\)，使得 \(x - k(\pi/2) = y\)（其中 \(0 \le y \le \pi/2\)），那么三角恒等式将取决于 \(k \pmod 4\) 的值：

Range Reduction mod \(\pi/2\)

To compute \(\sin(x)\) or \(\cos(x)\) efficiently, we reduce \(x\) to a range \([0, \pi/2]\). If we find an integer \(k\) such that \(x - k(\pi/2) = y\) (where \(0 \le y \le \pi/2\)), the trigonometric identity depends on \(k \pmod 4\):

from math import *

def reduced_sin(x, k):
     match k % 4:
        case 0: return sin(x)
        case 1: return cos(x)
        case 2: return -sin(x)
        case 3: return -cos(x)

朴素的归约区间方法

朴素的方法使用标准的浮点数除法来寻找 \(k\)。例如，当 \(x = 500\) 时：

Naive Range Reduction

A naive approach uses standard floating-point division to find \(k\). For example, with \(x = 500\):

def naive_sin(x):
    k = floor(x / (pi/2))
    y = x % (pi/2)
    return reduced_sin(y, k)

When \(x = 500\), this method yields an error on the order of \(1.7 \times 10^{-14}\), which is insufficient for high-precision applications.

更好的归约区间方法：Cody 和 Waite 法

Cody 和 Waite 方法通过将 \(\pi/2\) 表示为两个常数 \(C_1\) 和 \(C_2\) 的和来提高精度，这两个常数组合起来提供的精度超过了单个 64 位浮点数所能存储的范围。

Better Range Reduction: The Cody and Waite Method

The Cody and Waite method improves accuracy by representing \(\pi/2\) as the sum of two constants, \(C_1\) and \(C_2\), which together provide more precision than a single 64-bit float can store.

def Cody_Waite_sin(x):
    C1 = 1686629713 / 2**30
    C2 = 4701928774853425 / 2**86

    k = floor(x / (pi/2))
    y = (x - k*C1) - k*C2
    return reduced_sin(y, k)

为什么它有效

秘诀在于运算的顺序。虽然从数学上讲 \(C_1 + C_2\) 接近 \(\pi/2\)，但如果计算 y = x - k*(C1 + C2)，将会导致与朴素方法相同的精度损失。通过计算 (x - k*C1) - k*C2，我们保留了原本会在浮点数舍入中丢失的有效数字。

Why It Works

The trick lies in the order of operations. While \(C_1 + C_2\) is mathematically close to \(\pi/2\), computing y = x - k*(C1 + C2) would result in the same precision loss as the naive method. By calculating (x - k*C1) - k*C2, we preserve significant digits that would otherwise be lost to floating-point rounding.

测试精度

你可以将这些方法的准确度与系统内置的 math.sin（它使用了高度优化的硬件级归约区间）进行比较：

Testing Precision

You can compare the accuracy of these methods against the system's built-in math.sin (which uses highly optimized, hardware-level range reduction):

def compare(x):
    y0 = naive_sin(x) 
    y1 = Cody_Waite_sin(x)
    y2 = sin(x)
    print("Naive error:     ", y2 - y0)
    print("Cody Waite error:", y2 - y1)

虽然 Cody 和 Waite 算法不适用于极其巨大的 \(x\) 值，但它是一篇极好的入门教程，展示了如何利用巧妙的技术从浮点算术中榨取最大的精度。

While the Cody and Waite algorithm is not suitable for extremely large values of \(x\), it serves as an excellent introduction to the clever techniques used to squeeze maximum precision out of floating-point arithmetic.

注：\(C_1\) 的分子 \(n_1\) 为 \(\lfloor 2^{30} \pi/2 \rfloor\)。\(C_2\) 的分子 \(n_2\) 是方程 \(2^{86-30} n_1 + n_2 = \lfloor 2^{86} \pi/2 \rfloor\) 的解。

来源：John D. Cook 的 A simple range reduction method。

---

Note: The numerator \(n_1\) of \(C_1\) is \(\lfloor 2^{30} \pi/2 \rfloor\). The numerator \(n_2\) of \(C_2\) is the solution to \(2^{86-30} n_1 + n_2 = \lfloor 2^{86} \pi/2 \rfloor\).

Source: A simple range reduction method by John D. Cook.