Re: A problem with range objects and floating point numbers

[s.jo]

Daniel Sebald wrote
> [snip]
> 
> OK.  I created a changeset that does basically that:
> 
> https://savannah.gnu.org/patch/index.php?8548
> 
> The change uses the existing algorithm for integer elements because
> 
> output(i) = a + inc * i
> 
> has exact results when inc and a are integers.
> 
> However, for other cases the routine
> 
> output(i) = (a (N-i) + b (i)) / N-i * scale
> 
> is used because it has better floating point behavior (several Internet 
> discussion boards cover this).
> 
> John, this is actually a pretty straightforward changeset.  I think more 
> than anything it is a case of organizing in such a way that all 
> variations of limit/increment as far as int/non-int can be addressed.
> 
> Dan

Dear Dan,

I have looked into the changeset
https://savannah.gnu.org/patch/index.php?8548 
and your testscripts there over a couple of weeks.

I am happy to see the changeset, which is the best for the integer-bounded 
decimal ranges such as (-1000:0.1,1000). 
The range implies the equally spaced row vector.

I am not good at hg and git system. I could not compile it from the patches.
I should wait for next release. Will this patch is applied to 3.8.3 or so?

Well, based on your testscprits revised in the changset page, 
I was trying to discuss the decimal range computation in IEEE base-2 
floating-point format.
I applied your integerizing range alorithm to the non-integer-bounded 
decimal ranges such as (-1000.1:0.1,1000.1). But I failed to generalize it.
I found that there is no perfect solution for decimal range computation
without using decimal (base-10) floating-point arithmetic library 
or architecture.
But I believe that the new integerizing range algorithm is the best solution
with the integer-bounded decimal ranges, which are frequently used in 
many test script or applications.

In the following testscript, I summarized 8 algorithms, 4 range examples,
and 3 test methods applied to both Octave 3.8.2 and Matlab 2012b:

8 algorithms (rows in result matrix):
  A1. colon operator (colon.m or built-in function)
  A2. explicit incremental operation (colon_incr.m)
  A3. explicit incremental operation scaled for integer step
(colon_incr_int.m)
  A4. colon operator with explicit decimation (colon_dec.m)
  A5. linspace (built-in function)
  A6. explicit interpolation (linspace_intp.m)
  A7. explicit interpolation scaled for integer step (linspace_intp_int.m)
  A8. linspace operator with explicit decimation (linspace_dec.m)

4 range examples (columns in result matrix):
  R1. (-1000:0.1:1000)
  R2. (-1000.1:0.1:1000.1)
  R3. (-2:0.1:0.1)
  R4. (-0.3:0.1:0.3)

3 test methods (blocks below):
  T1. Count sin(range*pi)<eps (sin_eps_result)
  T2. Measure norm to see how close to the IEEE float-point equivalent
(norm_result)
  T3. Count mismatch numbers with the IEEE float-point equivalent
(count_result)

%% test result in octave 3.8.2 %%
cell2mat(sin_eps_result)
ans =

        1055        1920           3           1
        2001        1415           3           0
        2001        2001           3           0
        2001        2001           3           1
        2001        2001           2           1
        2001        1313           3           1
        2001        1553           3           1
        2001        2001           3           1

cell2mat(norm_result)
ans =

  1.0004e-11  7.9346e-12  3.8858e-16  6.2063e-17
  1.1963e-11  9.1503e-12  4.3885e-16  1.3311e-16
  7.2125e-12  7.2125e-12  4.7429e-16  1.2117e-16
           0           0           0           0
   1.711e-16   5.386e-12  3.1465e-17  4.3885e-17
  8.4071e-12  1.1133e-11  4.2093e-16  1.9626e-17
           0  9.8931e-12  3.3422e-16  4.3885e-17
           0           0           0           0

cell2mat(count_result)
ans =

       14002       11399          10           5
       10821       10208           9           5
        7192        7192           7           6
           0           0           0           0
           5        3976           3           4
       12096       11444           8           2
           0       14016           4           4
           0           0           0           0

%% test result in matlab 2012b %%
cell2mat(sin_eps_result)
ans =

        2001        2001           3           1
        2001        1415           3           0
        2001        2001           3           0
        2001        2001           3           1
        2001        1143           3           1
        2001        1313           3           1
        2001        1553           3           1
        2001        2001           3           1


cell2mat(norm_result)
ans =

   6.3738e-12   4.5688e-12   3.0022e-16   5.5511e-17
    1.196e-11   9.1474e-12   4.3885e-16   1.3311e-16
   7.2107e-12   7.2107e-12   4.7429e-16   1.2117e-16
            0            0            0            0
   7.5997e-12   1.4415e-11    5.446e-16   2.7756e-17
   8.4071e-12   1.1133e-11   4.2093e-16   1.9626e-17
            0   9.8931e-12   3.3422e-16   4.3885e-17
            0            0            0            0


cell2mat(count_result)
ans =

        7592        6462           6           4
       10820       10208           9           5
        7190        7190           7           6
           0           0           0           0
        9024       11538           9           1
       12096       11444           8           2
           0       14016           4           4
           0           0           0           0

----

In discussion, the decimal range with known decimal digits can be
computed perfectly through decimal string conversion, and converted
into base-2 floating-point format.
Integerizing range algorithms (A3 and A7) outperform both in Octave and
Matlab,
close to the explicit decimation in the integer-bounded ranges.
But, for non-integer-bounded ranges, Matlab's built-in colon operator 
outperforms in the non-integer-bounded ranges, besides the explicit
decimation.

Note that the explicit decimation with known decimal digits is the best
solution. But, without using the knowledge of the number of decimal digits,
the integerizing range algorithms (A3 and A7) seem to be good choices.
Matlab's colon operator still seems to have a treatment even in the 
non-integer-bounded ranges.

--jo

ps. the test codes are listed as below. 
(sorry for the long scripts, 
and omitted the subfunctions for readability)

----

%% rangetestscripts_review.m

% example list of test functions:
% N=1000
% NN=20*N+1;
% assume step is 0.1.
% 1. colon : -N:0.1:N
% 2. colon_incr : for ii=1:NN;range_test(ii)=(-N+(ii-1)*0.1);end
% 3. colon_incr_int : for ii=1:NN;range_test(ii)=(-N/0.1+(ii-1))*0.1;end
% 4. colon_dec : (str2num(sprintf('%.1f ',-N:0.1:N)));
% 5. linspace : linspace(-N,N,20*N+1)
% 6. linspace_intp : for
ii=1:NN;range_test(ii)=(-N*(NN-ii)/(NN-1)+N*(ii-1)/(NN-1));end
% 7. linspace_intp_int : for
ii=1:NN;range_test(ii)=(-N*(NN-ii)+N*(ii-1))/(NN-1);end
% 8. linspace_dec : (str2num(sprintf('%.1f ',linspace(-N,N,20*N+1))));

%% Generate ranges:
% 1. Range_bipolar_test of (-1000:0.1:1000)
% 1.1. colon operator (colon.m)
% 1.2. explicit incremental operation (colon_incr.m)
% 1.3. explicit incremental operation scaled for integer step
(colon_incr_int.m)
% 1.4. colon operator with explicit decimation (colon_dec.m)
% 1.5. linspace 
% 1.6. explicit interpolation (linspace_intp.m)
% 1.7. explicit interpolation scaled for integer step (linspace_intp_int.m)
% 1.8. linspace operator with explicit decimation (linspace_dec.m)
% 2. Range_bipolar_test of (-1000.1:0.1:1000.1)
% 2.x. subtests
% 3. Range_bipolar_test of (-2:0.1:0.1)
% 3.x. subtests
% 4. Range_bipolar_test of (-2:0.1:2)
% 4.x. subtests

range_test={};
% 1. Range_bipolar_test of (-1000:0.1:1000)
N=1000;
NN=round(2*N/0.1)+1;
range_test{1,1}=colon(-N,0.1,N);
range_test{2,1}=colon_incr(-N,0.1,N);
range_test{3,1}=colon_incr_int(-N,0.1,N);
range_test{4,1}=colon_dec(-N,0.1,N,1);
range_test{5,1}=linspace(-N,N,NN);
range_test{6,1}=linspace_intp(-N,N,NN);
range_test{7,1}=linspace_intp_int(-N,N,NN);
range_test{8,1}=linspace_dec(-N,N,NN,1);
% 2. Range_bipolar_test of (-1000.1:0.1:1000.1)
N=1000.1;
NN=round(2*N/0.1)+1;
range_test{1,2}=colon(-N,0.1,N);
range_test{2,2}=colon_incr(-N,0.1,N);
range_test{3,2}=colon_incr_int(-N,0.1,N);
range_test{4,2}=colon_dec(-N,0.1,N,1);
range_test{5,2}=linspace(-N,N,NN);
range_test{6,2}=linspace_intp(-N,N,NN);
range_test{7,2}=linspace_intp_int(-N,N,NN);
range_test{8,2}=linspace_dec(-N,N,NN,1);
% 3. Range_bipolar_test of (-2:0.1:0.1)
A = -2;
B = .1;
S = .1;
N = 22;
decimal_digits=1;
range_test{1,3}=colon(A,S,B);
range_test{2,3}=colon_incr(A,S,B);
range_test{3,3}=colon_incr_int(A,S,B);
range_test{4,3}=colon_dec(A,S,B,decimal_digits);
range_test{5,3}=linspace(A,B,N);
range_test{6,3}=linspace_intp(A,B,N);
range_test{7,3}=linspace_intp_int(A,B,N);
range_test{8,3}=linspace_dec(A,B,N,decimal_digits);
% 4. Range_bipolar_test of (-2:0.1:2)
A = -2;
B = 2;
S = .1;
N = round((B-A)/S)+1;
decimal_digits=1;
range_test{1,4}=colon(A,S,B);
range_test{2,4}=colon_incr(A,S,B);
range_test{3,4}=colon_incr_int(A,S,B);
range_test{4,4}=colon_dec(A,S,B,decimal_digits);
range_test{5,4}=linspace(A,B,N);
range_test{6,4}=linspace_intp(A,B,N);
range_test{7,4}=linspace_intp_int(A,B,N);
range_test{8,4}=linspace_dec(A,B,N,decimal_digits);


%% result
ss = size(range_test);
num_subtests = ss(1);
num_cases = ss(2);
sin_eps_result={};
norm_result={};
count_result={};
for ii=1:num_subtests
for jj=1:num_cases
        % Count sin(k*pi)<eps 

sin_eps_result{ii,jj}=sum(abs(sin(range_test{ii,jj}*pi))<eps(range_test{ii,jj}*pi));
        %
        range_test_benchmark = range_test{4,jj};
        % Measure how close to the IEEE float-point equivalent
        norm_result{ii,jj}=norm(range_test{ii,jj} - range_test_benchmark);
        % Count mismatch // Measure how close to the IEEE float equivalent
        count_result{ii,jj}=sum(range_test{ii,jj} ~= range_test_benchmark);
end
end

% display
cell2mat(sin_eps_result)
cell2mat(norm_result)
cell2mat(count_result)


----




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