Mathematica 3.0 LogZero: an expert's response

*To*: mathgroup at smc.vnet.net*Subject*: [mg24515] Mathematica 3.0 LogZero: an expert's response*From*: William Boyd <william at kamaz.demon.co.uk>*Date*: Thu, 20 Jul 2000 03:01:50 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

Well ... no sooner had I posted the query below I got the following blistering fast response from Wolfram support, which I'm very happy to share in case my first post got past the moderator (so okay, no more posts ... but errors are always instructive, so why not keep it in?) Very grateful to Wolfram: thank you. <<response follows>> This calculation is fundamentally flawed for reasons unrelated to Mathematica. The inputs 1. and 0.0015723270440251572 are ordinary low-precision (about 16 digits) machine numbers. It not useful to do high-precision calculations with low-precision inputs. There is in fact no reason to use SetPrecision in any calculation of this type. A deeper problem is that the arguments appear to be specifically chosen so that the argument 1 - r*s/(n*(1 - r)) of the logarithm is near zero, which means that the value of r*s/(n*(1 - r)), which is subtracted from 1, is very close to 1. If r*s/(n*(1 - r)) is equal to 1 to within, say, 16 digits, then obviously more than 16 digits will be required to get a reliable result. To get correct results for this calculation in Mathematica or in any other system it is necessary to either rewrite this program or choose inputs carefully to insure that the results are not artifacts of elementary numerical error. For example, with the machine number 1. as the first argument in your "result" function, the other parameters must be chosen so that r*s/(n*(1 - r)) always differs from 1 at least as many digits before the 16th decimal place as there are digits in the output of the Floor function. To get 5 reliable digits in the output of the Floor function, there must be at least 5 digits of precision in the argument of the Floor function, which means that there must be at least 5 reliable digits in the argument of the logarithm, which means that r*s/(n*(1 - r)) must be calculated to sufficient precision to leave at least 5 digits when it is subtracted from 1. If this were my calculation I might start with a definition such as In[1]:= result[n_, s_, r_] := With[{e = Log[1 - r*s/(n*(1 - r))]/Log[1 - r]}, {"precision" -> Precision[e], Floor[e] + 1}] which would report to me both the result and the precision of the critical intermediate result so that I could see for myself if the result was reliable. For example: In[2]:= result[1, 635, 0.001572327044025157200000] Out[2]= {precision -> 6, 24410} is reliable because the precision of the argument of Floor is greater than the number of digits in 24410. The result from In[3]:= result[1, 635, 0.00157232704402515722] Out[3]= {precision -> 2, 25011} is not reliable, because the precision (2 digits) of the argument of Floor is less than the number of digits reported in the final result. Dave Withoff Wolfram Research <<original request>> 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? All I'm doing is calculating the number of terms needed to reach a given sum in the geometric series: following is typical: the values are the first few terms of a denary limit sequence tending to 1/636, bar the last ordinate which is last digit failsafed: I hope I am working to 50 digit precision (I do find The Book hard work at times) - even if I'm not it should surely be better than this?. The first Log you see in 'result' is the one tending to zero. So we are doing LogZero / LogOne, which is a test of course, but not such a hard test. <<start cut and paste of Mathematica 3.0 notebook: formula is slightly different from standard as I am modelling with first term of series n*(1-r) >> 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] 19022 result[1.,635,0.0015723270440251] 19667 result[1.,635,0.00157232704402515] 20971 result[1.,635,0.001572327044025157] 23347 result[1.,635,0.0015723270440251572] 23347 result[1.,635,0.00157232704402515722] \[Infinity] << finish >> On the face of it, this isn't impressive. Compare my results for the same formula with the same ordinates using the Extended type (19/20 sig figs) in Delphi 3, which gives {19022, 19665, 20979, 23163, 24409, 25008}, while a binary divide algorithm in Delphi using a tried and trusted integer power raise of mine (no logs) gives {19022, 19665, 20979, 23163, 24410, 25014}. Comments appreciated. -- William Boyd