Mercurial > octave-libtiff
changeset 30396:bbf1293bd255 stable
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
* scripts/specfun/primes.m: Use double internally. Add BISTs for integer input.
author | Arun Giridhar <arungiridhar@gmail.com> |
---|---|
date | Mon, 29 Nov 2021 12:36:04 -0500 |
parents | f3f3e3793fb5 |
children | 7e50c0369940 d4d83344d653 |
files | scripts/specfun/primes.m |
diffstat | 1 files changed, 14 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/primes.m Mon Nov 29 10:03:47 2021 -0800 +++ b/scripts/specfun/primes.m Mon Nov 29 12:36:04 2021 -0500 @@ -61,6 +61,16 @@ error ("primes: N must be finite (not +Inf or NaN)"); endif + cls = class (n); # if n is not double, store its class + n = double (n); # and use only double for internal use. + # This conversion is needed for both calculation speed (twice as fast as + # integer) and also for the accuracy of the sieve calculation when given + # integer input, to avoid unwanted rounding in the sieve lengths. + + if (n > flintmax ()) + warning ("primes: input exceeds flintmax. Results may be inaccurate."); + endif + 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, ... @@ -102,15 +112,16 @@ p = sort ([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]); endif - if (! isa (n, "double")) - p = cast (p, class (n)); - endif + # cast back to the type of the input + p = cast (p, cls); endfunction %!assert (size (primes (350)), [1, 70]) %!assert (primes (357)(end), 353) +%!assert (primes (uint64 (358))(end), uint64 (353)) +%!assert (primes (int32 (1e6))(end), int32 (999983)) %!assert (class (primes (single (10))), "single") %!assert (class (primes (uint8 (10))), "uint8") %!assert (primes (-Inf), zeros (1,0))