为什么将GNU科学库的文件trigc中的PI分为三部分

Question

在下面的代码中，为什么Pi被分成三个常数P1, P2和P3?有相关的数学理论吗?如果是为了提高r的计算精度，我在更高的精度上运行代码，但没有比Pi有任何提高。(代码来自gsl/specfunc/trigc: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=%.16enfmod  =%.16enmod2pi=%.16en",x,fmod(x,2*M_PI), mod2pi(x));
  return 0;
}

与使用Magma在线计算器的多精度结果的比较

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个前导位不会导致精度的损失。