Re: Mathematica 3.0: reliability close to LogZero?
- To: mathgroup at smc.vnet.net
- Subject: [mg24547] Re: Mathematica 3.0: reliability close to LogZero?
- From: William Boyd <william at kamaz.demon.co.uk>
- Date: Mon, 24 Jul 2000 03:04:26 -0400 (EDT)
- References: <8l6ab4$sml@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <8l6ab4$sml at smc.vnet.net>, William Boyd
<william at kamaz.demon.co.uk> writes
>Could I ask an expert to cast a jaded eye over the following: is there
>something wrong with Log close to zero in Mathematica 3.0?
I would like to thank my two respondents for their very helpful replies
and advice. I am very grateful.
There are still some points of interest if anyone would care to bear
with me.
Briefly I was worried about computing something near LogZero / LogOne.
Expecting trouble, I found it:
result[n_Real,s_Integer,r_Real]:=
Floor[SetPrecision[Log[ 1 - r*s/(n*(1 - r))] / Log[1 - r],50]]+ 1
result[1.,635,0.001572327044025]
result[1.,635,0.0015723270440251]
result[1.,635,0.00157232704402515]
result[1.,635,0.001572327044025157]
result[1.,635,0.0015723270440251572]
result[1.,635,0.00157232704402515722]
19022
19667
20971
23347
23347
\[Infinity]
Well, as both my respondents pointed out, this is entirely misconceived.
I've only entered 1. to one significant figure and the call to
SetPrecision is pointless given that the inputs are (mostly) just
machine precision. Better would have been:
betterResult[n_Real,s_Integer,r_Real]:=
Floor[Log[ 1 - r*s/(n*(1 - r))] / Log[1 - r]]+ 1
result[1.0000000000000000000,635,0.0015723270440250000000]
result[1.0000000000000000000,635,0.0015723270440251000000]
result[1.0000000000000000000,635,0.0015723270440251500000]
result[1.0000000000000000000,635,0.0015723270440251570000]
result[1.0000000000000000000,635,0.0015723270440251572000]
result[1.0000000000000000000,635,0.0015723270440251572200]
19022
19665
20979
23163
24410
25011
where I've entered 20 significant digits. The output is now quite
plausible and in line with an independent 'binary divide and seek
better' machine precision computation.
This still doesn't answer the question: how much faith can I have in the
output?
My first respondent, Dave Withoff of Wolfram Research, suggested
reporting the precision of the calculation at the same time:
firstSuggestion[n_Real, s_Integer, r_Real] :=With[
{e = Log[1 - r*s/(n*(1 - r))]/Log[1 - r]},
{"precision" -> Precision[e], Floor[e] + 1}
]
firstSuggestion[1.0000000000000000000,635,0.0015723270440250000000]
firstSuggestion[1.0000000000000000000,635,0.0015723270440251000000]
firstSuggestion[1.0000000000000000000,635,0.0015723270440251500000]
firstSuggestion[1.0000000000000000000,635,0.0015723270440251570000]
firstSuggestion[1.0000000000000000000,635,0.0015723270440251572000]
firstSuggestion[1.0000000000000000000,635,0.0015723270440251572200]
{precision -> 7, 19022}
{precision -> 7, 19665}
{precision -> 6, 20979}
{precision -> 5, 23163}
{precision -> 4, 24410}
{precision -> 4, 25011}
so I should distrust the last two results as insufficiently precise.
My second respondent, Daniel Lichtblau also of Wolfram research,
suggested forcing precision in the following way:
secondSuggestion[n_Real, s_Integer, r_Real] := With[
{n20=SetPrecision[n, 20], r20=SetPrecision[r,20]},
Floor[Log[1-r20*s/(n20*(1-r20))]/Log[1-r20]] + 1
]
secondSuggestion[1.,635,0.001572327044025]
secondSuggestion[1.,635,0.0015723270440251]
secondSuggestion[1.,635,0.00157232704402515]
secondSuggestion[1.,635,0.001572327044025157]
secondSuggestion[1.,635,0.0015723270440251572]
secondSuggestion[1.,635,0.00157232704402515722]
19023
19665
20974
23522
24410
25011
and, for me, this is lightbulb idea: well, I was getting there ...
honest ;-)
But what is fascinating is how different that fourth result is (23522
against 23163 - that's 1.5 % different). Yet on the face of it both are
calculating to 20 significant figures.
So easy to check; as I mentioned in my original post this is just a seek
for the smallest integer n that gets a given partial sum (635 in this
case) in a geometric series. Using each idea, I trust the following is
correct:
firstSum[solution_Integer, n_Real, r_Real]:= With[
{e = (n * (1 - r) * ( 1 - (1 - r)^solution))/r},
{"precision" -> Precision[e], Floor[ e]}
]
secondSum[solution_Integer, n_Real, r_Real]:= With[
{n20 = SetPrecision[n,20], r20 = SetPrecision[r,20]},
Floor[(n20 * (1 - r20) * ( 1 - (1 - r20)^solution))/r20]
]
and this confirms the result 25011:
firstSum[25010,1.0000000000000000000,0.0015723270440251572200]
firstSum[25011,1.0000000000000000000,0.0015723270440251572200]
secondSum[25010, 1., 0.00157232704402515722]
secondSum[25011, 1., 0.00157232704402515722]
{precision -> 19,634}
{precision -> 19,635}
634
635
but advances me no further with 23163 versus 23522:
firstSum[23162,1.0000000000000000000,0.0015723270440251570000]
firstSum[23163,1.0000000000000000000,0.0015723270440251570000]
secondSum[23521, 1., 0.001572327044025157]
secondSum[23522, 1., 0.001572327044025157]
{precision -> 19, 634}
{precision ->19, 635}
634
635
Comments?
These figures are, by the way, not entirely unconnected with the number
of pages in the latest 'Harry Potter' saga (636 in the British edition).
My researches in this matter are perhaps a little too arcane for
detailed discussion in a newsgroup but, put in simple terms an
intelligent lawyer might reasonably be expected to understand, I'm at
pains to discover the maximum overwhelming feeling of boredom I can
allow myself in reading 'Harry Potter IV' consistent with eventually
finishing the b*****d sometime in the hopefully not too close future.
So is fundamentally important research of huge relevance parents the
world over and you no sneer because best thing I ever done. And okay I
no violate premises of newsgroup by giving address of Web page where to
find 'HairyPooter IV¾ ReadTime Calculator and Bookmark' but worth
searching for on Web I can tell you (but not yet 'cos I'm still
struggling with the hard maths ... oh dear).
I am indeed most grateful for my two responses: I learnt a lot about
Mathematica as a result.