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