Re: Pi day
Patrick Scheibe
imagine you would take two good rational approximations for Pi:
Rationalize[Pi, 10^-#] & /@ {11, 22}
gives
{833719/265381, 1299139324288/413528890451}
what you see is that the first approximation uses 6 digits and the last
one 13 at most.
Taking the numerators as x-values and the denominators as y-values and
you can think of this as a line between the points {833719,265381} and
{1299139324288,413528890451}.
All points "near" this line seem to be good approximations to Pi.
Furthermore, the numbers in the last points contain to much digits to
get a rational expression where every digit is different.
Now you could say, I'm going from the last point downwards on the line
and check the points near the line. If a point fulfills the constraint
that every digit in x and y value is unique, you probably have a good
solution.
I don't know what you have so far but I got (zeros included, since I
don't want to solve existing exercises):
8954632071/2850347916
which is 7.299378879110918*^-10 away from Pi.
Cheers
Patrick
On Sat, 2010-03-13 at 07:55 -0500, Tom wrote:
> Hello, I am a high school math teacher and the following puzzle was
> posed by a few math teachers I am in contact with.
>
> Create a fraction whose numerator has the digits 1 - 9 (used once)
> and whose denominator has the digits 1 - 9 (used one) .
>
> Which fraction has a value closest to the value of pi?
>
> I've worked on some "brute force" checks and managed to check all
> possible fractions with 2,3,4,5 and 6 digits. But after that, there
> are just too many possibilities.
>
> I don't have the programming ability to implement something elegant in
> Mathematica.
>
> Is there anyone who could suggest an approach to find the solution to
> this?
>
> Sincerely,
>
> Tom
>
DrMajorBob
(correct me if I'm wrong, someone!), but it fails in a peculiar way.
Here's an example with 3-digit numerator and denominator:
n = 3;
ge0 = {0, 1};
eq1 = {1, 0};
numerator = Sum[j v[i, j] 10^(i - 1), {i, n}, {j, n}];
denominator = Sum[j v[i + n, j] 10^(i - 1), {i, n}, {j, n}];
ratio = numerator/denominator;
variables = Flatten@Array[v, {2 n, n}];
c = Coefficient[Pi denominator - numerator, #] & /@ variables;
bounds = ConstantArray[{0, 1}, 2 n^2];
m1 = SparseArray[{row_, col_} /; row == Quotient[col, n, 1] + 1 ->
1, {2 n, 2 n^2}];
m2 = SparseArray[{row_, col_} /;
Mod[row - col, n] == 0 &&
Quotient[col, n^2, 1] == Quotient[row, n, 1] -> 1, {2 n, 2 n^2}];
m = ArrayFlatten[{{m1}, {m2}}];
b = ConstantArray[{1, 0}, Length@m];
"add a column (the new objective) and two inequalities, to minimize \
absolute error";
variables = Flatten@{objective, variables};
MatrixForm[m = ArrayFlatten@{{1, {c, -c}}, {0, m}}];
c = Flatten@{1, 0 c};
PrependTo[bounds, {-10^10, 10^10}];
b = Join[{{0, 1}, {0, 1}}, b];
Dimensions /@ {c, m, b, bounds}
{{19}, {14, 19}, {14, 2}, {19, 2}}
LinearProgramming[c, m, b, bounds]
(fail... with the errors N::meprec and LinearProgramming:lpsiter1... too
much extra precision and too many iterations needed)
In the first error message, Mathematica says it's reaching
$MaxExtraPrecision = 50 while evaluating
1/25 (1 - 100/37 (-(1/25) + 107/(100 \[Pi]))) +
41/37 (-(1/25) + 107/(100 \[Pi])) - 107/(100 \[Pi])
But that's exactly zero:
1/25 (1 - 100/37 (-(1/25) + 107/(100 \[Pi]))) +
41/37 (-(1/25) + 107/(100 \[Pi])) - 107/(100 \[Pi]) // Simplify
0
The code before "add a column..." produces an LP that works, either in the
form
LinearProgramming[c, m, b, bounds]
or
LinearProgramming[-c, m, b, bounds]
Those are the WORST approximations of the desired kind, above and below
Pi... but LP worked perfectly and gave Integer results.
The later code seems to produce the right matrices.
But LP runs into trouble with computing the Real value of an exact zero.
What's going on?
Bobby
On Sat, 13 Mar 2010 06:55:23 -0600, Tom <tideta...@gmail.com> wrote:
> Hello, I am a high school math teacher and the following puzzle was
> posed by a few math teachers I am in contact with.
>
> Create a fraction whose numerator has the digits 1 - 9 (used once)
> and whose denominator has the digits 1 - 9 (used one) .
>
> Which fraction has a value closest to the value of pi?
>
> I've worked on some "brute force" checks and managed to check all
> possible fractions with 2,3,4,5 and 6 digits. But after that, there
> are just too many possibilities.
>
> I don't have the programming ability to implement something elegant in
> Mathematica.
>
> Is there anyone who could suggest an approach to find the solution to
> this?
>
> Sincerely,
>
> Tom
>
da...@wolfram.com
> posed by a few math teachers I am in contact with.
>
> Create a fraction whose numerator has the digits 1 - 9 (used once)
> and whose denominator has the digits 1 - 9 (used one) .
>
> Which fraction has a value closest to the value of pi?
>
> I've worked on some "brute force" checks and managed to check all
> possible fractions with 2,3,4,5 and 6 digits. But after that, there
> are just too many possibilities.
>
> I don't have the programming ability to implement something elegant in
> Mathematica.
>
> Is there anyone who could suggest an approach to find the solution to
> this?
>
> Sincerely,
>
> Tom
Here is a method that is reasonably effective. It is a "greedy" approach
and I cannot offhand guarantee that is gives the best result (though I
believe it does).
The idea is to begin with all one digit (numerator,denominator) pairs.
Filter out the ones that cannot possibly be improved via further digits in
a way that will beat the best pairs. For the keepers, augment numerators
and denominators with all available single digits. Iterate this process
until we are out of digits. It is the pruning that keeps this from
becoming intractable as we progress over the digit pool. (Eventual
exhaustion of said pool also helps.)
In brief, the routines do as follows.
frac[] makes an integer pair into a rational.
closest[] takes a list of pairs and a target value, and finds the one
whose corresponding fraction is closest to that target.
intervalize[] takes a pair and creates the interval comprised of fractions
where first denominator, and then numerator, are increased by one. This
interval is a crude lower/upper bound on all values that can be attained
by adding allowable digits to the given numerator and denominator.
siftPairs[] takes pairs and a target. it first finds the pair that best
approximates the target. It then removes all others that cannot possibly
get within range via allowable augmentation to numerator and denominator.
successors[] takes a numerator/denominator pair and finds the set of allowed
pairs. These are defined as pairs with numerator and denominator one digit
longer, and neither one repeating digits.
With that intro, here is the code.
frac[{n_Integer, d_Integer}] := n/d
closest[pairs_, target_] :=
pairs[[Ordering[Abs[Map[frac, pairs, 1] - target], 1]]][[1]]
intervalize[{num_Integer, den_Integer}] :=
Interval[{num/(den + 1), (num + 1)/den}]
siftPairs[pairs_, target_] := Module[
{interval, num, den},
{num, den} = closest[pairs, target];
interval = intervalize[{num, den}];
Select[pairs, (IntervalIntersection[intervalize[#], interval] =!=
Interval[]) &]
]
successors[{num_Integer, den_Integer}] := Module[
{nd, dd, cnd, cdd, nums, dens},
{nd, dd} = IntegerDigits[{num, den}];
cnd = Complement[Range[9], nd];
cdd = Complement[Range[9], dd];
nums = Map[10*num + # &, cnd];
dens = Map[10*den + # &, cdd];
Flatten[Outer[List, nums, dens], 1]
]
We now form our initial pairs of one digit numerators and denominators. We
then nest our successorship function eight times (thus exhausting
available digits). In each step we sift to remove pairs that are not
contenders.
initpairs = Flatten[Outer[List, Range[9], Range[9]], 1];
In[377]:= Timing[
candidates =
Nest[Flatten[Map[successors, siftPairs[#, Pi]], 1] &, initpairs,
8];]
Out[377]= {16.864, Null}
In[378]:= Length[candidates]
Out[378]= 2576
In[380]:= best = closest[candidates, Pi]
Out[380]= {429751836, 136794258}
In[382]:= bestfrac = frac[best]
Out[382]= 23875102/7599681
In[386]:= N[bestfrac - Pi, 5]
Out[386]= 1.0186*10^-10
So we get within around 10^(-10) of our quarry.
It is instructive to print lengths of our pair lists at intermediate
steps. It goes over 70,000 before receding to the eventual 2600 or so.
Possibly a more aggressive sifter would so better (but would be more
trouble to code, or for that matter to figure out how to construct).
Given that we work from most significant digits downward, it seems
plausible that our eventual result will be optimal. But it is not obvious
to me that it is forced to be so. I need to give this more thought.
Daniel Lichtblau
Wolfram Research
Ray Koopman
I get a little better fit by walking the list.
s = FromDigits/@Permutations@Range@9;
i = j = 1; While[s[[i]]/s[[j]] < Pi, i++];
hi = s[[ i ]]/s[[j]]; ijhi = { i ,j};
lo = s[[i-1]]/s[[j]]; ijlo = {i-1,j};
k = Length@s; While[s[[-1]]/s[[k]] < Pi, k--];
While[j < k,
j++; While[s[[i]]/s[[j]] < Pi, i++];
If[s[[ i ]]/s[[j]] < hi, hi = s[[ i ]]/s[[j]]; ijhi = { i ,j}];
If[s[[i-1]]/s[[j]] > lo, lo = s[[i-1]]/s[[j]]; ijlo = {i-1,j}]];
{N[hi-Pi], s[[ijhi]]}
{N[lo-Pi], s[[ijlo]]}
{1.0185541299279066*^-10, {429751836, 136794258}}
{-8.499689840846258*^-11, {467895213, 148935672}}
Mark McClure
> In[380]:= best = closest[candidates, Pi]
> Out[380]= {429751836, 136794258}
>
> In[382]:= bestfrac = frac[best]
> Out[382]= 23875102/7599681
>
> In[386]:= N[bestfrac - Pi, 5]
> Out[386]= 1.0186*10^-10
>
> So we get within around 10^(-10) of our quarry.
I've got that beat by a bit:
N[Pi - 467895213/148935672]
8.49969*10^-11
In fact, this is provably the best possible. My code is only slightly
more clever than brute force. Of course, there are 9! = 362880 such
integers. Sort them from highest to lowest, say
{a(1),a(2),...a(362880)}. Then, on each sublist of the form
{a(k),a(k+1),...a(362880)}, find the fractions of the form a(k)/a(n)
and a(k)/a(n+1) that are just smaller and just larger than pi. (This
can be done quickly using a binary search.) Now, you've only got
around 600,000 fractions to check. Here's the code that runs in a few
minutes.
findPiInterval[n_Integer, rest_List] := Module[{a, b, c, len},
Which[
n/First[rest] > Pi, {1},
n/Last[rest] < Pi, {Length[rest]},
True,
list = rest;
a = 1; b = len = Length[list];
While[b - a > 1,
c = Floor[(a + b)/2];
If[n/list[[c]] < Pi, a = c, b = c]];
{a, b}]];
findPiInterval[all_List] := findPiInterval[First[all], Rest[all]] + 1;
ints1 = Reverse[FromDigits /@ Permutations[Range[9]]];
len = Length[ints1];
ints = Prepend[ints1, 0];
Monitor[close = Table[{k, findPiInterval[
ints = Rest[ints]]}, {k, 1, len - 1}], {k, len}]; // Timing
{235.83, Null}
close = Flatten[close /. {n_Integer, pos_List} :>
({ints1[[n]], #} & /@ ints1[[pos + n - 1]]), 1];
close = Append[#, #[[1]]/#[[2]]] & /@ close;
lo = Last[SortBy[Select[close, Last[#] < Pi &], Last]];
hi = First[SortBy[Select[close, Last[#] > Pi &], Last]];
{lo, hi}
{{467895213, 148935672, 51988357/16548408},
{429751836, 136794258, 23875102/7599681}}
N[Pi - Last[lo]]
8.49969*10^-11
Mark McClure
DrMajorBob
computer:
Timing[s = (Permutations@Range@9).(10^Range[8, 0, -1]);
j = 1; i = Position[s, x_ /; x > Pi s[[j]], 1, 1][[1, 1]];
hi = s[[i]]/s[[j]]; ijhi = {i, j};
lo = s[[i - 1]]/s[[j]]; ijlo = {i - 1, j};
k = -1 + Position[s, x_ /; s[[-1]] < Pi x, 1, 1][[1, 1]];
While[j < k, j++; While[s[[i]]/s[[j]] < Pi, i++];
If[s[[i]]/s[[j]] < hi, hi = s[[i]]/s[[j]]; ijhi = {i, j}];
If[s[[i - 1]]/s[[j]] > lo, lo = s[[i - 1]]/s[[j]];
ijlo = {i - 1, j}]];
{{N[hi - Pi], s[[ijhi]]}, {N[lo - Pi], s[[ijlo]]}}]
{4.90821, {{1.01855*10^-10, {429751836,
136794258}}, {-8.49969*10^-11, {467895213, 148935672}}}}
Position replaced the first two While loops, and Dot replaced Map in the
first line.
Using Position in the main loop is too slow, because it doesn't start at
the current position in the list. Using Drop[s, i - 1] and adding i - 1
would do that, but it's also too slow, presumably due to the expense of
allocating the Drop result.
Here's the same thing in base 11:
Timing[s = (Permutations@Range@10).(11^Range[9, 0, -1]);
j = 1; i = Position[s, x_ /; x > Pi s[[j]], 1, 1][[1, 1]];
hi = s[[i]]/s[[j]]; ijhi = {i, j};
lo = s[[i - 1]]/s[[j]]; ijlo = {i - 1, j};
k = -1 + Position[s, x_ /; s[[-1]] < Pi x, 1, 1][[1, 1]];
While[j < k, j++; While[s[[i]]/s[[j]] < Pi, i++];
If[s[[i]]/s[[j]] < hi, hi = s[[i]]/s[[j]]; ijhi = {i, j}];
If[s[[i - 1]]/s[[j]] > lo, lo = s[[i - 1]]/s[[j]];
ijlo = {i - 1, j}]];
{{N[hi - Pi], BaseForm[s[[ijhi]], 11]}, {N[lo - Pi],
BaseForm[s[[ijlo]], 11]}}]
{55.9061,{{4.2637*10^-12,{Subscript[a139786245, 11],Subscript[32498a5761,
11]}},{-8.34888*10^-13,{Subscript[a675231948, 11],Subscript[3415a68729,
11]}}}}
Base 12 seemed to steal every byte of memory on my 4GB iMac with 64-bit
Mathematic.
Permutations choked on base 13, immediately.
Bobby
>> In[380]:= best = closest[candidates, Pi]
>> Out[380]= {429751836, 136794258}
>>
>> In[382]:= bestfrac = frac[best]
>> Out[382]= 23875102/7599681
>>
>> In[386]:= N[bestfrac - Pi, 5]
>> Out[386]= 1.0186*10^-10
>>
>> So we get within around 10^(-10) of our quarry.
>>
>> It is instructive to print lengths of our pair lists at intermediate
>> steps. It goes over 70,000 before receding to the eventual 2600 or so.
>> Possibly a more aggressive sifter would so better (but would be more
>> trouble to code, or for that matter to figure out how to construct).
>>
>> Given that we work from most significant digits downward, it seems
>> plausible that our eventual result will be optimal. But it is not
>> obvious
>> to me that it is forced to be so. I need to give this more thought.
>>
>> Daniel Lichtblau
>> Wolfram Research
>
> I get a little better fit by walking the list.
>
> s = FromDigits/@Permutations@Range@9;
> i = j = 1; While[s[[i]]/s[[j]] < Pi, i++];
> hi = s[[ i ]]/s[[j]]; ijhi = { i ,j};
> lo = s[[i-1]]/s[[j]]; ijlo = {i-1,j};
> k = Length@s; While[s[[-1]]/s[[k]] < Pi, k--];
> While[j < k,
> j++; While[s[[i]]/s[[j]] < Pi, i++];
> If[s[[ i ]]/s[[j]] < hi, hi = s[[ i ]]/s[[j]]; ijhi = { i ,j}];
> If[s[[i-1]]/s[[j]] > lo, lo = s[[i-1]]/s[[j]]; ijlo = {i-1,j}]];
> {N[hi-Pi], s[[ijhi]]}
> {N[lo-Pi], s[[ijlo]]}
>
> {1.0185541299279066*^-10, {429751836, 136794258}}
> {-8.499689840846258*^-11, {467895213, 148935672}}
>
brien colwell
P/Q, then the smallest step forward (seems to be) either increasing P by
the smallest amount, or decreasing Q by the smallest amount and then
decreasing P to ensure strictly ascending.
The iteration would be start with the smallest number (123456789 /
987654321) and take the smallest steps forward until it crosses Pi.
I can't think of a clever way to iterate the numbers of interest in
ascending order, so I would probably just throw all 9! (400k) in a list and
sort it.
brien colwell
nums = Sort[FromDigits /@ Permutations[{1, 2, 3, 4, 5, 6, 7, 8, 9}]];
NIndex = 1;
DIndex = -1;
pv = 0;
v = 0;
ChooseNextIndexes[] := Module[{i = 1},
While[v < (nums[[NIndex - i]]/nums[[DIndex - 1]]),
i += 1];
If[nums[[NIndex + 1]]/nums[[DIndex]] <
nums[[NIndex - (i - 1)]]/nums[[DIndex - 1]],
{NIndex + 1, DIndex},
{NIndex - (i - 1), DIndex - 1}
]
];
(* This is a guess to speed up the search
*)
While[nums[[NIndex + 1]]/nums[[DIndex - 1]] < \[Pi],
NIndex += 1;
DIndex -= 1;
];
pv = 0;
pIndexes = {NIndex, DIndex};
While[(v = nums[[NIndex]]/nums[[DIndex]]) < \[Pi],
pv = v;
pIndexes = {NIndex, DIndex};
{NIndex, DIndex} = ChooseNextIndexes[];
];
If[N[\[Pi] - pv] < N[v - \[Pi]],
ToString[nums[[pIndexes[[1]]]]] <> " / " <>
ToString[nums[[pIndexes[[2]]]]],
ToString[nums[[NIndex]]] <> " / " <> ToString[nums[[DIndex]]]
]
N[If[N[\[Pi] - pv] < N[v - \[Pi]],
pv,
v
]]
Ray Koopman
concerned about the algorithm. 'hi' is meant to be the smallest value
of the form s[[i]]/s[[j]] that is greater than Pi, and 'lo' is meant
to be the largest value of the form s[[i]]/s[[j]] that is less than
Pi. 'hi' is OK, but I wonder about 'lo'.
Specifically, potential values of 'lo' with j > k are never checked,
where k is the largest value of j for which s[[-1]]/s[[j]] > Pi. When
the While[j < k, ...] loop exits, we have j = k, s[[i]]/s[[j]] > Pi,
and s[[i-1]]/s[[j]] < Pi. s[[-1]]/s[[k+1]] < Pi, but might it not
also be > the current 'lo'?
When I posted the code I thought that no further checking would be
necessary, but I am no longer convinced. However, I have not been
able to prove that checking is necessary, either. Mark McClure's
results indicate that checking is not necessary in this particular
case, but what about the general case, for targets other than Pi
or other 's'-lists?
----- DrMajorBob <btr...@austin.rr.com> wrote:
> A couple of tweaks improved the timing from 6.4 seconds to 4.9, on my
> computer:
>
> Timing[s = (Permutations@Range@9).(10^Range[8, 0, -1]);
> j = 1; i = Position[s, x_ /; x > Pi s[[j]], 1, 1][[1, 1]];
> hi = s[[i]]/s[[j]]; ijhi = {i, j};
> lo = s[[i - 1]]/s[[j]]; ijlo = {i - 1, j};
> k = -1 + Position[s, x_ /; s[[-1]] < Pi x, 1, 1][[1, 1]];
> While[j < k, j++; While[s[[i]]/s[[j]] < Pi, i++];
> If[s[[i]]/s[[j]] < hi, hi = s[[i]]/s[[j]]; ijhi = {i, j}];
> If[s[[i - 1]]/s[[j]] > lo, lo = s[[i - 1]]/s[[j]];
> ijlo = {i - 1, j}]];
> {{N[hi - Pi], s[[ijhi]]}, {N[lo - Pi], s[[ijlo]]}}]
>
> {4.90821, {{1.01855*10^-10, {429751836,
> 136794258}}, {-8.49969*10^-11, {467895213, 148935672}}}}
>
> Position replaced the first two While loops, and Dot replaced Map in the
> first line.
>
> Using Position in the main loop is too slow, because it doesn't start at
> the current position in the list. Using Drop[s, i - 1] and adding i - 1
> would do that, but it's also too slow, presumably due to the expense of
> allocating the Drop result.
>
> Here's the same thing in base 11:
>
> Timing[s = (Permutations@Range@10).(11^Range[9, 0, -1]);
> j = 1; i = Position[s, x_ /; x > Pi s[[j]], 1, 1][[1, 1]];
> hi = s[[i]]/s[[j]]; ijhi = {i, j};
> lo = s[[i - 1]]/s[[j]]; ijlo = {i - 1, j};
> k = -1 + Position[s, x_ /; s[[-1]] < Pi x, 1, 1][[1, 1]];
> While[j < k, j++; While[s[[i]]/s[[j]] < Pi, i++];
> If[s[[i]]/s[[j]] < hi, hi = s[[i]]/s[[j]]; ijhi = {i, j}];
> If[s[[i - 1]]/s[[j]] > lo, lo = s[[i - 1]]/s[[j]];
> ijlo = {i - 1, j}]];
> {{N[hi - Pi], BaseForm[s[[ijhi]], 11]}, {N[lo - Pi],
> BaseForm[s[[ijlo]], 11]}}]
>
> {55.9061,{{4.2637*10^-12,{Subscript[a139786245, 11],Subscript[32498a5761,
> 11]}},{-8.34888*10^-13,{Subscript[a675231948, 11],Subscript[3415a68729,
> 11]}}}}
>
> Base 12 seemed to steal every byte of memory on my 4GB iMac with 64-bit
> Mathematic.
>
> Permutations choked on base 13, immediately.
>
> Bobby
>
> --
> DrMaj...@yahoo.com
Daniel Lichtblau
>> Out[380]= {429751836, 136794258}
>>
>> In[382]:= bestfrac = frac[best]
>> Out[382]= 23875102/7599681
>>
>> In[386]:= N[bestfrac - Pi, 5]
>> Out[386]= 1.0186*10^-10
>>
>> So we get within around 10^(-10) of our quarry.
>
>
> N[Pi - 467895213/148935672]
> 8.49969*10^-11
> {{467895213, 148935672, 51988357/16548408},
> {429751836, 136794258, 23875102/7599681}}
>
> N[Pi - Last[lo]]
> 8.49969*10^-11
>
> Mark McClure
It took some time but I finally found the problems in my code. One is
that the ordering requires numbers (exact or approximate), and i was
using exact Pi. THe other is that I was forgetting an interval for a
candidate can intersect the negative of the interval of the best thus
far, and still be required for retaining for the next round.
Here is the corrected version. I use 20 digits precision for
approximating Pi. In some problems (but not this one), higher precision
could be required.
frac[{n_Integer, d_Integer}] := n/d
closest[pairs_, target_] :=
pairs[[Ordering[Abs[Map[frac, pairs, 1] - target], 1]]][[1]]
intervalize[{num_Integer, den_Integer}] :=
Interval[{num/(den + 1), (num + 1)/den}]
siftPairs[pairs_, target_] := Module[
{interval, num, den, neginterval, ntarget = N[target, 20]},
{num, den} = closest[pairs, ntarget];
interval = intervalize[{num, den}] - ntarget;
interval = IntervalUnion[interval, -interval];
Select[pairs,
IntervalIntersection[intervalize[#] - ntarget, interval] =!=
Interval[] &]]
successors[{num_Integer, den_Integer}] :=
Module[{nd, dd, cnd, cdd, nums, dens}, {nd, dd} =
IntegerDigits[{num, den}];
cnd = Complement[Range[9], nd];
cdd = Complement[Range[9], dd];
nums = Map[10*num + # &, cnd];
dens = Map[10*den + # &, cdd];
Flatten[Outer[List, nums, dens], 1]]
initpairs = Flatten[Outer[List, Range[9], Range[9]], 1];
In[308]:= Timing[
candidates =
Nest[Flatten[Map[successors, siftPairs[#, Pi]], 1] &, initpairs,
8];]
Out[308]= {9.57455, Null}
In[310]:= best = closest[candidates, N[Pi, 20]]
Out[310]= {467895213, 148935672}
In[311]:= bestfrac = frac[best]
Out[311]= 51988357/16548408
In[312]:= N[bestfrac - Pi, 5]
Out[312]= -8.4997*10^-11
One can prove this method captures the best value (even without knowing
that value from two other responses), assuming there are no further bugs
in the code.
A slightly faster variant replaces intervalize[] and siftPairs[] with
the code below.
intervalize2[{num_Integer, den_Integer}] := {num/(den + 1), (num + 1)/
den}
siftPairs2[pairs_, target_] := Module[
{imax, num, den, ntarget = N[target, 20]},
{num, den} = closest[pairs, ntarget];
imax = Max[Abs[intervalize2[{num, den}] - ntarget]];
Select[pairs, With[{int = ntarget - intervalize2[#]},
(Min[Abs[nint]] < imax || Apply[Times, int] < 0)] &]]
In[324]:= Timing[
candidates2 =
Nest[Flatten[Map[successors, siftPairs2[#, Pi]], 1] &, initpairs,
8];]
Out[324]= {4.13437, Null}
In[326]:= best = closest[candidates2, N[Pi, 20]]
Out[326]= {467895213, 148935672}
Daniel Lichtblau
Wolfram Research
Daniel Lichtblau
> Hello, I am a high school math teacher and the following puzzle was
> posed by a few math teachers I am in contact with.
>
> Create a fraction whose numerator has the digits 1 - 9 (used once)
> and whose denominator has the digits 1 - 9 (used one) .
>
> Which fraction has a value closest to the value of pi?
>
> I've worked on some "brute force" checks and managed to check all
> possible fractions with 2,3,4,5 and 6 digits. But after that, there
> are just too many possibilities.
>
> I don't have the programming ability to implement something elegant in
> Mathematica.
>
> Is there anyone who could suggest an approach to find the solution to
> this?
>
> Sincerely,
>
> Tom
>
Thought I would post a different method. We set this up as an integer
programming (ILP) problem. For each digit we have nine variables; they
will have values of zero or one, and exactly one of them will be one.
For example, if the sixth variable for the third numerator digit is one,
than that digit is to be regarded as six. We furthermore restrict so
that exactly one sixth variable in the numerator is one (that is,
exactly one digit must be equal to six).
We enforce these constraints by insisting that all these variables be
between zero and one inclusive, and that appropriate families of nine
sum to one. Conveniently, this amounts to having row sums and column
sums of one, when we lay the variables out in a grid (separate grids for
numerator and denominator variables).
We formulate an objective function that amounts to minimizing the
absolute value of numerator/denominator - Pi (by minimizing some new
variable constrained to be at least as large as that value and also its
negative). We suitably rationalize that Pi and clear denominators, so we
have an integer program.
We can last insist that all variables be integer valued, and call
Minimize. There is a problem with this. It works, but not very fast
(took around 13 hours overnight to run to completion, only for me to
discover I had a typo in a constraint and so the solution was invalid).
The bottleneck is that the internal ILP code uses branch-and-bound, and
has no good way to realize which variables to branch on. But we do. We
know it is important to settle the variables first that correspond to
high digits in numerator and denominator. So I'll solve relaxed (that
is, ordinary) LPs in a loop, and whenever we get an integer solution we
have a new minimizer. Else we branch on the highest variable taking a
non-integral value.
For (my programming) convenience, I set up with explicit variables
and use Minimize. One can bypass some run time overhead by using
LinearProgramming directly.
Here is the code.
ratApproximate[n_, target_] :=
Module[{nvars, dvars, vars, num, den, numer, denom, approx, pnum,
pden, cn1, cd1, cn2, cd2, obj, val, constraints, program, stack,
vals, min, best = Infinity, soln = {}, badvar, counter = 1,
avars},
nvars = Array[num, {n, n}];
dvars = Array[den, {n, n}];
vars = Flatten[Reverse[Riffle[nvars, dvars]]];
numer = Expand[(n + 1)^Range[0, n - 1].nvars.Range[n]];
denom = Expand[(n + 1)^Range[0, n - 1].dvars.Range[n]];
cn1 = Map[Total[#] == 1 &, nvars];
cn2 = Map[Total[#] == 1 &, Transpose[nvars]];
cd1 = Map[Total[#] == 1 &, dvars];
cd2 = Map[Total[#] == 1 &, Transpose[dvars]];
approx = Rationalize[N[target, 20], 10^(-16)];
pnum = Numerator[approx];
pden = Denominator[approx];
val = pden*numer - pnum*denom;
constraints =
Join[cn1, cn2, cd1, cd2,
Map[0 <= # <= 1 &, vars], {obj >= val, obj >= -val}];
stack = {constraints, {}};
avars = Append[vars, obj];
While[stack =!= {},
counter++;
{constraints, stack} = stack;
If[best =!= Infinity,
constraints = Join[constraints, {obj <= best - 1}]];
Quiet[vals = Minimize[{obj, constraints}, avars]];
If[Head[vals] == Minimize || ! FreeQ[vals, Indeterminate],
Continue[]];
{min, vals} = vals;
badvar =
Position[vars /. vals, (a_ /; ! IntegerQ[a]), {1}, 1,
Heads -> False];
If[badvar == {}, best = min; soln = vals; Continue[]];
badvar = vars[[badvar[[1, 1]]]];
stack = {Append[constraints, badvar == 0], stack};
stack = {Append[constraints, badvar == 1], stack};];
{counter, best, {numer, denom} /. soln}]
This is still not blindingly fast. Finds the pair {391625847,124658379}
without too much trouble. Then takes many iterations to get first
{35642934,{428913756,136527489}} and then
{1090368,{429751836,136794258}}. That latter appears after around 5000
trips through the While loop. {467895213,148935672} shows up before 9000
iterations. Alas, after two hours it still has not become convinced it
is finished. Large search space for the relaxed subproblems, I guess.
Daniel Lichtblau
Wolfram Research
Daniel Lichtblau
> The code I gave can certainly be speeded up, but right now I'm more
> concerned about the algorithm. 'hi' is meant to be the smallest value
> of the form s[[i]]/s[[j]] that is greater than Pi, and 'lo' is meant
> to be the largest value of the form s[[i]]/s[[j]] that is less than
> Pi. 'hi' is OK, but I wonder about 'lo'.
I have not looked carefully at the code. The basic idea is smarter than
I thought (more correctly, my understanding was stupider than necessary).
A variant would do an iteration over valid denominators. For each, find
the integer numerator that gives the closest approximation to pi. Then
walk in both directions (that is, increment and decrement) until you
find a valid numerator, and see if either gives a result better than teh
best found thus far.
Or only look at valid numerators (that is, those formed by permutations
of the available digits), and do a binary search to get to the one
giving a value closest to pi. For this one would need to sort the
available values. not to onerous when there are only 9! or so.
Daniel
dh
I did not take the trouble to go through all proposed solutions, so it
could be that somebody else already got the same idea.
a) Consider a sorted list of all permutations of 0..9, call it: per.
b )consider a function quot[x1_,x2_]:=per[[x1]]/per[[x2]]
c) consider a 2 dim grid of indices (1..Length[per]) X (1..Length[per])
d) consider the surface defined by f evaluated on this grid. Note that
this surface is monotonically ascending from the upper right to the
lower left (low raw-high column to high row-low column index)We are
looking for the contour line f==Pi. As this line may lay between grid
points, we are looking for a line on the grid with maximal f and the
condition f<=Pi. Note that along this line the column indeces are
increasing (not necessarily strict) with row indeces. This makes it
simple to calculate.
e) Next we check which grid point on this "contour line" is closest to
Pi from below.
On my machine the code to create and sort the possible numbers takes:
13 sec.
The code for the search takes 67 sec. If Pi is not rationalized the time
is: 77 sec
The result is: {9735046128 / 3098761425} what is 2.4 10^-12 below Pi.
The approximation from above should be obvious.
Here is the code:
===========================
per = Sort[FromDigits /@ Permutations[{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}]];
p = Rationalize[Pi, 10^-20];
len = Length[per]; Print["len=", len];
ind[x_] := Module[{i1 = 1, i2 = Length[per], i},
While[i2 - i1 > 1,
i = Floor[(i1 + i2)/2];
If[x < per[[i]], i2 = i, i1 = i]; Print[{i1, i2, x < per[[i]]}];
];
i1
];
quot[i1_, i2_] := per[[i1]]/per[[i2]];
col = row = 1;
max = 0; row0 = 1; col0 = 1;
While [quot[row + 1, col] < p, ++row]; Print["row0=", row];
While[row < len && col < len,
If[quot[row + 1, col] < p, ++row, ++col];
If[(t = quot[row, col]) > max, max = t; row0 = row; col0 = col;];
];
Print["Error,num,denum=", {max - Pi // N, per[[row0]], per[[col0]]}]
================================================
Daniel
--
Daniel Huber
Metrohm Ltd.
Oberdorfstr. 68
CH-9100 Herisau
Tel. +41 71 353 8585, Fax +41 71 353 8907
E-Mail:<mailto:d...@metrohm.com>
Internet:<http://www.metrohm.com>
Ray Koopman
> The code I gave can certainly be speeded up, but right now I'm more
> concerned about the algorithm. 'hi' is meant to be the smallest value
> of the form s[[i]]/s[[j]] that is greater than Pi, and 'lo' is meant
> to be the largest value of the form s[[i]]/s[[j]] that is less than
> Pi. 'hi' is OK, but I wonder about 'lo'.
>
> Specifically, potential values of 'lo' with j > k are never checked,
> where k is the largest value of j for which s[[-1]]/s[[j]] > Pi. When
> the While[j < k, ...] loop exits, we have j = k, s[[i]]/s[[j]] > Pi,
> and s[[i-1]]/s[[j]] < Pi. s[[-1]]/s[[k+1]] < Pi, but might it not
> also be > the current 'lo'?
One more check is necessary. If j > k then s[[i]]/s[[j]] < Pi for
all i, so we need to check the largest possible new value for 'lo'.
Here is the revised code, compiled.
hiloPi = Compile[{}, Module[{s,i,j,k,hi,lo,ijhi,ijlo},
s = (Permutations@Range@9).(10^Range[8,0,-1]);
i = j = 1; While[s[[i]]/s[[j]] < Pi, i++];
hi = s[[ i ]]/s[[j]]; ijhi = { i ,j};
lo = s[[i-1]]/s[[j]]; ijlo = {i-1,j};
k = Length@s; While[s[[-1]]/s[[k]] < Pi, k--];
While[j < k, j++; While[s[[i]]/s[[j]] < Pi, i++];
If[s[[ i ]]/s[[j]] < hi, hi = s[[ i ]]/s[[j]]; ijhi = { i ,j}];
If[s[[i-1]]/s[[j]] > lo, lo = s[[i-1]]/s[[j]]; ijlo = {i-1,j}]];
If[s[[-1]]/s[[k+1]] > lo, lo = s[[-1]]/s[[k+1]]; ijlo = {-1,k+1}];
{s[[ijhi]],s[[ijlo]]}]]
Timing@hiloPi[]
N[Divide@@#-Pi]&/@%[[2]]
{1.23 Second,{{429751836,136794258},{467895213,148935672}}}
{1.0185541299279066*^-10, -8.499689840846258*^-11}