Re: Mathematica 3.0: reliability close to LogZero?
- To: mathgroup at smc.vnet.net
- Subject: [mg24531] Re: [mg24514] Mathematica 3.0: reliability close to LogZero?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Mon, 24 Jul 2000 03:04:10 -0400 (EDT)
- References: <200007200701.DAA28925@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
William Boyd wrote:
>
> 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.
>
> sincerely,
>
> --
> William Boyd
I have not tested this on version 3, so I cannot say whether there is a
version-dependent issue. Regardless, your code, as written, cannot do
what you expect. The problem is that machine arithmetic is used before
SetPrecision, and there is no hope for correcting catastrophic
cancellation in machine arithmetic.
Rewriting a bit to remove the inexact 1., I get the following behavior.
In[57]:= result[s_Integer, r_Real] :=
Floor[SetPrecision[Log[1-r*s/(1-r)] / Log[1-r], 50]] + 1
In[58]:= result[635, 0.001572327044025157]
Out[58]= Indeterminate
In[59]:= result[635, 0.00157232704402515722]
Out[59]= 25011
I believe this is appropriate. The first one had machine numbers for
input, hence could not handle the cancellation and corresponding loss of
precision. The second time we have enough digits to handle this.
To force higher precision throughout, you might do something like the
following (I restore the use of n here).
In[63]:= result[n_Real, s_Integer, r_Real] := With[
{n50=SetPrecision[n, 50], r50=SetPrecision[r,50]},
Floor[Log[1-r50*s/(n50*(1-r50))]/Log[1-r50]] + 1
]
We see that this handles the bad case above that previously gave
Indeterminate.
In[64]:= result[1., 635, 0.001572327044025157]
Out[64]= 23522
Daniel Lichtblau
Wolfram Research
- References:
- Mathematica 3.0: reliability close to LogZero?
- From: William Boyd <william@kamaz.demon.co.uk>
- Mathematica 3.0: reliability close to LogZero?