[JZOJ5527]. 【清华冬令营2018模拟】 Silly

分析

一到大整数运算就贼恶心。
注意pollard_rho，死循环，一般是miller_rabin打错。
因为小质数十分麻烦，先对1e5以内的质数模一遍。
miller_rabin的时候搞8次就够了，随机的数不能是0。
大整数乘法取模最好用log的正确方法。类似快速幂

LL mul(LL x,LL y,LL mo)
{
    if (!y) return 0;
    LL s=mul(x,y>>1,mo); s=(s+s)%mo;
    if (y&1) s=(s+x)%mo;
    return s;
}

代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
using namespace std;
typedef long long ll;
typedef double db;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
#define cmax(a,b) (a=(a>b)?a:b)
#define cmin(a,b) (a=(a<b)?a:b)
#define ls tr[x][0]
#define rs tr[x][1]
const int N=1e5+5,M=2e5+5;
ll n,K,mo,ans,tmp,d[N],c[N],pri[N],p1,p2,Phi,Pri[N],Pd[N],sn;
int i,j,T,TT;
ll ran(ll x)
{
    return rand()%x;
}
ll ch(ll x,ll y,ll mo)
{
    return (__int128)x*y%mo;
}
ll gcd(ll x,ll y)
{
    if (!y) return x;
    return gcd(y,x%y);
}
ll ksm(ll x,ll y,ll mo)
{
    ll ret=1;
    while (y)
    {
        if (y&1) ret=ch(ret,x,mo);
        y>>=1;
        x=ch(x,x,mo);
    }
    return ret;
}
int test(ll x,ll a)
{
    ll ci=0,i,y,z,t;
    y=x-1;
    while (y%2==0)
        y/=2,ci++;
    z=ksm(a,y,x);
    while (ci--)
    {
        t=ch(z,z,x);
        if (t==1&&z!=1&&z!=x-1) return 0;
        z=t;
    }
    return z==1;
}
int miller_rabin(ll x)
{
    if (x==2) return 1;
    if (x%2==0) return 0;
    int i;
    fo(i,1,10) 
        if (!test(x,ran(x-1)+1))
            return 0;
    return 1;
}
ll rho(ll x)
{
    ll c=ran(x-1)+1,a=ran(x-1)+1,i=1,k=2;
    ll b=a,d;
    while (1)
    {
        i++;
        a=(ch(a,a,x)+c)%x;
        if (a==b) return 0;
        d=gcd(x,abs(a-b));
        if (d!=1&&d!=x) return d;
        if (i==k)
        {
            k<<=1;
            b=a;
        }
    }
}
void get(ll x)
{
    if (miller_rabin(x))
    {
        pri[++pri[0]]=x;
        return ;
    }
    int y=0;
    while (!y)
        y=rho(x);
    get(y);
    get(x/y);
}
ll phi(ll x)
{
    return x-x/d[i];
}
ll f(ll p,ll k)
{
    p1=1;
    p2=ksm(p,k,1e18);
    Phi=ch(Phi,phi(p2),1e18);
    ll ret=0,tp=0;
    int i;
    fo(i,0,k)
    {
        tp=ch(phi(p2),p2,mo);
        ret=ret+ch(ksm(phi(p1),K,mo),tp,mo);
        p1*=p;
        p2/=p;
    }
    return ret;
}
void clear()
{
    fd(i,d[0],0) d[i]=c[i]=0;
    fd(i,pri[0],0) pri[i]=0;
}
void work()
{
    clear();
    if (TT==188)
    {
        i=0;
    }
    scanf("%lld %lld %lld",&n,&K,&mo);
    sn=n;
    fo(i,1,Pri[0])
    {
        while (sn%Pri[i]==0)
        {
            sn/=Pri[i];
            pri[++pri[0]]=Pri[i];
        }
    }
    if (sn!=1)
        get(sn);
    sort(pri+1,pri+1+pri[0]);
    i=1;
    while (i<=pri[0])
    {
        j=i;
        while (pri[j+1]==pri[i]) j++;
        d[++d[0]]=pri[i];
        c[d[0]]=j-i+1;
        i=j+1;
    }
    ans=1;
    Phi=1;
    fo(i,1,d[0])
    {
        tmp=f(d[i],c[i]);
        ans=ch(ans,tmp,mo);
    }
    Phi=ksm(Phi,K,mo);
    ans=(ans+Phi)%mo;
    ans=ch(ans,n,mo);
    ans=ch(ans,ksm(2,mo-2,mo),mo);
    printf("%lld\n",ans);
}
int main()
{
    freopen("silly.in","r",stdin);
    freopen("silly.out","w",stdout);
    fo(i,2,100000)
    {
        if (!Pd[i]) Pri[++Pri[0]]=i;
        fo(j,1,Pri[0])
        {
            if (1ll*i*Pri[j]>100000) break;
            Pd[i*Pri[j]]=1;
            if (i%Pri[j]==0) break;
        }
    }
    scanf("%d",&T);
    srand(T);
    fo(TT,1,T)
        work();
}