BZOJ2082 POI2010 Divine divisor

最新推荐文章于 2018-07-13 13:02:00 发布

Kanosword

最新推荐文章于 2018-07-13 13:02:00 发布

阅读量920

点赞数

分类专栏： POI 2010 文章标签：数论 POI

本文链接：https://blog.csdn.net/kanosword/article/details/52717273

版权

POI 2010 专栏收录该内容

11 篇文章 2 订阅

订阅专栏

POI2010 题解整理

题目大意

给定一个n，找到n中出现次数最多的那些因子，并输出个数。（即要求找出满足 $n\%a^x=0$ ， $x$ 最大的a的个数）。

Input

第一行给出一个整数 $m(1\le m\le600)$

第二行是由空格分开的 $m$ 个数 $ai(2\le a_i\le 10^{18})$ ， $n$ 就是 $a_i$ 的乘积 $（n=a_1\times a_2\times a_3\times\cdots\times a_m）$ 。

Sample Input

3
4 3 4

1
6

Output

第一行为最大的 $k$ 。

第二行为非凡因子 $d$ 的个数。

Sample Output

4
1

1
3

感谢Matrix67大犇对Miller_Rabin的讲解，非常详细，一看就懂。

$Solution$ ：

总的来说这题的思路其实非常简单，而代码只是比较复杂而已。如果每个 $a_i$ 的数据范围小一点（至少 $\sqrt{a_i}$ 可以被筛表覆盖），我们显然是考虑将所有的数全部进行分解，将所有因子统计归类，然后全部扫一遍就可以做了。此时有一个最基本不过的结论：

如果 $\sqrt{a_i}$ 内的素数都无法整除 $a_i$ ，那么这个数一定是素数。

当范围提升到 $10^{18}$ 后，筛出 $10^6$ 以上的素数便显得有些吃力无用了。根据上面那个结论，我们同样有以下结论：

用 $10^6$ 内的素数去试除 $a_i$ ，则最后剩下的值一定只包含 $10^6$ 以上数量级的素数，并且最多只会剩下两个素数。

那么我们有以下这个算法：

将所有 $a_i$ 在 $10^6$ 数量级以下的质因子全部筛出，最后剩下的质因子 $p,q$ 只会符合三种情况： $a i = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ p * q p 2 P p, q \neq 1, p \neq q, g c d (p, q) = 1 p = q p \neq 1, q = 1$ $a_i= \begin{cases} p*q & p,q\neq 1,p\neq q,gcd(p,q)=1 \\[2ex] p^2 & p=q \\[2ex] P & p\neq 1,q=1 \end{cases}$

接下来就是考虑如何完成这个算法了。

首先就是如何判断这个数是一个大素数，还是两个素数相乘。这里就需要采用经典的Miller_Rabin算法进行判断。由于这个范围已经达到 $10^{18}$ ，所以还要考虑用快速乘来避免long long溢出。

接下来重点处理p*q的情况，这个可以进行两两gcd拆分，此时还会出现无论如何也无法将这两个拆开的情况，这个时候就把这个整体当素数维护，统计个数时算两份就可以了。

最后的计算结果是：

C 1 n + C 2 n + C 3 n + \dots + C n n = 2 n - 1

$C^{1}_{n}+C^{2}_{n}+C^{3}_{n}+\cdots+C^{n}_{n}=2^n-1$ 由于n可能比较大，所以还需要高精度计算。

#include <bits/stdc++.h>
#define M 1000005
#define N 605
using namespace std;
typedef long long ll;
//预处理与调试部分
inline void Rd(ll &ans){
    ans=0;char c;
    while(c=getchar(),c<48);
    do ans=(ans<<3)+(ans<<1)+(c^48);
    while(c=getchar(),c>47);
}
long long res[N];
bool vis[M];
map<ll,int>Mp;
int Map[M],mtop=0,ex=0;
void init(){
    for(int i=2;i<M;i++)if(!vis[i]){
        Map[++mtop]=i;
        for(int j=i+i;j<M;j+=i)vis[j]=1;
    }
}
void debug(){
    map<ll,int>::iterator it=Mp.begin();
    for(;it!=Mp.end();++it)//if(it->second)
        cout<<"first="<<it->first<<','<<"second="<<it->second<<endl;
}
//素性测试：Miller_Rabin
const int prime[]={2,3,7,61,21251};
ll mod(ll m,ll n,ll P){
    ll ans=m+n;
    if(ans>=P)ans-=P;
    return ans;
}//加法，带取模
ll mul(ll m,ll n,ll P){
    if(m<n)swap(m,n);
    ll ans=0;
    while(n){
        if(n&1)ans=mod(ans,m,P);
        n>>=1;m=mod(m,m,P);
    }
    return ans;
}//快速乘 
ll pow(ll m,ll n,ll P){
    ll ans=1;
    while(n){
        if(n&1)ans=mul(ans,m,P);
        n>>=1;m=mul(m,m,P);
    }return ans;
}//快速幂
bool Miller_Rabin(ll pri,ll n){//只有a^d%n==1 或者 a^{d*2^i}%n==n-1
    ll d=n-1;
    while(d%2==0)d>>=1;
    ll t=pow(pri,d,n);
    while(d!=n-1&&t!=1&&t!=n-1)
        t=mul(t,t,n),d<<=1;//当t==1时，上面的所有结果必然都是t==1
    return (t==n-1)||(d&1);
}
bool judge(ll n){
    for(int i=0;i<5;i++)
        if(!Miller_Rabin(1LL*prime[i],n))return false;
    return true;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}

set<ll>Twoprime;
//高精度部分
struct BigInt{
    static const int P=10000;
    int len,num[1005];
    BigInt(){
        len=1;
        memset(num,0,sizeof(num));
    }
    BigInt operator * (const BigInt &A)const{
        BigInt B;
        B.len=len+A.len-1;
        for(int i=0;i<len;i++)
            for(int j=0;j<A.len;j++){
                B.num[i+j]+=num[i]*A.num[j];
                if(B.num[i+j]>=P){
                    B.num[i+j+1]+=B.num[i+j]/P;
                    B.num[i+j]%=P;
                }
            }
        if(B.num[B.len])B.len++;
        return B;
    }
    void Print(){
        printf("%d",num[len-1]);
        for(int i=len-2;i>=0;i--)printf("%04d",num[i]);
        puts("");
    }
};
BigInt Pow(BigInt n,ll m){
    BigInt Ans;
    Ans.num[0]=1;
    while(m){
        if(m&1)Ans=Ans*n;
        m>>=1;n=n*n;
    }
    return Ans;
}

int main(){
    init();int n;
    scanf("%d",&n);
    for(int i=1;i<=n;i++)Rd(res[i]);
    for(int i=1;i<=n;i++){
        for(int j=1;Map[j]<=res[i]&&j<=mtop;j++)
            while(res[i]%Map[j]==0){
                res[i]/=Map[j];
                Mp[Map[j]]++;
            }
        if(res[i]==1)continue;
        ll p=sqrt(res[i])+0.01;
        if(p*p==res[i]){res[i]=1,Mp[p]+=2;continue;}
        if(judge(res[i]))Mp[res[i]]++,res[i]=1;
    }
    for(int i=1;i<=n;i++)
        for(int j=1;res[i]!=1&&j<=n;j++){
            if(i==j||res[j]==1)continue;
            ll c=gcd(res[i],res[j]);
            if(c!=1&&c!=res[i])Mp[c]+=0,Mp[res[i]/c]+=0;
        }
    for(int i=1;i<=n;i++)
        if(res[i]!=1){
            map<ll,int>::iterator it=Mp.begin();
            for(;it->first<=res[i]&&it!=Mp.end();++it)
                while(res[i]%it->first==0){
                    res[i]/=it->first;
                    it->second++;
                }
            if(res[i]!=1){
                Twoprime.insert(res[i]);
                Mp[res[i]]++;
                res[i]=1;
            }
        }
    ll ans=0,cnt=0;
    map<ll,int>::iterator it=Mp.begin();
    while(it!=Mp.end()){
        if(ans<it->second)ans=it->second,cnt=0;
        if(ans==it->second){
            if(Twoprime.find(it->first)!=Twoprime.end())cnt+=2;
            else cnt++;
        }++it;
    }
    BigInt Ans,Base;Base.num[0]=2;
    Ans=Pow(Base,cnt);Ans.num[0]--;
    cout<<ans<<endl;
    Ans.Print();
}