Re: A problem with range objects and floating point numbers

[s.jo]

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))&lt;eps(range_test1*pi))
&gt; sum(abs(sin(range_test2*pi))&lt;eps(range_test2*pi))
&gt; sum(abs(sin(range_test3*pi))&lt;eps(range_test3*pi))
&gt; sum(abs(sin(range_test4*pi))&lt;eps(range_test4*pi))

and here:

sum(abs(sin(range_test5*pi))&lt;eps(range_test5*pi))
sum(abs(sin(range_test6*pi))&lt;eps(range_test6*pi))


&gt; == 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))&lt;eps(range_test1*pi))
matrix_value, inc 0.314159
matrix_value, inc 0.314159
ans =  1034
octave:28&gt; sum(abs(sin(range_test2*pi))&lt;eps(range_test2*pi))
ans =  2001
octave:29&gt; sum(abs(sin(range_test3*pi))&lt;eps(range_test3*pi))
ans =  2001
octave:30&gt; sum(abs(sin(range_test4*pi))&lt;eps(range_test4*pi))
ans =  2001
octave:31&gt; sum(abs(sin(range_test5*pi))&lt;eps(range_test5*pi))
matrix_value, inc 0.314159
matrix_value, inc 0.314159
ans =  1034
octave:32&gt; sum(abs(sin(range_test6*pi))&lt;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 &lt; 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 &lt; 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&gt; 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
&lt;/quote>

Dan said,
> > range_test4=3D(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.
>
>

Here I want to find the largest digits for the perfect decimation, 
say '%.12f'. Roughly speaking, the number of decimation digits is 
not much related to step size.
You can try and see range_test4 with any number less than 12, and 
more than 13, such as '%.11f' and '%.13f'
With less than 12, there is no errors and no non-zero digit from 
10^-12 up-to 0.1. Why?
This is the difference between round-off and decimation.
It depends on from which base system number to which base system
we convert it.
In our example, we decimate base-2 numbers (into base-10 number).
For example, if we take the decimal (or base-10) range [... -0.2 -0.1 0 0.1
0.2 ....],
we may have base-2 representation equivalent to [... -0.199999 -0.1 0
0.111111 0.2 ...].
The floating-point numbers with errors from decimal number may look like
repeating decimals. http://en.wikipedia.org/wiki/Repeating_decimal
Thus, the difference between them can be bounded by a small number such
as+/-0.00...001 or +/-10^12 or so.
The exact number of decimation digits depends on step-size and epsilons
(floating-point space). Epsilon is also a function of numbers.
So I claim that we can choose '12'  as a decimation digit number of 
the base-2 double precision floating-point numbers for decimal numbers.
I believe that Matlab's colon operator employs such 'blind decimation' by
choosing a specific decimal digit.
How? It is not difficult. If we scale the step-size by a large number of 
multiple of 10, say 10^12, and take it back to the original scale,
we have the similar effect as '%.12f'. decimation.
Again, The range (A:s:B) can be decimated by
(A*10^12:s*10^12:B*10^12)/10^12,
whatever step-size s is. The scale factor 10^12 will be a design parameter.
Such multiplication and division with a large number causes a truncation
around a specific decimal position.

Dan expanded test cases up-to 6:
In my note pc (Octave 3.8.2/cygwin/intel-i5-M450), I have ...

+ N = 1000;
+ range_test1 = (-N:0.1:N);
+ range_test2 = [-N:0.1:N];
+ range_test3 = linspace (-N, N, 20 * N + 1);
+ range_test4 = (str2num (sprintf ('%.12f ', -N:0.1:N)));
+ range_test5 = (-N / 0.1:1:N / 0.1) * 0.1;
+ range_test6 = [-N / 0.1:1:N / 0.1] * 0.1;

+ sum (abs (sin (range_test1 * pi)) < eps (range_test1 * pi))
ans =3D       1055
+ sum (abs (sin (range_test2 * pi)) < eps (range_test2 * pi))
ans =3D       1168
+ sum (abs (sin (range_test3 * pi)) < eps (range_test3 * pi))
ans =3D       2001
+ sum (abs (sin (range_test4 * pi)) < eps (range_test4 * pi))
ans =3D       2001
+ sum (abs (sin (range_test5 * pi)) < eps (range_test5 * pi))
ans =3D       1055
+ sum (abs (sin (range_test6 * pi)) < eps (range_test6 * pi))
ans =3D       2001


I happened to think my system is based on a mobile processor. It may cause
test 2 to fail.

Replying to Dan's comments, I agree that Octave's colon operation may not
pass 
the sin(k*pi)<eps test. (tests 1,2,5)
It is natural because Octave follows IEEE floating-point operation rules.

Now I understand that the calling orders of 'range object' may be different
for each case.
That causes many variants: some case has integer start number and/or integer
step-size.
This can explain the difference from Matlab and each others. 
It is also natural for range object with floating-point numbers.

Last thing I want to point out is that Matlab's colon operation uses some
unknown decimation 
process. It can be verified by comparing the raw bit IEEE format between the
range operator's 
result and the arithmetic incremental routine. 
Matlab do not follow the IEEE floating-point operation rule in colon
operator; 
it looks close to the linspace function.

Matlab does not explain why and how they do this on the colon operation with
decimation.
I think that people using colon operator prefer to deal with decimal numbers
such as 
10, 5.2, [0 0.2 0.4 0.6], and so on. They use them in a short command line
to see some trivial cases.

The range operation with floating-point numbers happen to be directly coded
in user application.
Speaking the numerical algorithms, they are usually written by using 4
arithmetic operators and
some primitive functions. In such algorithm codes, it is hardly to see the
colon operators.

For the compatibility with Matlab, Octave can use the decimated colon
operator 
for people using decimal numbers.
Without loss of floating-point standard, Octave can apply the decimation 
only when the starting and ending of range are all integer. 

People always try to compare Octave and Matlab to get more confidence.

-- jo













--
View this message in context: 
http://octave.1599824.n4.nabble.com/A-problem-with-range-objects-and-floating-point-numbers-tp4664900p4666772.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.