[ICPC2022Macau]Pass the Ball!【双模术NTT】【分块】

分析：看到$\sum i*b_i$的式子以及每次传球向右移动一位，很容易想到HNOI2017项链的循环卷积。假设我们用FFT处理出来了每个环传0次，传1次，传2次，...，传环的大小n-1次对应的结果。那么对于每个询问k，只需要对于每个环，找出k模环的大小对应的答案即可。但是如果对每个环都找，时间复杂度是吃不消的。所以采用分块的方法，对于小于$\sqrt(n)$的环，把答案加进一个预处理数组中，对于大于$\sqrt(n)$的环，直接枚举所有询问，把结果加进去。最后再对每个询问枚举小于$\sqrt(n)$的环的累计答案。

#include<bits/stdc++.h>
using namespace std;

const int maxn = 102000;
const int mod1 = 998244353;
const int mod2 = 1004535809;

int n,m,num;
int p[maxn],vis[maxn],query[maxn];
vector<int> vec;
vector<int> ans1[500]; // xiaoyu sqrtn de daan
vector<int> ans2[500];

int res1[maxn],res2[maxn]; //zuizhongdaan

int A[maxn*16],B[maxn*16],C[maxn*16];
int r[maxn*16],len = 1,bit=0;
void build_P(){
    for(int i=0;i<9*num;i++) A[i] = B[i] = C[i] = 0;
    for(int i=0;i<num;i++) A[num-i-1] = vec[i];
    for(int i=0;i<num;i++) B[i] = B[num+i] = vec[i];
}

int fast_pow(int dt,int pw,int mod){
    int res=1,now=dt,p=1;
    while(p<=pw){
    if(p&pw) res = 1ll*res*now%mod;
    now=1ll*now*now%mod;
    p<<=1;
    }
    return res;
}

void FFT(int *a,int flag,int mod){
    for(int i=0;i<len;i++) if(i < r[i]) swap(a[i],a[r[i]]);
    for(int mid = 1;mid<len;mid<<=1){
    int nw = fast_pow(3,(mod-1)/(2*mid),mod);
    if(flag == -1) nw = fast_pow(nw,mod-2,mod);
    for(int i=0;i<len;i+=(mid<<1)){
        int tmp = 1;
        for(int j=0;j<mid;j++,tmp=1ll*tmp*nw%mod){
        int x = a[i+j],y=1ll*tmp*a[i+j+mid]%mod;
        a[i+j]=(x+y)%mod,a[i+j+mid]=(x-y+mod)%mod;
        }
    }
    }
    if(flag == -1){
    int dt = fast_pow(len,mod-2,mod);
    for(int i=0;i<len;i++){
        a[i] = 1ll*a[i]*dt%mod;
    }
    }
}

void fft(int mod){
    len = 1,bit = 0;
    while(len < 4*num){len<<=1,bit++;}
    for(int i=0;i<len;i++) r[i] = (r[i>>1]>>1)|((i&1)<<(bit-1));
    FFT(A,1,mod); FFT(B,1,mod);
    for(int i=0;i<len;i++) C[i] = 1ll*A[i]*B[i]%mod;
    FFT(C,-1,mod);
}

long long xxx,yyy;
void exgcd(int a,int b){
    if(!b){
    xxx=1;yyy=0;return;
    }
    exgcd(b,a%b);
    long long t = xxx; xxx = yyy; yyy = t-a/b*yyy;
}

long long fast_mul(long long dt,long long pw,long long mod){
    long long res=0,now=dt,p=1;
    while(p<=pw){
    if(p&pw) res = (res+now)%mod;
    now=(now+now)%mod;
    p<<=1;
    }
    return res;
}

long long merge(int b,int d){
    long long mod = 1ll*mod1*mod2;
    long long as;
    if(d-b > 0){
    as = fast_mul(fast_mul(xxx,(d-b),mod),mod1,mod)+b;
    as %= mod;
    as += mod; as%= mod;
    }else{
    long long xx = -xxx;
    as = fast_mul(fast_mul(xx,b-d,mod),mod1,mod)+b;
    as %= mod;
    as += mod; as%= mod;
    }
    return as;
}

int main(){
    //freopen("1.in","r",stdin);
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for(int i=1;i<=n;i++){
    cin >> p[i];
    }
    for(int i=1;i<=m;i++) cin >> query[i];
    
    for(int i=1;i*i<=n;i++){
    for(int j=0;j<i;j++) ans1[i].push_back(0),ans2[i].push_back(0);
    }
    
    for(int i=1;i<=n;i++){
    if(!vis[i]){
        int j = i;vec.clear();num = 0;
        while(!vis[j]){num++;vec.push_back(j);vis[j]=1;j = p[j];}
        build_P();
        fft(mod1);
        if(num*num <= n){
        for(int j=0;j<num;j++){
            ans1[num][j]+=C[num-1+j];
            ans1[num][j]%=mod1;
        }
        }else{
        for(int j=1;j<=m;j++){
            int dt = query[j]%num;
            res1[j] += C[num-1+dt];
            res1[j]%=mod1;
        }
        }
        build_P();
        fft(mod2);
        if(num*num <= n){
        for(int j=0;j<num;j++){
            ans2[num][j]+=C[num-1+j];
            ans2[num][j]%=mod2;
        }
        }else{
        for(int j=1;j<=m;j++){
            int dt = query[j]%num;
            res2[j] += C[num-1+dt];
            res2[j]%=mod2;
        }
        }
    }
    }
    for(int i=2;i*i<=n;i++){
    for(int j=1;j<=m;j++){
        res1[j] += ans1[i][query[j]%i];
        res1[j]%=mod1;
        res2[j] += ans2[i][query[j]%i];
        res2[j]%=mod2;
    }
    }
    exgcd(mod1,mod2);
    for(int i=1;i<=m;i++){
    //cout<<res1[i]<<" "<<res2[i]<<endl;
    cout<<merge(res1[i],res2[i])<<endl;
    }
    return 0;
}