Trial Division and MATLAB isprime Function

Apr. 27, 2023

Trial Division

试除法（trial division）是一种判定某个数$n$是否为素数的方式。想要判定某个数$n$是不是素数，我们只需要使用素数$2,3,5,\cdots$一个一个地去除$n$（因为任何一个合数都可以分解为素数之积，因此只需要尝试素数即可），直到其中有一个素数正好可以整除$n$，那么我们就可以说明$n$不是素数。这种尝试直到我们将要用来除的素数大于$\sqrt n$（即尝试的素数不超过$\sqrt n$，小于等于$\sqrt n$）时就停下来。


MATLAB isprime Function

试除法的思想是很简单的，MATLAB的自建函数isprime就是通过这种方式来判断一个数是否是素数：

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

function isp = isprime(X)
%ISPRIME True for prime numbers.
%   ISPRIME(X) is 1 for the elements of X that are prime, 0 otherwise.
%
%   Class support for input X:
%      float: double, single
%      integer: uint8, int8, uint16, int16, uint32, int32, uint64, int64
%   See also FACTOR, PRIMES.
%   Copyright 1984-2012 The MathWorks, Inc. 

isp = false(size(X));

if ~isempty(X)  
    X = X(:);
    if ~isreal(X) || any(X < 0) || any(floor(X) ~= X) || ...
            any(isinf(X))
        error(message('MATLAB:isprime:InputNotPosInt'));
    end
    
    n = max(X);
    if isinteger(X) || n <= flintmax(class(X))
        if (isa(X,'uint64') || isa(X,'int64')) && n > flintmax
            p = primes(2.^(nextpow2(n)/2));
        else
            p = primes(cast(sqrt(double(n)),class(X)));
        end
        for k = 1:numel(isp)
            Xk = X(k);
            isp(k) = (Xk>1) && all(rem(Xk, p(p<Xk)));
        end
    else
        fm = flintmax(class(X));
        p = primes(sqrt(fm));
        for k = 1:numel(isp)
            Xk = X(k);
            isp(k) = (Xk<fm) && (Xk>1) && all(rem(Xk, p(p<Xk)));
        end
    end
end

如果不考虑检查输入的步骤，isprime函数主要实现了以下的工作：

  • isprime函数接受标量，向量和矩阵输入。因此，第一步是将输入展开成一个列向量（X = X(:);），之后寻找到列向量中的最大值（n = max(X);）；
  • 之后对最大值开根号，并使用primes函数得到包含所有不超过该值的素数列表p（关于primes函数，可参考博客 [4]）；
    • 对于整型输入：p = primes(2.^(nextpow2(n)/2));，这是为了防止在计算的过程中将整型变量转换为浮点数；
    • 对于浮点类型的输入：p = primes(cast(sqrt(double(n)),class(X)));；
  • 之后，对输入中的所有元素进行（向量化的）试除法：isp(k) = (Xk>1) && all(rem(Xk, p(p<Xk)));，并将判断的结果保存在isp变量中；

例如对于：

1

primes([133,457,9,13,103])

我们可以简单地可视化上述的过程：

image-20230426201640597

但是，isprime函数中语句isp(k) = (Xk>1) && all(rem(Xk, p(p<Xk)))表示在对元素进行试除法时，尝试的是素数列表p中小于该元素的所有素数。例如，对于元素9，需要对2，3，5，7四个素数进行试除。而根据试除法的原理，我们只需要对素数列表p中小于等于$\sqrt{\text{9}}=3$的所有素数进行试除即可，因此这一步会增加计算时间！！！对于元素133，13和103的试除同样会面临这样的问题：

image-20230426202608380

上图中，浅紫色的元素表示没有必要进行试除的元素。

因此，我们可以把这试除的语句更改为：

1

isp(k) = (Xk>1) && all(rem(Xk, p(p<=sqrt(Xk))));

注：一定注意是<=，而不是<，<是不对的。例如对于元素4，开根号为2，若写为<，则不存在试除的元素，最终all(4,[])的值为1，因此就会判定4是素数，这显然是不正确的。

为了验证这一点，我复制了一份isprime函数，然后仅仅对试除部分做出了如上所述的修改，将其定义为helperIsprime函数，测试了一下两个函数“重复10次判断1:1e7内所有元素是否为素数”所花费的时间：

1
2
3
4
5
6
7
8
9
10
11
12
13

clc,clear,close all

tic
for i = 1:10
    result_helper = helperIsprime(1:1e7);
end
toc

tic
for i = 1:10
    result_builtin = isprime(1:1e7);
end
toc

做出修改的代码快了三分多钟：

1
2

Elapsed time is 598.448084 seconds.
Elapsed time is 805.830040 seconds.

并且判断的结果是一致的：

1
2
3

>> sum(result_builtin~=result_helper)
ans =
     0

可以预期的是，随着上限的增加，相比于原来的isprime函数，修改过的代码会节省更多的时间。


References

[1] Derbyshire J. Prime obsession: Bernhard Riemann and the greatest unsolved problem in mathematics[M]. Joseph Henry Press, 2003.

[2] Trial division - Wikipedia.

[3] isprime - MathWorks.

[4] Sieve of Eratosthenes and MATLAB primes Function - What a starry night~.