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

In code below, why is Pi divided in three constants P1, P2 and P3? Is there some related math theory? If it is for improve the compute accuracy for r, I run the code on higher precision but without any improve than just Pi.(The code from 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;

Solution

Test program in 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;
}

compared to multiprecision result using the Magma online calculator

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

with results

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

and

magma = 6.269073200367326567623794342882040802035079748091348034188201251009459335653510999632076033999854435

showing that indeed the higher effort leads to more precise results.

Why those constants

For some reason the developers decided to split the bits not directly of pi/4 but based on 10*pi/4=5/2*pi as can be seen in the next table where the top row are the bits of a long version of 5/2*pi while the next three are binary representations of the constants multiplied by 10.

111 11011010100111101000101001010101010011100001011110010110000011111010111110

111.1101101010011110100001
  0.00000000000000000000011001010101010011100001
  0.000000000000000000000000000000000000000000000111100101100000

A split based on pi/4 using 25 bits in each part is

0.1100100100001111110110101010001000100001011010001100001000110100110001001100

0.1100100100001111110110101
0.00000000000000000000000000100010001000010110100011
0.000000000000000000000000000000000000000000000000000000100011010011000100110

and would lead to constants

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

The idea being that integer multiples up to 2^27 of P1,P2,P3 are exact so that the successive reductions remove the leading identical bits without losing precision. Essentially, the input argument with 53 bits mantissa gets (virtually) extended to a 75 bit mantissa by filling with zeros, and then this number gets reduces exactly by multiples of 2*pi. Cancellation of up to 22 leading bits will not result in a loss of precision.