[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Uniform partition of an interval
From: |
Dirk Laurie |
Subject: |
Re: Uniform partition of an interval |
Date: |
Thu, 18 May 2000 10:18:00 +0200 (SAST) |
On Wed, 17 May 2000, John W. Eaton wrote:
> On 17-May-2000, address@hidden <address@hidden> 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
-----------------------------------------------------------------------