factor Pi into two singles A*B

gistfile1.txt

(* Mathematica source 
same method as the previous gist for double. This validates that the constants giving here are optimal: 
https://stackoverflow.com/questions/26692859/best-machine-optimized-polynomial-minimax-approximation-to-arctangent-on-1-1
*) 

(* round(Pi,48,RN) *)
s = 221069929750889

(* prime factors of 's' and nearby integers *)
f0 = FactorInteger[s];
f1 = FactorInteger[s + 1];
f2 = FactorInteger[s + 2];
f3 = FactorInteger[s + 3];
f4 = FactorInteger[s - 1];
f5 = FactorInteger[s - 2];
f6 = FactorInteger[s - 3];

(* bits need for largest prime *)
{N[Log2[f0[[Length[f0], 1]]]],
 N[Log2[f1[[Length[f1], 1]]]],
 N[Log2[f2[[Length[f2], 1]]]],
 N[Log2[f3[[Length[f3], 1]]]],
 N[Log2[f4[[Length[f4], 1]]]],
 N[Log2[f5[[Length[f5], 1]]]],
 N[Log2[f6[[Length[f6], 1]]]]}

(* output: {19.6988, 29.234, 17.9695, 45.6515, 36.3432, 39.0113, 30.5098} *) 

(* only f0 & f2 are candidates (24 or less) : dump the factors *)
f0
f2

(* output:
{{5867, 1}, {44279, 1}, {850973, 1}}
{{13, 1}, {61, 1}, {73, 1}, {14879, 1}, {256661, 1}}
*)

(* 
  442779*5867 requires more than 24 bits...f0 is out
  
  256661*61 fits & 13*73*14879 both fit: so f2 (s+2) is good
  
  // pi-(pi_a*pi_b) ~= 3.19744*10^-14 (1.0001*2^-45)
  const float pi_a = 15656321.f/8388608.f;
  const float pi_b = 14120171.f/8388608.f;
  
  then: fmaf(pi_a, pi_b, x) is ~= pi+x;
  
*)