YY模拟：和（CRT+倍增）

题意：
给定$m,T$,$T$次询问，每次询问给定$N$,求 ：

$$\sum_{i=1}^n i^k \pmod{m}$$

其中$n,m,k \le 1e18 ,T\le 3000$，且$m$的最大质因子$p$不超过$3e5$。

题解：

设

$$m=\prod_{i=1}^{s}p_i^{a_i}$$
我的做法是 $O(T\log^3 n + \max_{i=1}^s p_ia_i)$，满数据要跑70s+，标算据说少一个 $\log n$，不过并不会。

自然数幂和在模数较大的情况下的最优复杂度一般为$O(k)$（拉格朗日插值）。不过注意到$p_i \le 3e5$，这就是让我们用$\text{CRT}$来做了。

问题转化为

$$\sum_{i=1}^n i^k \pmod{p^e}$$

记$S=\lfloor \frac{n}{p}\rfloor$，答案转化为两部分：
1.

$$\sum_{i=1}^{n \text{%} p}(Sp+i)^k$$
二项式展开不在赘述，保留 $p$的次数小于$e$的项，为：
$$\sum_{z=0}^{e-1}\binom{k}{z}S^zp^z \sum_{i=1}^{n \text{%} p}i^{k-z}$$
预处理后面部分，做到 $O(e)$回答。

2.

$$\sum_{i=0}^{S-1}\sum_{j=1}^p (ip+j)^k$$

同样二项式展开，问题转化为求：

$$\sum_{z=0}^{e-1} \binom{k}{z} (\sum_{i=0}^{S-1}i^z)(\sum_{j=1}^p j^{k-z})$$

后面$j$的部分同样预处理，中间的部分$\sum_{i=0}^{S-1}i^z$要单独考虑。

注意在不存在逆元的前提下，插值等方法基本作废，那怎么办？

一种直观的想法是递归下去解决。 不过这个复杂度是爆炸的（粗略的分析了一下应该是$e!$），怎么办 ？观察到$z$很小， 我们有一种比较冷门的自然数幂和方法：倍增！

记$f_{n,k} = \sum\limits_{i=1}^n i^k \pmod{m}$，有：

$$f_{n,k} = \begin{cases} f_{n-1,k}+n^k & (n \equiv 1 \pmod{2}) \\ f_{l,k}+poly(l,k) & (n \equiv 0 \pmod{2}) \\ 0 & (n=0) \end{cases}$$

其中

$$l=\frac{n}{2}\\poly(l,k)= \sum_{i=1}^l (i+l)^k$$
同样二项式展开得
$$\begin{aligned} poly(l,k) &=& \sum_{i=1}^l\sum_{j=0}^{k}\binom{k}{j}i^kj^{k-j} \\&=& \sum_{j=0}^k \binom{k}{j}j^{k-j} (\sum_{i=1}^li^k) \end{aligned}$$

变为子问题， 记忆化这个过程，那么倍增共有$\log n$层，每层共有$k$个项，枚举$k$次， 时间复杂度为$O(k^2 \log n)$。

于是我们在$O(T \log^3 n + pe)$的时间内解决了。
TLE代码：

#include <bits/stdc++.h>
#include <tr1/unordered_map>
using namespace std;
typedef __int128 IL;
typedef long long LL;
typedef pair <LL,LL> pii;

const int N=3e5+50, L=100;

inline LL add(LL x,LL y,LL mod) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
inline LL mul(LL x,LL y,LL mod) {return (LL)((IL)x*y%mod);;}
inline LL power(LL a,LL b,LL mod,LL rs=1) {for(;b;b>>=1,a=mul(a,a,mod)) if(b&1) rs=mul(rs,a,mod);return rs;}
template <typename T> inline T gcd(T x,T y) {return (y ? gcd(y,x%y) : x);}

inline LL rd() {
    char ch=getchar(); LL i=0,f=1;
    while(!isdigit(ch)) {if(ch=='-')f=-1; ch=getchar();}
    while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=getchar();}
    return i*f;
}

LL m,K,T; 
LL mnpr[N],npr[N],pr[N],pt,C[L][L],R[L],f[L][N],gp[L][N],mxp,mxl;
vector <LL> fac;
vector <pii> factor;
tr1::unordered_map <LL,LL> g[L];

inline void sieve() {
    mnpr[1]=1;
    for(int i=2;i<N;++i) {
        if(!npr[i]) pr[++pt]=i, mnpr[i]=i;
        for(int j=1;j<=pt;j++) {
            LL k=i*pr[j];
            if(k>=N) break;
            mnpr[k]=pr[j]; npr[k]=1;
            if(!(i%pr[j])) break;
        }
    }
}

inline LL getC(LL n,LL a,LL mod) {
    vector <LL> vec;
    for(int i=1;i<=a;i++) vec.push_back(n-i+1);
    for(int i=1;i<=a;i++) {
        LL rs=i, t;
        for(int j=0;j<vec.size();++j)
            if((t=gcd(rs,vec[j]))!=1) vec[j]/=t, rs/=t;
    }
    LL rs=1;
    for(int j=0;j<vec.size();++j) rs=mul(rs,vec[j],mod);
    return rs;
}

inline void init() {
    for(int i=0;i<=mxl && i<=K;++i) {
        for(int j=1;j<=mxp;++j) {
            f[i][j]=(mnpr[j]==j) ? power(j,K-i,m) : mul(f[i][mnpr[j]],f[i][j/mnpr[j]],m);
            gp[i][j]=(mnpr[j]==j) ? power(j,i,m) : mul(gp[i][mnpr[j]],gp[i][j/mnpr[j]],m);
        }
        for(int j=1;j<=mxp;++j)
            f[i][j]=add(f[i][j],f[i][j-1],m), gp[i][j]=add(gp[i][j],gp[i][j-1],m);
    }
    for(int i=0;i<=mxl && i<=K;++i) R[i]=getC(K,i,m);
    for(int i=0;i<=mxl;++i)
        for(int j=0;j<=i;++j)
            C[i][j]=getC(i,j,m);
}

IL xl,yl;
inline void exgcd(IL a,IL b,IL &x,IL &y) {b ? (exgcd(b,a%b,y,x), y-=a/b*x) : (x=1, y=0);}  
inline void merge(IL &A,IL &B,IL C,IL D) {
    exgcd(A,C,xl,yl);
    xl*=(D-B);
    IL f1=A/gcd(A,C)*C, f2=A*xl+B;
    f2=(f2%f1+f1)%f1;
    A=f1; B=f2;
}

inline LL getrs(LL x,LL p,LL rs,LL mod) {
    LL tp=0, px=1, pw=1;
    for(int z=0;z<=K && z<factor[p].second; ++z) {
        LL v1=mul(R[z],mul(px,pw,mod),mod);
        LL v2=mul(v1,f[z][rs],mod);
        tp=add(tp,v2,mod);
        px=mul(px,x,mod); 
        pw=mul(pw,factor[p].first,mod);
    } return tp;
}

inline LL G(LL n,LL z) {
    if(n<=mxp) return gp[z][n];
    if(!n) return 0;
    if(g[z].count(n)) return g[z][n];
    if(n&1) return g[z][n]=add(G(n-1,z),power(n,z,m),m);
    LL s=0, l=n/2, pw=1;
    for(int j=0;j<=z;j++) {
        LL v1=mul(C[z][j],pw,m);
        LL v2=mul(v1,G(l,z-j),m);
        s=add(s,v2,m);
        pw=mul(pw,l,m);
    } return g[z][n]=add(s,G(l,z),m);
}

inline LL getF(LL n,LL p,LL mod) {
    LL rs=n%factor[p].first; 
    rs=rs ? getrs(n/factor[p].first,p,rs,mod) : 0;
    if(n<factor[p].first) return rs;  LL tp=1;
    for(int z=0;z<=K && z<factor[p].second; ++z) {
        LL v1=mul(R[z],tp,mod);
        LL v2=mul(G(n/factor[p].first-1,z)+(!z),f[z][factor[p].first],mod);
        rs=add(rs,mul(v1,v2,mod),mod);
        tp=mul(tp,factor[p].first,mod);
    }
    return rs;
}

inline void solve(LL n) {
    IL A=1, B=0;
    for(int i=0;i<factor.size();++i) {
        LL p=power(factor[i].first,factor[i].second,m+1);
        LL v=getF(n,i,p);
        merge(A,B,p,v);
    } 
    printf("%lld\n",(LL)B);
}

int main() {
    sieve();
    m=rd(), K=rd(), T=rd();
    LL tp=m;  
    for(int i=2;i<N && tp!=1; ++i) 
        while(!(tp%i)) fac.push_back(i), tp/=i;
    LL lst=0;
    for(int i=0;i<fac.size();++i) {
        if(lst!=fac[i]) factor.push_back(pii(fac[i],1));
        else ++factor.back().second;
        mxp=max(mxp,fac[i]); mxl=max(mxl,factor.back().second); lst=fac[i];
    }
    init();
    while(T--) solve(rd());
} 

std代码

#include<algorithm>
#include<iostream>
#include<stdio.h>
#include<map>

using namespace std;

typedef unsigned long long ull;
typedef long long ll;
typedef __int128 LL;

const int N=50005,M=70,L=2000005,G=300005;
const int Base=33333,SIZE=Base+5;

struct no{
    int p;int e;
    no(){p=e=0;}
    no(int p,int e):p(p),e(e){}
    ll num(){
        ll r=1;int i;
        for(i=1;i<=e;i++)r*=p;
        return r;
    }
};

ll m,k,q,pm,c[M][M],v[N],ret[N],now[N],rmod,mod,pl[L],*pp=pl,*f[G],p[G];
bool fl[L];
no di[M];
int n,s,pr[G],pc;

ll fpow(ll a,ll t){
    ll r=1;
    for(;t;t>>=1,a=(LL)a*a%m)if(t&1)r=(LL)r*a%m;
    return r;
}
void exgcd(ll a,ll b,ll&x,ll&y){
    if(b==0)return (void)(x=1,y=0);
    exgcd(b,a%b,y,x),y-=a/b*x;
}
void crt(){
    ll m1=rmod,m2=mod,r1,r2,k1,k2;int i;
    exgcd(m1,m2,k1,k2);
    for(i=1;i<=q;i++){
        r1=ret[i];r2=now[i];
        ret[i]=((LL)k1%m2*(r2-r1)%m2+m2)%m2*m1+r1;
    }
    rmod*=mod;
}
ll Inv(ll a){
    ll x,y;exgcd(a,mod,x,y);
    return (x%mod+mod)%mod;
}
ll inv_num[M],aa[M],bb[M];
int inv_cnt[M];
ll binom(no z,ll n,ll m){
    if(n>=0&&n<m)return 0;
    ll ret=1,r;
    int tim=0,i;
    for(i=1;i<=m;i++){
        r=((n-i+1)%mod+mod)%mod;
        if(r==0)return 0;
        while(r>0&&r%z.p==0)r/=z.p,tim++;
        tim-=inv_cnt[i];
        ret=(LL)ret*r%mod*inv_num[i]%mod;
    }
    while(tim)ret=(LL)ret*z.p%mod,tim--;
    return ret;
};
ll getval(no z,ll *a,ll b){
    int i;
    ll ret=0,ml=1;
    for(i=0;i<=z.e;i++){
        if(a[i])ret=(ret+(LL)a[i]*ml)%mod;
        ml=(LL)ml*b%mod*z.p%mod;
    }
    return ret;
}

ll li[M],st[SIZE][M],tmp[M];
void move(no z,int t){
    static ll tmp[M];
    int i,j;
    tmp[0]=1;
    for(i=1;i<=z.e;i++)tmp[i]=(LL)tmp[i-1]*t%mod;
    for(i=z.e;i>=0;i--)for(j=0;j<i;j++)li[i]=(li[i]+(LL)c[i][j]*tmp[i-j]%mod*li[j])%mod;
    for(i=0;i<=z.e;i++)li[i]=(li[i]+st[t][i]-(i==0))%mod;
}
void lmove(no z,ll n,ll m){
    int i,j;
    tmp[0]=1;
    for(i=1;i<=z.e;i++)tmp[i]=(LL)tmp[i-1]*n%mod;
    for(i=z.e;i>=0;i--){
        ll w=0;
        for(j=0;j<=i;j++)w=(w+(LL)c[i][j]*li[j]%mod*tmp[i-j]%mod*st[m-1][i-j])%mod;
        li[i]=w;
    }
}
void solve(no z,ll n){
    if(n==0){
        int i;
        for(i=0;i<=z.e;i++)li[i]=0;
        return;
    }
    solve(z,n/Base);
    if(n>=Base)lmove(z,n/Base,Base);
    if(n%Base)move(z,n%Base);
}

void ext(no z){
    int i,j,l;
    mod=z.num();
    for(i=1;i<=z.e;i++){
        ll r=i;int tim=0;
        while(r>0&&r%z.p==0)r/=z.p,tim++;
        inv_num[i]=Inv(i);
        inv_cnt[i]=tim;
    }
    for(i=0;i<=z.p;i++)f[i]=pp,pp+=z.e+1;
    for(i=0;i<=z.e;i++)bb[i]=binom(z,k,i);
    for(i=1;i<=z.p;i++){
        if(i==z.p){
            for(j=0;j<=z.e&&j<=k;j++)aa[j]=fpow(i,k-j)%mod;
        }else{
            aa[0]=p[i];
            ll w=Inv(i);
            for(j=1;j<=z.e;j++)aa[j]=(LL)aa[j-1]*w%mod;
        }
        for(j=0;j<=z.e;j++)f[i][j]=(f[i-1][j]+(LL)aa[j]*bb[j])%mod;
    }
    for(i=1;i<=q;i++){
        now[i]=getval(z,f[v[i]%z.p],v[i]/z.p);
        if(v[i]>=z.p){
            solve(z,v[i]/z.p-1);
            ll w=1;
            for(j=0;j<=z.e;j++){
                if(f[z.p][j])now[i]=(now[i]+(LL)f[z.p][j]*w%mod*(li[j]+(j==0)))%mod;
                w=(LL)w*z.p%mod;
            }
        }
    }
    crt();
    while(pp>pl)*--pp=0;
}

#include<assert.h>
int main(){
    int i,j;ll t;
    scanf("%lld%lld%lld",&m,&k,&q);
    assert(k>=2);assert(m<=1E18);
    for(i=1;i<=q;i++)scanf("%lld",v+i);
    for(i=0;i<M;i++)for(j=*c[i]=1;j<=i;j++)c[i][j]=(c[i-1][j-1]+c[i-1][j])%m;
    for(i=0;i<=Base;i++)for(j=st[i][0]=1;j<M;j++)st[i][j]=(LL)st[i][j-1]*i%m; 
    for(i=1;i<=Base;i++)for(j=0;j<M;j++)st[i][j]=(st[i-1][j]+st[i][j])%m;
    for(pm=t=m,i=2;i<=t;i++)if(t%i==0){
        pm=pm/i*(i-1);
        for(j=0;t%i==0;t/=i,j++);
        di[++s]=no(i,j);
        n=max(n,i);
    }
    p[1]=1;mod=m;
    for(i=2;i<=n;i++){
        if(!fl[i]){
            pr[++pc]=i,fl[i]=true;
            p[i]=fpow(i,k);
        }
        for(j=1;i*pr[j]<=n;j++){
            p[i*pr[j]]=(LL)p[i]*p[pr[j]]%m;
            fl[i*pr[j]]=true;
            if(i%pr[j]==0)break;
        }
    }
    rmod=1;
    for(i=1;i<=s;i++)ext(di[i]);
    for(i=1;i<=q;i++)printf("%lld\n",ret[i]);
    return 0;
}