MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Pi day

  • To: mathgroup at
  • Subject: [mg108424] Re: Pi day
  • From: Daniel Lichtblau <danl at>
  • 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,
   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 =!= {},
    {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],
    {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

  • Prev by Date: Re: Video compression
  • Next by Date: Re: Axeslabel containing capital n
  • Previous by thread: Re: Pi day
  • Next by thread: Re: Pi day