Re: A problem with range objects and floating point numbers

[Daniel J Sebald]

On 09/29/2014 04:07 AM, Daniel J Sebald wrote:
> On 09/28/2014 10:20 PM, s.jo wrote:

> 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))

Jo,

Here are a couple more cases for you to try:

NN=20*N+1;
range_test7=zeros(1,NN);
for (i=1:1:NN)
  range_test7(i)=(i-1-N/0.1)*0.1;
endfor
range_test8=zeros(1,NN); A=-N; B=N;
for (i=1:1:NN)
  range_test8(i)=(A*(NN-i) + B*(i-1))/(NN-1);
endfor

...

sum(abs(sin(range_test7*pi))<eps(range_test7*pi))
sum(abs(sin(range_test8*pi))<eps(range_test8*pi))

test7 should be the same as test6.  test7 is the scenario of

  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.

and that is what I'm seeing on my system.

test8 is the linear spacing (interpolation) algorithm I hinted at last 
night, and rather than try it in C code, I figured we could first try it 
as a script.  At first blush it looks promising.  Interestingly, the 
algorithm produces the benchmark test4 exactly.

octave:47> norm(range_test4-range_test8)
ans = 0

Without analysis, this is just a guess, but rounding floating point 
errors for one term may be counteracted by the rounding floating point 
errors for the second term.  I suspected such a thing when I saw the 
errors for the A + i*s approach build as i progressed (for i near the 
start there is zero error and for i near the end almost every term has 
an error).

JWE may want to experiment with the test8 linear-spacing algorithm.  It 
would be a lot of work to replace the dozen or so cases that linspace 
shows up in the C++ code, so would want to make sure it's a true 
phenomenon first.

Dan