CodeChef November Challenge2019 Winning Ways (3-FWT)

显然每个把每个数换成其因子个数-1，就能转为一个扩展的$\text{Nim}$游戏

每次操作$1,2,\cdots,k$堆的$\text{Nim}$游戏，其判定方法是：

将每个数二进制分解，对于每个二进制位上分别计算1的个数$\mod (k+1)$，如果均为0则先手必败

对于这道题$k=3$，我们考虑将其转为二进制之后的形式累成3进制，然后就能进行3进制按位不进位加法，即类异或

然后问题实际上有非常多的部分需要考虑

Part1 如何求因子个数

一个简单的方法是枚举$[1,\sqrt n]$判断是否整除，复杂度过高

对于$n=\prod p_i^{c_i}$($p_i$为质数)，其因子个数为$\prod (c_i+1)$

由这个式子对于$n$进行质因数分解，枚举$[1,\sqrt n]$中的质数，复杂度为$O(\pi(\sqrt n))=O(\frac{\sqrt n} {\log n})$，这个应该够了？

然后是一个常规套路型的分解方法：

先对于$[1,n^{\frac{1} {3} }]$的质数筛$n$，剩余的部分只有3种情况

1.$n$被筛成1了

2.$n$被筛到只剩一个质数，可以用$\text{Miller_Rabin}$算法快速判断，可以参考

3.$n$仍然是若干质数的乘积，此时质因子必然$>n^{\frac{1} {3} }$，因此最多只有两个

那么只需要判断$n$是否是完全平方数即可

总复杂度为$O(w\cdot \log n+\frac{n^{\frac{1} {3} }} {\log n})$，其中$w$为$\text{Miller_Rabin}$筛选次数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
const int N=2e5+10;
int notpri[N],pri[N],pc;

ll qpow(ll x,ll k=P-2,int P=::P) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}
int Trans(int x){ 
    static int buf[20],l=0;
    l=0;
    while(x) buf[++l]=x%2,x/=2;
    int s=0;
    drep(i,l,1) s=s*3+buf[i];
    cmax(ma,s);
    return s;
}

int Miller_Rabin(int x) {
    if(x<N) return !notpri[x];
    if(x==2) return 1;
    if(x<=1 || ~x&1) return 0;
    ll s=0,t=x-1;
    while(~t&1) s++,t>>=1;
    rep(i,1,4) {
        ll a=pri[rand()%pc+1],b=qpow(a,t,x),c;
        rep(j,1,s) {
            c=1ll*b*b%x;
            if(c==1 && b!=1 && b!=x-1) return 0;
            b=c;
        }
        if(b!=1) return 0;
    }
    return 1;
}
int CheckSqr(int n){
    int y=round(sqrt(n));
    return y*y==n;
}

int Count(int n){
    int res=1;
    for(int i=1;i<=pc && pri[i]*pri[i]*pri[i]<=n;++i) if(n%pri[i]==0) {
        int c=0;
        while(n%pri[i]==0) n/=pri[i],c++;
        res*=c+1;
    }
    if(n==1) return res;
    if(CheckSqr(n)) return res*3;
    if(Miller_Rabin(n)) return res*2;
    return res*4;
}

void Init(){
    rep(i,2,N-1) if(!notpri[i]) {
        pri[++pc]=i;
        for(int j=i+i;j<N;j+=i) notpri[j]=1;
    }
}

Part2 快速计算答案

$10^9$以内的数，最大因子个数为$1334$，这个数为$931170240$

转成二进制之后最多包含$11$位，三进制下最大为$3^{11}-1=177146$，令这个上界为$M$

一种非常暴力的方法就是直接枚举，$NM$计算每次选择一个数，复杂度为$O(NMK)$，应该可以通过$N\leq 12$的数据

一个比较浅显的优化可以用快速幂维护乘法，复杂度为$O(M^2\log K)$

由于是3进制类异或，接下来考虑用$\text{3-FWT}$优化乘法，可以参考

模数为$P=10^9+7$，不存在整数$3$阶单位根，因此要用类似拆系数$\text{FFT}$方法做

复杂度为$O(M\log M\log K)$，似乎已经比较小了，但是常数非常大，应该难以通过

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
int R;// R为上界
struct Cp{
    db x,y;
    Cp(){ }
    Cp(db x,db y):x(x),y(y) { }
    Cp operator + (const Cp t){ return Cp(x+t.x,y+t.y); }
    Cp operator - (const Cp t){ return Cp(x-t.x,y-t.y); }
    Cp operator * (const Cp t){ return Cp(x*t.x-y*t.y,x*t.y+y*t.x); }
} A[N],B[N],C[N],D[N];
Cp w[30];

int Add(int x,int y) {
    static int A[20],B[20];
    rep(i,0,19) A[i]=x%3,x/=3;
    rep(i,0,19) B[i]=y%3,y/=3;
    int ans=0;
    drep(i,19,0) ans=ans*3+((A[i]+B[i])%3);
    return ans;
}

void FWT(Cp *a,int f) {
    for(int i=1;i<R;i*=3) {
        for(int l=0;l<R;l+=i*3) {
            for(int j=l;j<l+i;++j) {
                static Cp t[3];
                if(f==1) {
                    t[0]=a[j]+a[j+i]+a[j+i*2];
                    t[1]=a[j]+w[1]*a[j+i]+w[2]*a[j+i*2];
                    t[2]=a[j]+w[2]*a[j+i]+w[1]*a[j+i*2];
                } else {
                    t[0]=a[j]+a[j+i]+a[j+i*2];
                    t[1]=a[j]+w[2]*a[j+i]+w[1]*a[j+i*2];
                    t[2]=a[j]+w[1]*a[j+i]+w[2]*a[j+i*2];
                }
                a[j]=t[0],a[j+i]=t[1],a[j+i*2]=t[2];
                if(f==-1) {
                    rep(d,0,2) {
                        a[j+i*d].x/=3,a[j+i*d].y/=3;
                    }
                }
            }
        }
    }
}

const int S=(1<<15)-1;

#define FWTs

struct Poly{ 
    int a[N];
    Poly operator * (const Poly __) const {
        Poly res;
        // 拆系数，任意模数3-FWT
#ifdef FWTs
        rep(i,0,R-1) A[i]=Cp((a[i]&S),(a[i]>>15)),B[i]=Cp((a[i]&S),0);
        rep(i,0,R-1) C[i]=Cp((__.a[i]&S),(__.a[i]>>15));
        FWT(A,1),FWT(B,1),FWT(C,1);
#define E(x) ((ll)(x+0.5))%P
        rep(i,0,R-1) A[i]=A[i]*C[i];
        rep(i,0,R-1) B[i]=B[i]*C[i];
        FWT(A,-1),FWT(B,-1);
        rep(i,0,R-1) {
            ll a=E(B[i].x),b=E(A[i].y),c=E(B[i].x-A[i].x);
            res.a[i]=(a+1ll*b*(S+1)+1ll*c*(S+1)*(S+1))%P;
        }
#else
        rep(i,0,R-1) res.a[i]=0;
        rep(i,0,R-1) if(a[i]) rep(j,0,R-1) if(__.a[j]) {
            int k=Add(i,j);
            res.a[k]=(res.a[k]+1ll*a[i]*__.a[j])%P;
        }
#endif
        return res;
    }
} x,res;



ll qpow(ll x,ll k=P-2) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}

int main(){
    rep(i,2,N-1) if(!notpri[i]) {
        pri[++pc]=i;
        for(int j=i+i;j<N;j+=i) notpri[j]=1;
    }
    w[0]=Cp(1,0);
    rep(i,1,29) w[i]=w[i-1]*Cp(cos(2*Pi/3),sin(2*Pi/3));
    // 复平面3阶单位根

    n=rd(),m=rd();
    rep(i,1,n) x.a[Count(rd())]++;
    for(R=1;R<=ma;R*=3);

    res.a[0]++;
    for(;m;m>>=1) {
        if(m&1) res=res*x;
        x=x*x;
    }
    int ans=0;
    rep(i,1,R-1) ans+=res.a[i],Mod1(ans);
    printf("%d\n",ans);
}

Further

对于形式幂级数多项式，我们知道$K$次幂的循环卷积可以直接 $\text{DFT}$一次，每一位快速幂，然后$\text{IDFT}$

同理的，如果你学习了$\text{K-FWT}$就知道这就是一个按$K$进制位，每一位分别进行循环卷积，因此也可以用类似的方法做

但是遇到一个非常大的问题就是无法找到模意义下的$3$阶单位根(指$3\not |(P-1)$)

如果用复平面单位根$\omega_n=cos(\frac{2\pi} {n})+sin(\frac{2\pi} {n})\cdot i$($i=\sqrt {-1})$，无法在计算时保证值域精度

这里由于$n=3$比较特殊，发现$\omega_3=cos(\frac{2\pi} {3})+sin(\frac{2\pi} {3})\cdot i=-\frac{1} {2}+\frac{\sqrt 3} {2}\cdot i$

而$3$在$\mod 10^9+7$下存在二次剩余，因此可以用一个模意义下的复数描述复平面单位根

应该是有通行的单位根求法，会根据$n$不同要用更复杂的高维复数描述，但是我并不会.jpg

总复杂度为$O(M(\log M+\log K))$，分别为进行$\text{3-FWT}$以及快速幂的复杂度

Code总览:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef long double db;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
int rd(){
    int s=0,f=0;
    while(!isdigit(IO=getchar())) if(IO=='-') f=1;
    do s=(s<<1)+(s<<3)+(IO^'0');
    while(isdigit(IO=getchar()));
    return f?-s:s;
}
bool Mbe;
const int N=2e5,P=1e9+7;
const int Quad3=82062379; // 3在Mod P意义下的二次剩余
const db Pi=acos((db)-1);

const int MaxX=931170240;


int n,m,ma,R;
int notpri[N],pri[N],pc;

ll qpow(ll x,ll k=P-2,int P=::P) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}
int Trans(int x){ 
    static int buf[20],l=0;
    l=0;
    while(x) buf[++l]=x%2,x/=2;
    int s=0;
    drep(i,l,1) s=s*3+buf[i];
    cmax(ma,s);
    return s;
}

int Miller_Rabin(int x) {
    if(x<N) return !notpri[x];
    if(x==2) return 1;
    if(x<=1 || ~x&1) return 0;
    ll s=0,t=x-1;
    while(~t&1) s++,t>>=1;
    rep(i,1,4) {
        ll a=pri[rand()%pc+1],b=qpow(a,t,x),c;
        rep(j,1,s) {
            c=1ll*b*b%x;
            if(c==1 && b!=1 && b!=x-1) return 0;
            b=c;
        }
        if(b!=1) return 0;
    }
    return 1;
}
int CheckSqr(int n){
    int y=round(sqrt(n));
    return y*y==n;
}

int Count(int n){
    int res=1;
    for(int i=1;i<=pc && pri[i]*pri[i]*pri[i]<=n;++i) if(n%pri[i]==0) {
        int c=0;
        while(n%pri[i]==0) n/=pri[i],c++;
        res*=c+1;
    }
    if(n==1) return res;
    if(CheckSqr(n)) return res*3;
    if(Miller_Rabin(n)) return res*2;
    return res*4;
}

struct Cp{
    int x,y;
    Cp(){ }
    Cp(int x,int y):x(x),y(y) { }
    Cp operator + (const Cp t){ 
        Cp res(x+t.x,y+t.y);
        Mod1(res.x),Mod1(res.y);
        return res;
    }
    Cp operator * (const Cp t){ return Cp((1ll*x*t.x+1ll*(P-y)*t.y)%P,(1ll*x*t.y+1ll*y*t.x)%P); }
} A[N]; // 模意义下 模拟复平面单位根
Cp w1,w2; // 3阶单位根及其平方

Cp qpow(Cp x,ll k=P-2) {
    Cp res(1,0);
    for(;k;k>>=1,x=x*x) if(k&1) res=res*x;
    return res;
}

// 下面是展开的FWT式子
void FWT() {
    for(int i=1;i<R;i*=3) {
        for(int l=0;l<R;l+=i*3) {
            for(int j=l;j<l+i;++j) {
                Cp a=A[j]+A[j+i]+A[j+i*2];
                Cp b=A[j]+w1*A[j+i]+w2*A[j+i*2];
                Cp c=A[j]+w2*A[j+i]+w1*A[j+i*2];
                A[j]=a,A[j+i]=b,A[j+i*2]=c;
            }
        }
    }
}

void IFWT() {
    for(int i=1;i<R;i*=3) {
        for(int l=0;l<R;l+=i*3) {
            for(int j=l;j<l+i;++j) {
                Cp a=A[j]+A[j+i]+A[j+i*2];
                Cp b=A[j]+w2*A[j+i]+w1*A[j+i*2];
                Cp c=A[j]+w1*A[j+i]+w2*A[j+i*2];
                A[j]=a,A[j+i]=b,A[j+i*2]=c;
            }
        }
    }
    ll base=qpow(R);
    rep(i,0,R-1) A[i].x=A[i].x*base%P;
}

int main(){
    rep(i,2,N-1) if(!notpri[i]) {
        pri[++pc]=i;
        for(int j=i+i;j<N;j+=i) notpri[j]=1;
    }
    w1=Cp(P-(P+1)/2,1ll*Quad3*(P+1)/2%P);
    w2=w1*w1;

    rep(kase,1,rd()) {
        n=rd(),m=rd();
        rep(i,1,n) A[Trans(Count(rd())-1)].x++;
        for(R=1;R<=ma;R*=3);
        FWT();
        rep(i,0,R-1) A[i]=qpow(A[i],m);
        IFWT();
        int ans=0;
        rep(i,1,R-1) ans+=A[i].x,Mod1(ans);
        rep(i,0,R-1) A[i].x=A[i].y=0;
        printf("%d\n",ans);
    }
}