Re: A problem with range objects and floating point numbers

[Daniel J Sebald]

On 09/28/2014 10:20 PM, s.jo wrote:
> Daniel Sebald wrote

[snip]

> From the result of the script, I can determine that pi multiplied by decimal
> points such 1.0*pi, 2.0*pi, ... 1000.0*pi are well generated from range
> operations.
> As a result, both (A:1:B)*d and [A:1:B]*d are poor than linspace in Octave.
> However, in Matlab, all results from range type, matrix-range type and
> linspace pass the 'abs(sin(k*pi))<eps' test for large k.
> See the complete script and their results (comparing Octave and Matlab) as
> follows:
>
> == sin test script ==
> version
> % generate ranges :
> %  1. range object, 2. matrix object, 3. linspace, 4. perfect decimation
> N=1000;
> range_test1=(-N:0.1:N);
> range_test2=[-N:0.1:N];
> range_test3=linspace(-N,N,20*N+1);


Nice coding for this line:

> range_test4=(str2num(sprintf('%.12f ',-N:0.1:N)));

but I think if you change this to '%.1f ', it pretty much retains only 
the 1/10 fractional digit, which is all we really want.  If there is a 
nonzero digit at 10^-12 or something and it gets converted back to a 
number, then that is not a correct benchmark.


> % count sin(k*pi)==0 : 0 or 1 expected due to floating point representation
> error
> sum(sin(range_test1*pi)==0)
> sum(sin(range_test2*pi)==0)
> sum(sin(range_test3*pi)==0)
> sum(sin(range_test4*pi)==0)

Could you add a couple more cases here?

range_test5=(-N/0.1:1:N/0.1)*0.1;
range_test6=[-N/0.1:1:N/0.1]*0.1;

> % count sin(k*pi)<eps : what we want is 2*N+1 in case of the best decimation
> sum(abs(sin(range_test1*pi))<eps(range_test1*pi))
> sum(abs(sin(range_test2*pi))<eps(range_test2*pi))
> sum(abs(sin(range_test3*pi))<eps(range_test3*pi))
> sum(abs(sin(range_test4*pi))<eps(range_test4*pi))

and here:

sum(abs(sin(range_test5*pi))<eps(range_test5*pi))
sum(abs(sin(range_test6*pi))<eps(range_test6*pi))


> == Octave 3.8.2 result ==
> + sum (abs (sin (range_test1 * pi))<  eps (range_test1 * pi))
> ans =       1055
> + sum (abs (sin (range_test2 * pi))<  eps (range_test2 * pi))
> ans =       1168

On my system, the above case is showing 2001, so it is only the first 
test1 case out of tolerance for me.


> + sum (abs (sin (range_test3 * pi))<  eps (range_test3 * pi))
> ans =       2001
> + sum (abs (sin (range_test4 * pi))<  eps (range_test4 * pi))
> ans =       2001


OK, on my system I printed out the matrix_value whenever the conversion 
was done from range to Matrix.  Here are my results, which should start 
to illustrate where the differences are:

octave:20> N=1000;
octave:21> range_test1=(-N:0.1:N);
octave:22> range_test2=[-N:0.1:N];
matrix_value, inc 0.1
octave:23> range_test3=linspace(-N,N,20*N+1);
octave:24> %range_test4=(str2num(sprintf('%.12f ',-N:0.1:N)));
octave:24> range_test4=(str2num(sprintf('%.1f ',-N:0.1:N)));
octave:25> range_test5=(-N/0.1:1:N/0.1)*0.1;
octave:26> range_test6=[-N/0.1:1:N/0.1]*0.1;
matrix_value, inc 1
octave:27> sum(abs(sin(range_test1*pi))<eps(range_test1*pi))
matrix_value, inc 0.314159
matrix_value, inc 0.314159
ans =  1034
octave:28> sum(abs(sin(range_test2*pi))<eps(range_test2*pi))
ans =  2001
octave:29> sum(abs(sin(range_test3*pi))<eps(range_test3*pi))
ans =  2001
octave:30> sum(abs(sin(range_test4*pi))<eps(range_test4*pi))
ans =  2001
octave:31> sum(abs(sin(range_test5*pi))<eps(range_test5*pi))
matrix_value, inc 0.314159
matrix_value, inc 0.314159
ans =  1034
octave:32> sum(abs(sin(range_test6*pi))<eps(range_test6*pi))
ans =  2001

Comments:

1) Note how the range class, i.e., (A:s:B), multiplication does not 
expand the range, it simply scales the constant s.  So (-N:0.1:N)*pi is 
actually changing the range to, I believe, (-N*pi:pi:N*pi).

2) Similarly, the 5th case in which (-N/0.1:1:N/0.1)*0.1 is exactly the 
same as (-N:0.1:N), such that (-N/0.1:1:N/0.1)*0.1 is also treated as 
(-N*pi:pi:N*pi).  Notice for test1 and test5 that the increment is 
0.1*pi = 0.314159 at the point a range-to-matrix_value conversion is done.

3) The test1 and test5 don't pass the test because this formula

          for (octave_idx_type i = 1; i < rng_nelem; i++)
            cache(i) = b + i * increment;

behaves particular to the numerical characterists of 'b' and of 
'increment'.  If b is an integer, say b=-1000, then that addition step 
may not be so susceptible to numerical effects.  However, if b=-1000pi, 
i.e., a fractional number, then the addition might be less accurate.

4) The same can probably be said for increment.  If that can be made to 
be strictly integer (i.e., fits in the mantissa of the float), then the 
multiplication is more accurate.  But really, I think the addition is 
the more troublesome of the two.  So this sort of comes back to 
representing the range slightly different as:

Given [A:s:B], if

  1) A is integer, and
  2) A/s is integer

then [A:s:B] can be implemented as

  [A/s:1:floor(B/s)] * s

where the range operation is integer-based.

5) I notice that in linspace() algorithms for the various class objects, 
the algorithm is similarly

  for (octave_idx_type i = 1; i < n-1; i++)
    retval(i) = x1 + i*delta;

There is a slightly different alternative for this.  For example,

retval(i) = (x1*(N-i) + x2*i)/N;

for i equal 0 to N, where N is the number of intervals.  Maybe that 
behaves better--don't know.  I will analyze the difference when I get 
the chance.

6) Jo, you are seeing test2 fail.  Could you please try test6 on your 
system:

octave:26> range_test6=[-N/0.1:1:N/0.1]*0.1;
matrix_value, inc 1
octave:32> sum(abs(sin(range_test6*pi))<eps(range_test6*pi))
ans =  2001

I'm thinking that one might pass on your system because this is the case 
where the range is comprised of integer start and integer increment, as 
described above.

7) Does all this matter?  I guess, simply because we tend to think in 
nice even number intervals and nice decimal fractions.  But for an 
increment which is an arbitrary float, it's probably not that critical.

Dan