Mercurial > octave
annotate scripts/specfun/isprime.m @ 33614:7f7d6bc5702b default tip @
maint: merge stable to default
author | Rik <rik@octave.org> |
---|---|
date | Mon, 20 May 2024 09:12:08 -0700 |
parents | 2e484f9f1f18 |
children |
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 ## |
32632
2e484f9f1f18
maint: update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
32159
diff
changeset
|
3 ## Copyright (C) 2000-2024 The Octave Project Developers |
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
27803
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:
23220
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:
23220
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:
23220
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 -*- | |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
27 ## @deftypefn {} {@var{tf} =} isprime (@var{x}) |
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
28 ## Return a logical array which is true where the elements of @var{x} are prime |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
29 ## numbers and false where they are not. |
11431
0d9640d755b1
Improve docstrings for all isXXX functions.
Rik <octave@nomad.inbox5.com>
parents:
10664
diff
changeset
|
30 ## |
19052
bb0c5e182c12
isprime.m: Update docstring to note that 1 is not prime.
Rik <rik@octave.org>
parents:
19050
diff
changeset
|
31 ## A prime number is conventionally defined as a positive integer greater than |
bb0c5e182c12
isprime.m: Update docstring to note that 1 is not prime.
Rik <rik@octave.org>
parents:
19050
diff
changeset
|
32 ## 1 (e.g., 2, 3, @dots{}) which is divisible only by itself and 1. Octave |
bb0c5e182c12
isprime.m: Update docstring to note that 1 is not prime.
Rik <rik@octave.org>
parents:
19050
diff
changeset
|
33 ## extends this definition to include both negative integers and complex |
bb0c5e182c12
isprime.m: Update docstring to note that 1 is not prime.
Rik <rik@octave.org>
parents:
19050
diff
changeset
|
34 ## values. A negative integer is prime if its positive counterpart is prime. |
bb0c5e182c12
isprime.m: Update docstring to note that 1 is not prime.
Rik <rik@octave.org>
parents:
19050
diff
changeset
|
35 ## This is equivalent to @code{isprime (abs (x))}. |
19597
db92e7e28e1f
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
19052
diff
changeset
|
36 ## |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
37 ## If @code{class (@var{x})} is complex, then primality is tested in the domain |
25143
13fd0610480f
doc: Use https whenever possible in @url entries.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
38 ## of Gaussian integers (@url{https://en.wikipedia.org/wiki/Gaussian_integer}). |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
39 ## Some non-complex integers are prime in the ordinary sense, but not in the |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
40 ## domain of Gaussian integers. For example, @math{5 = (1+2i)*(1-2i)} shows |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
41 ## that 5 is not prime because it has a factor other than itself and 1. |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
42 ## Exercise caution when testing complex and real values together in the same |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
43 ## matrix. |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
44 ## |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
45 ## Examples: |
5827 | 46 ## |
11431
0d9640d755b1
Improve docstrings for all isXXX functions.
Rik <octave@nomad.inbox5.com>
parents:
10664
diff
changeset
|
47 ## @example |
0d9640d755b1
Improve docstrings for all isXXX functions.
Rik <octave@nomad.inbox5.com>
parents:
10664
diff
changeset
|
48 ## @group |
0d9640d755b1
Improve docstrings for all isXXX functions.
Rik <octave@nomad.inbox5.com>
parents:
10664
diff
changeset
|
49 ## isprime (1:6) |
26600
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
50 ## @result{} 0 1 1 0 1 0 |
11431
0d9640d755b1
Improve docstrings for all isXXX functions.
Rik <octave@nomad.inbox5.com>
parents:
10664
diff
changeset
|
51 ## @end group |
0d9640d755b1
Improve docstrings for all isXXX functions.
Rik <octave@nomad.inbox5.com>
parents:
10664
diff
changeset
|
52 ## @end example |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
53 ## |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
54 ## @example |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
55 ## @group |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
56 ## isprime ([i, 2, 3, 5]) |
26600
f6730533820e
doc: clean up doc example blocks for accuracy and consistent formatting
Mike Miller <mtmiller@octave.org>
parents:
26376
diff
changeset
|
57 ## @result{} 0 0 1 0 |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
58 ## @end group |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
59 ## @end example |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
60 ## |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
61 ## Programming Note: @code{isprime} is suitable for all @var{x} |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
62 ## in the range abs(@var{x}) |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
63 ## @tex |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
64 ## $ < 2^{64}$. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
65 ## @end tex |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
66 ## @ifnottex |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
67 ## < 2^64. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
68 ## @end ifnottex |
32159
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
69 ## Cast inputs larger than @code{flintmax} to @code{uint64}. |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
70 ## |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
71 ## For larger inputs, use ‘sym’ if you have the Symbolic package installed |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
72 ## and loaded: |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
73 ## |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
74 ## @example |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
75 ## @group |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
76 ## isprime (sym ('58745389709258902525390450') + (0:4)) |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
77 ## @result{} 0 1 0 0 0 |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
78 ## @end group |
efcfafb7ad16
doc: Add note to isprime.m on inputs over 2^64.
Arun Giridhar <arungiridhar@gmail.com>
parents:
31706
diff
changeset
|
79 ## @end example |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
80 ## |
27803
e8823425a60b
doc: Change markup for "Matlab".
Markus Mützel <markus.muetzel@gmx.de>
parents:
26600
diff
changeset
|
81 ## Compatibility Note: @sc{matlab} does not extend the definition of prime |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
82 ## numbers and will produce an error if given negative or complex inputs. |
5827 | 83 ## @seealso{primes, factor, gcd, lcm} |
84 ## @end deftypefn | |
85 | |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
86 function tf = isprime (x) |
7125 | 87 |
28891
de5f2f9a64ff
maint: Use same coding style when checking for a minimum of 1 input.
Rik <rik@octave.org>
parents:
27923
diff
changeset
|
88 if (nargin < 1) |
19041
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
89 print_usage (); |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
90 elseif (any (fix (x) != x)) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
91 error ("isprime: X contains non-integer entries"); |
19041
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
92 endif |
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
93 |
19045
920a400929ca
* isprime.m: Return an empty array for empty input
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19041
diff
changeset
|
94 if (isempty (x)) |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
95 tf = x; |
19045
920a400929ca
* isprime.m: Return an empty array for empty input
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19041
diff
changeset
|
96 return; |
920a400929ca
* isprime.m: Return an empty array for empty input
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19041
diff
changeset
|
97 endif |
920a400929ca
* isprime.m: Return an empty array for empty input
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19041
diff
changeset
|
98 |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
99 if (iscomplex (x)) |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
100 tf = isgaussianprime (x); |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
101 return; |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
102 endif |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
103 |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
104 ## Code strategy is to quickly compare entries in x with small primes |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
105 ## using lookup(), then do direct division on larger numbers up to |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
106 ## a threshold, then call Miller-Rabin for numbers over the threshold. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
107 |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
108 x = abs (x); # handle negative entries |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
109 |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
110 ## Generate prime table of suitable length up to maxp. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
111 ## The value of maxp needs to be at least 37, |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
112 ## because of the method used by __isprimelarge__ below. |
30400
168da23530b4
maint: strip trailing spaces from code base.
Rik <rik@octave.org>
parents:
30261
diff
changeset
|
113 maxp = 37; |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
114 pr = [2 3 5 7 11 13 17 19 23 29 31 37]; |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
115 tf = lookup (pr, x, "b"); # quick search for table matches. |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
116 |
31050
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
117 THRESHOLD = 195e8; |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
118 ## FIXME: THRESHOLD is the input value at which Miller-Rabin |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
119 ## becomes more efficient than direct division. For smaller numbers, |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
120 ## use direct division. For larger numbers, use Miller-Rabin. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
121 ## |
31050
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
122 ## From numerical experiments in Jun 2022, this was observed: |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
123 ## THRESHOLD Division Miller-Rabin |
31050
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
124 ## 29e9 29.8196s 26.2484s (previous value) |
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
125 ## 20e9 26.7445s 26.0161s |
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
126 ## 10e9 20.9330s 25.3247s |
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
127 ## 19e9 26.5397s 26.8987s |
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
128 ## 195e8 26.5735s 26.4749s |
3a93a77dd4cf
Minor performance tweaks to isprime.m and __isprimelarge__.cc
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
129 ## which is close enough, so new threshold = 195e8. |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
130 ## |
30400
168da23530b4
maint: strip trailing spaces from code base.
Rik <rik@octave.org>
parents:
30261
diff
changeset
|
131 ## The test code was this: |
168da23530b4
maint: strip trailing spaces from code base.
Rik <rik@octave.org>
parents:
30261
diff
changeset
|
132 ## n = THRESHOLD - (1:1e7); tic; isprime(n); toc |
168da23530b4
maint: strip trailing spaces from code base.
Rik <rik@octave.org>
parents:
30261
diff
changeset
|
133 ## n = THRESHOLD + (1:1e7); tic; isprime(n); toc |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
134 ## |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
135 ## Two notes for future programmers: |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
136 ## |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
137 ## 1. Test and tune THRESHOLD periodically. Miller-Rabin is only CPU-limited, |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
138 ## while factorization by division is very memory-intensive. This is |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
139 ## plainly noticeable from the loudness of the computer fans when running |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
140 ## the division technique! CPU speed and RAM speed scale differently over |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
141 ## time, so test and tune THRESHOLD periodically. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
142 ## |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
143 ## 2. If you make improvements elsewhere in the code that favor one over |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
144 ## the other (not symmetric), you should also retune THRESHOLD afterwards. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
145 ## If the Miller-Rabin part is sped up, the optimum THRESHOLD will |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
146 ## decrease, and if factorization is sped up, it will increase. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
147 ## |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
148 |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
149 ## Process large entries that are still suitable for direct division |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
150 m = x (x > maxp & x <= THRESHOLD); |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
151 if ( ! isempty (m)) |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
152 ## Start by dividing through by the small primes until the remaining list |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
153 ## of entries is small (and most likely prime themselves). |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
154 pr2 = primes (sqrt (max (m))); |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
155 tf |= lookup (pr2, x, "b"); |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
156 for p = pr2 |
19041
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
157 m = m(rem (m, p) != 0); |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
158 if (numel (m) < numel (pr) / 10) |
19041
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
159 break; |
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
160 endif |
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
161 endfor |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
162 |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
163 ## Check the remaining list of possible primes against the remaining |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
164 ## prime factors which were not tested in the for loop. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
165 ## This is just an optimization to use arrayfun over for loop. |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
166 pr2 = pr2 (pr2 > p); |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
167 mm = arrayfun (@(x) all (rem (x, pr2)), m); |
19041
0976f9fccbbd
isprime.m: Verify that input is a real number (bug #43041)
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
168 m = m(mm); |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
169 |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
170 ## Add any remaining entries, which are truly prime, to the results. |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
171 if ( ! isempty (m)) |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
172 tf |= lookup (sort (m), x, "b"); |
10664
faff5367cc05
second isprime rewrite
Jaroslav Hajek <highegg@gmail.com>
parents:
10657
diff
changeset
|
173 endif |
7125 | 174 endif |
175 | |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
176 ## Process remaining entries (everything above THRESHOLD) with Miller-Rabin |
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
177 ii = (x(:)' > THRESHOLD); |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
178 tf(ii) = __isprimelarge__ (x(ii)); |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
179 |
5827 | 180 endfunction |
9984
d1cc2e0ddf55
isprime: produce logical result
John W. Eaton <jwe@octave.org>
parents:
7125
diff
changeset
|
181 |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
182 function tf = isgaussianprime (z) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
183 |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
184 ## Assume prime unless proven otherwise |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
185 tf = true (size (z)); |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
186 |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
187 x = real (z); |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
188 y = imag (z); |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
189 |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
190 ## If purely real or purely imaginary, ordinary prime test for |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
191 ## that complex part if that part is 3 mod 4. |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
192 xidx = y==0 & mod (x, 4) == 3; |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
193 yidx = x==0 & mod (y, 4) == 3; |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
194 |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
195 tf(xidx) &= isprime (x(xidx)); |
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
196 tf(yidx) &= isprime (y(yidx)); |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
197 |
19050
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
198 ## Otherwise, prime if x^2 + y^2 is prime |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
199 zidx = ! (xidx | yidx); # Skip entries that were already evaluated |
2f117c4b5cb0
isprime.m: Document negative and complex inputs.
Rik <rik@octave.org>
parents:
19046
diff
changeset
|
200 zabs = x(zidx).^2 + y(zidx).^2; |
30558
83aeaba707d8
doc: Use TF for output variable in documentation for isXXX functions in scripts/ directory.
Rik <rik@octave.org>
parents:
30400
diff
changeset
|
201 tf(zidx) &= isprime (zabs); |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
202 |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
203 endfunction |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
204 |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
205 |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
206 %!assert (isprime (3), true) |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
207 %!assert (isprime (4), false) |
30261
a49c635b179d
isprime.m: Speed up for large integers using Miller-Rabin test (bug #61312).
Arun Giridhar <arungiridhar@gmail.com>
parents:
29359
diff
changeset
|
208 %!assert (isprime (uint64 (18446744073709551557)), true) |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
209 %!assert (isprime (5i), false) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
210 %!assert (isprime (7i), true) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
211 %!assert (isprime ([1+2i, (2+3i)*(-1+2i)]), [true, false]) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
212 %!assert (isprime (-2), true) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
213 %!assert (isprime (complex (-2)), false) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
214 %!assert (isprime (2i), false) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
215 %!assert (isprime ([i, 2, 3, 5]), [false, false, true, false]) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
216 %!assert (isprime (0), false) |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
217 %!assert (isprime (magic (3)), logical ([0, 0, 0; 1, 1, 1; 0, 0, 1])) |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
218 |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
219 ## Test input validation |
28896
90fea9cc9caa
test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents:
28891
diff
changeset
|
220 %!error <Invalid call> isprime () |
19046
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
221 %!error <X contains non-integer entries> isprime (0.5i) |
89e275a4f6f6
Allow isprime to handle complex numbers (bug #43041)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
19045
diff
changeset
|
222 %!error <X contains non-integer entries> isprime (0.5) |