「理性愉悦」LOJ3069 集训队互测2019 整点计数

现场居然能过掉 (4) 个人，候选队实力恐怖如斯

面向数据范围不难考虑对 (f(x)) 是否是积性函数进行思考，而根据定义我们发现圆上的整点是平均分布在四个象限里面的，那么对 (f(i)/4) 打表

好耶，是积性函数！但是如何证明呢？

首先证明两个结论：

• 奇质数 (x) 可以被分解成 (a^2+b^2) 当且仅当 (xequiv 1 (mod 4))

  根据二次剩余的一些理论，得到此时存在一个 (a^2 equiv -1 (mod x))

  那么构造 (as_1+t_1equiv as_2+t_2mod p) 满足 (s_1,s_2,t_1,t_2le sqrt p)

  此时满足条件的的 (x=|s_1-s_2|,y=|t_1-t_2|)，证明带入即可，那么发现 (x^2+y^2 equiv 0mod p,x^2+y^2le p)

定义高斯整数 (x+yi[x,yinmathrm{Interger}]) 和其范数 (N(x+y_i)=x^2+y^2)

同时定义高斯质数为不能表示为若干个高斯质数乘积的数，这里使用上面的结论，发现高斯质数就是 (p equiv 3mod 4) 的质数

定义两个不同的高斯整数是相互伴随的当且仅当一个通过旋转 (90,180,270) 可以成为另一个

考虑这样一个性质：对于高斯整数，也存在唯一分解定律，这里分解出来的结果每个数可以替换为其伴随

这个性质考虑归奶证明，使用范数的乘积可以快速得到

那么这个函数是积性函数就完成了证明，同时不难发现质数处取值：

[g(p^k)=egin{cases}(2p+1)^k[pequiv 1mod 4]\1 [pequiv 3mod 4]end{cases} ]

剩下的部分不难使用 (min25) 筛进行求解，但是一个需要分开讨论

在进行 (min25) 筛的第一个 ( m{DP}) 的过程中，使用 (3 imes 3equiv 1) 的优美性质来进行 (pequiv 3mod 4) 的部分的转移

具体实现的时候分 (g_1(n,i),g_2(n,i)) 两个部分进行转移，一开始转移的都是个数，并不计算具体点值

剩下的都交给第二个 ( m{DP}) 啦，细节并不多，复杂度 (Theta(frac{n^{0.75}}{log n}))

#include<bits/stdc++.h>
using namespace std;
#define reg register
#define int long long
#define For(i,a,b) for(reg int i=a;i<=b;++i)
#define Down(i,a,b) for(reg int i=a;i>=b;--i) 
#define ull unsigned long long
#define rep(i,a,b) for(reg int i=a;i<=b;++i)
#define ll long long
inline int max(int x,int y){return x>y?x:y;}
inline int min(int x,int y){return x<y?x:y;}
inline void ckmin(int &x,int y){x=x<y?x:y; return ;}
inline void ckmax(int &x,int y){x=x>y?x:y; return ;}
inline int gcd(int n,int m){return m?gcd(m,n%m):n;}
inline int lcm(int x,int y){return x/gcd(x,y)*y;}
inline void swap(int &x,int &y){int t=x; x=y; y=t; return ;}
int mod;
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int del(int x,int y){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y){return x*y-x*y/mod*mod;}
inline void ckadd(int &x,int y){x=x+y>=mod?x+y-mod:x+y; return ;}
inline void ckdel(int &x,int y){x=x-y<0?x-y+mod:x-y; return ;}
inline void ckmul(int &x,int y){x=x*y-x*y/mod*mod; return ;}
inline int ksm(int x,int y){int res=1; for(;y;y>>=1,ckmul(x,x)) if(y&1) ckmul(res,x); return res;}
namespace yspm{
	inline int read(){
	    int res=0,f=1; char k;
	    while(!isdigit(k=getchar())) if(k=='-') f=-1;
	    while(isdigit(k)) res=res*10+k-'0',k=getchar();
	    return res*f;
    }
    char out[100];
    inline void print(int x){
        if(!x) return putchar('0'),putchar('
'),void();
        if(x<0) putchar('-'),x=-x; 
        int cnt=0; while(x) out[++cnt]=x%10,x/=10; 
        Down(i,cnt,1) putchar(out[i]+'0'); putchar('
'); return ;
    } 
    const int N=1e6+10;
    int id1[N],k,cnt,id2[N],R[N],tot,pri[N],s1[N],f[100],s3[N],g1[N],bl,n,g3[N];
    bool fl[N];
    inline int id(int x){return x<=bl?id1[x]:id2[n/x];}
    inline int calc(int n,int x){
        if(n<pri[x]) return 0; int res=add(mul(f[1],del(g1[id(n)],s1[x-1])),del(g3[id(n)],s3[x-1]));
        for(reg int j=x;pri[j]*pri[j]<=n;++j){
            int tim=1,prd=pri[j]*pri[j];
            while(prd<=n){
                ckadd(res,add(pri[j]%4==1?f[tim+1]:1,mul(pri[j]%4==1?f[tim]:1,calc(n*pri[j]/prd,j+1))));
                tim++; prd*=pri[j];
            }
        } return res;
    }
	signed main(){
	    n=read(); k=read(); mod=read(); bl=sqrt(n);
	    if(n==1) print(ksm(4,k)),exit(0);
	    rep(i,2,N-1){
	        if(!fl[i]){
	            pri[++cnt]=i; s1[cnt]=s1[cnt-1]; s3[cnt]=s3[cnt-1];
	            if(i%4==3) s1[cnt]=s1[cnt-1],s3[cnt]=s3[cnt-1]+1;
	            if(i%4==1) s1[cnt]=s1[cnt-1]+1,s3[cnt]=s3[cnt-1];
	        } for(reg int j=1;j<=cnt&&i*pri[j]<N;++j){
	            fl[i*pri[j]]=1; if(i%pri[j]==0) break;
	        }
	    }
	    for(reg int l=1,r;l<=n;l=R[tot]+1){
	        R[++tot]=n/(n/l),g1[tot]=((R[tot]-1)/4)%mod,g3[tot]=((R[tot]+1)/4)%mod;
	        if(R[tot]<=bl) id1[R[tot]]=tot; else id2[n/R[tot]]=tot;
	    }
	    rep(i,2,cnt){
	        if(pri[i]*pri[i]>n) break;
	        if(pri[i]%4==1){
	            for(reg int j=tot;pri[i]*pri[i]<=R[j];--j){
	                int x=id(R[j]/pri[i]);
	                ckdel(g1[j],del(g1[x],s1[i-1]));
	                ckdel(g3[j],del(g3[x],s3[i-1]));
	            }
	        }else{
	            for(reg int j=tot;pri[i]*pri[i]<=R[j];--j){
	                int x=id(R[j]/pri[i]);
	                ckdel(g1[j],del(g3[x],s3[i-1]));
	                ckdel(g3[j],del(g1[x],s1[i-1]));
	            }
	        }
	    } f[0]=1; rep(i,1,50) f[i]=ksm(i*2+1,k); print(mul(ksm(4,k),add(calc(n,1),2))); return 0;
    }
}signed main(){return yspm::main();}
//Use The Time To Enrich This Selfclosing Youth