数论之大素数检测与大整数分解(10^18)

大素数检测

大素数检测常用的方法为miller_rabin。网上讲这个方法的文章博客已经很多了，我在这里就只转载一篇别人的文章。

原文出处：http://www.cppblog.com/zoyi-zhang/archive/2008/09/23/62572.aspx

大素数的检验
费马小定理：a^(p-1) mod p = 1（p是素数＆＆a<p＆＆a>0）

 首先我们证明这样一个结论：如果p是一个素数的话，那么对任意一个小于p的正整数a，a, 2a, 3a, …, (p-1)a
除以p的余数正好是一个1到p-1的
排列。例如，5是素数，3, 6, 9, 12除以5的余数分别为3, 1, 4, 2，正好就是1到4这四个数。
反证法，假如结论不成立的话，那么就是说有两个小于p的正整数m和n使得na和ma除以p的余数相同。不妨假设n>m，
则p可以整除a(n-m)。但p是素
数，那么a和n-m中至少有一个含有因子p。这显然是不可能的，因为a和n-m都比p小。
用同余式表述，我们证明了：
(p-1)! ≡ a * 2a * 3a * … * (p-1)a (mod p)
    也即：
(p-1)! ≡ (p-1)! * a^(p-1) (mod p)
    两边同时除以(p-1)!，就得到了我们的最终结论：
1 ≡ a^(p-1) (mod p)

费马小定里的欧拉推广：a^φ(m) ≡ 1 (mod m)（其中φ(m)为与m互质的数的个数）

证明与费马小定理类似

但是费马小定理的逆命题并不正确，即，当满足a^(p-1) mod p = 1的数p不一定是素数，例如p=341，a=2，
而此时p=11＊13

后来，人们又发现了561, 645, 1105等数都表明a=2时Fermat小定理的逆命题不成立。虽然这样的数不多，
但不能忽视它们的存在。于是，人们把所
有能整除2^(n-1)-1的合数n叫做伪素数(pseudoprime)，意思就是告诉人们这个素数是假的。

不满足2^(n-1) mod n = 1的n一定不是素数；如果满足的话则多半是素数。这样，一个比试除法效率更高的
素性判断方法出现了：制作一张伪素数
表，记录某个范围内的所有伪素数，那么所有满足2^(n-1) mod n = 1且不在伪素数表中的n就是素数。之所以
这种方法更快，是因为我们可以使用
二分法快速计算2^(n-1) mod n 的值，这在计算机的帮助下变得非常容易；在计算机中也可以用二分查找有序
数列、Hash表开散列、构建Trie树等
方法使得查找伪素数表效率更高。

有人自然会关心这样一个问题：伪素数的个数到底有多少？换句话说，如果我只计算2^(n-1) mod n的值，事先
不准备伪素数表，那么素性判断出
错的概率有多少？研究这个问题是很有价值的，毕竟我们是OIer，不可能背一个长度上千的常量数组带上考场。
统计表明，在前10亿个自然数中共
有50847534个素数，而满足2^(n-1) mod n = 1的合数n有5597个。这样算下来，算法出错的可能性约为0.00011。
这个概率太高了，如果想免去建
立伪素数表的工作，我们需要改进素性判断的算法。

最简单的想法就是，我们刚才只考虑了a=2的情况。对于式子a^(n-1) mod n，取不同的a可能导致不同的结果。一
个合数可能在a=2时通过了测试，
但a=3时的计算结果却排除了素数的可能。于是，人们扩展了伪素数的定义，称满足 a^(n-1) mod n = 1的合数n叫
做以a为底的伪素数(pseudoprime to base a)
。前10亿个自然数中同时以2和3为底的伪素数只有1272个，这个数目不到刚才的1/4。这告诉我们如果同时验证a=2
和a=3两种情况，算法出错的概率
降到了0.000025。容易想到，选择用来测试的a越多，算法越准确。通常我们的做法是，随机选择

若干个小于待测数的正整数作为底数a进行若干次测试，只要有一次没有通过测试就立即把这个数扔回合数的世界。
这就是Fermat素性测试。

人们自然会想，如果考虑了所有小于n的底数 a，出错的概率是否就可以降到0呢？没想到的是，居然就有这样的合数，
它可以通过所有a的测试
（这个说法不准确，详见我在地核楼层的回复）。 Carmichael第一个发现这样极端的伪素数，他把它们称作Carmichael数。
你一定会以为这样的数一定很大。错。第一个Carmichael 数小得惊人，仅仅是一个三位数，561。前10亿个自然数中
Carmichael数也有600个之多。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的 Miller-Rabin素性测试算法。新的测
试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为
x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能
被p整除(此时x =p-1)。

这就是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。那么
我们需要计算的东西就变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。
如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足
a^(d * 2^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。这样，Fermat小定理加强为如下形式：
尽可能提取因子 2，把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) 
（注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）
Miller-Rabin 素性测试同样是不确定算法，我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素
数(strong pseudoprime)。
第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。

 Miller- Rabin算法的代码也非常简单：计算d和r的值（可以用位运算加速），然后二分计算a^d mod n的值
，最后把它平方r次。程序的代码比想像中的更简单，我写一份放在下边。虽然我已经转C了，但我相信还有很多
人看不懂C语言。我再写一次Pascal 吧。函数IsPrime返回对于特定的底数a，n是否是能通过测试。如果函数返回
False，那说明n不是素数；如果函数返回True，那么n极有可能是素数。注意这个代码的数据范围限制在longint，
你很可能需要把它们改成int64或高精度计算。
    我们下面来演示一下上面的定理如何应用在Fermat素性测试上。前面说过341可以通过以2为底的Fermat测试，
因为 2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于
1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，
面具后面真实的嘴脸显现了出来


  对于大数的素性判断，目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取，但当待测数不太大时，
选择测试底数就有一些技巧了。比如，如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足
够了。当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行
测试，所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内
唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在OI中非常实用。通常认为，Miller-Rabin素性测试
的正确率可以令人接受，随机选取 k个底数进行测试算法的失误率大概为4^(-k)。

同时，我将自己改的模板贴出来(参考网上多个模板)：

//Miller_Rabin素数测试法，质数返回1，否则返回0
//若用下面prm[]里的数进行测试，10^16以内只有唯一强伪素数46,856,248,255,981

/*#include <ctime>*/
const int MT = 5;      /*prm[]的元素个数，也即是需要测试的底的个数*/    
const int64 SPE = 46856248255981;             //强伪素数
typedef long long int64;
int64 prm[5] = {2, 3, 7, 61, 24251};

//要判断n是不是质数时，调用语句为miller_rabin(n, MT);

int64 mul_mod(int64 p, int64 q, int64 mod)
{
     int64 ret = 0;
     p %= mod;
     while (q){
         if (q & 1){
             ret += p;
             if (ret >= mod) ret -= mod;
         }
         p <<= 1; q >>= 1;
         if (p >= mod) p -= mod;
     }
     return ret;
}

int64 pow_mod(int64 p, int64 n, int64 mod)
{
    int64 ret = 1;
    p %= mod;
    while (n > 0){
        if (n & 1) ret = mul_mod(ret, p, mod);
        n >>= 1;
        p = mul_mod(p, p, mod);
    }
    return ret;
}

bool Wintess(int64 aa, int64 n, int64 mod)      /*是以aa为底的伪质数返回0,是合数返回1*/
{
    int t = 0;
    while ((n & 1) == 0){               /*n = d * 2^t; n = d;*/
        n >>= 1;
        ++ t;
    }

    int64 x = pow_mod(aa, n, mod), y;
    while (t --){
        y = mul_mod (x, x, mod);
        if (y == 1 && x != 1 && x != mod-1) return 1;
        x = y;
    }
    return (y != 1);
}

bool miller_rabin(int64 n, int tim)
{
    if (n == 2) return true;
    if (n == 1 || (n & 1) == 0 || n == SPE) return false;

    for (int i = 0; i < tim; ++ i){
        int64 aa = prm[i];           /*也可以用随机化的aa，在前面写srand(time(NULL));在这里写int64 aa=rand()%(n-1)+1;*/
        if (aa == n) return true;
        if (Wintess(aa, n-1, n)) return false;
    }
    return true;
}

大整数分解

一般我们都是使用试除法来进行分解质因数，但此方法对64位整数而言效率就太低了，所以我们需要新的方法。我在网上有看到两种大整数分解的方法，一种叫做整数分解费马方法，不过我觉得它并没有比试除法快到哪去......

另外一种就是Pollard rho整数分解法。这个分解方法网上也没看见太详细的讲解，一般都是直接给出了分解过程，所以我到现在也不太理解原理- -。

不过，还是转载一下别人写的东西吧，出处：http://echocx.diandian.com/post/2011-10-02/5471165

Pollard rho
原理：设n为待分解的大整数，用某种方法生成a和b，计算p=gcd(a-b,n),直到p不为1或a,b出现循环时为止，若p=n，则说明n是一个素数，否则p为n的一个约数。
算法步骤：选取一个小的随机数x1，迭代生成x[i] = x[i-1]^2+c，一般去c=1，若序列出现循环则退出，计算p=gcd(x[i-1]-x[i],n)，若p=1则返回上一步继续迭代，否则跳出迭代过程。若p=n，则n为素数，否则p为n的一个约数，并递归分解p和n/p。
可以在θ（sqrt(p)）的期望时间内找到n的一个小因子p。但对于因子很少，因子值很大的大整数n，该方法依然不是很有效。

只好弄了一个模板出来。

//Pollard分解质因数方法
//分解后的质因数全部存在
//时间复杂度一般为O(n^(1/4))

        /*调用miller_rabin*/
const int MT = 5;
int64 mul_mod(int64 p, int64 q, int64 mod)
int64 pow_mod(int64 p, int64 n, int64 mod)
bool miller_rabin(n, MT)
            /*end*/

typedef long long int64;
const int TIM = 240;       /*控制循环次数*/
int cnt;                /*分解后质因数的个数*/
int64 fac[];            /*存分解后的所有质因数*/

//要分解n，调用语句为cnt = 0; get_prime(n, TIM);

int64 gcd(int64 a, int64 b)
{
    return b ? gcd(b, a%b) : a;
}

int64 Pollard(int64 n, int c)
{
    srand(time(NULL));
    int64 i = 1, k = 2, x = rand()%n;
    int64 y = x;
    while (1){
        ++ i;
        x = (mul_mod(x, x, n) + c) % n;
        int64 d = gcd(y-x, n);
        if (d > 1 && d < n) return d;
        if (y == x) return n;
        if (i == k){
            y = x;
            k <<= 1;
        }
    }
}

void get_prime(int64 n, int c)
{
    if (n == 1) return;
    if (miller_rabin(n, MT)){
        fac[cnt++] = n;
        return;
    }

    int64 m = n;
    while (m == n) m = Pollard(m, c--);
    get_prime(m, c);
    get_prime(n/m, c);
}

Ps:因为pollard算法的时间复杂度和rand()生成的数有关，所以有可能因为控制循环次数的变量太大而导致TLE，比如有一道题2s的时限，有时候交是250msAC，有时候是TLE- -......

1、POJ 2191 Mersenne Composite Numbers

题意：给定n，对所有小于n的素数i，如果(2^1)-1是合数，就将其分解质因数并按题目要求格式输出。n <= 63。

解法：用快速幂求(2^1)-1，用miller_rabin判定是否质数，若不是，则用Pollard分解。模板题。

tag:math, miller_rabin, Pollard

Ps：在网上看见说，如果代码里面写了srand(time(NULL))，并且交POJ的G++就会错。不知道是不是真的，反正我下面这段代码交G++就RE，交C++就AC。

/*
 * Author:  Plumrain
 * Created Time:  2013-10-23 11:41
 * File Name: math-POJ-2191.cpp
 */
#include <iostream>
#include <cstdio>
#include <cstring>
#include<algorithm>
#include <vector>
#include <ctime>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
#define PB push_back
const int MT = 5;
const int TIM = 240; 
typedef vector<int> VI;
typedef long long int64;
const int64 SPE = 46856248255981;

VI cur;
int cnt;
int64 prm[5] = {2, 3, 7, 61, 24251};
int64 cha[19] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67};
int64 fac[100];

int64 gcd(int64 a, int64 b)
{
    return b ? gcd(b, a%b) : a;
}

int64 Mypow(int64 p, int64 n)
{
    int64 ret = 1;
    while (n){
        if (n & 1) ret *= p;
        n >>= 1;
        p *= p;
    }
    return ret;
}

int64 mul_mod(int64 p, int64 q, int64 mod)
{
     int64 ret = 0;
     p %= mod;
     while (q){
         if (q & 1){
             ret += p;
             if (ret >= mod) ret -= mod;
         }
         p <<= 1; q >>= 1;
         if (p >= mod) p -= mod;
     }
     return ret;
}

int64 pow_mod(int64 p, int64 n, int64 mod)
{
    int64 ret = 1;
    p %= mod;
    while (n > 0){
        if (n & 1) ret = mul_mod(ret, p, mod);
        n >>= 1;
        p = mul_mod(p, p, mod);
    }
    return ret;
}

bool Wintess(int64 aa, int64 n, int64 mod)
{
    int t = 0;
    while ((n & 1) == 0){
        n >>= 1;
        ++ t;
    }

    int64 x = pow_mod(aa, n, mod), y;
    while (t --){
        y = mul_mod (x, x, mod);
        if (y == 1 && x != 1 && x != mod-1) return 1;
        x = y;
    }
    return (y != 1);
}

bool miller_rabin(int64 n, int tim)
{
    if (n == 2) return true;
    if (n == 1 || (n & 1) == 0 || n == SPE) return false;

    for (int i = 0; i < tim; ++ i){
        int64 aa = prm[i];
        if (aa == n) return true;
        if (Wintess(aa, n-1, n)) return false;
    }
    return true;
}

int64 Pollard(int64 n, int c)
{
    srand(time(NULL));
    int64 i = 1, k = 2, x = rand()%n;
    int64 y = x;
    while (1){
        ++ i;
        x = (mul_mod(x, x, n) + c) % n;
        int64 d = gcd(y-x, n);
        if (d > 1 && d < n) return d;
        if (y == x) return n;
        if (i == k){
            y = x;
            k <<= 1;
        }
    }
}

void get_prime(int64 n, int c)
{
    if (n == 1) return;
    if (miller_rabin(n, MT)){
        fac[cnt++] = n;
        return;
    }

    int64 m = n;
    while (m == n) m = Pollard(m, c--);
    get_prime(m, c);
    get_prime(n/m, c);
}

bool cmp(int64 a, int64 b)
{
    return a < b;
}

int main()
{
    int n;
    while (scanf ("%d", &n) != EOF){
        cur.clear();
        for (int i = 0; cha[i] < n; ++ i)
            if (!miller_rabin(Mypow(2, cha[i])-1, MT))
                cur.PB (cha[i]);

        for (int i = 0; i < cur.size(); ++ i){
            cnt = 0;
            int64 tmp = Mypow(2, cur[i]) - 1;
            get_prime(tmp, TIM);
            sort(fac, fac+cnt, cmp);
            for (int j = 0; j < cnt; ++ j)
                printf ("%lld %c ", fac[j], j==cnt-1?'=':'*');
            printf ("%lld = ( 2 ^ %d ) - 1
", tmp, cur[i]); 
        }
    }
    return 0;
}

2、POJ 1811 Prime Test

题意：给定n，判断是不是素数，若不是，输出最小的质因子。

解法：模板题。

tag:math, miller_rabin, Pollard

/*
 * Author:  Plumrain
 * Created Time:  2013-10-24 09:16
 * File Name: math-POJ-1811.cpp
 */
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<ctime>

using namespace std;

const int MT = 5;
const int TIM = 240;
typedef long long int64;
const int64 SPE = 46856248255981;

int64 ans;
int64 prm[5] = {2, 3, 7, 61, 24251};

int64 min(int64 a, int64 b)
{
    return a < b ? a : b;
}

int64 gcd(int64 a, int64 b)
{
    return b ? gcd(b, a%b) : a;
}

int64 mul_mod(int64 p, int64 q, int64 mod)
{
     int64 ret = 0;
     p %= mod;
     while (q){
         if (q & 1){
             ret += p;
             if (ret >= mod) ret -= mod;
         }
         p <<= 1; q >>= 1;
         if (p >= mod) p -= mod;
     }
     return ret;
}

int64 pow_mod(int64 p, int64 n, int64 mod)
{
    int64 ret = 1;
    p %= mod;
    while (n > 0){
        if (n & 1) ret = mul_mod(ret, p, mod);
        n >>= 1;
        p = mul_mod(p, p, mod);
    }
    return ret;
}

bool Wintess(int64 aa, int64 n, int64 mod)
{
    int t = 0;
    while ((n & 1) == 0){
        n >>= 1;
        ++ t;
    }

    int64 x = pow_mod(aa, n, mod), y;
    while (t --){
        y = mul_mod (x, x, mod);
        if (y == 1 && x != 1 && x != mod-1) return 1;
        x = y;
    }
    return (y != 1);
}

bool miller_rabin(int64 n, int tim)
{
    if (n == 2) return true;
    if (n == 1 || (n & 1) == 0 || n == SPE) return false;

    for (int i = 0; i < tim; ++ i){
        int64 aa = prm[i];
        if (aa == n) return true;
        if (Wintess(aa, n-1, n)) return false;
    }
    return true;
}

int64 pollard(int64 n, int c)
{
    srand(time(NULL));
    int64 i = 1, k = 2, x = rand()%n;
    int64 y = x;
    while (1){
        ++ i;
        x = (mul_mod(x, x, n) + c) % n;
        int64 d = gcd(y-x, n);
        if (d > 1 && d < n) return d;
        if (y == x) return n;
        if (i == k){
            y = x;
            k <<= 1;
        }
    }
}

void get_prime(int64 n, int c)
{
    if (n == 1) return;
    if (miller_rabin(n, MT)){
        ans = min(ans, n);
        return;
    }

    int64 m = n;
    while (m == n) m = pollard(m, c--);
    get_prime(m, c);
    get_prime(n/m, c);
}


int main()
{
    int64 tmp = (int64)1 << 60;
    int T;
    scanf ("%d", &T);
    while (T--){
        int64 a;
        scanf ("%lld", &a);
        if (miller_rabin(a, MT))
            printf ("Prime
");
        else{
            ans = tmp;
            get_prime(a, TIM);            
            printf ("%lld
", ans);
        }
    }
    return 0;
}

3、POJ 2429 GCD & LCM Inverse

题意：对于两个正整数a, b，给出gcd(a, b)和lcm(a, b)，求a和b。多组数据的话，输出a + b最小的一组。

题解祥见http://www.cnblogs.com/plumrain/p/number_theory.html