On Thursday 29 January 2004 10:26 pm, J. Riley Bryant wrote:
> Let PR(x,y,f(n)) be the Product of f(n) for x<=n<=y.
>
> Probability = PR(1,8,(1-(n!-1)/n!)^1e6)
>
> = (1-(0/1)^1e6) *
> (1-(1/2)^1e6) *
> (1-(5/6)^1e6) *
> (1-(23/24)^1e6) *
> (1-(119/120)^1e6) *
> (1-(719/720)^1e6) *
> (1-(5039/5040)^1e6) *
> (1-(40319/40320)^1e6)
>
> And herein lies my problem. I have yet to find a
> calculator which is capable of handling any of
> the above given calculations (with the exception
> of the first and last two expressions. (in fact on
> a certain TI calculator all expressions evaluate
> to 1!)
>
> So I'm stuck with an approximation...
...
> So I propose that our probability is greater than
> 0.99999999988148537165638680164463 and less than 1.
> Kudos to anyone who can give a better approximation
> (or even a definite answer). I'd be extremely
> interested to see it.
OK. I'll do this rather formally, to be sure of not making
too many mistakes :-). This will make it look harder than
it really is, and therefore make me look cleverer than I
really am, which is an added bonus.
Define M := 10^6, for convenience.
Lemma 1: product(1-a_j) >= 1 - sum(a_j),
provided all a_j are >=0 and their sum <= 1.
Proof: Take logs of both sides and use the power series;
what we have to prove then is
sum(a+a^2/2+...) <= sum(a) + sum(a)^2/2 + ...
which is trivial.
Lemma 2: Each factor is <= 1.
Proof: Trivial.
Lemma 3: exp(-b/(a-1)) <= (1-1/a)^b <= exp(-b/(a-1/2)),
provided a,b >= 0.
Proof: Raising both sides to the power 1/b eliminates b
from the inequality, so we need only prove that
exp(-1/(a-1)) <= 1-1/a <= exp(-1/(a-1/2)).
Taking logs of both sides and using the power series,
this is equivalent to
1/(a-1/2) <= 1/a + 1/2a^2 + 1/3a^3 + ... <= 1/(a-1).
Expanding 1/(a-1/2) and 1/(a-1) by the binomial theorem,
these inequalities hold term by term and we're done.
Remark: I shan't actually use the first of these inequalities,
only the second. If we were satisfied with rather crude
bounds for the probability, using them both together
would be one way to get them.
Lemma 4: (1-1/n!)^M <= 10^-200 when n <= 6.
Proof: By lemma 3 it's enough to show that exp(-M/(n!-1/2)) <= 10^-200
when n <= 6. That's the same as n!-1/2 <= M/(200 log 10).
Since log 10 <= 5, this will be true provided n!-1/2 <= M/1000 = 1000,
which is certainly true for n <= 6.
Remark: In fact, (1-1/6!)^M < 10^-600.
For the next bit, unfortunately, we need to get our hands dirty with some
rather big numbers: it's time to calculate (1-1/7!)^M and (1-1/8!)^M
to 100 decimal places or so. One approach is to wheel out a high-precision
floating-point package, but I'm aiming for real, rigorous upper and lower
bounds, so here's something different.
Define A := 5039^5000, B := 5040^5000, C := 40319^40000, D := 40320^40000.
Then the numbers we want are (A/B)^200 and (C/D)^25.
Lemma 5A: A/10^18312 = 5269771432711962922797772729053779435700
3306925197638248587732614025743996470364
2682529403556003438342128396907785986976
1875748784748640171621187961146082636632
2143720298334934729715414441776354772118 + hA
where 0 <= hA < 1. (Call the big number A1.)
Lemma 5B: B/10^18312 = 1421288452744759435423898222790265050942
5405416550462194378268180305172687845769
7986109051946389287147480245974212005772
7984614307025265064636883921176404716321
55787621618757999727503295686887331079200 + hB
where 0 <= hB < 1. (Call the big number B1.)
Lemma 5C: C/10^184021 = 2455186319536720997394676499477484861209
6961077043822110290713894759829887860490
9497305157736347551149447104388193674108
5852569740617228041603426510294805581403
5712803327897508773278606167805672812767 + hC
where 0 <= hC < 1. (Call the big number C1.)
Lemma 5D: D/10^184021 = 6621212080452351946548268289624806230274
5889805059894616844352744147296627748863
8181291177918933843384166297558551040253
0232577723799885488948438644653963932640
9915717711584027952082406148988567516547 + hD
where 0 <= hD < 1. (Call the big number D1.)
Proofs: Get a computer to calculate A,B,C,D exactly and divide by
the appropriate power of 10.
Lemma 6: A1^200/(B1+1)^200 <= (A/B)^200 <= (A1+1)^200/B1^200;
C1^25/(D1+1)^25 <= (C/D)^25 <= (C1+1)^25/D1^25.
Proof: Trivial.
Lemma 7: (A/B)^200 lies between X and X + 10^-200, where
X = 0.0000000000000000000000000000000000000000
0000000000000000000000000000000000000000
0000006636059073569760123444030727811030
4962185423987277494669213425599007507328
4137452366231171911058932143962345736928.
Proof: Get a computer to calculate 10^200*A1^200 and (B1+1)^200
exactly for a lower bound, and 10^200*(A1+1)^200 and B1^200
for an upper bound.
Lemma 8: (C/D)^25 lies between Y and Y + 10^-200, where
Y = 0.0000000000169306611928046844874478276987
5571187874601424722391358424160013010721
3577983242881958747595256224286340810473
6494862762447195565026554299072497768865
6955670811963626008104251476655979008827.
Proof: Just like lemma 7.
Lemma 9: (1-(A/B)^200)*(1-(C/D)^200) = 1-Z + something
of absolute value < 3*10^-200, where
Z = round(X+Y-XY, 200dp)
= 0.0000000000169306611928046844874478276987
5571187874601424722391358424160013010721
3577989878941032205002511838270124200472
1765300425907875782610141581512192189174
2139024675906339250158229567980937981653.
Proof: Let h,k be the errors associated with X,Y, so that
(1-X-h)(1-Y-k) = 1-X-Y+XY-h-k+Yh+Xk+hk.
Everything here is >= 0, so the absolute value
of the error in Z is h+k-(Yh+Xk+hk) < h+k <= 2*10^-200.
Truncating XY to 200 places incurs at most a further
10^-200 error.
Theorem: The first 150 decimal places of the thing we're
trying to calculate are the same as those of 1-Z,
namely
0.99999999998306933880719531551255217230124428812125
39857527760864157583998698927864220101210589677949
97488161729875799527823469957409212421738985841848.
Proof: By lemmas 2 and 4, the contribution of the terms for
n<=6 is at most 6*10^-200. The error term in lemma 9
is at most 4*10^-200. So in fact at most the last
couple of digits of Z could be wrong.
Here are the first 2000 digits, calculated by doing floating-point
with 65536 bits of precision (which should be vastly more than
there's any need for). I expect they're all correct, but unlike
the 150-digit answer above it's at the mercy of a big pile of code
I didn't write and that might have subtle mistakes in it. The
150-digit answer depends only on the ability to add, subtract,
multiply and divide arbitrary-length integers correctly.
0.99999999998306933880719531551255217230124428812125
39857527760864157583998698927864220101210589677949
97488161729875799527823469957409212421738985841848
78078108257860975324093660749841770432019062018345
48747152344080723120184504811006034197073064335556
82141632534630176440969388177383530613588449981243
69149115750593004019539194569006349799715508094223
38183065349293201122429471638809686718098472840860
14444227009603306943152616078860259529981305126906
12401563602948391896507801877421182051111615759708
61590163194267433542182111431147018590713969698194
67955189833415128513109235102831893192723847679404
88955536430950612492907834496373190433905122541037
87992775405186468613075200206377931089523189278254
47543705167112641839907316917801572626452979319933
85012263419230925400115707038902158184960789139485
79246080583396118117027216223214009875508191592939
61497613988204038740805450897526546103308873084023
12305002335143312533441858960735777947327079069181
19378115072078712709568365487089460760152873317341
42331374758588112328351800454518853358512223560946
56567801110488274631481396515189636402013925894993
65274106970209580641754478099008330133011041906272
17675353426491330544667704019507505831937637714342
01449630767930931233263754273416990101382660945150
05127908688616635483985919910910153606727246723519
74481731343674185001201203668954763684165368100232
51896495965676992394612318823628327440342336960457
00476745682734100109926059092913728785645259662068
84562495156182835613986249168069542020515919785354
53837926039416965171040801174358814972868649267155
95436739345020205868522372252009235521795052068046
08581652036836723713547145973089520206999950651502
48369460712775136027718545730351753649484059156138
81739575818805033239936125653422551473737340511640
40908522486573109506495883076487483729465051004608
94285633013109108354746060391981465465718913908586
74419808443668479896095440484735011351446801312651
86673820779888844436742633756330948790093076249979
23463289968558271498088065117983280098344154145906
Oh, if you're distributing kudos, give some to Common Lisp
and in particular to the implementation called CLISP, which
has nice fast bignum operations and arbitrarily high-precision
floating-point. I could have used Perl for the calculations
above, but it would have been less pleasant.
--
g
