分解大质因数

参考博客
Description of the topic
In FZU ACM team, BroterJ and Silchen are good friends, and they often play some interesting games.
One day they play a game about GCD and LCM. firstly BrotherJ writes an integer A and Silchen writes an integer B on the paper. Then BrotherJ gives Silchen an integer X. Silchen will win if he can find two integers Y1 and Y2 that satisfy the following conditions:
• GCD(X, Y1) = A
• LCM(X, Y2) = B
• Fuction GCD(X, Y ) means greatest common divisor between X and Y .
• Fuction LCM(X, Y ) means lowest common multiple between X and Y .
BrotherJ loves Silchen so much that he wants Silchen to win the game. Now he wants to calculate how many number of X he can give to Silchen.

Input
Input is given from Standard Input in the following format:
A B
Constraints
1 ≤ A, B ≤ 10e18
Both A and B are integers.

Output
Print one integer denotes the number of X.

Sample input
3 12

Sample output
3
算法分析
假设c*A = X，d*x = B，我们的目标是求有多少个不同的x，换句话说就是求有多少个c或者d。c*d = B/A,c就是B/A的正约数。所以这道题就是求B/A的正约数个数。求正约数个数就是分解质因数。B/A的范围非常大。假设某个质数p>1e6，那么我可以断定B/A最多只包含两个大于1e6的质数。因为p*p*p>1e18。因此我们使用线筛求解1~1e7之间的质数，将B/A之间1e7以内的质因数分解求得结果ans。然后使用Miller-rabin算法判断剩余的数是一个质数，还是两个相同的质数，或者是两个不相同的质数。第一种情况答案就是2*ans，第二种情况答案是3*ans，第三种情况答案是4*ans。
代码中Miller-rabin的实现会误判。

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long LL;

LL a,b,kk;
LL prime[6] = {2, 3, 5, 233, 331};
LL gprime[10000005];
LL sum[100];
bool vis[10000005];
int cnt;
void primejudge()
{
    for(int i=2;i<10000005;i++)
    {
        if(!vis[i]) gprime[cnt++]=i;
        for(int k=0;k<cnt && i*gprime[k]<10000005;k++)
        {
            vis[i*gprime[k]]=true;
            if(i%gprime[k]==0) break;
        }
    }
}
LL qmul(LL x, LL y, LL mod)
{
    LL ret = 0;
    while(y) {
        if(y & 1)
            ret = (ret + x) % mod;
        x = x * 2 % mod;
        y >>= 1;
    }
    return ret;
}
LL qpow(LL a, LL n, LL mod)
{
    LL ret = 1;
    while(n)
    {
        if(n & 1)
            ret = qmul(ret, a, mod);
        a = qmul(a, a, mod);
        n >>= 1;
    }
    return ret;
}
bool Miller_Rabin(LL p)
{
    if(p < 2) return 0;
    if(p != 2 && p % 2 == 0) return 0;
    LL s = p - 1;
    while(! (s & 1)) s >>= 1;
    for(int i = 0; i < 5; ++i)
    {
        if(p == prime[i]) return 1;
        LL t = s, m = qpow(prime[i], s, p);
        while(t != p - 1 && m != 1 && m != p - 1)
        {
            m = qmul(m, m, p);
            t <<= 1;
        }
        if(m != p - 1 && !(t & 1)) return 0;
    }
    return 1;
}

LL Pollards_rho(LL x)
{
    int tot=1;
    LL ans=1;
    for(int i=0;gprime[i]*gprime[i]<=x && i<cnt ;i++)
    {
        if(x%gprime[i]==0)
        {
            while(x%gprime[i]==0)
            {
                x/=gprime[i];
                sum[tot]++;
            }
            tot++;
        }
        if(x==1) break;
    }
    for(int i=1;i<tot;i++) ans*=(sum[i]+1);
    if(x>1)
    {
        if(!Miller_Rabin(x))
        {
            LL tmp=sqrt(x);
            if(tmp*tmp==x) ans*=3;
            else ans*=4;
        }
        else ans*=2;
    }
    return ans;
}
int main()
{
    primejudge();
    scanf("%lld%lld",&a,&b);
    if (b%a)
    {
        printf("0\n");
        return 0;
    }
    kk=b/a;
    cout<<Pollards_rho(kk) << endl;
    return 0;
}