Re: A problem with range objects and floating point numbers

[Daniel J Sebald]

On 09/29/2014 12:51 PM, John W. Eaton wrote:
> Does it also bother you that the results of the following two
> expressions are not exactly bit for bit identical?
>
> linspace (-20, 1, 22) * .1
> linspace (-2, .1, 22)
>
> This is exactly the sort of thing you are comparing, and with binary
> representation of decimal numbers, I don't see a solution that will give
> you the results you seem to want in all cases. That's just the nature of
> doing binary floating point operations.
>
> I don't think Matlab gives the same results for these expressions
> either, so compatibility isn't an argument for changing Octave's
> behavior here.
>
> I plan to add an option to Octave (enabled by default when using the
> --traditional option) that will cause Octave's range operator to create
> a matrix immediately rather than storing just the base, limit, increment
> (and sometimes number of elements). That way the two expressions
>
> (-20:1:1) * .1
> [-20:1:1] * .1
>
> will produce the same values, because the range will be expanded
> immediately into a matrix in both cases.

I doubt such an option would be used that much.  So, it might be better 
to pick one way and stick with it.  I've no problem with the condensed 
(A:s:B) range class.  Ultimately it may work fine if it can be shown the 
underlying algorithm is numerically "nice".

The question I'm beginning to have is the numerical properties of

  result(i) = base + i * increment

which applies to linspace() as well as range and matrix_value classes. 
That is, the characteristics of the error residuals, in particular when 
that increment in the formula above is fractional.

My suspicion is that

linspace (-20, 1, 22) * .1

produces residuals that may in some sense be better than those produced by

linspace (-2, .1, 22)

Here are some results where I've added the third result which is more 
interpolate in nature:

octave:84> range_test_benchmark = str2num(sprintf('%.1f ',linspace(-2, 
.1, 22)));
octave:85> range_test10 = linspace (-20, 1, 22) * .1;
octave:86> norm(range_test10 - range_test_benchmark)
ans =    4.7429e-16
octave:87> range_test11 = linspace (-2, .1, 22);
octave:88> norm(range_test10 - range_test_benchmark)
ans =    4.7429e-16
octave:89> A = -2;
octave:90> B = .1;
octave:91> N = 22;
octave:92> for (i=1:1:N)
>   range_test12(i)=(A*(N-i) + B*(i-1))/(N-1);
> endfor
octave:93> norm(range_test12 - range_test_benchmark)
ans =    3.3422e-16
octave:94> [range_test10' range_test11' range_test12'] - 
range_test_benchmark'*[1 1 1]
ans =

   0.0000e+00   0.0000e+00   0.0000e+00
  -2.2204e-16   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   2.2204e-16
  -2.2204e-16   0.0000e+00  -2.2204e-16
   0.0000e+00   0.0000e+00   0.0000e+00
   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
  -2.2204e-16   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   1.1102e-16   1.1102e-16
   0.0000e+00   2.2204e-16   0.0000e+00
  -1.1102e-16   0.0000e+00   0.0000e+00
  -1.1102e-16   1.1102e-16   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   1.1102e-16   0.0000e+00
  -5.5511e-17   1.6653e-16   0.0000e+00
   0.0000e+00   5.5511e-17   0.0000e+00
   0.0000e+00   1.3878e-16   2.7756e-17
   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00

So the interpolating algorithm (third column) shows a smaller norm, but 
more than that we can see that the middle column (the result(i) = base + 
i*increment algorithm where increment is fractional 0.1) has a sort of 
skewed residual where it is zero at the base (A) end of the range and 
non-zero near the limit (B) end of the range.  Whether these small 
discrepancies are that important is up to the user, but someone with a 
numerical background might say that somehow the residuals of the middle 
column aren't so nice.

Dan