2

我目前正在尝试在 MATLAB 中编写一个程序来检查一个数字n是否为素数。对于初学者，我正在实施Fermat Primality Test。

费马指出，对于素数p和1 <= b < p：

b^(p-1) = 1  (mod p)

所以在 MATLAB 中p = 17，和b = 11

>> mod(b^(p-1),p)

或者

>> rem(b^(p-1),p)

我遇到的问题是，对于这种情况，MATLAB 返回0. 但是，如果p是素数，它应该返回1。我看不到我缺少什么，因此非常感谢任何帮助！

4

3 回答 3

9

@James 给出了正确的解释，我只是想扩展一点。

您在双精度浮点格式中看到，范围内的整数[2^52,2^53]是完全可表示的（那是因为我们有 52+1 位用于小数部分）。在下一个范围[2^53,2^54]中，可表示的整数是偶数（上一个范围乘以 2）。以此类推，每次我们走高时间距都会翻倍。

不幸的是，这个数字11^16（等于45949729863572161）不能用双精度精确表示。事实上，围绕该数字的可表示数字列表是：

45949729863572144
45949729863572152
45949729863572160
45949729863572168
45949729863572176

根据舍入模式，45949729863572161将近似为最接近的可表示数字，在这种情况下为45949729863572160。

要了解会发生什么，让我们尝试存储数字45949729863572100 + [44:76]并显示结果：

% build a cell array of strings containing the numbers, then convert to doubles
% (you could also enter the numbers as literals directly)
str = cellstr(num2str((44:76)', '459497298635721%d'));
num = str2double(str);

% print the original number, its stored value (in decimal and hex notations)
for i=1:numel(num)
    fprintf('%s %17.0f %bX\n', str{i}, num(i), num(i));
end

这是输出（带有一些注释）：

    actual           stored          stored in HEX
----------------------------------------------------
45949729863572144 45949729863572144 436467E125C16356    % exact representation
45949729863572145 45949729863572144 436467E125C16356
45949729863572146 45949729863572144 436467E125C16356
45949729863572147 45949729863572144 436467E125C16356
45949729863572148 45949729863572144 436467E125C16356
45949729863572149 45949729863572152 436467E125C16357
45949729863572150 45949729863572152 436467E125C16357
45949729863572151 45949729863572152 436467E125C16357
45949729863572152 45949729863572152 436467E125C16357    % exact representation
45949729863572153 45949729863572152 436467E125C16357
45949729863572154 45949729863572152 436467E125C16357
45949729863572155 45949729863572152 436467E125C16357
45949729863572156 45949729863572160 436467E125C16358
45949729863572157 45949729863572160 436467E125C16358
45949729863572158 45949729863572160 436467E125C16358
45949729863572159 45949729863572160 436467E125C16358
45949729863572160 45949729863572160 436467E125C16358    % exact representation
45949729863572161 45949729863572160 436467E125C16358
45949729863572162 45949729863572160 436467E125C16358
45949729863572163 45949729863572160 436467E125C16358
45949729863572164 45949729863572160 436467E125C16358
45949729863572165 45949729863572168 436467E125C16359
45949729863572166 45949729863572168 436467E125C16359
45949729863572167 45949729863572168 436467E125C16359
45949729863572168 45949729863572168 436467E125C16359    % exact representation
45949729863572169 45949729863572168 436467E125C16359
45949729863572170 45949729863572168 436467E125C16359
45949729863572171 45949729863572168 436467E125C16359
45949729863572172 45949729863572176 436467E125C1635A
45949729863572173 45949729863572176 436467E125C1635A
45949729863572174 45949729863572176 436467E125C1635A
45949729863572175 45949729863572176 436467E125C1635A
45949729863572176 45949729863572176 436467E125C1635A    % exact representation

如您所见，xxx44和之间不能有数字xxx52（因为它们的十六进制表示仅在最后一位不同）。介于两者之间的任何东西都必须近似为最接近的可表示数字。所以范围被二除，一半分配给下限，另一半分配给上限（注意中间有7个数字，所以中间的一个是特殊情况，分配给上/下限以交替方式限制）。

因此，输入 和 之间的任何数字45949729863572156（45949729863572164包括11^16）实际上将存储双值45949729863572160。

现在其他人建议使用bignum库来避免这些数值限制（来自 MathWorks的Symbolic Math Toolbox、John D'Errico的VPI或HPF，或 File Exchange 上可用的其他解决方案之一......）。例如：

>> b = sym(11);    % Symbolic Math Toolbox
>> b^16
ans =
45949729863572161
>> mod(b^16,17)
ans =
1

但是，在您的情况下，uint64能够准确存储这些数字：

>> b = uint64(11); p = uint64(17);
>> b^(p-1)
ans =
    45949729863572161
>> mod(b^(p-1),p)
ans =
                    1

请记住：

>> intmax('uint64')
ans =
 18446744073709551615
于 2014-01-21T04:07:10.197 回答
2

单个整数只能2^53由双精度数“连续”表示。11^16大于此，因此使用近似值。要执行此计算，您必须使用任意精度的整数数据结构。这是一个执行此操作的附加组件。

于 2014-01-21T00:48:44.080 回答
1

我相信这是由于浮点数的四舍五入。这是关于 的帮助部分mod。

EDU>> help mod

MOD    Modulus after division.
MOD(x,y) is x - n.*y where n = floor(x./y) if y ~= 0.  If y is not an
integer and the quotient x./y is within roundoff error of an integer,
then n is that integer.  The inputs x and y must be real arrays of the
same size, or real scalars.

The statement "x and y are congruent mod m" means mod(x,m) == mod(y,m).

By convention:
   MOD(x,0) is x.
   MOD(x,x) is 0.
   MOD(x,y), for x~=y and y~=0, has the same sign as y.

Note: REM(x,y), for x~=y and y~=0, has the same sign as x.
MOD(x,y) and REM(x,y) are equal if x and y have the same sign, but
differ by y if x and y have different signs.

See also rem.

Overloaded methods:
   sym/mod

Reference page in Help browser
   doc mod

在 MATLAB中实现 mod 似乎很奇怪floor(x./y)，但我有理由相信这就是原因。

编辑：我相信这可以帮助你。

于 2014-01-21T00:49:07.623 回答