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

[John W. Eaton]

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

"I can reproduce this problem with octave compiled with AVX/AVX2"

I tried building with CXXFLAGS="-g -O2 -mavx2" and I'm not able to reproduce
the problem.

Exactly what options are you using to build Octave?

In any case, the funny thing is that at comment #18, the values shown for x =
(-1:0.2:0.2) appear to be correct while the ones for [x] do not.  In the first
case, the range is stored as a special class containing the base, limit,
increment, final value, and number of elements.  The display of those values
has a special case in pr-output.cc and the elements are displayed in this
loop:


          while (col < num_elem)
            {
              octave_idx_type lim = (col + inc < num_elem ? col + inc
                                     : num_elem);

              pr_col_num_header (os, total_width, max_width, lim, col,
                                 extra_indent);

              os << std::setw (extra_indent) << "";

              for (octave_idx_type i = col; i < lim; i++)
                {
                  octave_quit ();

                  double val;
                  if (i == 0)
                    val = base;
                  else
                    val = base + i * increment;

                  if (i == num_elem - 1)
                    val = final_value;

                  os << "  ";

                  pr_float (os, fmt, val);
                }

              col += inc;
            }


For [x], an array of all the values is computed, stored in an NDArray object,
and display is handled by the NDArray display method which just prints the
elements of the array.  Computation of the elements is done in
range<T>::array_value in Range.h:


        // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment).
        retval(0) = m_base;

        if (m_reverse)
          for (octave_idx_type i = 1; i < nel - 1; i++)
            retval.xelem (i) = m_base - i * m_increment;
        else
          for (octave_idx_type i = 1; i < nel - 1; i++)
            retval.xelem (i) = m_base + i * m_increment;

        retval.xelem (nel - 1) = final_value ();


I *am* able to observe values that look more like the reported problem if I
make the following change to the loop in the octave_print_internal function
for range<double> objects:


diff --git a/libinterp/corefcn/pr-output.cc b/libinterp/corefcn/pr-output.cc
--- a/libinterp/corefcn/pr-output.cc
+++ b/libinterp/corefcn/pr-output.cc
@@ -2554,6 +2554,7 @@ octave_print_internal (std::ostream& os,
 
           pr_scale_header (os, fmt.scale_factor ());
 
+          double val = 0.0;
           octave_idx_type col = 0;
           while (col < num_elem)
             {
@@ -2569,11 +2570,10 @@ octave_print_internal (std::ostream& os,
                 {
                   octave_quit ();
 
-                  double val;
                   if (i == 0)
                     val = base;
                   else
-                    val = base + i * increment;
+                    val += increment;
 
                   if (i == num_elem - 1)
                     val = final_value;


Then I see:


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

 Columns 1 through 6:

                      -1                    -0.8     -0.6000000000000001    
-0.4000000000000001     -0.2000000000000001  -5.551115123125783e-17

 Column 7:

                     0.2


I'm compiling with -g -O2, not -mavx2, so maybe your compiler is "optimizing"
the loop in array_value to do the same kind of incrementing instead of
performing the multiplication that we asked for?


    _______________________________________________________

Reply to this item at:

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

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