为什么GNU Scientific Library trig.c 文件中的PI 被分成三部分？

Question

在下面的代码中，为什么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;

Answer 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 拆分位，如在下一个 table 中所示，其中第一行是位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 个前导位不会导致精度损失。

为什么GNU Scientific Library trig.c 文件中的PI 被分成三部分？

Why the PI in the file trig.c of the GNU Scientific Library be divided in three parts?

math

floating-point

pi

gnu

gsl

为什么这些常量