​​传送门​​

• 以下整点皆在第一象限讨论
  令 ，要求的就是
  考虑把整点 用复平面的高斯整数 表示，而 ，于是求 变成了求有多少个高斯整数使得它与它的共轭复数乘积为
• 定理：高斯整数可以写成一系列不能继续分解的高斯整数的乘积，并且为唯一表示（ 论文8.2）

费马平方和定理： 奇质数 可以表示为两个正整数 的平方和 ，当且
仅当 ，并且这样的表示若存在，在不计较排列顺序的情况下，是唯一的
证明：
首先 模 只可能为 故 不可能为 的质数
引理1：设 是质数，那么 ， 时成立，否则考虑乘法逆元不为本身的数
引理2：设质数 满足 ，那么存在 满足 ，注意到 ，因为 为 4 的倍数，所以我们可以将 进行配对，可以得到
由引理 2，对于 形的质数 ，我们知道 ，设 为高斯质数，有 或 ，而 并不是高斯整数，所以不成立，故我们可以令 ，有 ，由于 是质数，所以 ，又 是一个整数，所以 共轭，设 ，
下面来证明是唯一表示，设 ，若 为高斯质数，，或 ，又 ，故 是同种表示，否则令 ，有

• 记为的分解方式，容易发现是个积性函数，我们要做的就是求这个积性函数的前缀和，考虑在质数及质数次幂处的取值
• 对于的质数，的值是，容易发现分解方式为种
  对于的质数，存在唯一的分解方式
  容易发现的至少与的贡献均为 1
  那么我们可以知道（令）
  
  可以用简单计算

#include<bits/stdc++.h>
#define cs const
#define pb push_back
using namespace std;
typedef long long ll;
int Mod;
int add(int a, int b){ return a + b >= Mod ? a + b - Mod : a + b; }
int dec(int a, int b){ return a - b < 0 ? a - b + Mod : a - b; }
int mul(int a, int b){ return 1ll * a * b % Mod; }
void Add(int &a, int b){ a = add(a,b); }
void Dec(int &a, int b){ a = dec(a,b); }
void Mul(int &a, int b){ a = mul(a,b); }
int ksm(int a, ll b){ int as=1; for(;b;b>>=1,Mul(a,a)) if(b&1) Mul(as,a); return as; }
int fx(ll a){ return a >= Mod ? a % Mod : a; }
cs int N = 1e6 + 50;
ll n, k; int S, pw[N], prm[N], pc, Sm[2][N]; 
bool isp[N]; int f[2][N];
int A[N], B[N], nd; ll v[N];
void linear_sieve(int n){
  for(int i=2; i<=n; i++){
    if(!isp[i]){
      prm[++pc]=i;
      Sm[0][pc]=Sm[0][pc-1]+(i%4==1);
      Sm[1][pc]=Sm[1][pc-1]+(i%4==3);
    }
    for(int j=1; j<=pc; j++){
      if(i*prm[j]>n) break;
      isp[i*prm[j]]=true;
      if(i%prm[j]==0) break;
    }
  }
}
int idx(ll x){ return x<=S ? A[x] : B[n/x]; }
int work(ll n, int x){
  if(n<=1 || n<prm[x]) return 0;
  int ans = mul(pw[3],dec(f[0][idx(n)],Sm[0][x-1]));
  Add(ans, dec(f[1][idx(n)],Sm[1][x-1]));
  for(int i=x; i<=pc; i++){
    if((ll)prm[i]*prm[i]>n) break;
    for(ll t=prm[i],e=1; t*prm[i]<=n; t*=prm[i],++e){
      Add(ans, mul((prm[i]%4==1)?pw[2*e+1]:1,work(n/t,i+1)));
      Add(ans, (prm[i]%4==1)?pw[2*e+3]:1);
    }
  } return ans;
}
int main(){
  #ifdef FSYolanda
  freopen("1.in","r",stdin);
  #endif
  scanf("%lld%lld%lld",&n,&k,&Mod);
  linear_sieve(S=sqrt(n)); 
  for(int i=1; i<=120; i++) pw[i]=ksm(i,k);
  for(ll l=1,r,v;l<=n;l=r+1){
    v=n/l; r=n/v; if(v<=S) A[v]=++nd; 
    else B[n/v]=++nd; ::v[nd]=v;
  } for(int i=1; i<=nd; i++) // not involve 1 & 2 
  f[0][i]=fx((v[i]-1)/4), f[1][i]=fx((v[i]+1)/4);
  for(int i=1; i<=pc; i++) if(i>1)
  for(int j=1; j<=nd; j++){
    if((ll)prm[i]*prm[i]>v[j]) break;
    int k=idx(v[j]/prm[i]);
    if(prm[i]%4==1){
      Dec(f[0][j],dec(f[0][k],Sm[0][i-1]));
      Dec(f[1][j],dec(f[1][k],Sm[1][i-1]));
    } else{
      Dec(f[0][j],dec(f[1][k],Sm[1][i-1]));
      Dec(f[1][j],dec(f[0][k],Sm[0][i-1]));
    }
  } int ans = work(n,1) + 2;
  cout << mul(ans,pw[4]); return 0;
}