BZOJ2082[POI2010] Divine divisor

BZOJ2082[POI2010] Divine divisor

Description

Tz耍畸形，在寂寞的时候玩一个游戏，他随便找出一个n，然后算出n的所有因子，最后找出一个最大的k，即有一个因子d的k次方为n的因子，那个因子d就是非凡因子啦 。比如48他的非凡因子d就是2，k最大为4，因为16也是48的因子。一个整数的非凡因子可能不止一个，比如6就有3个：2，3，6（k最大是1）。

Input

有两行，第一行给出一个整数m(1<=m<=600)，第二行是由空格分开的m个数ai(2<=ai<=10^18(10的18次方))，n就是ai的乘积（$n=a1*a2*a3*....*am$）。

Output

也有两行（Tz真畸形），第一行为最大的k，第二行为非凡因子d的个数。

Sample Input

【样例输入1】

3

4 3 4 （$n=4*3*4=48$）

【样例输入2】

1

6

Sample Output

【样例输出1】

4 （k最大为4）

1 （一个非凡因子，2）

【样例输出2】

1

3

Solution：

题目的意思非常容易抽离出来：将此大整数分解质因数，求次数最多的质因子次数及个数。（表述略有不通）

刚开始拿到这道题时，甚至一直在想将它直接分解，以至于一直没有想法。这里要反思一下：有点缺乏举一反三的能力，不过后面想的还挺顺的，没有什么障碍。

于是与小数字质因数分解相同的，我们先造出$10^6$以内的素数表。然后的思路就比较顺畅了：剩下的数如果不是1，那么有三种情况：（设两个大于$10^6$的质数为p和q）

• p：本身是一个大质数
• $p^2$：完全平方数
• $p*q$：两个大质数的乘积

第两种很容易就能解决，然后我们对这些大质数两两取gcd，能除的就将它除掉，归入已确定的表中。剩下的无法被整除的第三种情况就对答案没有什么影响了（因为组成这个数的这两个数都肯定不在我们确定了的表中），唯一要确定的就是当这种质数对答案贡献时要判断到底是一个质数还是两个。朴素地判断复杂度为$O(\sqrt{n})$是肯定要超时的。于是就需要一个高效的素性测试：Miller-Rabin素性测试算法。

讲这个之前先提一句：第二个答案就是$2^n-1$，并且要高精，简单带过。

以下部分引用自：数论部分第一节：素数与素性测试 By Matrix67。

素数的五个性质：

1. 素数的个数无限多（不存在最大的素数）

证明：反证法，假设存在最大的素数$P$，那么我们可以构造一个新的数$2 * 3 * 5 * 7 * … * P + 1$（所有的素数乘起来加1）。显然这个数不能被任一素数整除（所有素数除它都余1），这说明我们找到了一个更大的素数。

2. 存在任意长的一段连续数，其中的所有数都是合数（相邻素数之间的间隔任意大）

证明：当$0<a<=n$时，$n!+a$能被$a$整除。长度为$n-1$的数列$n!+2, n!+3, n!+4, …, n!+n$中，所有的数都是合数。这个结论对所有大于1的整数$n$都成立，而$n$可以取到任意大。

3. 所有大于2的素数都可以唯一地表示成两个平方数之差。

证明：大于2的素数都是奇数。假设这个数是$2n+1$。由于$(n+1)^2=n^2+2n+1$，$(n+1)^2$和$n^2$就是我们要找的两个平方数。下面证明这个方案是唯一的。如果素数$p$能表示成$a^2-b^2$，则$p=a^2-b^2=(a+b)(a-b)$。由于$p$是素数，那么只可能$a+b=p$且$a-b=1$，这给出了$a$和$b$的唯一解。

4. 当n为大于2的整数时，$2^n+1$和$2^n-1$两个数中，如果其中一个数是素数，那么另一个数一定是合数。

证明：$2^n$不能被3整除。如果它被3除余1，那么$2^n-1$就能被3整除；如果被3除余2，那么$2^n+1$就能被3整除。总之，$2^n+1$和$2^n-1$中至少有一个是合数。

5. 如果$p$是素数，$a$是小于$p$的正整数，那么$a^{p-1}\equiv1(mod p)$。（费马小定理）

首先我们证明这样一个结论：如果$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)$

中间略过历史的发展，直接看到我们的MR算法：（文章写的真的很不错，可以好好看看）

这就是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素性测试同样是不确定算法。

如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341550071728320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46856248255981。

通常认为，Miller-Rabin素性测试的正确率可以令人接受，随机选取k个底数进行测试算法的失误率大概为$4^{-k}$。

处理longlong时直接乘可能会溢出，因此采用快速乘。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<iostream>
#include<map>
#include<vector>
#define M 605
#define Max 1000000
#define ll long long
using namespace std;
bool mark[Max+5];
int Mp[Max],sum[10005];
int tot;
int cnt[Max];//<10^6
ll K[M*4];
map<ll,int>S1,S2;//不确定的大质数;确定的大质数 
struct Prime{
    void Init(){
        for(int i=2;i<=Max;i++)if(!mark[i])
            for(int j=i+i;j<=Max;j+=i)
                mark[j]=true;
        for(int i=2;i<=Max;i++)
            if(!mark[i])Mp[++tot]=i;
    }
    void Change(ll x){
        for(int i=1;i<=tot;i++){
            while(x%Mp[i]==0){
                x/=Mp[i];
                cnt[i]++;
            }
            if(x==1)break;
        }
        if(x!=1)S1[x]++;
    }
    ll mul(ll x,ll y,ll mod) {//快速乘
        if(x<y)swap(x,y);
        ll ret=0;
        while(y) {
            if(y&1){
                ret+=x;
                if(ret>=mod)ret-=mod;
            }
            x+=x;
            if(x>=mod)x-=mod;
            y>>=1;
        }
        return ret;
    }
    ll Pow(ll x,ll y,ll mod){
        ll ret=1;
        while(y){
            if(y&1)ret=mul(ret,x,mod);
            x=mul(x,x,mod);
            y>>=1;
        }
        return ret;
    }
    bool check(ll u,ll t,ll now) {
        ll s=rand()%(now-2)+2;
        s=Pow(s,u,now);
        if(s==1||s==now-1)return 1;
        for(ll i=1;i<=t;i++){
            s=mul(s,s,now);
            if(s==now-1&&i<t)return 1;
        }
        return 0;
    }
    bool chk(ll now) {
        if(now<2)return 0;
        if(!(now&1))return now==2;
        ll u=now-1,t=0;
        while(!(u&1))u>>=1,t++;
        for(int c=1; c<=10; c++)
            if(!check(u,t,now))
                return 0;
        return 1;
    }
}Prime;
ll gcd(ll a,ll b){
    if(b==0)return a;
    else return gcd(b,a%b);
}
void WritePow(int x){
    sum[1]=1;
    int ed=1;
    for(int i=1;i<=x;i++){
        for(int i=1;i<=ed;i++)
            sum[i]*=2;
        for(int i=1;i<=ed;i++){
            if(sum[i]>=10){
                sum[i+1]+=sum[i]/10;
                sum[i]%=10;
                if(i==ed)ed++;
            }
        }
    }
    sum[1]-=1;
    int j=1;
    while(sum[j]<0){
        sum[j]+=10;sum[j+1]--;j++;
    }
    if(sum[ed]==0)ed--;
    for(int i=ed;i>=1;i--)
        printf("%d",sum[i]);
    puts("");
}
int main(){
    srand(time(NULL));
    Prime.Init();
    int n,sz=0;
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        ll x;
        cin>>x;
        Prime.Change(x);
    }
    //处理大质数 
    for(map<ll,int>::iterator it1=S1.begin();it1!=S1.end();it1++)
        for(map<ll,int>::iterator it2=S1.begin();it2!=S1.end();it2++){
            if(it1==it2)continue;
            ll a=it1->first,b=it2->first,d=gcd(a,b);
            if(d!=1)K[++sz]=d;
        }
    for(map<ll,int>::iterator it=S1.begin();it!=S1.end();){
        map<ll,int>::iterator it1=it++;
        ll a=it1->first;int b=it1->second;
        bool f=false;
        for(int i=1;i<=sz;i++)
            if(a%K[i]==0){
                if(a/K[i]!=1)S2[a/K[i]]+=b;
                S2[K[i]]+=b;
                f=true;
                break;
            }
        if(f)S1.erase(it1);
    }
    for(map<ll,int>::iterator it=S1.begin();it!=S1.end();){
        map<ll,int>::iterator it1=it++;
        ll a=it1->first,b=floor(sqrt(a));
        int c=it1->second;
        if(b*b==a){
            S1.erase(it1);
            S2[b]+=c*2;
        }
    }
    //统计结果 
    int ans1=0,ans2=0;
    for(int i=1;i<=tot;i++){
        if(cnt[i]>ans1){
            ans1=cnt[i];
            ans2=1;
        }else if(cnt[i]==ans1)ans2++;
    }
    for(map<ll,int>::iterator it=S2.begin();it!=S2.end();it++){
        int b=it->second;
        if(b>ans1){
            ans1=b;
            ans2=1;
        }else if(b==ans1)ans2++;
    }
    for(map<ll,int>::iterator it=S1.begin();it!=S1.end();it++){
        ll a=it->first;
        int b=it->second,f;
        if(Prime.chk(a))f=1;
        else f=2;
        if(b>ans1){
            ans1=b;
            ans2=f;
        }else if(b==ans1)ans2+=f;
    }
    printf("%d\n",ans1);
    WritePow(ans2);
    return 0;
}