Re: Pi day
- To: mathgroup at smc.vnet.net
- Subject: [mg108424] Re: Pi day
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 17 Mar 2010 04:40:05 -0500 (EST)
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
>
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