• 【2018北京集训6】Lcm DFT&FWT


    首先我们来看下此题的模数232792561。

    232792561=lcm(1,2,3.......20)+1。这个性质将在求值时用到。

    我们将n分解质因数,令$m$为$n$的素因子个数,设n=$Pi_{j=0}^{m-1} p_j^{b_j}$ ,其中$p_j$是素数且$p_0$至$p_{m-1}$从小到大排列。考虑到$n≤10^{18}$,则$m≤15$。

    我们用 $f[i][j]$ 表示当前$n$的因数$x$所表示的状态为$i$,且模$k$为$j$时的方案数。

    下面讲下如何用一个已知的因数$x$求出$i$。

    设$x=Pi_{j=0}^{m-1} p_j^{d_j}$ ,则$i=sum _{j=0}^{m-1} 2^{j} imes [d_j==b_j]$ 。(此处举个例子,令$n=12$,则当$x=6$时,$i=2$,当$x=12$时,$i=3$,如果还是不清楚的话可以去看我的代码)

    据此定义,则答案显然为$f[2^{m}-1][0]$。

    初始化:对于$forall n \% x==0$ 有 $f[i][x \% k]=1$

    下面考虑如何转移。

    显然,朴素的转移方法是会超时的(说白了就是只有10分),那么考虑如何加速转移。

    首先,我们现在$f[i][0...2^{m}-1][]$中完成同层的转移(详情见代码)

    假设我们去取出数量相同的$c$和b,使得$Pi f[b_l][c_i]$ 可以对$f[2^m-1][0]$产生贡献。

    考虑可对$f[2^m-1][0]$产生贡献的条件,显然有三:

    $b_1 or b_2 or .....b_p=2^m-1 $

    $b_1 e b_2 e ..... e b_p$

    $c_1 + c_2 + ..... c_p equiv 0 (mod k)$

    对于这两条限制,直接对每一行(既$f[i]$)先做一次DFT以满足性质3,然后对每一列(即$f[][j]$)做一次FWT以满足性质1和性质2(没错通过这种嵌套即可同时满足两个性质),然后将最后一列做一次IDFT,即可求出$f[2^m-1][0]$

    此题看题5分钟,懵逼4小时,写题1小时.....

     

     1 #include<bits/stdc++.h>
     2 #define M 20000000
     3 #define MOD 232792561
     4 #define L long long
     5 using namespace std;
     6 L pow_mod(L x,int k){
     7     L ans=1;
     8     while(k){
     9         if(k&1) ans=ans*x%MOD;
    10         x=x*x%MOD; k>>=1;
    11     }
    12     return ans;
    13 }
    14 L f[1<<15][20]={0}; 
    15 L n,m,nn,e[20]={0},d[20]={0},p[20]={0},a[M]={0}; int k;
    16 
    17 void DFT(L a[],L w[],int n,int on){
    18     L b[20]; memset(b,0,sizeof(b)); 
    19     for(int i=0;i<n;i++)
    20     for(int j=0;j<n;j++)
    21     b[i]=(b[i]+a[j]*w[i*j])%MOD;
    22     memcpy(a,b,160);
    23     if(on==-1){
    24         L inv=pow_mod(n,MOD-2);
    25         for(int i=0;i<n;i++) a[i]=a[i]*inv%MOD;
    26     }
    27 }
    28 
    29 void fenjie(L n){//分解质因数 
    30     for(int i=2;i<=100;i++) 
    31     if(n%i==0){
    32         m++; p[m]=i; e[m]=1; 
    33         while(n%i==0) n/=i,d[m]++,e[m]*=i;
    34     }
    35 }
    36 void coushu(int x,L ji){//凑出所有约数 
    37     if(x>m) {a[++nn]=ji; return;}
    38     coushu(x+1,ji);
    39     for(int i=1;i<=d[x];i++){
    40         ji*=p[x];
    41         coushu(x+1,ji);
    42     }
    43 }
    44 
    45 void add(int x,int nk){
    46     L b[20]; 
    47     memcpy(b,f[x],160);
    48     for(int i=0;i<k;i++)
    49     f[x][i]=(f[x][i]+b[(i-nk+k)%k])%MOD;
    50     f[x][nk]=(f[x][nk]+1)%MOD;
    51 }    
    52 L w1[500]={0},w2[500]={0}; 
    53 void sub(){
    54     memset(e,0,sizeof(e)); memset(d,0,sizeof(d));
    55     memset(p,0,sizeof(p)); memset(a,0,sizeof(a));
    56     memset(f,0,sizeof(f)); m=nn=0;
    57     cin>>n>>k;
    58     L w=pow_mod(71,(MOD-1)/k); w1[0]=w2[0]=1; 
    59     for(int i=1;i<500;i++){
    60         w1[i]=w1[i-1]*w%MOD;
    61         w2[i]=pow_mod(w1[i],MOD-2);
    62     }
    63     m=0; fenjie(n);
    64     coushu(0,1);
    65     for(int i=1;i<=nn;i++){
    66         int id=0;
    67         for(int j=1;j<=m;j++)
    68         if(a[i]%e[j]==0) id=id|(1<<(j-1));
    69         add(id,a[i]%k);
    70     }
    71     int mm=1<<m;
    72     for(int i=0;i<mm;i++) DFT(f[i],w1,k,1);
    73     for(int i=0;i<mm;i++)
    74     for(int j=0;j<k;j++) f[i][j]++;
    75     
    76     for(int i=1;i<mm;i<<=1)
    77     for(int j=0;j<mm;j++)
    78     if(i&j){
    79         for(int l=0;l<k;l++)
    80         f[j][l]=f[j][l]*f[i^j][l]%MOD;
    81     }    
    82     for(int i=1;i<mm;i<<=1)
    83     for(int j=0;j<mm;j++)
    84     if(i&j){
    85         for(int l=0;l<k;l++)
    86         f[j][l]=(f[j][l]-f[i^j][l]+MOD)%MOD;
    87     }    
    88     DFT(f[mm-1],w2,k,-1);
    89     printf("%lld
    ",f[mm-1][0]);
    90 }
    91 
    92 int main(){
    93     int cas; scanf("%d",&cas);
    94     while(cas--) sub();
    95 }
  • 相关阅读:
    yii之behaviors
    查看windows系统信息
    idm chrome扩展被阻止解决办法
    音乐乐理基础
    bootstrap4
    七牛上传整合CI
    提升上传速度
    卡漫绘图
    指针的操作
    定语从句八个易混淆
  • 原文地址:https://www.cnblogs.com/xiefengze1/p/8603141.html
Copyright © 2020-2023  润新知