CodeChef November Challenge2019 Winning Ways (3-FWT)
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);
}
}
