参考ARC064F
令$h(n)=egin{cases}n(n为奇数)\frac{n}{2}(n为偶数)end{cases}$，$f(n)$定义与ARC064F相同，答案即$sum_{d|n}h(d)f(d)$
考虑$f(n)$的转移，即$sum_{d|n}f(d)=k^{lceilfrac{n}{2} ceil}$，将后者记作$g(n)$，即$(f*1)(x)=g(x)$，那么$(g*mu)(x)=f(x)$，代入即得到$f(n)=sum_{d|n}mu(frac{n}{d})g(d)$
将之代入，原式即$sum_{d|n}h(d)sum_{p|d}mu(frac{d}{p})g(p)$
交换枚举顺序，即$sum_{p|n}g(p)sum_{d'|frac{n}{p}}mu(d')h(d'p)$
注意到当$h(d'p) e d'h(p)$当且仅当$d$为偶数且$p$为奇数，先来对$n$分类讨论：
1.若$n$为奇数，则$d$必然为奇数，即恒有$h(d'p)=d'h(p)$
2.若$n$为偶数，则对于奇因子$p$，对$d'$中2的幂次分类讨论——
(1)$4|d'$，根据莫比乌斯函数的定义有$mu(d')=0$
(2)$4
otmid d'$，这些因子一定可以被表示为$frac{n}{p}$的奇因子以及它们的2倍，那么考虑枚举这个奇因子，即
$$
sum_{d'|frac{n}{p},d'为奇数}mu(d')h(d'p)+mu(2d')h(2d'p)=sum_{d'|frac{n}{p},d'为奇数}(mu(d')+mu(2d'))d'p=0
$$
由此，我们知道了当$n$为偶数且$p$为奇数时后面的式子一定为0，因此不妨在$n$为偶数时强制$p$为偶数
综上，原式即$sum_{p|n,n为奇数或p为偶数}h(p)g(p)sum_{d'|frac{n}{p}}d'mu(d')$
对于后者，考虑对$frac{n}{p}$质因数分解，即$prod_{i=1}^{k}p_{i}^{a_{i}}$，那么$d'$必然是从中选择若干个素数，假设为集合$S$，注意到此时的贡献即$(-1)^{|S|}prod_{xin S}p_{x}=prod_{xin S}(-p_{x})$
考虑$prod_{i=1}^{k}(1-p_{i})$，对于每一项选$p_{i}$即加入$S$中，不选即可看作1，也即答案
再假设$n=prod_{i=1}^{k}p_{i}^{a_{i}}$，注意到其因子个数并不不多，暴力递归$p$每一个素数的幂次（特判$2$这个素因子），并在递归时维护后式即可
使用Pollard-Rho算法进行$o(n^{frac{1}{4}})$来质因数分解，再枚举因子以及$g(p)$的快速幂，复杂度为$o(sigma(n)log n)$
（关于Rollard-Rho算法，可以参考洛谷4718）
不难证明$sigma(n)le 2^{16}$（即$n<2 imes 3 imes... imes 53$），因此复杂度可过
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define ll long long 4 vector<ll>v; 5 int t,m,mod,ans,P[5]={3,5,7,11,13},tot[105]; 6 ll n,k,p[105]; 7 ll mul(ll x,ll y,ll n){ 8 return (__int128)x*y%n; 9 } 10 ll gcd(ll x,ll y){ 11 if (!y)return x; 12 return gcd(y,x%y); 13 } 14 ll pow(ll n,ll m,ll mod){ 15 ll s=n,ans=1; 16 while (m){ 17 if (m&1)ans=mul(ans,s,mod); 18 s=mul(s,s,mod); 19 m>>=1; 20 } 21 return ans; 22 } 23 bool check(ll n){ 24 if ((n==1)||(n==2))return 1; 25 if (n%2==0)return 0; 26 ll m=n-1,k=0; 27 while (m%2==0){ 28 k++; 29 m/=2; 30 } 31 for(int i=0;i<5;i++){ 32 if (n%P[i]==0)return n==P[i]; 33 ll s=pow(P[i],m,n); 34 for(int j=0;j<k;j++){ 35 if ((mul(s,s,n)==1)&&(s!=1)&&(s!=n-1))return 0; 36 s=mul(s,s,n); 37 } 38 if (s>1)return 0; 39 } 40 return 1; 41 } 42 ll get_fac(ll n){ 43 ll x=rand()%n,c=rand()%n,y=(mul(x,x,n)+c)%n,z=1,tot=0; 44 while (1){ 45 if ((x==y)||(tot%100==0)){ 46 if (gcd(z,n)>1)return gcd(z,n); 47 if (x==y)return 1; 48 z=1; 49 } 50 tot++; 51 z=mul(z,(x+n-y)%n,n); 52 if (!z)return gcd((x+n-y)%n,n); 53 x=(mul(x,x,n)+c)%n; 54 y=(mul(y,y,n)+c)%n; 55 y=(mul(y,y,n)+c)%n; 56 } 57 } 58 void calc(ll n){ 59 if (check(n)){ 60 v.push_back(n); 61 return; 62 } 63 while (1){ 64 ll p=get_fac(n); 65 if (p>1){ 66 calc(p); 67 calc(n/p); 68 return; 69 } 70 } 71 } 72 void dfs(int k,ll d,int s){ 73 if (k>p[0]){ 74 if ((n%2==0)&&(d&1))return; 75 if (d&1)ans=(ans+1LL*s*pow(m,(d+1)/2,mod)%mod*(d%mod))%mod; 76 else ans=(ans+1LL*s*pow(m,(d+1)/2,mod)%mod*(d/2%mod))%mod; 77 return; 78 } 79 for(int i=1;i<=tot[k];i++)d*=p[k]; 80 dfs(k+1,d,s); 81 s=1LL*s*(mod+1-p[k]%mod)%mod; 82 for(int i=1;i<=tot[k];i++){ 83 d/=p[k]; 84 dfs(k+1,d,s); 85 } 86 } 87 int main(){ 88 srand(time(0)); 89 scanf("%d",&t); 90 while (t--){ 91 scanf("%lld%lld%d",&n,&k,&mod); 92 m=k%mod; 93 v.clear(); 94 calc(n); 95 sort(v.begin(),v.end()); 96 ans=p[0]=0; 97 for(int i=0;i<v.size();i++){ 98 if (v[i]!=p[p[0]]){ 99 p[++p[0]]=v[i]; 100 tot[p[0]]=0; 101 } 102 tot[p[0]]++; 103 } 104 dfs(1,1,1); 105 printf("%d ",ans); 106 } 107 }