logo标签列表

小型数字的最快素数测试

mathprimessieve
7

我在业余时间玩项目欧拉，现在需要进行一些重构。我已经实现了Miller-Rabin算法和一些筛子。我之前听说过对于小型数字（即几百万以下），筛子实际上更快。有人对此有什么信息吗？谷歌并没有提供很多帮助。

在第10个问题上，我的试除法算法在计算到根号n时比米勒-拉宾算法更加出色。 - afkbowflexin
为什么不在Trie中记忆先前出现的质数呢？这是一项非常便宜的操作。 - Hamish Grubijan
为什么不试试呢？看看这个答案，其中包含一个简单的Java筛法，我在欧拉计划中使用了很多次：https://dev59.com/vXNA5IYBdhLWcg3wPLL-#1043247 - starblue
5个回答
9

大多数算法都可以通过牺牲空间来换取时间。换句话说，通过允许使用更多的内存，可以大大提高速度*a。

我并不真正了解米勒-拉宾算法，但是，除非它比单个移位/加法和内存提取更简单，否则它将被预先计算的筛子算法击败。

这里重要的是预先计算。从性能的角度来看，像这样的事情最好预先计算，因为前一百万个质数在不久的将来很可能不会改变 :-)

换句话说，使用以下代码创建筛子：

unsigned char primeTbl[] = {0,0,1,1,0,1,0,1,0,0,0,1};
#define isPrime(x) ((x < sizeof(primeTbl) ? primeTbl[x] : isPrimeFn(x))

在不将像a++这样的内容传递到宏中的所有典型警告下，这为您提供了最佳选择，对于“较小”的素数，具有令人惊叹的快速表查找功能，并针对超出范围的素数采用计算方法。

显然，您将编写一个程序来使用其他方法之一生成该查找表 - 您确实不希望手动键入所有内容。

但是，与所有优化问题一样，请“测量，不要猜测！”

*a 这种情况的经典案例是我曾经为嵌入式系统编写的一些三角函数。这是一个竞争性合同投标，该系统的存储空间比CPU强大一点。

我们实际上赢得了合同，因为我们的函数基准值将竞争者甩在了后面。

为什么？因为我们事先将值预先计算到另一台计算机上计算的查找表中。通过巧妙地使用缩减（将输入值降至90度以下）和三角属性（余弦只是正弦的相位移动，而其他三个象限与第一个象限相关），我们将查找表减少到180个条目（每半度一个）。

最好的解决方案是优雅且狡猾的:-)

值得一提的是，以下C代码将为您生成这样的表格，所有小于四百万的素数（283,000个）。

#include <stdio.h>

static unsigned char primeTbl[4000000];

int main (void) {
    int i, j;

    for (i = 0; i < sizeof(primeTbl); i++)
        primeTbl[i] = 1;

    primeTbl[0] = 0;
    primeTbl[1] = 0;
    for (i = 2; i < sizeof(primeTbl); i++)
        if (primeTbl[i])
            for (j = i + i; j < sizeof(primeTbl); j += i)
                primeTbl[j] = 0;

    printf ("static unsigned char primeTbl[] = {");
    for (i = 0; i < sizeof(primeTbl); i++) {
        if ((i % 50) == 0) {
            printf ("\n   ");
        }
        printf ("%d,", primeTbl[i]);
    }
    printf ("\n};\n");
    printf ("#define isPrime(x) "
        "((x < sizeof(primeTbl) ? primeTbl[x] : isPrimeFn(x))\n");

    return 0;
}

如果您将primeTbl表增加到1600万个条目，就会发现足以使素数计数超过一百万（前1,031,130个素数）。现在有一些方法可以减少存储空间，例如仅存储奇数并调整宏以处理此问题，或使用位掩码而不是无符号字符。如果内存可用，我更喜欢简单的算法。
6
+1表示同意，“the first million primes will be unlikely to change in the near future”意为前一百万个质数在不久的将来不太可能改变。哈哈，我不熟悉Project Euler的规则，也许这是不允许的？ - Mark Ransom
2
@Mark：在欧拉计划中没有正式的规则。 - You
@phkahler：这个程序可以处理到436,273,009，也就是一个22MB的列表。如果你从5开始而不是3，并且只存储一半的差值，你可以处理更大的数字，例如304,599,508,537；这将产生一个11GB的列表。 - Charles
@Mark：Project Euler的非正式规则是，你可以花费任意长的时间来编写程序，但它应该在最多一分钟内运行（或类似的时间限制）。然后就是自我监管——如果你能在纸上计算出答案，那可能算作0的运行时间，但如果你编写了一个需要一个小时才能编译的C++模板元程序，那可能不算。只要质数表在许多不同的问题之间共享（这将是这样），个人认为我可能会允许自己分摊生成它的成本。 - Steve Jessop
存在明显的灰色地带，例如，如果我的第一次尝试需要两分钟，但是一旦我看到答案，我意识到有一个优化可以使其更快，那么我是否“作弊”？这取决于我能否写出正确性优化的正式证明。我倾向于认为我应该能够在纸上草拟证明（不使用其他慢程序的结果），证明我的程序找到了正确的答案。因此，一个大的质数表可能应该在每个问题的时间限制内重新生成，而且我记得过去我已经这样做了。 - Steve Jessop
显示剩余2条评论
6
我建议采用分层方法。首先，确保没有小的质因数。试除法使用前20或30个质数可以起到作用，尽管如果你使用聪明的方法，可以通过使用gcds来减少需要的除法次数。这一步可以过滤掉约90%的合数。
接下来，测试数字是否是一个强概率素数（Miller-Rabin测试）基于2。这一步可以消除几乎所有剩余的合数，但是一些罕见的合数可能会通过。
最终证明的步骤取决于你想要达到的规模。如果你愿意在一个小范围内工作，请在允许的最大值上对2-伪素数列表进行二进制搜索。如果是2 ^ 32，则列表将仅有10,403个成员，因此查询应该只需要14个。
如果你想要达到2的64次方，现在可以通过检查数字是否为BPSW伪素数（感谢Jan Feitisma的工作）来实现。 （您还可以下载所有异常的3 GB列表，删除试除法将删除的那些，并编写基于磁盘的二进制搜索。）T. R. Nicely有一个很好的页面，解释如何以相对高效的方式实现此方法。
如果您需要更高的数字，请实现上述方法并将其用作Pocklington-style测试的子程序。这扩展了“较小”的定义；如果您需要更多关于这些方法的信息，请随时询问。
2
作为预计算概念的一种变体，您可以首先廉价地检查候选数字 p 是否可被 2、3、5、7 或 11 整除。如果不行，则声明如果 2p-1 = 1 (mod p)，则 p 为质数。这会失败，但在一亿以内有效，因为我已经测试过了（预计算）。
换句话说，所有基于 2 的小 Fermat 伪质数都能被 3、5、7 或 11 中的一个整除。
编辑：
如 @starblue 所指出的那样，以上是错误的。我的程序有一个错误。我能做的最好的是修改上面的内容：
如果候选数字 p 可被 2、3、5、7 或 11 整除，则声明它为合数；
否则，如果 p 是 {4181921、4469471、5256091、9006401、9863461} 中的一个，则声明它为合数；
否则，如果 p 对于基数 2 和 5 通过 Miller-Rabin 检验，则声明它为质数；
否则，声明它为合数。
我测试了小于10,000,000的整数。也许使用不同的底数对会更好。
请原谅我的错误。
编辑2:
嗯，看起来我要找的信息已经在Miller-Rabin算法的维基百科页面上了，标题为"确定性变体的测试"。
通过"="，我指的是同余。我应该在模p的部分加上括号。已更正。 - President James K. Polk
我还没有找到。也许我会在晚上运行它，看看是否能找到一个。 - President James K. Polk
这不就是使用2作为a的Miller-Rabin算法吗？如果是的话，有趣的是它可以处理到1亿，因为我可能错误地认为Miller-Rabin是概率性的，需要用多个不同的a运行以获得准确性... - afkbowflexin
@GregS：哦，我没看到关于质数除法的那部分。 - afkbowflexin
1
+1：新测试可达到14,709,241。2进制和11进制非常好；5个例外使其能够工作到63,388,033。2进制和733进制稍微好一些，可达到74,927,161。 - Charles
显示剩余6条评论
1
唯一的方法是对自己进行基准测试。当你完成后，将其写下来并在网上发布。
认真点，你已经完成了实现，为什么不自己计时呢？如果你担心可能错过了最快的算法，可以将你的最佳算法发布为一个新问题，看看是否有人可以做得更好。 - Mark Ransom
我可以做到这一点，但我的一些测试实现非常糟糕。我相信我编写的 Miller-Rabin 算法实现方式相当糟糕。实际上，我知道它很糟糕。我想知道最佳情况下的方案，这样我就可以在测试之前只专注于“正确”的实现，而不需要重构每一个“足够好”的实现。 - afkbowflexin
Java库中也有Miller-Rabin算法，适用于BigInteger。 - starblue
0
假设 n < 4669921，它将非常快速：
if ((n == 1) == (n & 1)) return n == 2;
return ((n & 1) & ((n < 6) * 42 + 0x208A2882) >> n % 30 && (n < 49 || (n % 7 && n % 11 && n % 13 && n % 17 && n % 19 && n % 23 && n % 29 && (n < 961 || (n % 31 && n % 37 && n % 41 && n % 43 && n % 47 && n % 53 && n % 59 && n % 61 && n % 67 && (n < 5041 || (n % 71 && n % 73 && n % 79 && n % 83 && n % 89 && n % 97 && n % 101 && n % 103 && n % 107 && (n < 11881 || (n % 109 && n % 113 && n % 127 && n % 131 && n % 137 && n % 139 && n % 149 && n % 151 && n % 157 && (n < 26569 || (n % 163 && n % 167 && n % 173 && n % 179 && n % 181 && n % 191 && n % 193 && n % 197 && n % 199 && (n < 44521 || (n % 211 && n % 223 && n % 227 && n % 229 && n % 233 && n % 239 && n % 241 && n % 251 && n % 257 && (n < 69169 || (n % 263 && n % 269 && n % 271 && n % 277 && n % 281 && n % 283 && n % 293 && n % 307 && n % 311 && (n < 97969 || (n % 313 && n % 317 && n % 331 && n % 337 && n % 347 && n % 349 && n % 353 && n % 359 && n % 367 && (n < 139129 || (n % 373 && n % 379 && n % 383 && n % 389 && n % 397 && n % 401 && n % 409 && n % 419 && n % 421 && (n < 185761 || (n % 431 && n % 433 && n % 439 && n % 443 && n % 449 && n % 457 && n % 461 && n % 463 && n % 467 && (n < 229441 || (n % 479 && n % 487 && n % 491 && n % 499 && n % 503 && n % 509 && n % 521 && n % 523 && n % 541 && (n < 299209 || (n % 547 && n % 557 && n % 563 && n % 569 && n % 571 && n % 577 && n % 587 && n % 593 && n % 599 && (n < 361201 || (n % 601 && n % 607 && n % 613 && n % 617 && n % 619 && n % 631 && n % 641 && n % 643 && n % 647 && (n < 426409 || (n % 653 && n % 659 && n % 661 && n % 673 && n % 677 && n % 683 && n % 691 && n % 701 && n % 709 && (n < 516961 || (n % 719 && n % 727 && n % 733 && n % 739 && n % 743 && n % 751 && n % 757 && n % 761 && n % 769 && (n < 597529 || (n % 773 && n % 787 && n % 797 && n % 809 && n % 811 && n % 821 && n % 823 && n % 827 && n % 829 && (n < 703921 || (n % 839 && n % 853 && n % 857 && n % 859 && n % 863 && n % 877 && n % 881 && n % 883 && n % 887 && (n < 822649 || (n % 907 && n % 911 && n % 919 && n % 929 && n % 937 && n % 941 && n % 947 && n % 953 && n % 967 && (n < 942841 || (n % 971 && n % 977 && n % 983 && n % 991 && n % 997 && n % 1009 && n % 1013 && n % 1019 && n % 1021 && (n < 1062961 || (n % 1031 && n % 1033 && n % 1039 && n % 1049 && n % 1051 && n % 1061 && n % 1063 && n % 1069 && n % 1087 && (n < 1190281 || (n % 1091 && n % 1093 && n % 1097 && n % 1103 && n % 1109 && n % 1117 && n % 1123 && n % 1129 && n % 1151 && (n < 1329409 || (n % 1153 && n % 1163 && n % 1171 && n % 1181 && n % 1187 && n % 1193 && n % 1201 && n % 1213 && n % 1217 && (n < 1495729 || (n % 1223 && n % 1229 && n % 1231 && n % 1237 && n % 1249 && n % 1259 && n % 1277 && n % 1279 && n % 1283 && (n < 1661521 || (n % 1289 && n % 1291 && n % 1297 && n % 1301 && n % 1303 && n % 1307 && n % 1319 && n % 1321 && n % 1327 && (n < 1852321 || (n % 1361 && n % 1367 && n % 1373 && n % 1381 && n % 1399 && n % 1409 && n % 1423 && n % 1427 && n % 1429 && (n < 2053489 || (n % 1433 && n % 1439 && n % 1447 && n % 1451 && n % 1453 && n % 1459 && n % 1471 && n % 1481 && n % 1483 && (n < 2211169 || (n % 1487 && n % 1489 && n % 1493 && n % 1499 && n % 1511 && n % 1523 && n % 1531 && n % 1543 && n % 1549 && (n < 2411809 || (n % 1553 && n % 1559 && n % 1567 && n % 1571 && n % 1579 && n % 1583 && n % 1597 && n % 1601 && n % 1607 && (n < 2588881 || (n % 1609 && n % 1613 && n % 1619 && n % 1621 && n % 1627 && n % 1637 && n % 1657 && n % 1663 && n % 1667 && (n < 2785561 || (n % 1669 && n % 1693 && n % 1697 && n % 1699 && n % 1709 && n % 1721 && n % 1723 && n % 1733 && n % 1741 && (n < 3052009 || (n % 1747 && n % 1753 && n % 1759 && n % 1777 && n % 1783 && n % 1787 && n % 1789 && n % 1801 && n % 1811 && (n < 3323329 || (n % 1823 && n % 1831 && n % 1847 && n % 1861 && n % 1867 && n % 1871 && n % 1873 && n % 1877 && n % 1879 && (n < 3568321 || (n % 1889 && n % 1901 && n % 1907 && n % 1913 && n % 1931 && n % 1933 && n % 1949 && n % 1951 && n % 1973 && (n < 3916441 || (n % 1979 && n % 1987 && n % 1993 && n % 1997 && n % 1999 && n % 2003 && n % 2011 && n % 2017 && n % 2027 && (n < 4116841 || (n % 2029 && n % 2039 && n % 2053 && n % 2063 && n % 2069 && n % 2081 && n % 2083 && n % 2087 && n % 2089 && (n < 4405801 || (n % 2099 && n % 2111 && n % 2113 && n % 2129 && n % 2131 && n % 2137 && n % 2141 && n % 2143 && n % 2153)))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))));
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接