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}\) 筛选次数

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)\) ，似乎已经比较小了，但是常数非常大，应该难以通过

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总览:

#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);
    }
}