JZOJ 5527 Silly

Silly

Description

Data Constraint

$N<=10^{18}$,$K<=10^9$,$Mo<=10^{18}$且$Mo≠2$且$Mo$为质数

Solution

原式等于

$$Answer=\sum_{d|n}\varphi(d)^K*\{\sum_{i=1}^{n\over d}i*n*[(i,n)=1]\}$$

$$={\varphi(n)^K \over 2}+n*\sum_{d|n}\varphi(d)^K*{{{n\over d}*\varphi({n\over d})}\over 2}$$

$$={\varphi(n)^K \over 2}+{n\over 2}*\sum_{d|n}\varphi(d)^K*{{{n\over d}*\varphi({n\over d})}}$$

设$f(n)=ID(n)*\varphi(n)$，由于$ID$和$\varphi$都是积性函数，所以$f$也是积性函数。
显然有$g(n)=\varphi(n)^K$也为积性函数。
设$h(n)=\sum_{d|n}\varphi(d)^K*{{{n\over d}*\varphi({n\over d})}}$，则有
$h=f*g$，因为$f$和$g$都是积性函数，所以$h$为积性函数。

设$n=\Pi p_i^{a_i}$，则$h(n)=\Pi h(p_i^{a_i})$
我们只需分解质因数计算每一个质数的贡献即可。
分解质因数要用$Pollar\ rho$。

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#pragma GCC optimize(3)

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)
#define random(x) rand()%(x)

using namespace std;
typedef unsigned long long ll;
const ll N=666,K=1e6,sj=5e5;

struct note{
    ll ba,index,pro; 
}re[N];

ll x[N],oo,n,k,mod,ans,s[N];
int T,ss[sj],from[K],num;
bool bz[K];

inline ll times(ll x,ll y,ll mo)
{
    x=(x>mo)?x%mo:x; y=(y>mo)?y%mo:y; 
    ll ltime=(ll)((long double)x*y/mo+1e-8)*mo;
    return (x*y-ltime+mo)%mo;
}

inline ll ksm(ll o,ll t,ll mo)
{
    ll y=1; o=(o>=mo)?(o%mo):o; t=(t>=mo-1)?t%(mo-1):t;
    for(;t;t>>=1,o=times(o,o,mo))
    if(t&1)y=times(y,o,mo);
    return y;
}

inline bool Miller_Rabin(ll o)
{
    if(!(o&1))return false;
    int y; ll u=o-1;
    for(y=0;!(u&1);y++)u>>=1;
    fo(i,1,3){
        x[0]=random((o-1))+1;
        x[0]=ksm(x[0],u,o);
        fo(l,1,y){
            x[l]=times(x[l-1],x[l-1],o);
            if(x[l]==1&&x[l-1]!=1&&x[l-1]!=(o-1))return false;
        }
        if(x[y]!=1)return false;
    } 
    return true;
}

inline ll ab(ll o)
{return o>=0?o:(-o);}

ll gcd(ll a,ll b)
{return (!b)?a:gcd(b,a%b);}

void find(ll o)
{
    if(o<=sj){
        for(;o>1;o/=from[o])s[++oo]=from[o];
        return;
    }
    if(Miller_Rabin(o)){
        s[++oo]=o; return;
    }
    ll z=random((o-1))+1,p=0;
    ll x1=random(o),x2=x1;
    for(int i=2,k=2;true;++i){
        x1=(times(x1,x1,o)+z)%o;
        p=gcd(ab(x2-x1),o);
        if(p>1&&p<o)break;
        if(x1==x2)x1=x2=random(o),z=random(o-1)+1;
        if(i==k)k<<=1,x2=x1;
    }
    find(o/p); find(p);
}

int main()
{
    freopen("silly.in","r",stdin);
    freopen("silly.out","w",stdout);
    srand(1684168);
    cin>>T; oo=0;
    fo(i,2,sj){
        if(!bz[i])ss[++oo]=i,from[i]=i;
        fo(j,1,oo){
            if(i*ss[j]>sj)break;
            from[i*ss[j]]=ss[j];
            bz[i*ss[j]]=true;
            if(i%ss[j]==0)break;
        }
    }
    fo(tt,1,T){
        scanf("%lld%lld%lld",&n,&k,&mod);
        oo=num=0; ll kk=n; 
        for(;!(kk&1);kk>>=1)s[++oo]=2; 
        find(kk);   sort(s+1,s+oo+1);
        ss[0]=0;
        fo(i,1,oo)if(s[i]==s[i-1])++re[num].index;
        else{
            re[++num].ba=s[i];
            re[num].index=1;
        }
        ll ans=1,o1=1;
        fo(i,1,num){
            ll op=1,da=0,o3=1;
            fo(l,1,re[i].index-1)op*=re[i].ba;
            o3=op*re[i].ba;
            op=op*(re[i].ba-1);
            o1=o1*op;
            da=times(o3,op,mod);
            ll o2=1,o4=o3,o5;
            fo(j,1,re[i].index){
                o2=j==1?re[i].ba-1:o2*re[i].ba;
                if(op==re[i].ba-1)op=1;else op=op/re[i].ba;
                o3/=re[i].ba;
                o5=times(op,o3,mod);
                da=(da+times(ksm(o2,k,mod),o5,mod))%mod;
            }
            ans=times(ans,da,mod);
        }
        ans=(ans+ksm(o1,k,mod))%mod;
        ll nm=ksm(2,mod-2,mod);
        ans=times(ans,nm,mod);
        ans=times(ans,n%mod,mod);
        printf("%lld\n",ans);
    }
}