matlab - 费马小定理在 MATLAB 中失败？

Question

我目前正尝试在 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 .我看不到我遗漏了什么，所以非常感谢任何帮助!

Answer 1

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

您在 double-precision floating-point format 中看到，[2^52,2^53] 范围内的整数是可精确表示的(那是因为小数部分有 52+1 位)。在下一个范围 [2^53,2^54] 中，可表示的整数是偶数(前一个范围乘以二)。以此类推，每次我们走得更高时，间距都会加倍。

不幸的是，数字 11^16(等于 45949729863572161)不能用 double 精确表示。事实上，围绕那个数字的可表示数字列表是:

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 之间不能有数字(因为它们的 HEX 表示仅在最后一位不同)。介于两者之间的任何东西都必须近似为最接近的可表示数字。所以范围一分为二，一半分配给下限，另一半分配给上限(注意中间有7个数字，所以中间一个是特例，分配给上限/下限以交替方式边界)。

因此，输入 45949729863572156 和 45949729863572164 之间的任何数字(包括 11^16)实际上将存储 double 值 45949729863572160.

---

现在其他人建议使用 bignum库以避免这些数字限制(来自 MathWorks 的 Symbolic Math Toolbox，John D'Errico 的 VPI 或 HPF，或文件交换中可用的其他解决方案之一......)。例如:

>> 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