Fast Trigonometric Functions using the RLIBM Approach

基本信息

• 论文ID: 2510.13426
• 标题: Fast Trigonometric Functions using the RLIBM Approach
• 作者: Sehyeok Park, Santosh Nagarakatte (Rutgers University)
• 分类: cs.PL (Programming Languages)
• 发表会议: International Workshop on Verification of Scientific Software (VSS 2025)
• 论文链接: https://arxiv.org/abs/2510.13426

摘要

本文描述了使用RLIBM方法开发三角函数多项式逼近的经验，该方法能够为多种表示和舍入模式产生正确舍入的结果。三角函数的关键挑战在于涉及π的范围缩减，它将32位浮点数域中的输入缩减到小域。π值中的任何舍入误差在范围缩减过程中会被放大，可能导致错误结果。作者描述了实现快速范围缩减技术的经验，这些技术在浮点和整数计算中都能维护大量π的位数。最终的三角函数实现既快速又能为所有输入产生正确舍入的结果，支持多达32位的多种表示，且只需单一实现。

研究背景与动机

核心问题

1. 正确舍入的挑战: 科学计算广泛使用数学库提供的基本函数，但为所有输入产生正确舍入结果极其困难（即"制表者困境"），主流数学库无法为所有输入产生正确结果。
2. 可移植性和可重现性问题: 缺乏正确舍入的数学库会导致应用程序在不同机器上产生完全不同的结果，影响可移植性和可重现性。
3. 多种表示格式的需求: 随着自定义格式（如bfloat16, tensorfloat32, FP8）的增加，需要一个能为多种表示和舍入模式提供正确结果的参考库。

现有方法的局限性

• Minimax多项式逼近: 传统方法生成最小化所有输入最大误差的多项式逼近，但当实值输出非常接近舍入边界时，自由度显著减小。
• 性能与正确性权衡: 现有库在性能（如Payne-Hanek实现）或正确性（如GCC的libm）方面做出权衡。

核心贡献

1. 高效范围缩减技术: 开发了结合浮点和整数运算的高效范围缩减算法，能够维护足够的π位数以产生正确结果。
2. 多表示单一实现: 实现了单一多项式逼近，能为10位到32位的多种表示和所有标准舍入模式产生正确舍入结果。
3. 性能优化: 整数基础的范围缩减相比浮点策略提升19%的性能，整体比主流库更快或性能相当。
4. 完整的三角函数库: 为sin、cos、tan函数提供了快速且正确的实现。

方法详解

RLIBM方法核心思想

RLIBM方法的关键洞察是直接逼近正确舍入结果，而非函数的实值。对于给定输入的正确舍入结果，存在一个实值区间，该区间内任何值都会舍入到正确结果。这提供了比minimax方法更大的自由度（对所有输入都是1 ULP）。

多表示支持机制

为支持多种表示，RLIBM项目提出生成(n+2)位表示的多项式逼近，使用round-to-odd舍入模式。这种方法的优势在于：

• round-to-odd结果保留了直接舍入到目标表示所需的所有信息
• 后续舍入到较低位宽表示能产生正确结果
• 避免了双重舍入错误

范围缩减算法

基本原理

三角函数的范围缩减将输入x∈-∞,∞映射到缩减输入x'∈-π/2^(t+1), π/2^(t+1)，其中：

x = x' + kπ/2^t
k = [2^t * x/π]
x' = π/2^t * r, 其中r = 2^t*x/π - k

浮点实现策略

小输入处理 (|x| < 2^30):

• 使用80位的256/π，分为两个double值存储
• 避免中间舍入误差
• 利用部分乘积精确计算k和分数部分r

大输入处理 (2^30 ≤ |x|):

• 版本1: 将256/π分为28位片段存储在double数组中，每片段使用截断模式生成
• 版本2: 使用53位精度片段，利用fused-multiply-add指令减少舍入误差

整数实现策略

小输入优化:

• 使用80位的256/π，分为两个40位整数P1和P0
• 通过位移操作识别整数k和分数位
• 避免浮点运算的精度损失

大输入处理:

• 使用192位的256/π，分为三个64位整数
• 计算128位部分乘积
• 通过位移操作提取相关位

输出补偿

利用三角恒等式进行输出补偿：

sin(x) = sin(k'π/2^t)cos(x') + cos(k'π/2^t)sin(x')
cos(x) = cos(k'π/2^t)cos(x') - sin(k'π/2^t)sin(x')

通过预计算表和周期性/对称性优化，将所需预计算值减少到512个。

实验设置

测试环境

• 硬件: 2.10GHz Intel Xeon(R) Silver 4310服务器，256GB RAM
• 操作系统: Ubuntu 24.04.1 LTS
• 测量工具: 性能计数器

对比库

• GLIBC: float和double libm
• Core-Math: 正确舍入库
• RLIBM实现: 多种范围缩减策略的变体

评价指标

• 正确性: 通过完全枚举验证所有输入的正确性
• 性能: 相对于其他库的加速比

实验结果

正确性验证

• RLIBM函数: 为10位到32位所有表示的所有输入产生正确舍入结果
• GLIBC float libm: 对32位float输入的sin、cos、tan有数千个错误结果
• GLIBC double libm: 比float版本更准确但仍有错误
• Core-Math: 仅对32位产生正确结果，对10-32位范围因双重舍入错误而失败

性能结果

范围缩减优化效果

混合方法（小输入用浮点，大输入用整数）相比其他策略：

• 比初始浮点方法(FP V1)快19%
• 比替代浮点方法(FP V2)有显著提升
• 比纯整数方法快4%

与其他库的比较

• 比Core-Math平均快10%
• 比GLIBC double函数平均快137%
• 性能提升主要归因于高效的范围缩减和整数运算的精度优势

技术创新点

1. 精度与性能的平衡

• 整数运算提供比64位double更高的精度（uint64_t和uint128_t）
• 减少了获得足够精度缩减输入所需的部分乘积数量

2. 混合范围缩减策略

• 小输入使用浮点运算（当256*x/π的整数部分足够小时）
• 大输入使用整数运算（提供更高精度和更简单的位操作）

3. 位操作优化

• 使用位移操作识别256*x/π中与缩减输入和k的低位相关的部分
• 避免了浮点运算中的舍入累积

相关工作

传统方法

• Minimax逼近: Remez算法等，但在舍入边界附近自由度有限
• Payne-Hanek算法: 经典范围缩减方法，但实现效率是挑战

正确舍入研究

• CR-LIBM: 早期正确舍入库，但性能较慢
• Core-Math: 现代正确舍入实现，但仅支持单一表示

RLIBM项目发展

• 从基本函数(e^x, log等)扩展到三角函数
• 多表示支持的创新方法

结论与讨论

主要结论

1. 可行性证明: 证明了为三角函数生成快速且正确的实现是可能的
2. 范围缩减关键性: 高效范围缩减与低次多项式逼近同等重要
3. 整数运算优势: 整数基础实现在大输入时显著优于浮点方法

局限性

1. 复杂性: 实现复杂度较高，需要精确的位操作和多种策略
2. 内存开销: 需要预计算表和多精度常数存储
3. 可扩展性: 扩展到更高精度表示需要重新设计

未来方向

1. GPU平台: 探索GPU平台的正确舍入库
2. 标准化: 参与IEEE-754标准委员会推动强制正确舍入
3. 主流集成: 与主流数学库开发者合作集成这些方法

深度评价

优点

1. 理论与实践结合: 将RLIBM理论成功应用到具有挑战性的三角函数
2. 全面的工程优化: 从算法到实现的全方位优化
3. 严格的验证: 通过完全枚举验证正确性
4. 实用价值: 解决了实际应用中的重要问题

不足

1. 实现复杂性: 多种策略的组合增加了实现和维护复杂性
2. 可读性: 大量位操作代码的可读性和可维护性有待提高
3. 理论分析: 缺乏对为什么整数方法更优的深入理论分析

影响力

1. 学术贡献: 为数值计算领域提供了新的正确舍入实现方法
2. 实用价值: 可直接应用于需要高精度数值计算的科学计算
3. 标准推动: 可能影响未来浮点标准的发展

适用场景

1. 科学计算: 需要高精度和可重现性的数值模拟
2. 金融计算: 要求精确结果的金融建模
3. 嵌入式系统: 需要支持多种浮点格式的系统
4. 参考实现: 作为其他库的正确性基准

参考文献

本文引用了数值分析、浮点运算和正确舍入领域的重要文献，包括：

• Muller的基本函数参考书
• MPFR高精度库
• Payne-Hanek范围缩减算法
• IEEE-754浮点标准相关研究

这篇论文在数值计算领域做出了重要贡献，成功地将理论方法转化为实用的高性能实现，为科学计算中的正确舍入问题提供了有效解决方案。