1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
|
########################################################################
##
## Copyright (C) 2000-2024 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## <https://www.gnu.org/licenses/>.
##
########################################################################
## -*- texinfo -*-
## @deftypefn {} {@var{tf} =} isprime (@var{x})
## Return a logical array which is true where the elements of @var{x} are prime
## numbers and false where they are not.
##
## A prime number is conventionally defined as a positive integer greater than
## 1 (e.g., 2, 3, @dots{}) which is divisible only by itself and 1. Octave
## extends this definition to include both negative integers and complex
## values. A negative integer is prime if its positive counterpart is prime.
## This is equivalent to @code{isprime (abs (x))}.
##
## If @code{class (@var{x})} is complex, then primality is tested in the domain
## of Gaussian integers (@url{https://en.wikipedia.org/wiki/Gaussian_integer}).
## Some non-complex integers are prime in the ordinary sense, but not in the
## domain of Gaussian integers. For example, @math{5 = (1+2i)*(1-2i)} shows
## that 5 is not prime because it has a factor other than itself and 1.
## Exercise caution when testing complex and real values together in the same
## matrix.
##
## Examples:
##
## @example
## @group
## isprime (1:6)
## @result{} 0 1 1 0 1 0
## @end group
## @end example
##
## @example
## @group
## isprime ([i, 2, 3, 5])
## @result{} 0 0 1 0
## @end group
## @end example
##
## Programming Note: @code{isprime} is suitable for all @var{x}
## in the range abs(@var{x})
## @tex
## $ < 2^{64}$.
## @end tex
## @ifnottex
## < 2^64.
## @end ifnottex
## Cast inputs larger than @code{flintmax} to @code{uint64}.
##
## For larger inputs, use ‘sym’ if you have the Symbolic package installed
## and loaded:
##
## @example
## @group
## isprime (sym ('58745389709258902525390450') + (0:4))
## @result{} 0 1 0 0 0
## @end group
## @end example
##
## Compatibility Note: @sc{matlab} does not extend the definition of prime
## numbers and will produce an error if given negative or complex inputs.
## @seealso{primes, factor, gcd, lcm}
## @end deftypefn
function tf = isprime (x)
if (nargin < 1)
print_usage ();
elseif (any (fix (x) != x))
error ("isprime: X contains non-integer entries");
endif
if (isempty (x))
tf = x;
return;
endif
if (iscomplex (x))
tf = isgaussianprime (x);
return;
endif
## Code strategy is to quickly compare entries in x with small primes
## using lookup(), then do direct division on larger numbers up to
## a threshold, then call Miller-Rabin for numbers over the threshold.
x = abs (x); # handle negative entries
## Generate prime table of suitable length up to maxp.
## The value of maxp needs to be at least 37,
## because of the method used by __isprimelarge__ below.
maxp = 37;
pr = [2 3 5 7 11 13 17 19 23 29 31 37];
tf = lookup (pr, x, "b"); # quick search for table matches.
THRESHOLD = 195e8;
## FIXME: THRESHOLD is the input value at which Miller-Rabin
## becomes more efficient than direct division. For smaller numbers,
## use direct division. For larger numbers, use Miller-Rabin.
##
## From numerical experiments in Jun 2022, this was observed:
## THRESHOLD Division Miller-Rabin
## 29e9 29.8196s 26.2484s (previous value)
## 20e9 26.7445s 26.0161s
## 10e9 20.9330s 25.3247s
## 19e9 26.5397s 26.8987s
## 195e8 26.5735s 26.4749s
## which is close enough, so new threshold = 195e8.
##
## The test code was this:
## n = THRESHOLD - (1:1e7); tic; isprime(n); toc
## n = THRESHOLD + (1:1e7); tic; isprime(n); toc
##
## Two notes for future programmers:
##
## 1. Test and tune THRESHOLD periodically. Miller-Rabin is only CPU-limited,
## while factorization by division is very memory-intensive. This is
## plainly noticeable from the loudness of the computer fans when running
## the division technique! CPU speed and RAM speed scale differently over
## time, so test and tune THRESHOLD periodically.
##
## 2. If you make improvements elsewhere in the code that favor one over
## the other (not symmetric), you should also retune THRESHOLD afterwards.
## If the Miller-Rabin part is sped up, the optimum THRESHOLD will
## decrease, and if factorization is sped up, it will increase.
##
## Process large entries that are still suitable for direct division
m = x (x > maxp & x <= THRESHOLD);
if ( ! isempty (m))
## Start by dividing through by the small primes until the remaining list
## of entries is small (and most likely prime themselves).
pr2 = primes (sqrt (max (m)));
tf |= lookup (pr2, x, "b");
for p = pr2
m = m(rem (m, p) != 0);
if (numel (m) < numel (pr) / 10)
break;
endif
endfor
## Check the remaining list of possible primes against the remaining
## prime factors which were not tested in the for loop.
## This is just an optimization to use arrayfun over for loop.
pr2 = pr2 (pr2 > p);
mm = arrayfun (@(x) all (rem (x, pr2)), m);
m = m(mm);
## Add any remaining entries, which are truly prime, to the results.
if ( ! isempty (m))
tf |= lookup (sort (m), x, "b");
endif
endif
## Process remaining entries (everything above THRESHOLD) with Miller-Rabin
ii = (x(:)' > THRESHOLD);
tf(ii) = __isprimelarge__ (x(ii));
endfunction
function tf = isgaussianprime (z)
## Assume prime unless proven otherwise
tf = true (size (z));
x = real (z);
y = imag (z);
## If purely real or purely imaginary, ordinary prime test for
## that complex part if that part is 3 mod 4.
xidx = y==0 & mod (x, 4) == 3;
yidx = x==0 & mod (y, 4) == 3;
tf(xidx) &= isprime (x(xidx));
tf(yidx) &= isprime (y(yidx));
## Otherwise, prime if x^2 + y^2 is prime
zidx = ! (xidx | yidx); # Skip entries that were already evaluated
zabs = x(zidx).^2 + y(zidx).^2;
tf(zidx) &= isprime (zabs);
endfunction
%!assert (isprime (3), true)
%!assert (isprime (4), false)
%!assert (isprime (uint64 (18446744073709551557)), true)
%!assert (isprime (5i), false)
%!assert (isprime (7i), true)
%!assert (isprime ([1+2i, (2+3i)*(-1+2i)]), [true, false])
%!assert (isprime (-2), true)
%!assert (isprime (complex (-2)), false)
%!assert (isprime (2i), false)
%!assert (isprime ([i, 2, 3, 5]), [false, false, true, false])
%!assert (isprime (0), false)
%!assert (isprime (magic (3)), logical ([0, 0, 0; 1, 1, 1; 0, 0, 1]))
## Test input validation
%!error <Invalid call> isprime ()
%!error <X contains non-integer entries> isprime (0.5i)
%!error <X contains non-integer entries> isprime (0.5)
|
