Mercurial > octave
annotate scripts/specfun/primes.m @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 28 Dec 2021 18:22:40 -0500 |
| parents | bbf1293bd255 |
| children | 9a722a4316b6 |
| rev | line source |
|---|---|
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 ######################################################################## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 ## |
|
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
30396
diff
changeset
|
3 ## Copyright (C) 2000-2022 The Octave Project Developers |
|
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 ## |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 ## See the file COPYRIGHT.md in the top-level directory of this |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 ## distribution or <https://octave.org/copyright/>. |
| 5827 | 7 ## |
| 8 ## This file is part of Octave. | |
| 9 ## | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
24515
diff
changeset
|
10 ## Octave is free software: you can redistribute it and/or modify it |
| 5827 | 11 ## under the terms of the GNU General Public License as published by |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
24515
diff
changeset
|
12 ## the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
13 ## (at your option) any later version. |
| 5827 | 14 ## |
| 15 ## Octave is distributed in the hope that it will be useful, but | |
| 16 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
18 ## GNU General Public License for more details. |
| 5827 | 19 ## |
| 20 ## You should have received a copy of the GNU General Public License | |
| 7016 | 21 ## along with Octave; see the file COPYING. If not, see |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
24515
diff
changeset
|
22 ## <https://www.gnu.org/licenses/>. |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 ## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 ######################################################################## |
| 5827 | 25 |
| 26 ## -*- texinfo -*- | |
|
24515
070e9b036141
primes.m: Correctly return output of same class as input (bug #52787).
Sahil Yadav <yadavsahil5198@gmail.com>
parents:
23220
diff
changeset
|
27 ## @deftypefn {} {@var{p} =} primes (@var{n}) |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
28 ## Return all primes up to @var{n}. |
| 5827 | 29 ## |
|
24515
070e9b036141
primes.m: Correctly return output of same class as input (bug #52787).
Sahil Yadav <yadavsahil5198@gmail.com>
parents:
23220
diff
changeset
|
30 ## The output data class (double, single, uint32, etc.@:) is the same as the |
|
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
31 ## input class of @var{n}. The algorithm used is the Sieve of Eratosthenes. |
|
9141
c1fff751b5a8
Update section 17.1 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
32 ## |
|
24515
070e9b036141
primes.m: Correctly return output of same class as input (bug #52787).
Sahil Yadav <yadavsahil5198@gmail.com>
parents:
23220
diff
changeset
|
33 ## Note: If you need a specific number of primes you can use the fact that the |
|
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
34 ## distance from one prime to the next is, on average, proportional to the |
|
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
35 ## logarithm of the prime. Integrating, one finds that there are about |
|
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
36 ## @math{k} primes less than |
|
9165
8c71a86c4bf4
Update section 17.5 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9141
diff
changeset
|
37 ## @tex |
|
8c71a86c4bf4
Update section 17.5 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9141
diff
changeset
|
38 ## $k \log (5 k)$. |
|
8c71a86c4bf4
Update section 17.5 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9141
diff
changeset
|
39 ## @end tex |
|
8c71a86c4bf4
Update section 17.5 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9141
diff
changeset
|
40 ## @ifnottex |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
41 ## k*log (5*k). |
|
9165
8c71a86c4bf4
Update section 17.5 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9141
diff
changeset
|
42 ## @end ifnottex |
| 19096 | 43 ## |
|
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
44 ## See also @code{list_primes} if you need a specific number @var{n} of primes. |
|
9165
8c71a86c4bf4
Update section 17.5 (Utility Functions) of arith.txi
Rik <rdrider0-list@yahoo.com>
parents:
9141
diff
changeset
|
45 ## @seealso{list_primes, isprime} |
| 5827 | 46 ## @end deftypefn |
| 47 | |
| 19096 | 48 function p = primes (n) |
| 5827 | 49 |
|
28891
de5f2f9a64ff
maint: Use same coding style when checking for a minimum of 1 input.
Rik <rik@octave.org>
parents:
28886
diff
changeset
|
50 if (nargin < 1) |
| 5827 | 51 print_usage (); |
| 52 endif | |
| 53 | |
| 29985 | 54 if (! (isscalar (n) && isreal (n))) |
| 55 error ("primes: N must be a real scalar"); | |
| 56 endif | |
| 57 if (ischar (n)) | |
| 58 n = double (n); | |
| 59 endif | |
| 60 if (! isfinite (n) && n != -Inf) | |
| 61 error ("primes: N must be finite (not +Inf or NaN)"); | |
| 5827 | 62 endif |
| 63 | |
|
30396
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
64 cls = class (n); # if n is not double, store its class |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
65 n = double (n); # and use only double for internal use. |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
66 # This conversion is needed for both calculation speed (twice as fast as |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
67 # integer) and also for the accuracy of the sieve calculation when given |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
68 # integer input, to avoid unwanted rounding in the sieve lengths. |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
69 |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
70 if (n > flintmax ()) |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
71 warning ("primes: input exceeds flintmax. Results may be inaccurate."); |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
72 endif |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
73 |
| 29985 | 74 if (n < 353) |
| 75 ## Lookup table of first 70 primes | |
| 76 a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, ... | |
| 77 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, ... | |
| 78 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ... | |
| 79 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ... | |
| 80 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ... | |
| 81 293, 307, 311, 313, 317, 331, 337, 347, 349]; | |
| 82 p = a(a <= n); | |
| 83 elseif (n < 100e3) | |
| 84 ## Classical Sieve algorithm | |
| 85 ## Fast, but memory scales as n/2. | |
| 86 len = floor ((n-1)/2); # length of the sieve | |
| 87 sieve = true (1, len); # assume every odd number is prime | |
| 88 for i = 1:(sqrt (n)-1)/2 # check up to sqrt (n) | |
| 89 if (sieve(i)) # if i is prime, eliminate multiples of i | |
| 90 sieve(3*i+1:2*i+1:len) = false; # do it | |
| 91 endif | |
| 92 endfor | |
| 93 p = [2, 1+2*find(sieve)]; # primes remaining after sieve | |
| 94 else | |
| 95 ## Sieve algorithm optimized for large n | |
| 96 ## Memory scales as n/3 or 1/6th less than classical Sieve | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
97 lenm = floor ((n+1)/6); # length of the 6n-1 sieve |
|
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
98 lenp = floor ((n-1)/6); # length of the 6n+1 sieve |
|
10657
c6833d31f34e
optimize primes and isprime
Jaroslav Hajek <highegg@gmail.com>
parents:
10549
diff
changeset
|
99 sievem = true (1, lenm); # assume every number of form 6n-1 is prime |
|
c6833d31f34e
optimize primes and isprime
Jaroslav Hajek <highegg@gmail.com>
parents:
10549
diff
changeset
|
100 sievep = true (1, lenp); # assume every number of form 6n+1 is prime |
| 5827 | 101 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
102 for i = 1:(sqrt (n)+1)/6 # check up to sqrt (n) |
| 5827 | 103 if (sievem(i)) # if i is prime, eliminate multiples of i |
|
10657
c6833d31f34e
optimize primes and isprime
Jaroslav Hajek <highegg@gmail.com>
parents:
10549
diff
changeset
|
104 sievem(7*i-1:6*i-1:lenm) = false; |
|
c6833d31f34e
optimize primes and isprime
Jaroslav Hajek <highegg@gmail.com>
parents:
10549
diff
changeset
|
105 sievep(5*i-1:6*i-1:lenp) = false; |
| 5827 | 106 endif # if i is prime, eliminate multiples of i |
| 107 if (sievep(i)) | |
|
10657
c6833d31f34e
optimize primes and isprime
Jaroslav Hajek <highegg@gmail.com>
parents:
10549
diff
changeset
|
108 sievep(7*i+1:6*i+1:lenp) = false; |
|
c6833d31f34e
optimize primes and isprime
Jaroslav Hajek <highegg@gmail.com>
parents:
10549
diff
changeset
|
109 sievem(5*i+1:6*i+1:lenm) = false; |
| 5827 | 110 endif |
| 111 endfor | |
| 19096 | 112 p = sort ([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]); |
| 113 endif | |
| 114 | |
|
30396
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
115 # cast back to the type of the input |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
116 p = cast (p, cls); |
| 5827 | 117 |
| 118 endfunction | |
|
12815
918610ea2f34
codesprint: new tests for specfun directory
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
119 |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
120 |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
121 %!assert (size (primes (350)), [1, 70]) |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
122 %!assert (primes (357)(end), 353) |
|
30396
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
123 %!assert (primes (uint64 (358))(end), uint64 (353)) |
|
bbf1293bd255
primes.m: Convert input to double internally for accuracy and speed (bug #61300).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29985
diff
changeset
|
124 %!assert (primes (int32 (1e6))(end), int32 (999983)) |
|
24515
070e9b036141
primes.m: Correctly return output of same class as input (bug #52787).
Sahil Yadav <yadavsahil5198@gmail.com>
parents:
23220
diff
changeset
|
125 %!assert (class (primes (single (10))), "single") |
|
070e9b036141
primes.m: Correctly return output of same class as input (bug #52787).
Sahil Yadav <yadavsahil5198@gmail.com>
parents:
23220
diff
changeset
|
126 %!assert (class (primes (uint8 (10))), "uint8") |
| 29985 | 127 %!assert (primes (-Inf), zeros (1,0)) |
|
12815
918610ea2f34
codesprint: new tests for specfun directory
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
128 |
|
28886
d8318c12d903
test: remove unnecessary BIST tests in m-files checking for excessive number of inputs.
Rik <rik@octave.org>
parents:
27985
diff
changeset
|
129 %!error <Invalid call> primes () |
| 29985 | 130 %!error <N must be a real scalar> primes (ones (2,2)) |
| 131 %!error <N must be a real scalar> primes (5i) | |
| 132 %!error <N must be finite> primes (Inf) | |
| 133 %!error <N must be finite> primes (NaN) |