传送门
- 以下整点皆在第一象限讨论
令,要求的就是
考虑把整点用复平面的高斯整数
表示,而
,于是求
变成了求有多少个高斯整数使得它与它的共轭复数乘积为
- 定理:高斯整数可以写成一系列不能继续分解的高斯整数的乘积,并且为唯一表示(
论文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;
}