参考ARC064F

令$h(n)=\begin{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^{\lceil\frac{n}{2}\rceil}$，将后者记作$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)\ne 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\not\mid 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_{x\in S}p_{x}=\prod_{x\in 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\times 3\times…\times 53$），因此复杂度可过

1 #include
2 using namespace std;
3 #define ll long long
4 vectorv;
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;j1)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\n",ans);
106 }
107 }