一文说尽32bit素数检测（幂模运算，米勒拉宾算法，素数检测）

#include "Prime1.h" bool DemoPrimeQ_U32(unsigned int P) { /* if ((P == 0) || (P == 1)) return false; bool isPrime = true; unsigned int k; unsigned int LIMIT = (unsigned int)sqrt(P); for (k = 2; k <= LIMIT; k++) { if (P % k == 0) { isPrime = false; break; } } return isPrime; */ /* 以上这些注释掉的代码是基本算法，先读懂它。 下面的代码使用素数2和3两个 筛子，把速度加快约3倍。 */ switch (P) { case 0: case 1: case 4: return false; case 2: case 3: case 5: return true; default: break; } /* 用2和3两个筛子，那就是模6步进。先筛2和3，如果整除，那就是合数，如果不 整除，那么其后的含有因子2和3的数都没有必要进行测试。6*n、6*n + 1、 6*n + 2、6*n + 3、6*n + 4、6*n + 5，n从0开始循环，包括了所有数， 因为 只有6*n + 1和6*n + 5不含有因子2和3。但是n不能从0开始循环， 否则遇到的 第一个数就是6*0 + 1 = 1，1整除任何数，程序会失败。 n必须从1开始循环， 这就会漏掉素数5，所以先把5也筛掉。 */ if (P % 2 == 0) return false; if (P % 3 == 0) return false; if (P % 5 == 0) return false; /* 如果P就是2，3或者5，前面的switch语句直接判定结果返回。 */ bool isPrime = true; /* 先假定是素数，一旦测试为合数，则置false。 */ unsigned int LIMIT = (unsigned int)sqrt(P) / 6; unsigned int k; unsigned int n; /* 下边的代码，一看就懂，但是要思考以下几个问题。 1.循环里的最大的数6*LIMIT + 5确保大于等于根号P的向下取整吗？ 因为至少 要循环到根号P的向下取整这个数才行。 2.循环里的最大的数6*LIMIT + 5确保严格小于P吗？因为万一某个6*k + 1或者 6*k + 5恰好等于P，if语句成立判定为合数，可这是错的。 自己整除自己不能 作为依据，不能因为11整除11判定11是合数吧？所以必须确保6*LIMIT + 5严格 小于P。 3.如果P小于等于35，那么LIMIT = 0， 不满足循环的初始条件循环不会执行， 此时这个代码正确吗？ */ for (k = 1; k <= LIMIT; k++) { n = 6 * k; if (P % (n + 1) == 0) { isPrime = false; break; } if (P % (n + 5) == 0) { isPrime = false; break; } } return isPrime; /* 6是怎样被判定为合数的？第一个if (6 % 2 == 0)就判定了。 7是怎样被判定为素数的？试除法试除到根号7的向下取整就可以， 也就是说测 试到2就行，第一个if (7 % 2 == 0)不成立，按理说就可以结束了。 但是程序 是继续运行，先假定了isPrime为true（素数），计算出LIMIT = 0，for循环不 满足初始条件，跳过，没有代码改变过isPrime的初始假定，7被判定为素数。 31是怎样被判定为素数的？同理，试除到根号31的向下取整5就可以，前面三个 if刚好筛到5，下面的for循环又不满足初始条件，跳过，31被判定为素数。 3613是怎样被判定为素数的？2，3，5都不是因子，LIMIT = 10，k从1开始循环 到10，测试以下这些数：7，11，13，17，19，23，25，29，31，35，37，41， 43，47，49，53，55，59，61，65，全部不满足， 65大于根号3613的向下取整 值60，3613被判定为素数。 4087怎样被判定为合数的？2，3，5都不是因子，LIMIT = 10，k从1开始循环到 10，测试以下这些数：7，11，13，17，19，23，25，29，31，35，37，41， 43，47，49，53，55，59，61，65，当测试到61时满足条件，isPrime被改写为 false，终止循环，判定合数返回。 */ } unsigned int PowerMod_U32(unsigned int x, unsigned int N, unsigned int m) { if (m == 0) { std::cout << "除数为零，程序强行终止！\n"; exit(0); } if (m == 1) return 0; /* 模m不能为0。 模m是1，余数必定为0。 */ unsigned int k; unsigned long Max; unsigned int MaskCode = 1; unsigned long long X = x; unsigned long long Remainder = 1; /* 如果N = 0，对任意x，x^N = 1，模m是结果1。 注意m = 1被前面的if语句事先 滤除了，任何数模1结果是0。如果N较小，直接使用循环法更快。 */ if (N == 0) Remainder = 1; else if (N <= 12) { for (k = 0; k < N; k++) { Remainder *= X;//Remainder初始值为1，每次乘以X，取余，循环N次。 Remainder %= m; /* unsigned int的乘法需要用 unsigned long long来处理中间结果才不 会溢出，然后把乘法的结果求余数，余数必定小于m， 而m不会超出32 bit，无论循环多少次都不会溢出。 如果Remainder被设置为unsigned int类型，那么Remainder *= X只有低32bit乘法结果，高32bit被编译 器丢弃了。 */ } } else { _BitScanReverse(&Max, N); /* _BitScanReverse这个函数取得N的二进制里最高有效位的次数，基于下标0 计算，这是cpu硬件完成。 MaskCode初始化为1，每次循环左移一次， MaskCode的二进制只有一个1， 位与运算&用来测试N的二进制对应位是否为1。 示例：N = 18 =（10010）二进制，计算出Max = 4。 k if测试 每次循环后result(基于mod m) X(基于mod m) 0 假 1 X^2 1 真 X^2 X^4 2 假 X^2 X^8 3 假 X^2 X^16 4 真 X^18 X^32 注意到没，在这个例子里，明明计算到X^16(mod m)就够用了，可是平白无 故地多计算了一个X^32。其次最后一次循环的if测试必定成立， 多执行了 这次测试。最后，MaskCode <<= 1和k++，怎么都感觉都有一个是多余的。 */ for (k = 0; k <= Max; k++) { if (N & MaskCode) { Remainder *= X; Remainder %= m; } MaskCode <<= 1; X *= X; X %= m; } } return (unsigned int)Remainder; //余数必定小于模m，而m又必定小于2^32，所以强制转换不会丢失数据。 } /* 下面这个幂模函数改正了上一个存在的各种不足，提速15%。 */ void PowerMod_U32(unsigned int x, unsigned int N, \ unsigned int m, unsigned int& R) { if (m <= 1) { if (m == 0) { std::cout << "除数为零，程序强行终止！\n"; exit(0); } else { R = 0; return; } } if (N <= 1) { if (N == 0) { R = 1; return; } else { R = x % m; return; } } unsigned int k; unsigned long Max; unsigned long long Product; R = 1; if (N <= 12) { for (k = 0; k < N; k++) { Product = (unsigned long long)R * x; _udiv64(Product, m, &R); } } else { _BitScanReverse(&Max, N); Max = (1 << (Max - 1)); for (k = 1; k <= Max; k <<= 1) { if (N & k) { Product = (unsigned long long)R * x; _udiv64(Product, m, &R); } Product = (unsigned long long)x * x; _udiv64(Product, m, &x); } Product = (unsigned long long)R * x; _udiv64(Product, m, &R); /* 这段代码和上一段代码有何不同？ 1.循环变量k兼做位与测试，省去了一条MaskCode <<= 1的操作。 2.在上一段代码中Remainder是64bit变量，求模运算符%是64bit除以64bit 的运算；在这段代码中R是32bit变量，_udiv64是64bit除以32bit的运算， 速度稍快。 3.在上一段代码中，最后一次循环if测试必定为真， 里面的运算必定要执 行。但是if下面的两条语句被多余计算了（因为没有再下一次循环了）。 */ } } bool StdMillerRabin_U32(unsigned int Witness, unsigned int Composite) { if (Witness < 2) { std::cout << "Witness必须大于等于2！\n"; exit(0); } if (Composite <= 2) { if (Composite == 2) return false; else return true; } if ((Composite & 1) == 0) return true; /* 进行规则检查。 1.Witness必须大于等于2。 2.Composite = 0、1视同合数（素数从2开始定义）， 2是唯一偶素数，直接判 定返回。如果Composite是偶数，直接判定返回。 */ if (Composite == Witness) return false;//无法断言 else { if (Composite % Witness == 0) return true; } /* 上面的if-else代码需要好好理解，这是烧脑所在。 1.如果Composite和Witness相等（第一个if部分），这是不允许的。 米勒拉宾 算法的起始条件要求这两个数不能相同。在这种情况下，直接返回false，表示 无法断言Composite是合数。米勒拉宾算法是一种概率型的算法，在这种情况下 需要更换另外一些不同的证据数来重新检测。 2.如果Composite和Witness不相等（else部分）， 并且Witness整除Composite （这意味着Composite严格大，因为它们不相同的else假设）， 此时Composite 必定是合数，因为Composite有一个大于等于2的因子并且这个因子不是本身。 3.如果Composite和Witness不相等（else部分）， 又假设else里面的if语句不 成立，即Witness不整除Composite，那正是这个代码正常的流程会向下继续。 */ /* 如果没有这一组if-else代码会发生什么问题呢？假定要检测61，而证据数也是 61，在接下来的一系列运算中求得的余数都是零， 满足米勒拉宾算法的结论性 条件，会给出61是合数的荒谬结论， 这个错误的根源就在于证据数和被检测数 相同。再如，用证据数7来检测49，显然7整除49， 49可以直接被判定为合数而 无需继续。 */ unsigned int C_