[Octave-bug-tracker] [bug #64361] sinc() return NaN on sinc(-1:0.20:1-0.20) expression

[Markus Mützel]

Follow-up Comment #23, bug #64361 (project octave):

I can reproduce on Windows with "-O2 -march=x86-64-v3". That's with GCC 13.1.0
on an AMD Ryzen 7 5800X if that should matter.
I'm seeing the precision issue when I compile like that. The range and the
array have the same value though:

>> format longg
>> x = (-1:0.2:1-0.2)
x =

 Columns 1 through 6:

                      -1                    -0.8                    -0.6      
             -0.4                    -0.2   5.551115123125783e-17

 Columns 7 through 10:

      0.2000000000000001      0.4000000000000001      0.6000000000000001      
              0.8

>> x = [x]
x =

 Columns 1 through 6:

                      -1                    -0.8                    -0.6      
             -0.4                    -0.2   5.551115123125783e-17

 Columns 7 through 10:

      0.2000000000000001      0.4000000000000001      0.6000000000000001      
              0.8

>> x - [x]
ans =

 Columns 1 through 6:

                       0                       0                       0      
                0                       0                       0

 Columns 7 through 10:

                       0                       0                       0      
                0



The following might be interesting:

>> x = -1:0.2:1-0.2;
>> sin(x)./x
ans =

   0.8415   0.8967   0.9411   0.9735   0.9933   1.0000   0.9933   0.9735  
0.9411   0.8967

>> sin(x(:))./x(:)
ans =

   0.8415
   0.8967
   0.9411
   0.9735
   0.9933
   1.0000
   0.9933
   0.9735
   0.9411
   0.8967

>> sin([x])./[x]
ans =

   0.8415   0.8967   0.9411   0.9735   0.9933   1.0000   0.9933   0.9735  
0.9411   0.8967

>> sin(x(1:10))./x(1:10)
ans =

   0.8415   0.8967   0.9411   0.9735   0.9933      NaN   0.9933   0.9735  
0.9411   0.8967

>> x - x(1:10)
ans =

 Columns 1 through 7:

            0            0            0  -5.5511e-17            0   5.5511e-17
 -1.1102e-16

 Columns 8 through 10:

  -5.5511e-17            0            0


So, the result differs depending on whether I'm using indices or not. But
`x(:)` or `[x]` is fine here...

I agree with you that this very much smells like the compiler using different
instructions for something that we think should be equivalent.

Should we really mess with compiler optimizations and redesign our code
because of that? Or should we just accept that using limited precision
floating point numbers is tricky sometimes?
In the latter case, the patch from comment #3 could deal with that in the sinc
function.



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?64361>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/