1. 程式人生 > 實用技巧 >BZOJ-4802 尤拉函式(Pollard-Rho演算法模板)

BZOJ-4802 尤拉函式(Pollard-Rho演算法模板)

題目描述

  給出 \(n(1\leq n\leq 10^{18})\),求 \(\varphi(n)\)

分析

  根據 \(\varphi(n)=n\times \displaystyle\prod\limits_{i=1}^{k}(1-\frac{1}{p_i})\)\(\text{Pollard-Rho}\) 分解即可。

程式碼

#include<bits/stdc++.h>
using namespace std;
#define int128 __int128
long long quick_pow(long long a,long long b,long long mod)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
            ans=(int128)ans*a%mod;
        a=(int128)a*a%mod;
        b>>=1;
    }
    return ans;
}
long long max_factor,n;
bool Miller_Rabin(long long p)
{
    if(p<2)
        return 0;
    if(p==2)
        return 1;
    if(p==3)
        return 1;
    long long d=p-1,r=0;
    while(d%2==0)
    {
        r++;
        d>>=1;
    }
    for(long long k=0;k<10;k++)
    {
        long long a=rand()%(p-2)+2;
        long long x=quick_pow(a,d,p);
        if(x==1||x==p-1)
            continue;
        for(int i=0;i<r-1;i++)
        {
            x=(int128)x*x%p;
            if(x==p-1)
                break;
        }
        if(x!=p-1)
            return 0;
    }
    return 1;
}
long long f(long long x,long long c,long long n)
{
    return ((int128)x*x+c)%n;
}
long long Pollard_Rho(long long n)
{
    long long s=0,t=0;
    long long c=rand()%(n-1)+1;
    int step=0,goal=1;
    long long val=1;
    for(goal=1; ;goal<<=1,s=t,val=1)
    {
        for(step=1;step<=goal;step++)
        {
            t=f(t,c,n);
            val=(int128)val*abs(t-s)%n;
            if(step%127==0)
            {
                long long d=__gcd(val,n);
                if(d>1)
                    return d;
            }
        }
        long long d=__gcd(val,n);
        if(d>1)
            return d;
    }
}
long long factor[300010];
int cnt;
void fac(long long x)
{
    if(x<=max_factor||x<2)
        return ;
    if(Miller_Rabin(x))
    {
        factor[++cnt]=x;
        //max_factor=max(max_factor,x);
        return ;
    }
    long long p=x;
    while(p>=x)
        p=Pollard_Rho(x);
    while(x%p==0)
        x=x/p;
    fac(x);fac(p);
}
int main()
{
    int T=1;
    //cin>>T;
    while(T--)
    {
        srand((unsigned)time(NULL));
        max_factor=0;
        scanf("%lld",&n);
        fac(n);
        sort(factor+1,factor+1+cnt);
        cnt=unique(factor+1,factor+1+cnt)-factor-1;
        long long ans=n;
        for(int i=1;i<=cnt;i++)
            ans=ans/factor[i]*(factor[i]-1);
        cout<<ans<<endl;
        /*if(max_factor==n)
            puts("Prime");
        else
            printf("%lld\n",max_factor);
         */
    }
    return 0;
}