费马小定理在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给出了正确的解释,我只是想扩展一点.

您可以看到双精度浮点格式,范围内的整数[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(因为它们的十六进制表示仅通过一个在最后一位不同).中间的任何内容都必须近似为最接近的可表示数字.所以范围除以2,一半分配到下限,另一半分配到上限(注意中间有7个数字,所以中间的一个是特殊情况,并被分配到上/下以交替的方式界限).

因此,在45949729863572156和之间输入任意数字45949729863572164(包括11^16)将实际存储双值45949729863572160.

现在,其他人使用建议BIGNUM库,以避免这些数值限制(的符号数学工具箱从MathWorks公司,VPI或HPF由约翰D'ERRICO,或其他解决方案的一个对文件交换可用...).例如:

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