Mercurial > hg > octave-nkf

comparison scripts/specfun/isprime.m @ 19242:2f117c4b5cb0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
isprime.m: Document negative and complex inputs. Simplify Gaussian prime test in 3rd case where x^2 + y^2 is prime. * isprime.m: Document negative and complex inputs. Simplify Gaussian prime test in 3rd case where x^2 + y^2 is prime.
author Rik <rik@octave.org>
date Mon, 25 Aug 2014 21:30:23 -0700
parents 89e275a4f6f6
children bb0c5e182c12
comparison
equal deleted inserted replaced
19241:c573d9c70ae5 19242:2f117c4b5cb0
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {} isprime (@var{x}) 21 ## @deftypefn {Function File} {} isprime (@var{x})
22 ## Return a logical array which is true where the elements of @var{x} are 22 ## Return a logical array which is true where the elements of @var{x} are
23 ## prime numbers and false where they are not. 23 ## prime numbers and false where they are not.
24 ## 24 ##
25 ## @code{isprime} is appropriate if the maximum value in @var{x} is not too 25 ## A prime number is conventionally defined as a positive integer (1, 2, 3,
26 ## large (< 1e15). For larger values special purpose factorization code 26 ## @dots{}) which is divisible only by itself and 1. Octave extends this
27 ## should be used. 27 ## definition to include both negative integers and complex values. A
28 ## negative integer is prime if its positive counterpart is prime. This is
29 ## equivalent to @code{isprime (abs (x))}.
30 ##
31 ## If @code{class (@var{x})} is complex, then primality is tested in the domain
32 ## of Gaussian integers (@url{http://en.wikipedia.org/wiki/Gaussian_integer}).
33 ## Some non-complex integers are prime in the ordinary sense, but not in the
34 ## domain of Gaussian integers. For example, @math{5 = (1+2i)*(1-2i)} shows
35 ## that 5 is not prime because it has a factor other than itself and 1.
36 ## Exercise caution when testing complex and real values together in the same
37 ## matrix.
38 ##
39 ## Examples:
28 ## 40 ##
29 ## @example 41 ## @example
30 ## @group 42 ## @group
31 ## isprime (1:6) 43 ## isprime (1:6)
32 ## @result{} [0, 1, 1, 0, 1, 0] 44 ## @result{} [0, 1, 1, 0, 1, 0]
33 ## @end group 45 ## @end group
34 ## @end example 46 ## @end example
35 ## 47 ##
36 ## If @code{class (@var{x})} is complex, then primality will be tested
37 ## in the domain of Gaussian integers. This means that some real prime
38 ## numbers will fail the primality test in this case. For example, 5 =
39 ## (1+2i)*(1-2i) and 2 = i*(1-i)^2.
40 ##
41 ## @example 48 ## @example
42 ## @group 49 ## @group
43 ## isprime ([i, 2, 3, 5]) 50 ## isprime ([i, 2, 3, 5])
44 ## @result{} [0, 0, 1, 0] 51 ## @result{} [0, 0, 1, 0]
45 ## @end group 52 ## @end group
46 ## @end example 53 ## @end example
54 ##
55 ## Programming Note: @code{isprime} is appropriate if the maximum value in
56 ## @var{x} is not too large (< 1e15). For larger values special purpose
57 ## factorization code should be used.
58 ##
59 ## Compatibility Note: @var{matlab} does not extend the definition of prime
60 ## numbers and will produce an error if given negative or complex inputs.
47 ## @seealso{primes, factor, gcd, lcm} 61 ## @seealso{primes, factor, gcd, lcm}
48 ## @end deftypefn 62 ## @end deftypefn
49 63
50 function t = isprime (x) 64 function t = isprime (x)
51 65
63 if (iscomplex (x)) 77 if (iscomplex (x))
64 t = isgaussianprime (x); 78 t = isgaussianprime (x);
65 return; 79 return;
66 endif 80 endif
67 81
68 ## Handle negative entries 82 ## Code strategy is to build a table with the list of possible primes
69 x = abs (x); 83 ## and then quickly compare entries in x with the table of primes using
84 ## lookup(). The table size is limited to save memory and computation
85 ## time during its creation. All entries larger than the maximum in the
86 ## table are checked by straightforward division.
70 87
88 x = abs (x); # handle negative entries
71 maxn = max (x(:)); 89 maxn = max (x(:));
72 ## generate prime table of suitable length. 90 ## generate prime table of suitable length.
73 maxp = min (maxn, max (sqrt (maxn), 1e7)); # FIXME: threshold not optimized. 91 ## 1e7 threshold requires ~0.15 seconds of computation, 1e8 requires 1.8.
92 maxp = min (maxn, max (sqrt (maxn), 1e7));
74 pr = primes (maxp); 93 pr = primes (maxp);
75 ## quick search for table matches. 94 t = lookup (pr, x, "b"); # quick search for table matches.
76 t = lookup (pr, x, "b"); 95
77 ## take the rest. 96 ## process any remaining large entries
78 m = x(x > maxp); 97 m = x(x > maxp);
79 if (! isempty (m)) 98 if (! isempty (m))
80 ## there are still possible primes. filter them out by division.
81 if (maxn <= intmax ("uint32")) 99 if (maxn <= intmax ("uint32"))
82 m = uint32 (m); 100 m = uint32 (m);
83 elseif (maxn <= intmax ("uint64")) 101 elseif (maxn <= intmax ("uint64"))
84 m = uint64 (m); 102 m = uint64 (m);
85 else 103 else
86 warning ("isprime: X contains integers too large to be tested"); 104 warning ("isprime: X contains integers too large to be tested");
87 endif 105 endif
106
107 ## Start by dividing through by the small primes until the remaining
108 ## list of entries is small (and most likely prime themselves).
88 pr = cast (pr(pr <= sqrt (maxn)), class (m)); 109 pr = cast (pr(pr <= sqrt (maxn)), class (m));
89 for p = pr 110 for p = pr
90 m = m(rem (m, p) != 0); 111 m = m(rem (m, p) != 0);
91 if (length (m) < length (pr) / 10) 112 if (numel (m) < numel (pr) / 10)
92 break; 113 break;
93 endif 114 endif
94 endfor 115 endfor
116
117 ## Check the remaining list of possible primes against the
118 ## remaining prime factors which were not tested in the for loop.
119 ## This is just an optimization to use arrayfun over for loo
95 pr = pr(pr > p); 120 pr = pr(pr > p);
96 mm = arrayfun (@(x) all (rem (x, pr)), m); 121 mm = arrayfun (@(x) all (rem (x, pr)), m);
97 m = m(mm); 122 m = m(mm);
123
124 ## Add any remaining entries, which are truly prime, to the results.
98 if (! isempty (m)) 125 if (! isempty (m))
99 m = cast (sort (m), class (x)); 126 m = cast (sort (m), class (x));
100 t |= lookup (m, x, "b"); 127 t |= lookup (m, x, "b");
101 endif 128 endif
102 endif 129 endif
103 130
104 endfunction 131 endfunction
105 132
106 function t = isgaussianprime (z) 133 function t = isgaussianprime (z)
107 ## Assumed prime unless proven otherwise 134 ## Assume prime unless proven otherwise
108 t = true (size (z)); 135 t = true (size (z));
109 136
110 x = real (z); 137 x = real (z);
111 y = imag (z); 138 y = imag (z);
112 139
116 yidx = x==0 & mod (y, 4) == 3; 143 yidx = x==0 & mod (y, 4) == 3;
117 144
118 t(xidx) &= isprime (x(xidx)); 145 t(xidx) &= isprime (x(xidx));
119 t(yidx) &= isprime (y(yidx)); 146 t(yidx) &= isprime (y(yidx));
120 147
121 ## Otherwise, prime if x^2 + y^2 is prime and NOT 3 mod 4. 148 ## Otherwise, prime if x^2 + y^2 is prime
122 zabs = x.^2 + y.^2; 149 zidx = ! (xidx | yidx); # Skip entries that were already evaluated
123 zidx = mod (zabs, 4) != 3 & !xidx & !yidx; 150 zabs = x(zidx).^2 + y(zidx).^2;
124 151 t(zidx) &= isprime (zabs);
125 t(zidx) &= isprime (zabs (zidx));
126 endfunction 152 endfunction
127
128 153
129 154
130 %!assert (isprime (3), true) 155 %!assert (isprime (3), true)
131 %!assert (isprime (4), false) 156 %!assert (isprime (4), false)
132 %!assert (isprime (5i), false) 157 %!assert (isprime (5i), false)
142 %% Test input validation 167 %% Test input validation
143 %!error isprime () 168 %!error isprime ()
144 %!error isprime (1, 2) 169 %!error isprime (1, 2)
145 %!error <X contains non-integer entries> isprime (0.5i) 170 %!error <X contains non-integer entries> isprime (0.5i)
146 %!error <X contains non-integer entries> isprime (0.5) 171 %!error <X contains non-integer entries> isprime (0.5)
172