math - 为什么GNU科学库的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 进行拆分，如下所示如下表所示，其中顶行是长版本 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 个前导位不会导致精度损失。