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