-
以下整点皆在第一象限讨论
令 f ( x ) = ∑ i ≥ 1 ∑ j ≥ 0 [ i 2 + j 2 = x 2 ] f(x)=\sum_{i\ge 1}\sum_{j\ge 0}[i^2+j^2=x^2] f(x)=∑i≥1∑j≥0[i2+j2=x2],要求的就是 4 k ∑ i = 1 N f ( i ) k 4^k\sum_{i=1}^Nf(i)^k 4k∑i=1Nf(i)k
考虑把整点 ( x , y ) (x,y) (x,y) 用复平面的高斯整数 x + y i x+yi x+yi 表示,而 ( x + y i ) ( x − y i ) = x 2 + y 2 (x+yi)(x-yi)=x^2+y^2 (x+yi)(x−yi)=x2+y2,于是求 f ( x ) f(x) f(x) 变成了求有多少个高斯整数使得它与它的共轭复数乘积为 x 2 x^2 x2 -
定理:高斯整数可以写成一系列不能继续分解的高斯整数的乘积,并且为唯一表示( x y x xyx xyx 论文8.2)
费马平方和定理: 奇质数 p p p 可以表示为两个正整数 a , b a,b a,b 的平方和 p = a 2 + b 2 p=a^2 +b^2 p=a2+b2 ,当且
仅当 p ≡ 1 ( m o d 4 ) p \equiv1 (mod\ 4) p≡1(mod 4) ,并且这样的表示若存在,在不计较排列顺序的情况下,是唯一的
证明:
首先 a 2 a^2 a2 模 4 4 4 只可能为 1 1 1 故 a 2 + b 2 a^2+b^2 a2+b2 不可能为 4 k + 3 4k+3 4k+3 的质数
引理1:设 p p p 是质数,那么 ( p − 1 ) ! ≡ − 1 ( m o d p ) (p-1)!\equiv -1(mod\ p) (p−1)!≡−1(mod p), p = 2 p=2 p=2 时成立,否则考虑乘法逆元不为本身的数
引理2:设质数 p p p 满足 p ≡ 1 ( m o d 4 ) p\equiv 1(mod\ 4) p≡1(mod 4),那么存在 x x x 满足 p ∣ x 2 + 1 p|x^2+1 p∣x2+1,注意到 a b = ( p − a ) ( p − b ) ab = (p-a)(p-b) ab=(p−a)(p−b),因为 p − 1 p-1 p−1 为 4 的倍数,所以我们可以将 1 ∗ 2 , ( p − 1 ) ∗ ( p − 2 ) … 1*2,(p-1)*(p-2)\dots 1∗2,(p−1)∗(p−2)… 进行配对,可以得到 ( p − 1 2 ) ! 2 + 1 ≡ ( p − 1 ) ! + 1 ( m o d p ) (\frac{p-1}{2})!^2+1\equiv (p-1)!+1(mod\ p) (2p−1)!2+1≡(p−1)!+1(mod p)
由引理 2,对于 4 k + 1 4k+1 4k+1 形的质数 p p p,我们知道 p ∣ ( x + i ) ( x − i ) p|(x+i)(x-i) p∣(x+i)(x−i),设 p p p 为高斯质数,有 p ∣ x + i p|x+i p∣x+i 或 p ∣ x − i p|x-i p∣x−i,而 x + i p \frac{x+i}{p} px+i 并不是高斯整数,所以不成立,故我们可以令 p = u v p=uv p=uv,有 N ( p ) = N ( u v ) = N ( u ) N ( v ) N(p)=N(uv)=N(u)N(v) N(p)=N(uv)=N(u)N(v),由于 p p p 是质数,所以 N ( u ) , N ( v ) = p N(u),N(v)=p N(u),N(v)=p,又 u v uv uv 是一个整数,所以 u , v u,v u,v 共轭,设 u = a + b i , v = a − b i u=a+bi,v=a-bi u=a+bi,v=a−bi, p = ( a + b i ) ( a − b i ) = a 2 + b 2 p=(a+bi)(a-bi)=a^2+b^2 p=(a+bi)(a−bi)=a2+b2
下面来证明是唯一表示,设 p = a 2 + b 2 = c 2 + d 2 p=a^2+b^2=c^2+d^2 p=a2+b2=c2+d2,若 a + b i a+bi a+bi 为高斯质数, a + b i ∣ c + d i a+bi|c+di a+bi∣c+di,或 a + b i ∣ c − d i a+bi|c-di a+bi∣c−di,又 N ( a + b i ) = N ( c ± d i ) = p N(a+bi)=N(c\pm di)=p N(a+bi)=N(c±di)=p,故 a 2 + b 2 = c 2 + d 2 a^2+b^2=c^2+d^2 a2+b2=c2+d2 是同种表示,否则令 a + b i = u v a+bi=uv a+bi=uv,有 ( a + b i ) ( a − b i ) = ( u v ) ( u v ‾ ) = ( u v ‾ ) ( u ‾ v ) (a+bi)(a-bi)=(uv)(\overline {uv})=(u\overline v)(\overline uv) (a+bi)(a−bi)=(uv)(uv)=(uv)(uv)
- 记 w a y ( p 2 ) = f ( p ) way(p^2)=f(p) way(p2)=f(p) 为 p 2 p^2 p2 的分解方式,容易发现 f ( p ) f(p) f(p) 是个积性函数,我们要做的就是求这个积性函数的前缀和,考虑在质数及质数次幂处的取值
- 对于
4
k
+
1
4k+1
4k+1 的质数
p
p
p,
p
2
k
p^{2k}
p2k 的值是
(
a
+
b
i
)
2
k
(
a
−
b
i
)
2
k
(a+bi)^{2k}(a-bi)^{2k}
(a+bi)2k(a−bi)2k,容易发现分解方式为
2
k
+
1
2k+1
2k+1 种
对于 4 k + 3 4k+3 4k+3 的质数 p p p, p 2 k p^{2k} p2k 存在唯一的分解方式 p k ∗ p k p^k*p^k pk∗pk
容易发现 4 k + 3 4k+3 4k+3 的至少与 2 2 2 的贡献均为 1
那么我们可以知道(令 n = ∏ p i a i n=\prod p_i^{a_i} n=∏piai)
f ( n ) = ∏ ( 1 + 2 ∗ a i ∗ [ p i ≡ 1 ( m o d 4 ) ] ) f(n)=\prod(1+2*a_i*[p_i\equiv 1(mod\ 4)]) f(n)=∏(1+2∗ai∗[pi≡1(mod 4)])
可以用 m i n 25 min25 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;
}