bugGNU Octave - Bugs: bug #61300, integer range might exceed upper...

 
 

bug #61300: integer range might exceed upper limit

Submitter:  Arun Giridhar <arungiridhar>
Submitted:  Wed 06 Oct 2021 03:13:19 PM UTC
   
 
Category:  Interpreter Severity:  4 - Important
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 02 Dec 2021 05:08:58 PM UTC, comment #44: 

Thanks for testing.

Closing as fixed.

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Dec 2021 03:26:46 PM UTC, comment #43: 

At my end this patch for #61300 looks stable, and so does the effect on nchoosek (pending the patch for bug #61132 being merged with stable).

Arun Giridhar <arungiridhar>
Group Member
Tue 30 Nov 2021 10:20:29 PM UTC, comment #42: 

I think the problems are local to nchoosek and not to Range, becuase nchoosek's heavy use of ranges was done with floating point and not integers, which round differently.

It's not a regression as such. It's more that the fix to Range in this report exposes a previously hidden bug in nchoosek.


    119     ## For a ~25% performance boost, multiply values pairwise so there
    120     ## are fewer elements in do/until loop which is the slow part.
    121     ## Since Odd*Even is guaranteed to be Even, also take out a factor
    122     ## of 2 from numerator and denominator.
    123     if (rem (k, 2))  # k is odd
    124       numer = [(v-k+1:v-(k+1)/2) .* (v-1:-1:v-(k-1)/2) / 2, v];
    125       denom = [(1:k/2) .* (k-1:-1:(k+1)/2) / 2, k];
    126     else             # k is even
    127       numer = (v-k+1:v-k/2) .* (v:-1:v-k/2+1) / 2;
    128       denom = (1:k/2) .* (k:-1:k/2+1) / 2;
    129     endif


Arun Giridhar <arungiridhar>
Group Member
Tue 30 Nov 2021 10:13:34 PM UTC, comment #41: 

I am seeing a regression with nchoosek. It looks like it was because of this change because it affects stable too.

Before this patch, nchoosek only failed for unsigned types as described in bug #61565 but DID work for signed integers. Now that fails too with a length mismatch error:


>> nchoosek (int8 (10), int8 (5))
error: product: nonconformant arguments (op1 is 1x3, op2 is 1x2)
error: called from
    nchoosek at line 125 column 13


We can certainly patch nchoosek (much easier than patching Range) but I wanted to raise it here because it might be a symptom of something deeper.

Arun Giridhar <arungiridhar>
Group Member
Tue 30 Nov 2021 04:06:26 PM UTC, comment #40: 

@Markus: Thanks for making this change.  I made a few small additional changes here:

https://hg.savannah.gnu.org/hgweb/octave/rev/597275db9c7f

John W. Eaton <jwe>
Group administrator
Tue 30 Nov 2021 03:13:39 PM UTC, comment #39: 

I pushed a patch that moves to change to specializations for the octave_int<T> types here:
https://hg.savannah.gnu.org/hgweb/octave/rev/20cefb3b0da6

Marking as ready for test.

Re comment #38: Afaict, only instantiations of the range<T> template for double, float and octave_int<T> are available from .m file code. However, a programmer could instantiate that template with other types in their own code, e.g. for .oct files.

Markus Mützel <mmuetzel>
Group administrator
Mon 29 Nov 2021 08:08:01 PM UTC, comment #38: 

@Markus in comment #36: I think the denominator in Range.h should only say m_increment not m_increment.value().

If I understand this class architecture, the base version will work for any type T, but the integer-specialized ones will call foo.value() while the float version and double version will process NaN etc first. So if we are passing raw C++ integer values (not Octave integer values), which one is called? The base version? I take it that is not possible from an end-user M-file?

Arun Giridhar <arungiridhar>
Group Member
Mon 29 Nov 2021 07:21:14 PM UTC, comment #37: 

Pushed your patch with minor changes to the stable branch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/bbf1293bd255

Markus Mützel <mmuetzel>
Group administrator
Mon 29 Nov 2021 07:04:50 PM UTC, comment #36: 

Attaching a patch that moves the previous changes to specializations for octave_int<T> types.

(file #52386)

Markus Mützel <mmuetzel>
Group administrator
Mon 29 Nov 2021 05:38:41 PM UTC, comment #35: 

Updated primes.m patch attached with the warning about flintmax.

(file #52385)

Arun Giridhar <arungiridhar>
Group Member
Mon 29 Nov 2021 05:34:43 PM UTC, comment #34: 

I think we have to specialize for each octave_int<T> type but we can do it with a


template <typename T>
xinit (octave_int<T>& base, ...)
{
  ...
}


template function to avoid at least some code duplication.

John W. Eaton <jwe>
Group administrator
Mon 29 Nov 2021 05:08:07 PM UTC, comment #33: 

Re comment #31: Maybe it would make sense to add a warning if the input is larger than flintmax. Like you wrote, it is very unlikely that anyone will ever see this warning (unless accidentally). But that additional check wouldn't cost much either, I'd guess.

Re comment #32: Afaict, we are currently instantiating the range<T> template with float, double and octave_int<T> types. There is already a specialization of the init method for float and double. I guess, we can also specialize for all octave_int<T> types.
Is it possible to specialize a template with a template class in C++11? Or would we need to repeat the same specialization for each octave_int<T> type?

Markus Mützel <mmuetzel>
Group administrator
Mon 29 Nov 2021 04:48:36 PM UTC, comment #32: 

By using the T::value() method in the range<T>init() method, you are restricting the generic implementation of this class to only work for types that have value() methods.  That works for octave_int<T> but is there a better solution?  Does the original code work properly for normal integer types (int8_t, uint8_t, etc.)?  Maybe we could specialize the init function for all the octave_int<T> types we care about?

John W. Eaton <jwe>
Group administrator
Mon 29 Nov 2021 04:23:09 PM UTC, comment #31: 

@Markus: Thank you for rearranging the BISTs in range.tst. They are much more compact now.

Re primes.m, since Octave's generation of primes is limited to all primes less than a given input, limiting it to flintmax is not a problem at all since memory and time will preclude outputs that big. If it were a different case, like "generate all primes from A to B, where A and B can be around 1e17" then it would not work but Octave has no plans for such a function (nor Matlab). Typically for high arbitrary ranges, such as (int64 (10)^17 + (1:100)), isprime() is a better function. For more serious applications involving prime numbers in general, I use the external tool primesieve (https://github.com/kimwalisch/primesieve) which uses very complicated sieves with heavy CPU-specific optimization. That is outside the scope of Octave though.

Performance using integer ranges inside primes is about half the speed of using double:

octave:1> n = 1e9; tic; p1 = primes(n); toc, tic; p2 = primes(uint64(n)); toc, assert (isequal(p1,p2))
Elapsed time is 3.78581 seconds.
Elapsed time is 7.02282 seconds.


Performance with the internal cast to double is closer to equal:

octave:2> n = 1e9; tic; p1 = primes(n); toc, tic; p2 = primes(uint64(n)); toc, assert (isequal(p1,p2))
Elapsed time is 3.7974 seconds.
Elapsed time is 3.9996 seconds.


Arun Giridhar <arungiridhar>
Group Member
Mon 29 Nov 2021 03:33:09 PM UTC, comment #30: 

Looking at the changes for `primes.m`: Could the cast to double cause issues with input values larger than flintmax?
Does `primes` even work for numbers as large?

Should we instead try to do the "right thing" when the input is an integer type?

Markus Mützel <mmuetzel>
Group administrator
Mon 29 Nov 2021 03:28:16 PM UTC, comment #29: 

I slightly re-arrainged the BISTs for the integer range bug and pushed them together with the actual fix to the stable branch:
https://hg.savannah.gnu.org/hgweb/octave/rev/5793f0355bfa

Initially I thought this would be just one single bug. Instead this became one report about two different issues. Normally, we try to avoid that. Sorry for mixing this.

I'll keep the status "confirmed" for the original issue in `primes.m`.

Markus Mützel <mmuetzel>
Group administrator
Sun 28 Nov 2021 11:05:49 PM UTC, comment #28: 

I updated the BISTs inside primes.m -- the ones I uploaded an hour ago had a type mismatch between double and integer. Please use the attached version instead.

(file #52376)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 10:46:12 PM UTC, comment #27: 

Here is the patch to range.tst. Added 72 new tests for this bug.


octave:1> test ../test/range.tst
warning: imaginary part of complex colon arguments is ignored
warning: called from
    __test__ at line 2 column 1
    test at line 677 column 11

PASSES 455 out of 455 tests



(file #52375)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 09:56:44 PM UTC, comment #26: 

Here is the patch for primes.m (attached). I also added BISTs to primes.m. In addition to the answer being more accurate, it also boosts the performance by about 2X.

Next I will update BISTs for range.tst. It will take me maybe an hour to test everything first.

(file #52371)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 09:35:06 PM UTC, comment #25: 

Here is the patch for the integer range (attached).

I left it 0 instead of octave_idx_type (0), which seemed like an unnecessary cast but if any compilers complain I can add it at that point.

Don't forget the independent second patch for primes.m -- I'll be sending one in about 15 minutes.

(file #52370)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 08:16:26 PM UTC, comment #24: 

That change looks good to me.

Are compilers intelligent enough to deduce the correct type for the `0` in that expression?

One way to find out would probably be to push that change and check if the CI is fine with it...

Could you please prepare a patch including a commit message?
It would also be good to add some BISTs for that change to make sure it doesn't break again accidentally in the future.

Markus Mützel <mmuetzel>
Group administrator
Sun 28 Nov 2021 07:56:37 PM UTC, comment #23: 

Yes, you are right: the problem in the previous test (comment #21) was the values were too close to intmax. Reducing the values seems to work better:


octave:1> uint64(0) : uint64(1e18) : uint64(1e19)
ans =

 Columns 1 through 6:

                     0   1000000000000000000   2000000000000000000   3000000000000000000   4000000000000000000   5000000000000000000

 Columns 7 through 11:

   6000000000000000000   7000000000000000000   8000000000000000000   9000000000000000000  10000000000000000000

octave:2> int64(0) : int64(1e17) : int64(1e18)
ans =

 Columns 1 through 6:

                    0   100000000000000000   200000000000000000   300000000000000000   400000000000000000   500000000000000000

 Columns 7 through 11:

   600000000000000000   700000000000000000   800000000000000000   900000000000000000  1000000000000000000


Converged patch:


diff -r 3c4639e42855 liboctave/array/Range.h
--- a/liboctave/array/Range.h   Fri Nov 26 16:22:19 2021 -0800
+++ b/liboctave/array/Range.h   Sun Nov 28 14:54:56 2021 -0500
@@ -325,8 +325,9 @@
       m_numel = ((m_increment == T (0)
                   || (m_limit > m_base && m_increment < T (0))
                   || (m_limit < m_base && m_increment > T (0)))
-                 ? T (0)
-                 : (m_limit - m_base + m_increment) / m_increment);
+                 ? 0
+                 : (m_limit.value () - m_base.value () + m_increment.value ())
+                   / m_increment.value ());

       m_final = m_base + (m_numel - 1) * m_increment;
     }


Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 07:49:09 PM UTC, comment #22: 

IIUC, an array can't have more elements than the maximum of `octave_idx_type`. `octave_idx_type` is either `int64_t` or `int32_t` depending on the platform.
The same should probably be valid for ranges. So, the assignment in comment #20 should be save.

IIUC, in your examples in comment #21, the sum `m_limit.value () - m_base.value () + m_increment.value ()` overflows. Would that be different if the result would be cast to `int64_t`?

We might want to detect or avoid that overflow. But that is a different issue imho.

Markus Mützel <mmuetzel>
Group administrator
Sun 28 Nov 2021 07:37:41 PM UTC, comment #21: 

Here is the problem with not specifying "int64" for numer and denom:


octave:2> uint64(0) : uint64(1e18) : intmax("uint64")
ans = [](1x0)

octave:3> int64(0) : int64(1e18) : intmax("int64")
error: invalid range


How to work around this?

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 07:27:54 PM UTC, comment #20: 

Yes, there are specializations of init for float and double. They are in Range.cc and they handle NaN values etc:


  template <>
  void
  range<double>::init (void)
  {
    return xinit (m_base, m_limit, m_increment, m_final, m_numel);
  }

  template <>
  void
  range<float>::init (void)
  {
    return xinit (m_base, m_limit, m_increment, m_final, m_numel);
  }


Not sure of the consequences yet of casting to int64_t or not casting. So far, make check performs correctly, so no regressions at least. Here is a version that avoids numer and denom entirely:


    void init (void)
    {
      m_numel = ((m_increment == T (0)
                  || (m_limit > m_base && m_increment < T (0))
                  || (m_limit < m_base && m_increment > T (0)))
                 ? 0
                 : (m_limit.value () - m_base.value () + m_increment.value ())
                   / m_increment.value ());

      m_final = m_base + (m_numel - 1) * m_increment;
    }


But now I have no idea what the conversion is going to be from foo.value() to octave_idx_type -- is that a safe conversion internally?

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 07:14:23 PM UTC, comment #19: 

Is it ok to cast to a (signed) `int64_t` type? Wouldn't that limit the range for (unsigned) `uint64_t` values?

Do we need that cast at all? Couldn't we just use the value with the type that the `value ()` function returns directly?

Do we already have a specialization of the `init ()` function for `float` and `double`?

Markus Mützel <mmuetzel>
Group administrator
Sun 28 Nov 2021 07:04:02 PM UTC, comment #18: 

Since we no longer need floor() for the calculation, I removed cmath from the included header files. Updated patch attached. Please verify.


diff -r 3c4639e42855 liboctave/array/Range.h
--- a/liboctave/array/Range.h        Fri Nov 26 16:22:19 2021 -0800
+++ b/liboctave/array/Range.h        Sun Nov 28 14:01:47 2021 -0500
@@ -322,11 +322,14 @@

     void init (void)
     {
+      int64_t numer = m_limit.value () - m_base.value () + m_increment.value ();
+      int64_t denom = m_increment.value ();
+
       m_numel = ((m_increment == T (0)
                   || (m_limit > m_base && m_increment < T (0))
                   || (m_limit < m_base && m_increment > T (0)))
-                 ? T (0)
-                 : (m_limit - m_base + m_increment) / m_increment);
+                 ? 0
+                 : (numer / denom));

       m_final = m_base + (m_numel - 1) * m_increment;
     }



(file #52363)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 07:00:38 PM UTC, comment #17: 

I see that using int64_t instead of uint64_t for numer and denom does allow for both descending and ascending ranges.

Make check passes without regressions.

Updated patch attached. Please verify.


(file #52362)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 06:52:58 PM UTC, comment #16: 

OK, using value() works for ascending ranges but fails for descending ranges, giving an empty range instead. How to fix this behavior?

Updated code inside Range.h:

    void init (void)
    {
      uint64_t numer = m_limit.value () - m_base.value () + m_increment.value ();
      uint64_t denom = m_increment.value ();

      m_numel = ((m_increment == T (0)
                  || (m_limit > m_base && m_increment < T (0))
                  || (m_limit < m_base && m_increment > T (0)))
                 ? 0
                 : (numer / denom));

      m_final = m_base + (m_numel - 1) * m_increment;
    }


Wrong result for descending ranges:

>> int64 (100) : int64 (-6) : 0
ans = [](1x0)


Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 06:39:05 PM UTC, comment #15: 

@Markus, thank you! I had not known about the value function of the octave_int type. Let me look into using that.

I had the same question about whether casting to double was the right thing to do because of values above flintmax, but I realized that the number of elements in a range is not going to exceed flintmax for decades to come, and cautiously used the casting. Let me see how to achieve that without casting.

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 06:14:49 PM UTC, comment #14: 

I'm not sure if this is the right way to fix this. IIUC, the whole point of introducing integer type ranges was to correctly support ranges above flintmax.
Would that still be working correctly when we cast to `double`?

Could we instead use the `value ()` function of the octave_int type to get the underlying C++ type value? IIUC, C++ arithmetic would do the right thing in this case.

Markus Mützel <mmuetzel>
Group administrator
Sun 28 Nov 2021 05:52:34 PM UTC, comment #13: 

Improved patch attached for the integer range bug in Range.h. This one gets rid of the if-check and uses the standard technique of converting to double and taking the floor. Since m_numel is of type octave_idx_type, I chose to cast it to that type instead of type T.

All range tests from comment #12 pass with both integers and floating points. Make check passes with no regressions and no new failures:

  PASS                            17028
  FAIL                                0
  XFAIL (reported bug)               27
  SKIP (missing feature)             43
  SKIP (run-time condition)          25


Please verify the behavior and the performance.

This is independent of the integer division bug in primes.m (patch in comment #9).

(file #52361)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 11:24:48 AM UTC, comment #12: 

Regarding the integer range bug outside of primes.m (not the integer division bug inside primes.m), I have isolated one location in Range.h that triggers it. Very ironically, it comes down to integer division again (!) but this time in C++. You can try the attached WIP patch to see what happens.

The place that causes the problem of m_final exceeding m_limit is this statement inside range<T>::init, specifically the division inside the conditional statement:


      m_numel = ((m_increment == T (0)
                  || (m_limit > m_base && m_increment < T (0))
                  || (m_limit < m_base && m_increment > T (0)))
                 ? T (0)
                 : (m_limit - m_base + m_increment) / m_increment);

      m_final = m_base + (m_numel - 1) * m_increment;


When the template class T is one of the Octave integer types, that division causes rounding rather than flooring, and m_numel sometimes becomes one unit bigger than it should be, which causes m_final to exceed its value as well.

I demonstrate this bug by adding extra cout statements and by manually decrementing m_numel and m_final with this check right afterwards. (This code is also in the attached WIP patch).


      if ( (m_increment > T (0) && m_final > m_limit)
           || (m_increment < T (0) && m_final < m_limit) )
        {
          std::cout << "Out of limits, old == ";
          std::cout << m_base << ':' << m_increment << ':' << m_limit
                    << " Final=" << m_final << " Numel=" << m_numel
                    << "\tAdjusting..., new == ";

          m_final -= m_increment;
          m_numel--;

          std::cout << m_base << ':' << m_increment << ':' << m_limit
                    << " Final=" << m_final << " Numel=" << m_numel << '\n';
        }


The output is as follows:

  octave:9> uint64 (0) : uint64 (6) : uint64 (100)
  Out of limits, old == 0:6:100 Final=102 Numel=18        Adjusting..., new == 0:6:100 Final=96 Numel=17
  ans =
     0   6  12  18  24  30  36  42  48  54  60  66  72  78  84  90  96

  octave:10> uint64 (1) : uint64 (6) : uint64 (100)
  Out of limits, old == 1:6:100 Final=103 Numel=18        Adjusting..., new == 1:6:100 Final=97 Numel=17
  ans =
     1   7  13  19  25  31  37  43  49  55  61  67  73  79  85  91  97

  octave:11> uint64 (2) : uint64 (6) : uint64 (100)
  ans =
     2   8  14  20  26  32  38  44  50  56  62  68  74  80  86  92  98

  octave:12> uint64 (3) : uint64 (6) : uint64 (100)
  ans =
     3   9  15  21  27  33  39  45  51  57  63  69  75  81  87  93  99

  octave:13> uint64 (4) : uint64 (6) : uint64 (100)
  ans =
      4   10   16   22   28   34   40   46   52   58   64   70   76   82   88   94  100

  octave:14> uint64 (5) : uint64 (6) : uint64 (100)
  Out of limits, old == 5:6:100 Final=101 Numel=17        Adjusting..., new == 5:6:100 Final=95 Numel=16
  ans =
     5  11  17  23  29  35  41  47  53  59  65  71  77  83  89  95


and for decreasing ranges:


octave:17> int64 (100) : int64 (-6) : 0
Out of limits, old == 100:-6:0 Final=-2 Numel=18        Adjusting..., new == 100:-6:0 Final=4 Numel=17
ans =
  100   94   88   82   76   70   64   58   52   46   40   34   28   22   16   10    4

octave:18> int64 (99) : int64 (-6) : 0
Out of limits, old == 99:-6:0 Final=-3 Numel=18        Adjusting..., new == 99:-6:0 Final=3 Numel=17
ans =
  99  93  87  81  75  69  63  57  51  45  39  33  27  21  15   9   3

octave:19> int64 (98) : int64 (-6) : 0
ans =
  98  92  86  80  74  68  62  56  50  44  38  32  26  20  14   8   2

octave:20> int64 (97) : int64 (-6) : 0
ans =
  97  91  85  79  73  67  61  55  49  43  37  31  25  19  13   7   1

octave:21> int64 (96) : int64 (-6) : 0
ans =
  96  90  84  78  72  66  60  54  48  42  36  30  24  18  12   6   0

octave:22> int64 (95) : int64 (-6) : 0
Out of limits, old == 95:-6:0 Final=-1 Numel=17        Adjusting..., new == 95:-6:0 Final=5 Numel=16
ans =
  95  89  83  77  71  65  59  53  47  41  35  29  23  17  11   5


If you work through the arithmetic for each of the above 12 cases, you will see that all the cases with overruns are the cases where the expression (m_limit - m_base + m_increment) / m_increment is rounded upwards instead of being floored, and the other cases are where the value is rounded downwards, therefore the same as being floored.

That if-check in the WIP patch stops m_final from exceeding m_limit, but that is not a good solution. The real solution would be to find a way to calculate "(m_limit - m_base + m_increment) / m_increment)" while ensuring it is always floored, not rounded. Is there a way to do that?

(file #52360)

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 10:40:09 AM UTC, comment #11: 

@Markus, there are two independent bugs here. The integer division bug inside primes is independent of the integer range bug. The incorrect result of the integer division used to calculate len, lenp and lenm (because of forced rounding) is then used to create integer ranges. Even after the integer range bug is fixed, primes will still return values larger than the input for certain integer inputs like 358 or 1e6.

I added my patch for the division bug to this thread since the original bug report had been about primes. Would it be better to create a new report so it can be tracked separately?

Arun Giridhar <arungiridhar>
Group Member
Sun 28 Nov 2021 07:42:47 AM UTC, comment #10: 

@Arun: Your proposed fix might be good enough for `prime`. But it is not all we need.
Ranges shouldn't exceed their limits. That is a big that should be fixed.
There might be other functions in Octave or in user code that are also affected by this issue.
If I understood what you wrote correctly, we wouldn't need any change in `prime` if integer ranges worked correctly.

Markus Mützel <mmuetzel>
Group administrator
Sun 28 Nov 2021 12:53:51 AM UTC, comment #9: 

All, we DO need to cast n to double inside primes.m for the len calculations, because of rounding behavior that has nothing to do with integer ranges. Try this:


max (primes (uint64 (358)))
max (primes (uint64 (1e6)))


In both cases you will see that the result exceeds the input, which only happens when the input is an integer but not a floating point. That happens because the variables len, lenm and lenp are computed as (n-1)/2, (n-1)/6 or (n+1)/6, and when n is an integer that division causes rounding not flooring, so for certain values it returns an extra prime number.

The simplest fix is to cast n to double when calculating len, lenp, and lenm. Such a diff is as follows and is also attached:


diff -r 0b3d917678a8 scripts/specfun/primes.m
--- a/scripts/specfun/primes.m  Fri Nov 26 20:53:48 2021 -0800
+++ b/scripts/specfun/primes.m  Sat Nov 27 19:41:51 2021 -0500
@@ -73,7 +73,13 @@
   elseif (n < 100e3)
     ## Classical Sieve algorithm
     ## Fast, but memory scales as n/2.
-    len = floor ((n-1)/2);        # length of the sieve
+    len = floor (double (n-1)/2);        # length of the sieve
+    ## The casting to double is required for integer division.
+    ## Without this, primes (uint64 (358)) returns 359 which exceeds n.
+    ## That is because of implicit rounding in the division above:
+    ## For n == uint64(358), len becomes 357/2 which is rounded UP
+    ## to 179, and 179 * 2 + 1 == 359 is prime, which is returned.
+
     sieve = true (1, len);        # assume every odd number is prime
     for i = 1:(sqrt (n)-1)/2      # check up to sqrt (n)
       if (sieve(i))               # if i is prime, eliminate multiples of i
@@ -84,8 +90,12 @@
   else
     ## Sieve algorithm optimized for large n
     ## Memory scales as n/3 or 1/6th less than classical Sieve
-    lenm = floor ((n+1)/6);       # length of the 6n-1 sieve
-    lenp = floor ((n-1)/6);       # length of the 6n+1 sieve
+    lenm = floor (double (n+1)/6);       # length of the 6n-1 sieve
+    lenp = floor (double (n-1)/6);       # length of the 6n+1 sieve
+    ## The casting to double is required for integer division.
+    ## Without it, primes (uint64 (1e6)) returns 1000003 which exceeds n.
+    ## That is because of implicit rounding in the division above.
+
     sievem = true (1, lenm);      # assume every number of form 6n-1 is prime
     sievep = true (1, lenp);      # assume every number of form 6n+1 is prime



(file #52357)

Arun Giridhar <arungiridhar>
Group Member
Wed 06 Oct 2021 07:10:44 PM UTC, comment #8: 

Thanks for this report.  I'll take another look at integer ranges soon, definitely before version 7 is released.

John W. Eaton <jwe>
Group administrator
Wed 06 Oct 2021 06:16:11 PM UTC, comment #7: 

It can also underflow, which can cause negative indices like here:


>> double(10) : double(-6) : int32(1)
ans =
   10    4   -2

>> double(10) : int32(-6) : int32(1)
ans =
   10    4   -2

>> int32(10) : int32(-6) : int32(1)
ans =
   10    4   -2


Arun Giridhar <arungiridhar>
Group Member
Wed 06 Oct 2021 05:42:09 PM UTC, comment #6: 

I can confirm on Windows with Octave 7.

Re-titling to the better describe the underlying issue and raising severity.

Markus Mützel <mmuetzel>
Group administrator
Wed 06 Oct 2021 05:37:33 PM UTC, comment #5: 

Note: depending on input value n, sometimes the output of primes(n) can exceed n because of the range violation:


>> n = double(1.5e9); p = primes(n); p(end-5:end)
ans =
     1499999909     1499999927     1499999929     1499999933     1499999951     1499999957

>> n = uint64(1.5e9); p = primes(n); p(end-5:end)
ans =
  1499999927  1499999929  1499999933  1499999951  1499999957  1500000001


That last value 1500000001 is more than n.

Arun Giridhar <arungiridhar>
Group Member
Wed 06 Oct 2021 05:17:35 PM UTC, comment #4: 

Yes, this is a regression, on version 6.0.1 (via octave-online.net) I have


format long;n = uint64(1.5e9), lenp = floor((n+1)/6), a=(8:7:lenp)(end-4:end)
n = 1500000000
lenp = 250000000
a =

   249999968   249999975   249999982   249999989   249999996


Michael Leitner <mleitner>
Wed 06 Oct 2021 05:10:10 PM UTC, comment #3: 

I know now exactly why this slowdown happens. It's because for some reason the integer range is violated, causing the upper bound to exceed lenp, then the array is resized and rewritten with all that overhead. For i = 1, the range inside sievep becomes 8:7:lenp. Look at this experiment for that range:


>> n = double(1.5e9), lenp = floor((n+1)/6), (8:7:lenp)(end-4:end)
n =       1500000000
lenp =        250000000
ans =
               249999968               249999975               249999982               249999989               249999996

>> n = uint64(1.5e9), lenp = floor((n+1)/6), (8:7:lenp)(end-4:end)
n = 1500000000
lenp = 250000000
ans =
  249999975  249999982  249999989  249999996  250000003


That last value 250000003 should not be there. Somehow the mixed range with doubles and integers is causing inappropriate rounding, therefore the whole array is being resized and rewritten. And this can happen for both sievem and sievep, therefore twice the overhead.

One more experiment. Stick one or both these statements inside the sieve calculation, just before sievem and sievep are created:

lenp = double(lenp);
lenm = double(lenm);


You will see that using only one causes 50% of the slowdown, using both causes no slowdown at all, and using neither causes the full slowdown.

So this is an integer range bug, not necessarily a bug in primes.

Copying jwe for his input.

Arun Giridhar <arungiridhar>
Group Member
Wed 06 Oct 2021 04:07:03 PM UTC, comment #2: 

Wasn't there some recent work by jwe on integer ranges? Because without the patch, if you call primes with an integer class input, this ends up as the upper end of a range that is used for indexing, which perhaps is slow. But why it would be so very slow I do not see.

Michael Leitner <mleitner>
Wed 06 Oct 2021 04:02:49 PM UTC, comment #1: 

It looks like a regression. In 6.3.1 the difference is smaller.


octave:1> n = 1.5e9; tic; p1 = primes(n); toc, tic; p2 = primes(uint64(n)); toc
Elapsed time is 4.6902 seconds.
Elapsed time is 5.02894 seconds.


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 06 Oct 2021 03:13:19 PM UTC, original submission:  

The primes() function is slow for integer types by a factor of 4 to 5 times for large values of n:

>> n = 1.5e9; tic; p1 = primes(n); toc, tic; p2 = primes(uint64(n)); toc
Elapsed time is 5.46344 seconds.
Elapsed time is 25.0824 seconds.


The slowdown is entirely inside the sieve calculation loop, likely because lenm and lenp are integers like n but the loop counter i is double, and the mixed types are somehow causing a slowdown.

The attached patch fixes the slowdown:


diff -r 39a4ab124fd0 scripts/specfun/primes.m
--- a/scripts/specfun/primes.m        Wed Oct 06 10:09:48 2021 +0900
+++ b/scripts/specfun/primes.m        Wed Oct 06 10:54:29 2021 -0400
@@ -57,6 +57,10 @@
   if (ischar (n))
     n = double (n);
   endif
+
+  cls = class (n);     # if n is not double, store its class
+  n = double (n);      # and use only double for internal use
+
   if (! isfinite (n) && n != -Inf)
     error ("primes: N must be finite (not +Inf or NaN)");
   endif
@@ -102,9 +106,7 @@
     p = sort ([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]);
   endif

-  if (! isa (n, "double"))
-    p = cast (p, class (n));
-  endif
+  p = feval (cls, p);            # cast back to the type of the input

 endfunction


New performance with the patch is much closer to parity:

>> n = 1.5e9; tic; p1 = primes(n); toc, tic; p2 = primes(uint64(n)); toc
Elapsed time is 5.52101 seconds.
Elapsed time is 5.84542 seconds.


Thanks to mleitner for the discussions in bug #61129 in locating this performance bug in primes().

Arun Giridhar <arungiridhar>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #52385:  primes2.patch added by arungiridhar (2KiB - text/x-patch)
file #52376:  primes_fixed.patch added by arungiridhar (2KiB - text/x-patch)
file #52375:  range_tests.patch added by arungiridhar (4KiB - text/x-patch)
file #52371:  primes.patch added by arungiridhar (2KiB - text/x-patch)
file #52370:  integer_range.patch added by arungiridhar (981B - text/x-patch)
file #52363:  range_patch3.patch added by arungiridhar (737B - text/x-patch)
file #52362:  range_patch2.patch added by arungiridhar (848B - text/x-patch)
file #52361:  range_patch.patch added by arungiridhar (1KiB - text/x-patch)
file #52357:  primes_double.patch added by arungiridhar (2KiB - text/x-patch)
file #52056:  primes_patch1.patch added by arungiridhar (717B - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by arungiridhar (Integer range bug, paging jwe)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by nrjank (Updated the item)
  • -email is unavailable- added by arungiridhar (Submitted the item)
  • -email is unavailable- added by arungiridhar
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 24 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-12-02 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-11-30 mmuetzel StatusConfirmed Ready For Test
    2021-11-29 mmuetzel Attached File- Added bug61300-range-init-octave_int-specializations.patch, #52386
    2021-11-29 arungiridhar Attached File- Added primes2.patch, #52385
    2021-11-28 arungiridhar Attached File- Added primes_fixed.patch, #52376
    2021-11-28 arungiridhar Attached File- Added range_tests.patch, #52375
    2021-11-28 arungiridhar Attached File- Added primes.patch, #52371
    2021-11-28 arungiridhar Attached File- Added integer_range.patch, #52370
    2021-11-28 arungiridhar Attached File- Added range_patch3.patch, #52363
    2021-11-28 arungiridhar Attached File- Added range_patch2.patch, #52362
    2021-11-28 arungiridhar Attached File- Added range_patch.patch, #52361
    2021-11-28 arungiridhar Attached File- Added range_bug_isolation_WIP.patch, #52360
    2021-11-28 arungiridhar Attached File- Added primes_double.patch, #52357
    2021-10-06 mmuetzel CategoryOctave Function Interpreter
        Severity3 - Normal 4 - Important
        Item GroupPerformance Incorrect Result
        StatusPatch Submitted Confirmed
        Operating SystemGNU/Linux Any
        Summaryprimes() is 4 times slower for integer types than for double, patch attached integer range might exceed upper limit
    2021-10-06 arungiridhar Carbon-Copy- Added jwe
    2021-10-06 nrjank StatusNone Patch Submitted
    2021-10-06 arungiridhar Attached File- Added primes_patch1.patch, #52056
        Carbon-Copy- Added mleitner

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code