自动搬运

查看原文

来自洛谷，原作者为

	NaCly_Fish 北海虽赊，扶摇可接。

搬运于2025-08-24 22:39:05，当前版本为作者最后更新于2022-07-17 19:21:13，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

题目链接

今天比赛快结束时，听说有人阿克了这场，就进来看了看 F，这不是典题吗，，

首先这个模数是个合数，可以分解为 $24000001\times 30000001$，分别设为 $p_1$ 和 $p_2$。那么只要分别对 $p_1$ 和 $p_2$ 求解后用 crt 合并即可。

可以发现给出的 $n$ 除了 $13$ 都可以整除 $p_1-1$ 和 $p_2-1$，这些情况下 $\omega_n$ 在模 $p_1$ 和 $p_2$ 下都存在，直接用

$$\frac 1n\sum_{i=0}^{n-1}(\omega_n^i+\omega_n^{-i})^m \ \ \texttt{(1)} $$

计算，时间复杂度是 $\Theta(n \log p)$ 的。

还有另外一种做法是进行二项式展开，单位根反演即得

$$\sum_{k=0}^m \binom mk \frac 1n\sum_{i=0}^{n-1}\omega_n^{(2k-m)i}=\sum_{k=0}^m\binom mk[n|(2k-m)] \ \ \texttt{(2)} $$

用这种方法计算一组数据的时间复杂度是 $\Theta(p/n)$，需要 $\Theta(p)$ 预处理阶乘。

综合一下这两种算法，解方程 $n \log p = p/n$ 可以得到 $n=1000$ 时为临界数据，小于 $1000$ 的用方法 $\texttt{(1)}$，大于 $1000$ 的用方法 $\texttt{(2)}$ 即可。

对于 $n=13$ 这个特别数据可以发现，在模 $p=30000001$ 意义下有：

$$(x+x^{12})^p \equiv x^5+x^8 \pmod{(x^n-1)}$$ $$(x^5+x^8)^p\equiv 2+x^2+x^{11}\equiv (x+x^{12})^2 \pmod{(x^n-1)} $$

故直接将 $m$ 对 $p^2-1$ 取模即可。特别地，对于 $p=24000001$ 时循环节长度只有 $p-1$。

最后就是常数优化，主要有五点：

1. 对于题目中限定的 $n \neq 13$ 的数据，$m$ 为奇数时答案为零。

2. 由单位根的性质，若 $m$ 为偶数，在利用方法 $\texttt{(1)}$ 时 $i$ 只用取到一半，之后取二倍即可。

3. 对于 $n=1000$ 的情况，若使用方法 $\texttt{(1)}$，可以对快速幂做优化，做预处理。这里我选用 $8$ 进制快速幂。

4. 在使用方法 $\texttt{(2)}$ 时，可以发现答案也是对称的，也可以只算一半。

5. 使用无符号整数运算，效率会明显提升。

参考代码：

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#define ll unsigned long long
#define int unsigned int
#define N 30000003
#define P 720000054000001
#define i128 __int128_t
using namespace std;

const int pr[2] = {24000001,30000001};
const int g[2] = {23,14};
int ifac[N],sw[503][9][8];
int p,inv1k;

inline int power(int a,int t,int m){
    int res = 1;
    while(t){
        if(t&1) res = (ll)res*a%m;
        a = (ll)a*a%m;
        t >>= 1;
    }
    return res;
}

void init(int id){
    inv1k = power(1000,p-2,p);
    ifac[0] = 1,ifac[p-1] = p-1;
    for(int i=p-2;i;--i) ifac[i] = (ll)ifac[i+1]*(i+1)%p;
    int r = power(g[id],(p-1)/1000,p),w = 1,iw = 1;
    int ir = power(r,p-2,p);
    for(int i=0;i<500;++i){ 
        sw[i][0][0] = 1,sw[i][0][1] = w+iw;
        for(int k=2;k<8;++k) sw[i][0][k] = (ll)sw[i][0][k-1]*sw[i][0][1]%p;
        for(int j=1;j<9;++j){
            sw[i][j][0] = 1,sw[i][j][1] =  (ll)sw[i][j-1][7]*sw[i][j-1][1]%p;
            for(int k=2;k<8;++k) sw[i][j][k] = (ll)sw[i][j][k-1]*sw[i][j][1]%p;
        }
        w = (ll)w*r%p,iw = (ll)iw*ir%p;
    }
}

inline void multiply(const int *f,const int *g,int n,int *r){
    static int h[28];
    memset(h,0,sizeof(h));
    for(int i=0;i<n;++i)
    for(int j=0;j<n;++j)
        h[i+j] = (h[i+j]+(ll)f[i]*g[j])%p;
    for(int i=0;i<n;++i) r[i] = (h[i]+h[i+n])%p;
}

inline int solve(int n,ll m,int rt){
    int pw,t;
    ll res = 0;
    if(n==13){
        static int f[28],g[28];
        memset(f,0,sizeof(f));
        memset(g,0,sizeof(g));
        f[1] = f[n-1] = g[0] = 1;
        while(1){
            if(m&1) multiply(f,g,n,g);
            m >>= 1;
            if(m==0) break;
            multiply(f,f,n,f);
        }
        return g[0];
    }
    if(m&1) return 0;
    if(n==1000){
        for(int i=0;i<500;++i){
            pw = 1,t = m;
            for(int j=0;t;++j){
                if(t&7) pw = (ll)pw*sw[i][j][t&7]%p;
                t >>= 3;
            }
            res += pw;
        }
        return (res<<1)%p*inv1k%p;
    }
    if(n<=100){
        int r = power(rt,(p-1)/n,p);
        int ir = power(r,p-2,p),w = 1,iw = 1;
        n >>= 1;
        for(int i=0;i<n;++i){
            res += power(w+iw,m,p);
            w = (ll)w*r%p,iw = (ll)iw*ir%p;
        }
        return (res<<1)%p*power(n<<1,p-2,p)%p;
    }
    int st = (m>>1)%(n>>1),len = n>>1;
    for(int k=st;(k<<1)<m;k+=len) res += (ll)ifac[k]*ifac[m-k]%p;
    return ((res<<1)+(ll)ifac[m>>1]*ifac[m>>1])%p*power(ifac[m],p-2,p)%p;
}

inline ll crt(int x,int y){
    return ((i128)x*pr[1]*power(pr[1],pr[0]-2,pr[0])+(i128)y*pr[0]*power(pr[0],pr[1]-2,pr[1]))%P;
}

struct query{
    int n,m0,m1,x,y;
    ll m13;
    inline query(int _n=0,int _m0=0,int _m1=0,int _x=0,int _y=0):n(_n),m0(_m0),m1(_m1),x(_x),y(_y){}
}a[2503];

char str[10003];
int m0,m1,x,y,l,cnt;
ll m13;

signed main(){
    int n;
    p = pr[0];
    init(0);
    while(scanf("%d%s",&n,str)==2){
        l = strlen(str);
        m1 = m0 = m13 = 0;
        for(int i=0;i<l;++i) m0 = (m0*10+str[i]-'0')%(pr[0]-1);
        if(n==13){
            for(int i=0;i<l;++i) 
                m13 = (m13*10+str[i]-'0')%((ll)pr[1]*pr[1]-1);
        }else{
            for(int i=0;i<l;++i)
                m1 = (m1*10+str[i]-'0')%(pr[1]-1);
        }
        x = solve(n,m0,g[0]);
        a[++cnt] = (query(n,m0,m1,x));
        a[cnt].m13 = m13;
    }
    p = pr[1];
    init(1);
    for(int i=1;i<=cnt;++i){
        x = a[i].x;
        y = solve(a[i].n,a[i].n==13?a[i].m13:a[i].m1,g[1]);
        printf("%lld\n",crt(x,y));
    }
    return 0;
} 

后记：我没看比赛奖评，但这个出题人咋跑那么快啊？更神必的是代码量达到了惊人的 $29.49\text {KB}$，不知道是牛逼科技还是打表优化常数。