「2019 集训队互测 Day 1」整点计数（min25）

传送门

• 以下整点皆在第一象限讨论
  令 $f(x)={\sum }_{i\ge 1}{\sum }_{j\ge 0}[i^2+j^2=x^2]$，要求的就是 $4^k{\sum }_{i=1}^{N}f(i{)}^k$
  考虑把整点 $(x,y)$ 用复平面的高斯整数 $x+yi$ 表示，而 $(x+yi)(x-yi)=x^2+y^2$，于是求 $f(x)$ 变成了求有多少个高斯整数使得它与它的共轭复数乘积为 $x^2$

• 定理：高斯整数可以写成一系列不能继续分解的高斯整数的乘积，并且为唯一表示（$xyx$ 论文8.2）

费马平方和定理： 奇质数 $p$ 可以表示为两个正整数 $a,b$ 的平方和 $p=a^2+b^2$ ，当且
仅当$p\equiv 1(mod4)$ ，并且这样的表示若存在，在不计较排列顺序的情况下，是唯一的
证明：
首先 $a^2$ 模 $4$ 只可能为 $1$ 故 $a^2+b^2$ 不可能为 $4k+3$ 的质数
引理1：设 $p$ 是质数，那么 $(p-1)!\equiv -1(modp)$，$p=2$ 时成立，否则考虑乘法逆元不为本身的数
引理2：设质数 $p$ 满足 $p\equiv 1(mod4)$，那么存在 $x$ 满足 $p|x^2+1$，注意到 $ab=(p-a)(p-b)$，因为 $p-1$ 为 4 的倍数，所以我们可以将 $1\ast 2,(p-1)\ast (p-2)\dots$ 进行配对，可以得到 $(\frac{p-1}{2}){!}^2+1\equiv (p-1)!+1(modp)$
由引理 2，对于 $4k+1$ 形的质数 $p$，我们知道 $p|(x+i)(x-i)$，设 $p$ 为高斯质数，有 $p|x+i$ 或 $p|x-i$，而 $\frac{x+i}{p}$ 并不是高斯整数，所以不成立，故我们可以令 $p=uv$，有 $N(p)=N(uv)=N(u)N(v)$，由于 $p$ 是质数，所以 $N(u),N(v)=p$，又 $uv$ 是一个整数，所以 $u,v$ 共轭，设 $u=a+bi,v=a-bi$，$p=(a+bi)(a-bi)=a^2+b^2$
下面来证明是唯一表示，设 $p=a^2+b^2=c^2+d^2$，若 $a+bi$ 为高斯质数，$a+bi|c+di$，或 $a+bi|c-di$，又 $N(a+bi)=N(c\pm di)=p$，故 $a^2+b^2=c^2+d^2$ 是同种表示，否则令 $a+bi=uv$，有 $(a+bi)(a-bi)=(uv)(\overline{uv})=(u\overline{v})(\overline{u}v)$

• 记 $way(p^2)=f(p)$ 为 $p^2$ 的分解方式，容易发现 $f(p)$ 是个积性函数，我们要做的就是求这个积性函数的前缀和，考虑在质数及质数次幂处的取值
• 对于 $4k+1$ 的质数 $p$，$p^{2k}$ 的值是 $(a+bi{)}^{2k}(a-bi{)}^{2k}$，容易发现分解方式为 $2k+1$ 种
  对于 $4k+3$ 的质数 $p$，$p^{2k}$ 存在唯一的分解方式 $p^k\ast p^k$
  容易发现 $4k+3$ 的至少与 $2$ 的贡献均为 1
  那么我们可以知道（令 $n=\prod p_i^{a_i}$）
  $$f(n)=\prod (1+2\ast a_i\ast [p_i\equiv 1(mod4)])$$
  可以用 $min25$ 简单计算

#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;
}