[Octave-bug-tracker] [bug #61300] integer range might exceed upper limit

[Arun Giridhar]

Follow-up Comment #9, bug #61300 (project octave):

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

Additional Item Attachment:

File name: primes_double.patch            Size:1 KB
    <https://file.savannah.gnu.org/file/primes_double.patch?file_id=52357>



    _______________________________________________________

Reply to this item at:

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

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