Miller-Rabin

费马小定理：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)。














#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <algorithm>
typedef long long ll;
#define Time 15 //随机算法判定次数，Time越大，判错概率越小
using namespace std;
ll n,ans,factor[10001];//质因数分解结果（刚返回时是无序的）
ll tol;//质因数的个数，数组下标从0开始
//****************************************************************
// Miller_Rabin 算法进行素数测试
//速度快，而且可以判断 <2^63的数
//****************************************************************
long long mult_mod(ll a,ll b,ll c)//计算 (a*b)%c.   a,b都是ll的数，直接相乘可能溢出的
{
    a%=c;//                           利用二分思想减少相乘的时间
    b%=c;
    ll ret=0;
    while(b)
    {
        if(b&1)
        {
            ret+=a;
            ret%=c;
        }
        a<<=1;
        if(a>=c)a%=c;
        b>>=1;
    }
    return ret;
}
ll pow_mod(ll a,ll x,ll n)//a^x%n
{
    if(x==1)return a%n;
    a%=n;
    ll tmp=a;
    ll ret=1;
    while(x)
    {
        if(x&1) ret=mult_mod(ret,tmp,n);
        tmp=mult_mod(tmp,tmp,n);
        x>>=1;
    }
    return ret;
}
//以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
//一定是合数返回true,不一定返回false
//二次探测
bool check(ll a,ll n,ll x,ll t)
{
    ll ret=pow_mod(a,x,n);
    ll last=ret;
    for(int i=1; i<=t; i++)
    {
        ret=mult_mod(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1) return true;//合数
        last=ret;
    }
    if(ret!=1) return true;
    return false;
}

// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
bool Miller_Rabin(ll n)
{
    if(n<2)return false;
    if(n==2||n==3||n==5||n==7)return true;
    if(n==1||(n%2==0)||(n%3==0)||(n%5==0)||(n%7==0)) return false;//偶数
    ll x=n-1;
    ll t=0;
    while((x&1)==0)
    {
        x>>=1;
        t++;
    }
    for(int i=0; i<Time; i++)
    {
        ll a=rand()%(n-1)+1;//rand()需要stdlib.h头文件
        if(check(a,n,x,t))
            return false;//合数
    }
    return true;
}
//************************************************
//pollard_rho 算法进行质因数分解
//************************************************
ll gcd(ll a,ll b)
{
    if(a==0)return 1;
    if(a<0) return gcd(-a,b);
    while(b)
    {
        long long t=a%b;
        a=b;
        b=t;
    }
    return a;
}
ll Pollard_rho(ll x,ll c)
{
    ll i=1,k=2;
    ll x0=rand()%x;
    ll y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        long long d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k)
        {
            y=x0;
            k+=k;
        }
    }
}
//对n进行素因子分解
void findfac(ll n)
{
    if(Miller_Rabin(n))//素数
    {
        factor[tol++]=n;
        return;
    }
    ll p=n;
    while(p>=n) p=Pollard_rho(p,rand()%(n-1)+1);
    findfac(p);//递归调用
    findfac(n/p);
}
int main()
{
    int T=1;
    //srand(time(NULL));加上RE不懂
    //scanf("%d",&T);
    while(T--)
    {
        n=341;
        //scanf("%lld",&n);//(n>=2)
        /*if(n==1)
        {
            printf("1\n");
            continue;
        }*/
        if(Miller_Rabin(n))
        {
            printf("Prime\n");
            continue;
        }
        tol=0;
        findfac(n);//对n分解质因子
        ll ans=factor[0];
        for(int i=1; i<tol; i++)
            if(factor[i]<ans)
                ans=factor[i];
        for(int i=0;i<tol;i++)
        {
            printf("%lld\n",factor[i]);
        }
    }
    return 0;
}