写一下这道题的类欧部分。

$$y\equiv kx\pmod p$$

这个方程在 $x,y\in[1,n\bmod p]$ 的范围内的整数解的个数便是题解所说的 $c_k$。记 $m=n\bmod p$。

也就是说 $kx-y=tp$ 其中 $t$ 为整数。不难发现，有多少个合法的 $t$ 就有多少个合法的整数解。

对于一个 $t$ 其合法的充分必要条件是 $kx-tp\in[1,m]$。注意到 $kx\equiv 0\pmod p$ 不成立。所以一个 $t$ 合法当且仅当 $kx-m\le tp\le kx$。也就是说 $t$ 满足 $\lceil\frac{kx-m}{p}\rceil\le t\le \lfloor\frac{kx}{p}\rfloor$。注意到对于固定的 $x$ 这个区间左右区间至多差 1。

所以我们将对于固定的 $x$ 的 $t$ 的数量表示为 $\lfloor\frac{kx}{p}\rfloor-\lceil\frac{kx-m}{p}\rceil+1$。 所以 $c_k=\sum\limits_{x=1}^m \lfloor\frac{kx}{p}\rfloor-\lceil\frac{kx-m}{p}\rceil+1$。我们将这个难看的上取整去掉。

$$\lceil\frac{a}{p}\rceil=\lfloor\frac{a}{p}\rfloor+1-[p|a] $$

带入上式 $c_k=\sum\limits_{x=1}^m \lfloor\frac{kx}{p}\rfloor-\sum\limits_{x=1}^m\lfloor\frac{kx-m}{p}\rfloor+[k^{-1}m\bmod p\le m]$。

很显然这个是 $\sum\limits_{x=0}^n\lfloor\frac{ax+b}{c}\rfloor$ 的类欧几里得的板子了。

记住：不要写unordered_map的记忆化，不写记忆化也可以随便跑。

提供一个这题用的类欧的板子。

long long f(int a,int b,int c,int n) {
	if(!a)return 1ll*b/c*(n+1);
//	unsigned long long has=gethas(a,b,c,n);
//	if(mem.find(has)!=mem.end())return mem[has];
	if(a>=c||b>=c)return /*mem[has]=*/1ll*(1ll*n*(n+1)/2)
	*(a/c)+1ll*(n+1)*(b/c)+f(a%c,b%c,c,n);long long m=(1ll*a*n+b)/c;
	return /*mem[has]=*/1ll*n*m-f(c,c-b-1,a,m-1);
}

注意不要爆int或long long 附上一份正解代码对拍使用

```
#include<bits/stdc++.h>
using namespace std;
const int N=5e6+5;
typedef long long ll;
int mod;ll n;
int down(const int&x){
    if(x>=mod)return x-mod;
    return x;
}
struct Modint{
    int v;
    Modint(int x=0){v=x;}
    Modint operator+(const Modint&x)const{return down(v+x.v);}
    Modint operator-(const Modint&x)const{return down(v+mod-x.v);}
    Modint operator*(const Modint&x)const{return 1ll*v*x.v%mod;}
    Modint operator=(const Modint&x){return v=x.v;}
    void operator+=(const Modint&x){v=down(v+x.v);}
    void operator-=(const Modint&x){v=down(v+mod-x.v);}
    void operator*=(const Modint&x){v=1ll*v*x.v%mod;}
};
struct Complex{
    Modint re,im;
    Complex(Modint x=0,Modint y=0){re=x,im=y;}
    Complex(int x,int y){re=x,im=y;}
    Complex operator+(const Complex&x)const{return Complex(re+x.re,im+x.im);}
    Complex operator-(const Complex&x)const{return Complex(re-x.re,im-x.im);}
    Complex operator*(const Complex&x)const{return Complex(re*x.re-im*x.im,re*x.im+im*x.re);}
    Complex operator=(const Complex&x){re=x.re,im=x.im;return *this;}
    void operator+=(const Complex&x){re+=x.re,im+=x.im;}
    void operator-=(const Complex&x){re-=x.re,im-=x.im;}
    void operator*=(const Complex&x){*this=*this*x;}
    int Im()const{return im.v;}
    int Re()const{return re.v;}
}ans;
typedef long long ll;
ll calc(int a,int b,int c,int n){
    if(!a)return 1ll*b/c*(n+1);
    if(a>=c||b>=c)return 1ll*n*(n+1)/2*(a/c)+1ll*(n+1)*(b/c)+calc(a%c,b%c,c,n);
    ll m=(1ll*a*n+b)/c;return 1ll*n*m-calc(c,c-b-1,a,m-1);
}
namespace Find_Yuangeng{
    const int MAXN=N;
    int t,p,cnt,tot,ctans,fc[MAXN],ans[MAXN],pri[MAXN],rt[MAXN],q[MAXN],phi[MAXN];
    void init () {
        phi[1]=1;
        for (int i=2;i<=MAXN-5;i++) {
            if (!q[i]) {pri[++tot]=i,phi[i]=i-1;}
            for (int j=1;j<=tot&&pri[j]*i<=MAXN-10;j++) {
                q[i*pri[j]]=1;
                if (i%pri[j]==0) {
                    phi[i*pri[j]]=phi[i]*pri[j];
                    break;
                }
                phi[i*pri[j]]=phi[i]*(pri[j]-1);
            }
        }
        return;
    }
    int qpow (int a,int b,int p) {
        int res=1;
        while (b) {
            if (b&1) {res=(1ll*res*a)%p;}
            a=(1ll*a*a)%p;
            b>>=1;
        }
        return res;
    }
    bool chk (int x,int p) {
        if (qpow(x,phi[p],p)!=1) {return 0;}
        for (int i=1;i<=cnt;i++) {
            if (qpow(x,phi[p]/fc[i],p)==1) {return 0;}
        }
        return 1;
    }
    void proc (int p) {
        for (int i=2;i*i<=p;i++) {
            if (p%i==0) {
                fc[++cnt]=i;
                while (p%i==0) {p/=i;}
            }
        }
        if (p>1) {fc[++cnt]=p;}
        return;
    }
    int findrt (int p) {
        ctans=cnt=0;
        proc(phi[p]);
        for (int i=1;i<p;i++) {
            if (chk(i,p)) {return i;}
        }
        return 0;
    }
}
int g;
int val[N],pw[N];
int qpow(int a,int b){
    a%=mod;
    if(!a)return 0;
    return pw[1ll*val[a]*b%(mod-1)];
}
int inv(int x){
    return qpow(x,mod-2);
}
ll solve(int k,int m){
    ll ans=((inv(k)*m%mod)<=m);ans+=m;
    ans+=calc(k,0,mod,m);
    ans-=calc(k,mod-m,mod,m);
    return ans%(mod-1);
}
int main(){
    Find_Yuangeng::init();
    scanf("%d%lld",&mod,&n);
    if(mod%4==1&&n>=mod){
        puts("0 0");
        return 0;
    }
    g=Find_Yuangeng::findrt(mod);
    val[1]=0;pw[0]=1;
    for(int i=1;i<mod-1;i++){
        pw[i]=1ll*pw[i-1]*g%mod;
		val[pw[i]]=i;
    }
    ll lim=n/mod,mul=1ll,m=n-lim*mod;
    for(int i=1;i<=m;i++)mul=4ll*mul*i*i%mod;
    mul=qpow(mul,lim%(mod-1));int t=lim*m%4;
    ans.re=mul,ans.im=0;
    for(int i=0;i<t;i++)ans*=Complex(0,1);
    for(int i=1;i<=m;i++)ans*=Complex(0,1ll*i*i%mod);
    if(m){
        ll nw_mul=1;
        for(int i=1;i<=m;i++)nw_mul=1ll*nw_mul*i%mod;
        nw_mul=qpow(nw_mul,m);
        ans*=Complex(nw_mul,0);
        for(int i=1;i<=m;i++)ans*=Complex(1,1);
        //(1+i)^m
        for(int i=2;i<mod-1;i++)if(inv(i)>i){
            int x=inv(i),val=solve(i,m);
            ans*=Complex(qpow(x+i,val),0);
            for(int i=0;i<val%4;i++)ans*=Complex(0,1);
        }
        int val=solve(mod-1,m);
        for(int i=1;i<=val;i++)ans*=Complex(1,mod-1);
    }
    printf("%d %d\n",ans.Re(),ans.Im());
    return 0;
}
```