Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save Marc-B-Reynolds/948b550fbcba3726da8bdb680b50e440 to your computer and use it in GitHub Desktop.
Save Marc-B-Reynolds/948b550fbcba3726da8bdb680b50e440 to your computer and use it in GitHub Desktop.
factor Pi into two singles A*B
(* 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;
*)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment