[JZOJ5527] Silly

Description

求$\sum\limits_{i=1}^{n}φ^K(\gcd(i,n))*lcm(i,n)\pmod mo$

T组询问
$N\leq 10^{18},K\leq 10^9,T\leq 300,mo\leq 10^9+7$且mo为奇质数

Solution

化式子

$$=\sum\limits_{d|n}φ^K(d)\sum\limits_{d|i}^{n}[\gcd(id,n)=d]{i*n\over d}$$
$$=\sum\limits_{d|n}φ^K(d)*n\sum\limits_{i=1}^{n \over d}[\gcd(i,{n\over d})=1]i$$
我们知道，小于等于n且与n互质的数的和是 $φ(n)*n\over 2$，特别的，当n=1时和为1

$$=n*\sum\limits_{d|n}φ^K(d)*{φ({n\over d})*{n\over d}\over 2}$$

1的情况要特判
n的约数个数大约是n^1/4级别的，这样会被卡

不考虑1的情况，后面的式子是个狄利克雷卷积的形式

根据狄利克雷卷积的性质，两个积性函数卷积起来还是积性函数

那么对n分解质因子的时候，对每个质因子分开来做，最后再乘在一起，再加上1的情况即可

Code

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <cstdlib>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define N 1000005
#define LL long long
using namespace std;
LL mo,n,m,phi[N],pr[N][3],ans,ny,p1[N],d[5]={2,3,5,11,10007};
int l;
LL mul(LL x,LL y,LL mo)
{
    if(x>mo) x-=mo;
    if(y>mo) y-=mo;
    LL s=x*y-(LL)((long double)x*y/mo+1e-8)*mo;
    return (s+mo)%mo;
}
LL ksd(LL x,LL y,LL mo)
{
    if(y==0||x==0) return 0;
    if(y==1) return x;
    LL s=ksd(x,y/2,mo);
    return(y&1)?((s+s)%mo+x)%mo:(s+s)%mo;
}
LL ksm(LL k,LL n)
{
    LL s=1;
    k=k%mo;
    for(;n;k=k*k%mo,n>>=1) if(n&1) s=s*k%mo;
    return s;
}
void dfs(int k,LL sv)
{
    if(k>l)
    {
        ans=(ans+sv)%mo;
        return;
    }
    LL v1=1,s1=0;
    fo(i,0,pr[k][1])
    {       
        s1=(s1+ksm(v1-v1/pr[k][0],m)%mo*((pr[k][2]/v1-pr[k][2]/v1/pr[k][0])%mo)%mo*((pr[k][2]/v1)%mo)%mo)%mo;
        v1=v1*pr[k][0];
    }
    dfs(k+1,sv*s1%mo);
}
LL rd(LL k)
{
    return (rand()*rand()%k+rand())%k;
}
LL km(LL k,LL n,LL mo)
{
    LL s=1;
    for(;n;k=mul(k,k,mo),n>>=1) if(n&1) s=mul(s,k,mo);
    return s;
}
bool pd1(LL k)
{
    fo(i,0,4)
    {
        if(k==d[i]) return 1;
        if(k%d[i]==0) return 0;
        LL t,n=k-1;
        while((n&1)==0) n>>=1;
        t=km(d[i],n,k);
        while(n<k-1&&t!=1&&t!=k-1) t=mul(t,t,k),n<<=1;
        if(!((n&1)||t==k-1)) return 0;
    }   
    return 1;
}
LL gcd(LL x,LL y)
{
    if(y==0) return x;
    else return gcd(y,x%y);
}
LL get(LL k)
{
    LL c=rd(1e9)+1,r,x,y,i=1,n=2;
    x=y=rd(k-1)+1;
    while(1)
    {
        i++;
        x=(mul(x,x,k)+c)%k;
        if(x==y) return k;
        r=gcd(abs(x-y),k);
        if(1<r&&r<k) return r;
        if(i==n) y=x,n<<=1;
    }
}
void fd(LL k)
{
    if(k==1) return;
    if(pd1(k))
    {
        p1[++p1[0]]=k;
        return;
    }
    LL p=k;
    while(p==k&&p!=1) 
        p=get(k);
    fd(k/p),fd(p);
}
int main()
{
    int t;
    srand(1657479);
    cin>>t; 
    while(t--)
    {
        cin>>n>>m>>mo;
        ny=ksm(2,mo-2)%mo;
        p1[0]=0;
        fd(n);
        sort(p1+1,p1+p1[0]+1);
        l=ans=0;
        pr[++l][0]=p1[1],pr[l][1]=1,pr[l][2]=p1[1];
        LL phi=p1[1]-1;
        fo(i,2,p1[0])
        {
            if(p1[i]!=p1[i-1])
            {
                pr[++l][0]=pr[l][2]=p1[i];
                pr[l][1]=1;
                phi*=(LL)p1[i]-1;
            }
            else pr[l][1]++,pr[l][2]*=p1[i],phi*=(LL)p1[i];
        }
        dfs(1,1);
        printf("%lld\n",(ans+ksm(phi,m))*ny%mo*(LL)(n%mo)%mo);
    }
}