Re: A problem with range objects and floating point numbers

[Daniel J Sebald]

On 09/26/2014 02:43 AM, s.jo wrote:
> Daniel Sebald wrote
[snip]
> > To be factorable doesn't necessarily require that the floating point
> > step size be an exact number system equivalent, i.e., that .1 equal
> > exactly 1/10.  It just requires the machine arithmetic to be as expected
> > when producing an operation result, i.e.,
> >
> > octave-cli:1>  rem(-2,.1) == 0
> > ans =  1
> > octave-cli:2>  rem(0,.1) == 0
> > ans =  1
> > octave-cli:3>  [-20:1:0]*.1
> > ans =
> >
> >    Columns 1 through 5:
> >
> >     -2.00000  -1.90000  -1.80000  -1.70000  -1.60000
> >
> >    Columns 6 through 10:
> >
> >     -1.50000  -1.40000  -1.30000  -1.20000  -1.10000
> >
> >    Columns 11 through 15:
> >
> >     -1.00000  -0.90000  -0.80000  -0.70000  -0.60000
> >
> >    Columns 16 through 20:
> >
> >     -0.50000  -0.40000  -0.30000  -0.20000  -0.10000
> >
> >    Column 21:
> >
> >      0.00000
> >
> >
> > Jo, what do you get for the above scripts?
> >
> > Dan
>
> Dan and others,
> Sorry for late reply.
>
> Your integerized stepping seems to be a decimation process.
> You used step-size 1. Any bigger integers such as 10 or 10^17 are allowed in
> the decimation process.

By decimation, I take it you mean somehow converting the number to true 
decimal representation, similar to the decimal arithmetic toolbox Oliver 
alluded to.  I don't think that's the case.


> Comparing the results of colon operators with floating-point numbers,
> Matlab seems to use such decimation process, but Octave does not.
> I agree that this decimation feature of colon operator is far from any IEEE
> standard.
>
> Also I point out that the results of linspace function and colon operator in
> Matlab are different.
> This may imply decimation processes are different.
> By inspection, the linspace result is much close to decimal number.
> I guess that the decimation process of Matlab's colon operator is simpler
> for speed.

No, I doubt that.  The linspace and colon operator are probably two 
slightly different algorithms.  For linspace the spacing is exactly an 
integer, by definition.  That's not necessarily so for the colon 
operator where the spacing can result in a fraction of the step size at 
the last interval.

Looking at the code a bit, I see that the range (colon) operator does 
attempt to use integer multiplication as much as possible, e.g., from 
Range::matrix_value:

      cache.resize (1, rng_nelem);
      double b = rng_base;
      double increment = rng_inc;
      if (rng_nelem > 0)
        {
          // The first element must always be *exactly* the base.
          // E.g, -0 would otherwise become +0 in the loop (-0 + 
0*increment).
          cache(0) = b;
          for (octave_idx_type i = 1; i < rng_nelem; i++)
            cache(i) = b + i * increment;
        }

So, the accumulation of addition errors doesn't seem to be an issue here.


> So far, I learned that integer range can be generated 'safely' with colon
> operator.
> Regarding to floating-point number range, we may want to use the function
> "linspace" for the best decimation.
>
> I don't know who is responsible for the colon operator source code.
> But I think that if there is a document for comparing colon operator with
> linspace function,
> that will be helpful.
>
> There is another question: I found that (-20:1:1)*.1 and [-20:1:1]*0.1 have
> different results.
>
>
> The range with [...] is better decimated than the case of (...).

Interesting.  You may have it right there.  And in fact, it appears the 
conversation in the bug report

https://savannah.gnu.org/bugs/?42589

turns to this issue.  John describes how a "Range" class retains the 
description of the range rather than converting it to a matrix class. 
That actually may be a good idea because it could possibly allow 
retaining double precision if needed somehow (e.g., machines where the 
single precision is markedly different from "double").  Whereas 
converting the range to a matrix might make the result single precision 
depending on compiler settings or whatnot.  I'm just speaking in the 
abstract without any example.

So, in your example, the difference is (-20:1:1) is a Range description 
(base, increment, limit), where [-20:1:1] is a matrix_value 
implementation of that range (all the actual 21 values).

Let's try to discern which of these might be the more accurate result:

> > t_1 = (-20:1:1)*.1;
> > t_2 = [-20:1:1]*.1;
> > t_f = [-2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 ...
       -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1];
> > t_d = double(zeros(1,21));
> > t_d(1) = -2; t_d(2) = -1.9; t_d(3) = -1.8; t_d(4) = -1.7; t_d(5) = -1.6;
> > t_d(6) = -1.5; t_d(7) = -1.4; t_d(8) = -1.3; t_d(9) = -1.2; t_d(10) = -1.1;
> > t_d(11) = -1; t_d(12) = -0.9; t_d(13) = -0.8; t_d(14) = -0.7; t_d(15) = -0.6;
> > t_d(16) = -0.5; t_d(17) = -0.4; t_d(18) = -0.3; t_d(19) = -0.2; t_d(20) = -0.1;
> > t_d(21) = 0; t_d(22) = 0.1;
> > max(abs(t_d - t_f))
ans = 0
> > norm(abs(t_d - t_f))
ans = 0
> > max(abs(t_d - t_1))
ans =    2.2204e-16
> > norm(abs(t_d - t_1))
ans =    4.3088e-16
> > max(abs(t_d - t_2))
ans =    2.2204e-16
> > norm(abs(t_d - t_2))
ans =    4.7429e-16
> > max(abs(t_1 - t_2))
ans =    2.2204e-16
> > [abs(t_d-t_1)' abs(t_d-t_2)']
ans =

   0.0000e+00   0.0000e+00
   0.0000e+00   2.2204e-16
   0.0000e+00   0.0000e+00
   0.0000e+00   2.2204e-16
   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00
   0.0000e+00   2.2204e-16
   2.2204e-16   0.0000e+00
   0.0000e+00   2.2204e-16
   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00
   1.1102e-16   0.0000e+00
   2.2204e-16   0.0000e+00
   0.0000e+00   1.1102e-16
   1.1102e-16   1.1102e-16
   0.0000e+00   0.0000e+00
   1.1102e-16   0.0000e+00
   1.6653e-16   5.5511e-17
   5.5511e-17   0.0000e+00
   1.3878e-16   0.0000e+00
   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00

From the norm(abs(t_d - t_f)) = 0 result, I'd say my machine/compiler 
is using double precision by default.  Also, there is a difference 
between (A:1:B)*d and [A:1:B]*d but no implementation seems better than 
the other.  If anything (A:1:B)*d is marginally better than [A:1:B]*d, 
but statistically speaking the norm difference of 4.3405e-17 is probably 
not conclusive.  However, if I were to guess, I'd say (A:1:B)*d retains 
accuracy better because it doesn't cast to a IEEE float representation 
and then multiply.  The computations may be all done with more bits in 
the ALU.  But we're talking the very last bit here, I believe.


> More interestingly, if I transpose the range array  such as (-20:1:1)'*.1,
> then I have the same result in [-20:1:1]'*.1.

The transpose likely forces the Range description into a matrix_value 
implementation.


> Is this a bug? Did someone say that there is a patch for this?
> I am useing Ocatve 3.8.2 on cygwin.

OK, Cygwin.  So there's where there could be some differences in 
compiler settings that make your computations single precision and 
therefore showing larger errors in the arithmetic than what I'm seeing. 
 Try the script lines above on your machine and see what you get.  You 
probably should have the latest code with the Range patch applied.  In 
any case, though, if you can identify a single/double precision issue, 
it might suggest that the Octave code isn't retaining precision...or it 
could simply be that on your system Matlab is using double precision at 
compilation whereas Octave is not.

Dan