From help-octave-request at bevo dot che dot wisc dot edu Thu May 18 03:18:05 2000 Subject: Re: Uniform partition of an interval From: Dirk Laurie To: john at arrows dot demon dot co dot uk (John Smith) Cc: jwe at bevo dot che dot wisc dot edu (John W. Eaton), emil@gollum.fri.uni-lj.si, help-octave@bevo.che.wisc.edu (Octave help) Date: Thu, 18 May 2000 10:18:00 +0200 (SAST) On Wed, 17 May 2000, John W. Eaton wrote: > On 17-May-2000, emil at gollum dot fri dot uni-lj dot si wrote: > > | Hello! > | > | What is the MAIN reason that 1.8:0.05:1.9 produces [1.8000 1.8500] > | and not [1.8000 1.8500 1.9000]? > | I am using 2.0.14 version of Octave. > | Thank you for your answer. > | Best regards, > | Emil Zagar > > I'd guess that the MAIN reason is that there is a bug in the way > Octave is trying (very hard) to compute the correct number of elements > for ranges. If you're in a debugging mood, the code to look at is in > the Range::nelem_internal and related functions in liboctave/Range.cc. > > jwe It is obscene to write code that intimately explores the delicate secrets of floating-point arithmetic on a particular implementation. The decent programmer uses the colon range operator only on integers, or in cases where the total range is not close to a multiple of the step size. I.e. '1.8+(0:2)*0.05' or '1.8:0.05:1.91'. But we live in an age where seemliness is not part of the popular ethos, and people will write things like 'y=1.8:0.05:1.9'. We should agree what that should do. Intuitively one feels that a:h:b with h>0 should be equivalent to: y=[]; x=a; while x<=b, y=[y x]; x += h; end And indeed, if I run the above in Octave on my i686 machine, I get [1.8000 1.8500]. Yet it is unsatisfactory, because with pencil and paper, or on a decimal machine, or on some binary machines, I would have got [1.8000 1.8500 1.9000]. One can get round the problem by saying it should be equivalent to: r=(b-a)/h; n=round(r); if h*abs(n-r)>max(a,b)*eps, n=floor(r); end y=a+h*(0:n); But doing so would treat one case of a pervasive problem: the non-intuitiveness of floating-point comparison. A good cure should work in other places too. I think Octave should borrow an idea from the grandfather of interactive matrix languages, namely APL. This language has a built-in variable which in Octave we would call 'comparison_tolerance'. Then we could write: > comparison_tolerance = 0; > 1.8+0.1 == 1.9 ans = 0 > 1.8:0.05:1.9 ans = 1.8000 1.8500 > comparison_tolerance = eps; > 1.8+0.1 == 1.9 ans = 1 > 1.8:0.05:1.9 ans = 1.8000 1.8500 1.9000 There should also be a built-in variable 'absolute_comparison'. The definition of comparison operators would be e.g.: if absolute_comparison, equal = x+eps>=y && y+eps>=x; else equal = x*+eps*abs(x)>=y && y+eps*abs(y)>=x; end Just imagine this: > n=5; > comparison_tolerance = eps*n; > absolute_comparison = 1; > A=rand(n); A*inv(A)==eye(n) ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 The overheads to someone not requiring tolerant comparison would be minimal: test a boolean variable equal to comparison_tolerance>0 before making the floating-point comparison. The APL implementation with which I worked had only relative comparison, and restricted comparison_tolerance to be in the range 2^(-53) to 2^(-24), obviously in order to allow efficient implementation on that machine. Dirk Laurie ----------------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.che.wisc.edu/octave/octave.html How to fund new projects: http://www.che.wisc.edu/octave/funding.html Subscription information: http://www.che.wisc.edu/octave/archive.html -----------------------------------------------------------------------