Mercurial > octave
changeset 29985:3eb6483241d0
primes.m: Overhaul function.
Add support for char inputs for Matlab compatibility. Do stricter
input validiton to weed out NaN and Inf values.
* NEWS: Announce support for char inputs in Matlab Compatibility section.
* primes.m: Cast char inputs to double during input validation. Check
for non-finite values +Inf and NaN and throw an error. Emit an error
for any complex input. Re-order if/elseif/else tree to have most common
branches first. Add new BIST tests.
| author | Rik <rik@octave.org> |
|---|---|
| date | Tue, 17 Aug 2021 21:11:45 -0700 |
| parents | c633d34960b4 |
| children | 3f598ae945f2 |
| files | NEWS scripts/specfun/primes.m |
| diffstat | 2 files changed, 38 insertions(+), 27 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Tue Aug 17 16:35:23 2021 -0700 +++ b/NEWS Tue Aug 17 21:11:45 2021 -0700 @@ -183,6 +183,8 @@ - The function `factor` now supports uint64 inputs larger than `flintmax`. +- The function `primes` now supports char inputs. + - The functions `quantile` and `prctile` now permit operating on dimensions greater than `ndims (x)`.
--- a/scripts/specfun/primes.m Tue Aug 17 16:35:23 2021 -0700 +++ b/scripts/specfun/primes.m Tue Aug 17 21:11:45 2021 -0700 @@ -51,16 +51,39 @@ print_usage (); endif - if (! (isnumeric (n) && isscalar (n))) - error ("primes: N must be a numeric scalar"); + if (! (isscalar (n) && isreal (n))) + error ("primes: N must be a real scalar"); + endif + if (ischar (n)) + n = double (n); + endif + if (! isfinite (n) && n != -Inf) + error ("primes: N must be finite (not +Inf or NaN)"); endif - if (n > 100e3) - ## Optimization: 1/6 less memory, and much faster (asymptotically) - ## 100K happens to be the cross-over point for Paul's machine; - ## below this the more direct code below is faster. At the limit - ## of memory in Paul's machine, this saves .7 seconds out of 7 for - ## n = 3e6. Hardly worthwhile, but Dirk reports better numbers. + if (n < 353) + ## Lookup table of first 70 primes + a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, ... + 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, ... + 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ... + 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ... + 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ... + 293, 307, 311, 313, 317, 331, 337, 347, 349]; + p = a(a <= n); + elseif (n < 100e3) + ## Classical Sieve algorithm + ## Fast, but memory scales as n/2. + len = floor ((n-1)/2); # length of the sieve + 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 + sieve(3*i+1:2*i+1:len) = false; # do it + endif + endfor + p = [2, 1+2*find(sieve)]; # primes remaining after sieve + 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 sievem = true (1, lenm); # assume every number of form 6n-1 is prime @@ -77,23 +100,6 @@ endif endfor p = sort ([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]); - elseif (n > 352) # nothing magical about 352; must be > 2 - len = floor ((n-1)/2); # length of the sieve - 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 - sieve(3*i+1:2*i+1:len) = false; # do it - endif - endfor - p = [2, 1+2*find(sieve)]; # primes remaining after sieve - else - a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, ... - 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, ... - 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ... - 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ... - 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ... - 293, 307, 311, 313, 317, 331, 337, 347, 349]; - p = a(a <= n); endif if (! isa (n, "double")) @@ -107,7 +113,10 @@ %!assert (primes (357)(end), 353) %!assert (class (primes (single (10))), "single") %!assert (class (primes (uint8 (10))), "uint8") +%!assert (primes (-Inf), zeros (1,0)) %!error <Invalid call> primes () -%!error <N must be a numeric scalar> primes ("1") -%!error <N must be a numeric scalar> primes (ones (2,2)) +%!error <N must be a real scalar> primes (ones (2,2)) +%!error <N must be a real scalar> primes (5i) +%!error <N must be finite> primes (Inf) +%!error <N must be finite> primes (NaN)