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