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.