uoj #86 mx的组合数 FFT 原根

最新推荐文章于 2022-08-08 11:07:35 发布

原创最新推荐文章于 2022-08-08 11:07:35 发布 · 750 阅读

0 ·

CC 4.0 BY-SA版权

FFT 同时被 2 个专栏收录

5 篇文章

订阅专栏

原根

1 篇文章

订阅专栏

本文介绍了一种利用快速傅里叶变换解决模数运算问题的方法，通过将乘法转换为加法卷积的形式，实现了高效的计算。文章详细讲解了如何使用Lucas定理、快速幂算法及离散对数等技术来处理特定范围内的组合数计算。

把l,r的条件看成前缀相减，那么就是求 $C_x^n\%p(1\le x\le bound)$
设n的p进制形式为 $a_1a_2a_3....a_{cnt}$
bound的p进制形式为 $b_1b_2b_3.....b_{cnt}$
由lucas定理原式为 $\prod C_{x_i}^{a_i}$ 。
设 $g[i][j]$ 表示由低到高前i位的xi任意取值且乘积为j的方案数
$f[i][j]$ 表示前i位组成的数小于n的前i位组成的数且乘积为j的方案数。

转移时枚举 $x=[0,p-1]$ ，设 $z=C_x^{a_{i+1}}$
$g[i+1][jz\%p]+=g[i][j]$
$f[i+1][jz\%p]+=g[i][j](x<b_{i+1})$
$f[i+1][jz\%p]+=f[i][j](x=b_{i+1})$
设G为p的原根，那么可以把x转成 $G^{ind[x]}$ 的形式。
这样就可以把乘法转成加法卷积的形式。
注意ind[0]是不存在的，因此0需要在计算完其他后用总的减掉。

#include <bits/stdc++.h>
using namespace std;
#define ll long long 
#define mod 998244353
#define P 31000
#ifdef ONLINE_JUDGE
typedef __int128_t sll;
#endif
#ifndef ONLINE_JUDGE
typedef long long sll;
#endif
int p,pg,cnt,cnt1,len;
sll n,l,r;
int ind[P],f[110][P],g[110][P],a1[110],b1[110],pw[P];
int jc[P],njc[P],C[110][P],ans[P];
int a[65537],b[65537],c[65537];
int qpow(int x,int y,int MOD)
{
    int ret=1;
    while(y)
    {
        if(y&1)ret=(ll)ret*x%MOD;
        x=(ll)x*x%MOD;y>>=1;
    }
    return ret;
}
int check(int x)
{
    for(int i=1,now=x;i<=p-2;i++)
    {
        if(now==1)return 0;
        now=now*x%p;
    }
    return 1;
}
void init()
{
    for(pg=1;;pg++)
        if(check(pg))break;
    pw[0]=1;
    for(int i=1;i<p;i++)
        pw[i]=pw[i-1]*pg%p;
    for(int i=1,now=pg;i<=p-1;i++)
        ind[now]=i,now=now*pg%p;
    while(n)a1[++cnt]=n%p,n/=p;

    jc[0]=njc[0]=1;
    for(int i=1;i<=p-1;i++)
        jc[i]=jc[i-1]*i%p;
    njc[p-1]=qpow(jc[p-1],p-2,p);
    for(int i=p-2;i>=1;i--)
        njc[i]=njc[i+1]*(i+1)%p;

    for(int i=1;i<=100;i++)
        for(int j=a1[i];j<p;j++)
            C[i][j]=jc[j]*njc[j-a1[i]]%p*njc[a1[i]]%p;
}
void NTT(int *a,int len,int type)
{
    for(int i=0,t=0;i<len;i++)
    {
        if(i>t)swap(a[i],a[t]);
        for(int j=len>>1;(t^=j)<j;j>>=1);
    }
    for(int i=2;i<=len;i<<=1)
    {
        int wn=qpow(3,(mod-1)/i,mod);
        for(int j=0;j<len;j+=i)
        {
            int w=1,t;
            for(int k=0;k<i>>1;k++,w=(ll)w*wn%mod)
            {
                t=(ll)a[j+k+(i>>1)]*w%mod;
                a[j+k+(i>>1)]=(a[j+k]-t+mod)%mod;
                a[j+k]=(a[j+k]+t)%mod;
            }
        }
    }
    if(type==-1)
    {
        for(int i=1;i<len>>1;i++)swap(a[i],a[len-i]);
        int t=qpow(len,mod-2,mod);
        for(int i=0;i<len;i++)a[i]=(ll)a[i]*t%mod;
    }
}
void calc(sll x)
{
    cnt1=0;memset(b1,0,sizeof(b1));
    while(x)b1[++cnt1]=x%p,x/=p;
    cnt1=max(cnt,cnt1);
    for(len=1;len<p<<1;len<<=1);
    g[0][1]=1;f[0][1]=1;

    for(int i=1;i<=cnt1;i++)
    {
        memset(g[i],0,sizeof(g[i]));
        memset(f[i],0,sizeof(f[i]));
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        for(int j=a1[i];j<p;j++)
            if(C[i][j])a[ind[C[i][j]]]++;
        for(int j=1;j<p;j++)
            (b[ind[j]]+=g[i-1][j])%=mod;
        NTT(a,len,1);NTT(b,len,1);
        for(int j=0;j<len;j++)
            c[j]=(ll)a[j]*b[j]%mod;
        NTT(c,len,-1);
        for(int j=0;j<len;j++)
            (g[i][pw[j%(p-1)]]+=c[j])%=mod;
        memset(a,0,sizeof(a));
        for(int j=a1[i];j<b1[i];j++)
            if(C[i][j])a[ind[C[i][j]]]++;
        NTT(a,len,1);
        for(int j=0;j<len;j++)
            c[j]=(ll)a[j]*b[j]%mod;
        NTT(c,len,-1);
        for(int j=0;j<len;j++)
            (f[i][pw[j%(p-1)]]+=c[j])%=mod;
        if(C[i][b1[i]])for(int j=1;j<p;j++)
            (f[i][C[i][b1[i]]*j%p]+=f[i-1][j])%=mod;
    }
}
template<typename T> inline T &Read(T &x)
{
    static char c;
    while (!isdigit(c = getchar()));
    x = c - '0';
    while (isdigit(c = getchar())) (x *= 10) += c - '0';
    return x;
}
int main()
{
    //freopen("tt.in","r",stdin);
    Read(p);Read(n);Read(l);Read(r);
    ans[0]=(r-l+1)%mod;
    init();
    calc(r);
    for(int i=1;i<p;i++)
        ans[i]=f[cnt1][i];
    calc(l-1);
    for(int i=1;i<p;i++)
        ans[i]=(ans[i]-f[cnt1][i]+mod)%mod;
    for(int i=1;i<p;i++)
        ans[0]=(ans[0]-ans[i]+mod)%mod;
    for(int i=0;i<p;i++)
        printf("%d\n",ans[i]);
    return 0;
}