为什么 GNU 科学库的文件 trig.c 中的 PI 分为三部分？答案

【问题标题】：Why the PI in the file trig.c of the GNU Scientific Library be divided in three parts?为什么 GNU 科学库的文件 trig.c 中的 PI 分为三部分？
【发布时间】：2016-09-24 10:37:27
【问题描述】：

在下面的代码中，为什么 Pi 被分成三个常数 P1、P2 和 P3？有没有相关的数学理论？如果是为了提高 r 的计算精度，我会以更高的精度运行代码，但除了 Pi 没有任何改进。（来自 gsl/specfunc/trig.c:576 的代码）

  const double P1 = 4 * 7.85398125648498535156e-01;
  const double P2 = 4 * 3.77489470793079817668e-08;
  const double P3 = 4 * 2.69515142907905952645e-15;
  const double TwoPi = 2*(P1 + P2 + P3);

  const double y = 2*floor(theta/TwoPi);

  double r = ((theta - y*P1) - y*P2) - y*P3;

【问题讨论】：

它应该比fmod(theta, M_2PI)更准确地减少theta到1e+8的幅度。

标签： math floating-point gnu gsl pi

【解决方案1】：

C 语言测试程序

#include<math.h>
#include<stdio.h>


double mod2pi(double theta) {
  const double P1 = 4 * 7.85398125648498535156e-01;
  const double P2 = 4 * 3.77489470793079817668e-08;
  const double P3 = 4 * 2.69515142907905952645e-15;
  const double TwoPi = 2*(P1 + P2 + P3);

  const double y = 2*floor(theta/TwoPi);

  return ((theta - y*P1) - y*P2) - y*P3;
}

int main() {
  double x = 1.234e+7;

  printf("x=%.16e\nfmod  =%.16e\nmod2pi=%.16e\n",x,fmod(x,2*M_PI), mod2pi(x));

  return 0;
}

与使用Magma online calculator 的多精度结果相比

RR := RealField(100);
pi := Pi(RR);
x := 1.234e+7;
n := 2*Floor(x/(2*pi));
"magma =",RR!x-n*pi;

有结果

x=1.2340000000000000e+07
fmod  =6.2690732008483607e+00
mod2pi=6.2690732003673268e+00

和

magma = 6.269073200367326567623794342882040802035079748091348034188201251009459335653510999632076033999854435

表明确实更高的努力会导致更精确的结果。

为什么使用这些常量

出于某种原因，开发人员决定不直接拆分pi/4 的位，而是基于10*pi/4=5/2*pi，如下表所示，其中第一行是5/2*pi 的长版本的位，而接下来的三个是常量乘以10 的二进制表示。

111 11011010100111101000101001010101010011100001011110010110000011111010111110

111.1101101010011110100001
  0.00000000000000000000011001010101010011100001
  0.000000000000000000000000000000000000000000000111100101100000

基于pi/4 在每个部分中使用 25 位的拆分是

0.1100100100001111110110101010001000100001011010001100001000110100110001001100

0.1100100100001111110110101
0.00000000000000000000000000100010001000010110100011
0.000000000000000000000000000000000000000000000000000000100011010011000100110

并且会导致常量

const double P1 = 4 * 7.85398155450820922852e-01;
const double P2 = 4 * 7.94662735614792836714e-09;
const double P3 = 4 * 3.06161646971842959369e-17;

这个想法是P1,P2,P3 的整数倍数达到2^27 是精确的，因此连续减少会删除前导相同位而不会丢失精度。本质上，带有 53 位尾数的输入参数（实际上）通过填充零扩展为 75 位尾数，然后这个数字精确地减少了 2*pi 的倍数。取消最多 22 个前导位不会导致精度损失。

【讨论】：

这种三角函数的参数缩减风格被称为 Cody-Waite 缩减，自 1970 年代后期以来一直在使用。参见：William J. Cody 和 William M. Waite。《基本功能软件手册》。普伦蒂斯霍尔，1980 年。
最近的后续工作：[1] Jean-Michel Muller 和 Peter Kornerup。 “扩展 Cody 和 Waite 范围缩小方法的范围”。 2005.draft. [2] 西尔维·博尔多、马克·道马斯和李仁仓。 “使用融合乘法加法正式验证了参数减少。” IEEE Transactions on Computers 58.8 (2009): 1139-1145。 preprint
我傻了，我没有意识到它和fmod一样的功能，我现在明白了代码，感谢LutzL的回答。也感谢 njuffa 的推荐，我正在阅读这篇文章。这是我在 stackoverflow 上的第一个问题，我非常感谢这个平台帮助我解决了我的问题。
我能找到的最早的印刷参考文献是：J. Y. Wang。 “具有大输入参数的周期函数的评估。” ACM SIGNUM 通讯 13.4 (1978): 7-8.
@Star 当您说pi 时，大概您指的是最接近数学π 的double 数字？如果是这样，您看到的结果是这台机器 PI 与数学 π 之间差异的结果。这里使用 Cody-Waite 缩减的要点是通过数学 π 有效缩减，在一定范围内输入theta。