Why $\pi$ // Rationalize does not give a rational multiple of $\pi$ in my example?
I want to show that without hypotheses, the problem of determining whether $x=a/\pi$ came from a rational number is unsolvable. But, spoiler alert, the OP's number 2.47...
is probably not a rational multiple of Pi
or it is not precise enough to determine whether it is.
First, a floating-point number $x$ comes with an uncertainty $dx$. If $x$ comes from rounding a rational number $y$ to the nearest floating-point number, then we have $$x-dx \le y \le x+dx \,.$$ In Mathematica, $dx$ is given roughly by
dx = x * 10.^-Precision[x]
I say "roughly" because $dx$ will be a power of $2$ and the left and right uncertainties can be asymmetric at the boundaries where the floating-point exponent changes. For instance, if $y$ rounds to $1$ in the current binary64 standard for machine precision (double precision), then we have $1-2^{-54} \le y \le 1+2^{-53}$.
Second, the rational numbers are dense in the real numbers, which means there are infinitely many rational numbers $\tilde y$ such that $$x-dx \le \tilde y \le x+dx \,. \tag{1}$$
The problem is this: Given that we know $x$ (and therefore $dx$), find the rational number $y$ from among the infinitely many $\tilde y$ satisfying (1). It is impossible, since too much information was lost when $y$ was rounded.
A solvable problem is this: Find the rational number $\tilde y$ satisfying (1) with the smallest denominator. This problem is solved in Mathematica by any of these calls:
yy = Rationalize[x, dx]
yy = Rationalize[x, 0] (* uses the precision of x to determine dx *)
yy = Last@Convergents[x] (* uses the continued fraction expansion of x *)
Example.
Let's take a simpler version of the OP's second problem, without the Pi
, which is a red herring anyway.
brat = (12343556574747/50000); (* == b/Pi *)
x = N[brat];
dx = x * 10.^-Precision[x];
yy = Rationalize[x, dx]
Rationalize[x, 0]
Last@Convergents[x]
(* 1048955437722 / 4249 *)
There is no way to infer from x
that the original denominator was 50000
. It should be clear that yy
is not the closest fraction, just the one with the smallest denominator within the uncertainty dx
. In fact brat
is closer to x
than yy
. One might think there is still hope to recover brat
from x
, but there is not.
A continued fraction convergent of x
gives the rational number closest to x
that has a denominator smaller than the denominator of the next convergent. We will see from the convergents of x
that there is always a better choice of a fraction than brat
.
Let's look at the number x
. It is a fraction of the form $x_0/2^n$, which we can get using SetPrecision
:
xx = SetPrecision[x, Infinity]
(* 8283620594510023 / 33554432 <-- x0/2^25 *)
Here are some of the convergents related to our problem. Since we have converted x
to an infinite-precision fraction, we get more convergents:
Convergents[xx][[6 ;; 11]]
(*
{ 268348919935 / 1087,
1048955437722 / 4249, <-- smallest at mach. prec.
1317304357657 / 5336,
2366259795379 / 9585,
10782343539173 / 43676, <-- best with denom. < 140612
34713290412898 / 140613}
*)
So if you knew that the denominator was between, say, 10000
and 100000
, the best approximation is not brat
, but 10782343539173/43676
. What we hoped to recover turns out to be a less than optimal solution. I have no idea how to choose the correct solution from among all the less than optimal ones. I don't think it can be done.
Finally, we can come to a way to solve this problem, assuming $x$ is a good enough approximation of $y$. In general rational numbers are harder to approximate by rational numbers (other than themselves) than irrational numbers. We can use this to determine when a number is probably a rational number. (I am not an expert in this field: I do not know whether there is a way to assign a probability to the result.)
Rationalize[x]
solves the problem if the denominator of $y=p/q$ satisfies the criterion given in the docs:
The left-hand side of the inequality is $dx = |y-x|$. The exponent $2$ shows up in various Diophantine approximation theorems. The most elementary is
$$\left| {p' \over q'} - {p \over q} \right| \ge {1 \over q'q} \ge {1 \over \mathop{\text{min}}(q',q)^2}\,.$$
From this we can infer on condition x
should satisfy: Since the a small error in the denominator of x
implies the error in approximating y = p/q
is at least around 1/q^2
, we should have
Accuracy[x] >= 2 Log10[q]
Then there are Dirichlet's theorem and its generalizations: For an irrational $y$, there are infinitely many coprime integers $p'$ and $q'>0$ such that $$\left| {p' \over q'} - y \right| \le {c \over (q')^2} \,,$$ where the constant $c = 1$ is Dirchlet's theorem. See Diophantine approximation for more.
For a rational number, the error scaled by square of the denominator should be large, except once, when the error is zero. If $x \approx p/q$, then one might expect the scaled error $(q')^2\left| {p' / q'} - x \right|$ will be small only once, provided $x$ is a good enough approximation. (Since we're changing $y$ to $x$, this needs proof, which I haven't done, but it seems to hold computationally.) On the other hand, an irrational number will have several peaks in the scaled error.
Here's a look at some known irrational numbers:
Block[{$MaxExtraPrecision = 500},
convE = Convergents[E, 100];
scerrE = (E - convE) Denominator[convE]^2;
convPi = Convergents[Pi, 100];
scerrPi = (Pi - convPi) Denominator[convPi]^2;
conv2 = Convergents[10 Sqrt[2], 100];
scerr2 = (10 Sqrt[2] - conv2) Denominator[conv2]^2;
ListLinePlot[RealExponent@{scerrE, scerrPi, scerr2}
, PlotRange -> All, Frame -> True,
PlotLegends -> {E, Pi, 10 Sqrt[2]}]
]
Here's an approximation of a rational number, similar to the problem in the OP. There's a shallow "peak" in 5th convergent, similar in magnitude to the irrational numbers above. This is followed by deep spike, followed by a large jump. It is so large because it is scaled by q^2
, where q
is the denominator of the 11th convergent; therefore 1/q
will be less than the rounding error in x = N[y]
, and the jump in scaled error will be greater than Accuracy[x]
.
x = N[12346/33333 Pi, 30]/Pi;
conv = Convergents[SetPrecision[x, Infinity]];
scerr = (x - conv) Denominator[conv]^2;
ListLinePlot[RealExponent@scerr
, PlotRange -> All, Frame -> True]
There's a disturbing issue in the behavior of Convergents[x]
. It doesn't find the approximation we seek. Apparently we need to increase the precision of x
.
Convergents[x]
Convergents[SetPrecision[x, Infinity]]
(*
{0, 1/2, 1/3, 3/8, 10/27, 1023/2762, 1033/2789, 2056/5551, 5145/13891}
{0, 1/2, 1/3, 3/8, 10/27, 1023/2762, 1033/2789, 2056/5551, 5145/13891,
12346/33333,
192070822042862116379943873/518572550717213909387062135,
...10 fractions omitted...,
939034248967552886981567682091/2535301200456458802993406410752}
*)
The deep spike is at the desired fraction 12346/33333
. This should be typical when $x$ is a good approximation to $y$.
If we put the ideas above together we can write a function findRationalMultiples[x, m]
that will return a rational multiple of m
in the form {p * m / q}
, if x
seems to be a rational multiple of m
. It returns an empty list {}
if x
appears not to be a rational multiple of m
. It may also return a list of rational multiples, depending on the parameters to PeakDetect[]
. If x
is a rational multiple of m
, then the last element of the list is probably the correct answer; but such cases should be examined carefully.
The call findRationalMultiples[{x, m}, sm, sh, th]
passes the smoothing, sharpness, and threshold parameters to PeakDetect[]
.
The call findRationalMultiples[{x, m}, sm, sh, th, True]
returns an Association
containing auxiliary data helpful in analyzing rational approximations of x/m
.
ClearAll[findRationalMultiples];
findRationalMultiples[xm_, smoothing_ : 1, sharpness_ : 5,
t_ : Automatic, verbose_ : False] :=
Module[{x, m, conv, scalederror, data, threshold, res},
{x, m} = Replace[xm, {y_Real :> {y, 1}}];
x = x/m;
res = <||>;
res["Convergents"] = Convergents[x];
res["ScaledError"] = Abs@(x - res["Convergents"])*
Denominator[res["Convergents"]]^2;
data = -RealExponent@res["ScaledError"];
threshold =
t /. {Automatic -> -Max@data/3, Scaled[s_] :> -s*Max@data};
res["Peaks"] = PeakDetect[data, smoothing, sharpness, threshold];
res["Answers"] =
m*Extract[res["Convergents"], Position[res["Peaks"], 1]];
If[verbose,
res["SmoothData"] = GaussianFilter[data, {3, smoothing}];
res,
res["Answers"]
]
];
Example: A rational multiple of Pi
. You need enough digits
to make a good approximation, and we set the working precision with SetPrecision
a bit higher so that Convergents[]
doesn't quit too soon.
digits = 15. + Log10[33000^2];
x = SetPrecision[
Round[29/33000 Pi, 2^-Ceiling[Log2[10^digits]]],
100];
findRationalMultiples[{x, Pi}]
(* {(29 π)/33000} *)
Rationalize
fails on this example:
Rationalize[x/Pi]
(* 0.000878787878787878787878773416025526... *)
OP's main question. The function findRationalMultiples
fails to find a rational multiple of Pi
. We asked for the verbose association, so that we can examine the details.
x = 2.475188114422334324055850325497184848859183731753159368668810119536078273196602710715450734911249987;
frm = findRationalMultiples[{x, Pi}, 1, 5, Automatic, True];
frm@"Answers"
(* {} *)
A plot of the scaled error shows several peaks, not just one as we would expect for a rational multiple. Consequently, I think it is likely that the OP's x
is not a rational multiple of Pi
or x
is a bad approximation of the real number it represents.
ListLinePlot[RealExponent@frm["ScaledError"][[;; 88]], PlotRange -> All]
Numerical caveat: As with most, if not all, numerical routines, it is possible to construct pathological examples that fool the routine. For instance, y = Round[Pi, 1/2^b]
will have the same initial convergents as Pi
until the convergents differ from Pi
by around 1/2^b
or less. But y
is a rational number and findRationalMultiples[{N@y, 1}]
works up to b == 27
, returns multiple candidates up to b == 36
, and returns no solutions for b >= 37
.
Appendix.
Well, why not?, the post is already so long.... Here's a slightly different approach to findRationalMultiples
using Mathematica's notion of accuracy to perform the role of scalederror
. It also returns an association, with slightly different keys and through Reap[findRationalMultiples[{x, m}]]
; it uses Sow
with the tag findRationalMultiples
if you want to use Reap
with tags. Generally, because we increase the precision before calling Convergents
, we get at least one, and often several, convergents at the end of the Convergents[]
list whose error is zero at the scaled accuracy; we drop all but the one with the least denominator.
The Accuracy[]
for comparing the difference of each convergent p/q
from x
is set to -Log10[1/q^2]
plus the parameter scale
. The default value 4
corresponds roughly to the constant $c$ used in Rationalize[x]
(see docs or quote above); however, the behavior here is somewhat different than Rationalize
. If x
is found by rounding an unpathological rational number, then findRationalMultiples[{x, 1}, s]
should find a single good candidate at both smaller and larger scales. (The chances of finding more good candidates increases as the scale s
decreases, and vice versa.) At small enough scales, findRationalMultiples[{x, 1}, s]
will fail on any x
probably.
[TBD: Just realized that trim
needs some adjustment; sometimes the good candidate is adjacent to the bad candidates, but I have to run. Drats.]
ClearAll[findRationalMultiples];
findRationalMultiples::mult = "Multiple good candidates found at scale ``; returning first. Use Reap for more information.";
findRationalMultiples::noans = "No good candidates found at scale ``; returning first convergent with zero error. Use Reap for more information.";
findRationalMultiples[xm_, scale_ : 4] :=
Module[{x, m, conv, trim, zeros, res},
{x, m} = Replace[xm, {y_Real :> {y, 1}}];
x = x/m;
res = <||>;
res["Convergents"] = Convergents[SetPrecision[x, Infinity]];
res["Errors"] = MapThread[SetAccuracy,
{SetPrecision[x, N@Precision[x]] - res["Convergents"],
scale + 2 Max[5, -Log10[x], #] & /@
Log10@N@Denominator@res["Convergents"]}];
trim = LengthWhile[Reverse@res["Errors"], # == 0 &];
(* compare first two zeros to be dropped
IF the accuracy jump > Precision[x] or Log10[denominator] >= Precision[x]
THEN first one is a good candidate, so decrement trim *)
If[trim > 1 &&
Accuracy[res["Errors"][[-trim]]] < Accuracy[res["Errors"][[1 - trim]]] - Precision[x],
trim--
];
zeros = Position[Drop[res["Errors"], 1 - trim], z_ /; z == 0];
res["Candidates"] = m*Extract[res["Convergents"], zeros];
res["Answer"] = Replace[res["Candidates"], {a_, ___} :> a];
Sow[res, findRationalMultiples];
Switch[Length@res["Candidates"],
0 | 1, Message[findRationalMultiples::noans, scale],
2, Null,
_, Message[findRationalMultiples::mult, scale]
];
res["Answer"]
];
My example from above has stable results at different scales:
digits = 15. + Log10[33000^2];
x = SetPrecision[
Round[29/33000 Pi, 2^-Ceiling[Log2[10^digits]]],
100];
findRationalMultiples[{x, Pi}]
findRationalMultiples[{x, Pi}, 2]
findRationalMultiples[{x, Pi}, 8]
(*
(29 π)/33000
(29 π)/33000
(29 π)/33000
*)
Rationalization of the OP's number is unstable with change of scale. The first trial below reproduces one of the other posts' answers, but the instability suggests it's doubtful that x
came from rounding.
x = 2.475188114422334324055850325497184848859183731753159368668810119536078273196602710715450734911249987;
findRationalMultiples[{x, Pi}]
(* (528518980406819854 π)/670814204566564617 *)
findRationalMultiples[{x, Pi}, 2]
findRationalMultiples::mult : Multiple good candidates found at scale 2; returning first. Use Reap for more information. >>
(* (105618056127 π)/134054016857 *)
findRationalMultiples[{x, Pi}, 8]
findRationalMultiples::noans : No good candidates found at scale 8; returning first convergent with zero error. Use Reap for more information.
(* (127783177086971763037570942964097163555065759781296 π) /
162186739686443508647264007143241471430512176744729 *)
b = (12343556574747/50000) \[Pi];
N[b, 20]/Pi // Rationalize[#, 0] &
(* 12343556574747/50000 *)
Rationalize[(2.\
4751881144223343240558503254971848488591837317531593686688101195360782\
73196602710715450734911249987/Pi), #] & /@ (10^Range[-10, -5])
(* {86111/109295, 24592/31213, 12257/15557, 11711/14864, 8123/10310, \
26/33} *)
Expansion into a continuous fraction can be used.
In your example
a = 2.4751881144223343240558503254971848488591837317531593686688101195\
36078273196602710715450734911249987;
ContinuedFraction[a/\[Pi]]
{0,1,3,1,2,2,472,1,2,3,1,1,1,1,1,2,3,3,2,4,2,2,2,4,5,243,1,2,6,2,3,1,1,2,1,16,1,24388,1,4,2,2,1,1,1,1,20,2,1,1,6,6,1,3,1,8,1,1,2,5,4,13,1,1,6,3,2,300,6,103,5,2,3,3,2,1,3,1,3,16,2,1,4,1,1,3,34,1}
There is a large element 24388. Dropping all elements starting with it gives
b = FromContinuedFraction@{0, 1, 3, 1, 2, 2, 472, 1, 2, 3, 1, 1, 1, 1,
1, 2, 3, 3, 2, 4, 2, 2, 2, 4, 5, 243, 1, 2, 6, 2, 3, 1, 1, 2, 1,
17}
$$ \frac{528518980406819854}{670814204566564617} $$ with error
a/\[Pi] - b///N
$$\text{9.11146$\grave{ }$*${}^{\wedge}$-41}$$
If all digits in your number are correct it is not a rational number with a small denominatior. If only 40 digits or so are true then b
looks like the best candidate. But if only a few digits are true one can take segments up to the previous relatively large numbers, say
FromContinuedFraction@{0, 1, 3, 1, 2, 2}
$$\frac{26}{33}$$