MathGroup Archive 2010

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

Search the Archive

Re: Pi day

  • To: mathgroup at smc.vnet.net
  • Subject: [mg108402] Re: Pi day
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Tue, 16 Mar 2010 04:47:31 -0500 (EST)

Mark McClure wrote:
> On Mon, Mar 15, 2010 at 1:05 AM,  <danl at wolfram.com> wrote:
>> 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
> [...]
> {{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


  • Prev by Date: Re: Video compression
  • Next by Date: Re: Pi day
  • Previous by thread: Re: Pi day
  • Next by thread: Re: Pi day