【2018北京集训（六）】Lcm

时间 2019-12-05

标签 2018北京集训（六） lcm 繁體版

原文原文链接

Portal --> 出错啦qwq（好吧实际上是没有）ios

Description

　　给定两个正整数\(n,k\)，选择一些互不相同的正整数，知足这些数的最小公倍数刚好为\(n\)，而且这些数的和为\(k\)的倍数数组

　　求选择的方案数对\(232792561\)取模spa

　　数据范围：多组数据，组数\(T<=10,n<=10^{18},k<=20\)，且\(n\)的全部质因子不大于\(100\)code

Solution

　　这题。。好神仙啊qwq敲爆脑子都想不出来系列qwqblog

　　注意到\(n<=10^{18}\)意味着\(n\)至多有\(15\)个不一样的质因子（前\(15\)个质数乘一下就知道了），而且\(k\)的范围也是比较小的ip

　　而后还有一点就是这个模数比较有趣，它\(=lcm(1,2,...,20)+1\)，也就是说它能够对长度在\(1..20\)的区间DFTstring

　　那无论别的咱们能够先考虑一个最朴素的dpit

　　记\(f[i][j]\)表示当前所选的数状态为\(i\)，和\(\%k=j\)的方案数，其中这个“所选数的状态”具体计算方式是：假设写成分解质因数的形式后\(n=\sum p_i^{mi_i}\)，其中\(p_i\)为\(n\)的每个不一样的质因子，咱们的状态是一个二进制数，而且二进制的每一位对应一个\(n\)的质因子，对于二进制的第\(i\)位（假设对应\(p_i\)），若是说当前所选的数中存在一个数\(x\)知足\(p_i^{mi_i}|x\)，那么这位为\(1\)，不然为\(0\)io

　　这样一来，咱们最后的答案就应该是\(f[(1<<Cnt)-1][0]\)，其中\(Cnt\)为\(n\)的不一样质因子个数class

　　至于求解。。分解质因数和预处理因数能够暴力求解，可是后面的dp直接暴力转移的话稳稳的T啊（光枚举子集就要\(3^{Cnt}\)了还要再乘个\(k\))，因此如今的问题是咱们要怎么比较快速地转移

　　观察这个dp两维的转移，咱们会发现第一维的转移能够写成一个集合或卷积的形式（或起来等于某个数就把值累计进去），而第二维的转移则直接枚举因数什么的大力转移就行了

　　而后这里有一个颇有意思的想法：首先咱们发现这两维若是能够分开处理就会比较好一点，而后注意到这个模数颇有趣，容许咱们进行区间长度\(\in[1,20]\)的DFT，而DFT操做以后，dp的这两维在某种意义上就是独立的了，换句话来讲，若是说咱们先在同行内转移\(f\)数组，再对转之后每行的不彻底计算的\(f\)数组作一次DFT，这行的\(f\)的第二维就不会互相影响了，也就是说咱们能够经过这种方式实现第一维和第二维转移的分离

　　这样咱们就能够先处理\(f[][i]\)，直观一点来讲就是先在同行转移\(f\)数组，具体一点就是枚举\(n\)的因数，而后考虑加进去的贡献，也就是假设当前枚举到的因数是\(a[i]\)，\(a[i]\)对应的状态是\(st\)，那么咱们能够对于全部的\(j\in [0,k)\)转移：
\[ f[st][(j+a[i]\%k)]+=f'[st][j] \]
　　之因此在后面的\(f\)打了个\('\)是由于咱们这里要累加进去的是用\(a[i]\)更新前的\(f[st]\)的版本

　　这样咱们就获得了整合了同行数据以后的\(f\)数组（为了防止弄混在接下来的描述中咱们仍是将这个不彻底转移的\(f\)数组记为\(f1\)好了），而后对于每一行DFT一下，接下来就是考虑第一维的转移了，更加直观一点就是。。同列的\(f\)进行转移

　　接下来为了让描述变得更加简洁，咱们将第二维省去（由于反正转移的时候都是同列转移）

　　如今咱们要作的事情就是挑出若干个\(f1[i]\)，知足\(i_1\ or\ i_2\ or\ ...\ i_m=st\)而后将这堆\(f1[i]\)的值累加到\(f[st]\)去

　　咱们记\(F(S)=\sum\limits_{i\subseteq S}f[i]\)，若是说咱们知道了\(F(S)\)的取值，那么只要大力容斥一下就能够得出\(f\)数组了，具体一点的话就是：
\[ f[T]=\sum\limits_{S\subseteq T}(-1)^{|T|-|S|}F(S) \]
　　咱们将\(F(S)\)进行一下转化，会发现：
\[ \begin{aligned} F(S)&=\sum\limits_{j\subseteq S}f[j]\\ &=\prod\limits_{i\subseteq S}(1+f1[i]) \end{aligned} \]
　　具体的话就是由于对于\(f\)来讲，每个\(f[j]\)都是由若干个\(f1[i]\)组成的，咱们考虑将第二个等号后面的连乘的括号拆掉，展开以后会发现囊括了\(i\)的全部组合方式（为了方便理解能够本身用比较小的规模模拟一下），而后由于咱们限制了\(i\subseteq S\)因此能够保证任意的组合方式都是知足组合后\(\subseteq S\)的

　　而后因为\(f1\)是已知的，因此咱们如今就要想如何快速求\(F(S)\)

　　显然枚举子集不现实

　　若是咱们直接从小到大枚举状态，每次将这个状态对应的\(f1\)值转移的话，可能会出现重复计算的状况（好比说\(010\)能够转移到\(011\)和\(110\)，而\(011\)和\(110\)又均可以转移到\(111\)，这样若是按照这种转移方式的话，\(111\)这里\(010\)的\(f1\)值就会被重复计算），因此这里考虑一种很玄学的按顺序转移方式：

　　咱们考虑二进制一位一位转移，从低位枚举到高位，每枚举一位就把全部的这位为\(0\)的状态的值累加到将这位修改成\(1\)后的状态里面去，具体一点的话就是：

　　咱们以\(Cnt=3\)为例子，状态转移大概是这样：

　　其中橙黄色的箭头表示转移的方向，而后位数是从最右边开始数的

　　这样的话显然全部的状态都能转移到它能转移到的位置，如今咱们来讲明一下为何这样不会重复计算：首先能够确认的一点是，这种重复计算的状况只会出如今转移到一个被修改了两个及以上位的状态的时候，而咱们这样的枚举方式每次只转移到修改了一位的地方，而且是按照数位从低到高，状态从小到大的顺序转移，好比说刚才提到的\(010\)在\(111\)中被重复计算的状况，放在这个转移方式中就应该是：在转移第一位的时候，\(010\)会先转移到\(011\)，而后在转移第三位的时候再由累加了\(010\)的\(011\)转移到\(111\)，每次都是从只修改一位的地方转移过来，不会出现重复计算的状况~~（但其实上面的证实也不算。。特别严谨。。仍是感性理解既视感qwq）~~

　　而后咱们就能够在\(O(2^{Cnt})\)的时间内求得\(F(S)\)，接下来直接大力容斥一下就能够得出\(f[(1<<Cnt)-1][0]\)的值最后IDFT一下就能够得出答案啦

　　代码大概长这个样子

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
const int MOD=232792561,G=71,N=110;
const int p[26]={0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
ll mi[N],rec_mx[N],Tmp[N],have[N];
ll a[500010];
int f[(1<<15)+10][20];
ll n,K,T,Cnt,all;
void add(int &x,int y){x=(1LL*x+y+MOD)%MOD;}
bool in(int st,int x){return st>>(x-1)&1;}
int St(int x){return 1<<x-1;}
int ksm(int x,int y){
    int ret=1,base=x;
    for (;y;y>>=1,base=1LL*base*base%MOD)
        if (y&1) ret=1LL*ret*base%MOD;
    return ret;
}
namespace DFT{/*{{{*/
    int tmp[110],W[1010][2];
    int inv;
    void init(){
        int rt=ksm(G,(MOD-1)/K);
        W[0][0]=1;
        for (int i=1;i<=1000;++i)
            W[i][0]=1LL*W[i-1][0]*rt%MOD;
        for (int i=0;i<=1000;++i)
            W[i][1]=ksm(W[i][0],MOD-2);
    }
    void dft(int *a,int n,int op){
        for (int i=0;i<n;++i) tmp[i]=0;
        for (int i=0;i<n;++i) 
            for (int j=0;j<n;++j)
                tmp[i]=(1LL*tmp[i]+1LL*a[j]*W[i*j][op==-1]%MOD)%MOD;
        if (op==-1){
            inv=ksm(n,MOD-2);
            for (int i=0;i<n;++i) tmp[i]=1LL*tmp[i]*inv%MOD;
        }
        for (int i=0;i<n;++i) a[i]=tmp[i];
    }
}/*}}}*/
void Div(ll x){
    for (int i=1;i<=25;++i)
        if (x%p[i]==0){
            ++Cnt; rec_mx[Cnt]=1; mi[Cnt]=0; have[Cnt]=p[i];
            while (x%p[i]==0) 
                x/=p[i],++mi[Cnt],rec_mx[Cnt]*=p[i];
        }
}
void dfs(int now,ll prod){
    if (now>Cnt){
        a[++a[0]]=prod;return;
    }
    ll tmp=prod;
    for (int i=0;i<=mi[now];++i){
        dfs(now+1,tmp);
        tmp*=have[now];
    }
}
void prework(){
    a[0]=0; Cnt=0;
    Div(n);
    dfs(1,1);
}
void update(int st,int r){
    for (int i=0;i<K;++i) Tmp[i]=f[st][i];
    for (int i=0;i<K;++i) add(f[st][(i+r)%K],Tmp[i]);
}
void solve(){
    int st,op;
    memset(f,0,sizeof(f));
    all=(1<<Cnt)-1;
    for (int i=0;i<=all;++i) f[i][0]=1;//init
    for (int i=1;i<=a[0];++i){
        st=0;
        for (int j=1;j<=Cnt;++j)
            if (a[i]%rec_mx[j]==0) st|=St(j);
        update(st,a[i]%K);
    }
    for (int i=0;i<=all;++i)
        DFT::dft(f[i],K,1);

    for (int i=1;i<=Cnt;++i)
        for (int j=0;j<=all;++j)
            if (in(j,i))
                for (int k=0;k<K;++k)
                    f[j][k]=1LL*f[j][k]*f[j^St(i)][k]%MOD;
    for (int i=0;i<all;++i){
        op=1;
        for (int j=1;j<=Cnt;++j)
            if (!in(i,j)) op*=-1;
        for (int j=0;j<K;++j)
            add(f[all][j],op*f[i][j]);
    }
    DFT::dft(f[all],K,-1);
    printf("%d\n",f[all][0]);
}

int main(){
#ifndef ONLINE_JUDGE
    freopen("a.in","r",stdin);
#endif
    scanf("%d",&T);
    for (int o=1;o<=T;++o){
        scanf("%lld%lld",&n,&K);
        DFT::init();
        prework();
        solve();
    }
}