[氧化镍]七星连珠

最新推荐文章于 2023-11-24 22:27:48 发布

OneInDark

最新推荐文章于 2023-11-24 22:27:48 发布

阅读量159

点赞数

分类专栏： C++ 数学

本文链接：https://blog.csdn.net/qq_42101694/article/details/118853720

版权

C++ 同时被 2 个专栏收录

547 篇文章 3 订阅

订阅专栏

数学

280 篇文章 1 订阅

订阅专栏

题目

数据范围与提示
$N\le 50$ 且 $k\in\{2,3\}$ 且 $a_i<k^7$ 。

思路

首先这个每一行、每一列选一个，很让人联想到行列式。但是行列式是选出的元素作 “乘法”（也可以是自定义的乘法），这里是 $\rm xor$ ，怎么办？

当然是使用 $\rm fwt$ 啦！矩阵每个元素都是一个 $\rm fwt$ 数组，就可以直接乘起来了。恰好加法也支持。可是行列式里有一个 $(- 1)$ 的系数——可能某个值可以用两种方式得到，一个的系数是 $1$ ，另一个的系数是 $- 1$ ，相加就消失了！

一个骚操作是，给每个元素一个 随机权值，那么刚好抵消为 $0$ 的概率就很小。

那么行列式又怎么求？对 $\rm fwt$ 数组的每一位分开考虑。于是这就是一个高斯消元。

时间复杂度 $\mathcal O(n^4+n^3\cdot V\log V)$ 。至于 $k = 3$ ，显然是为了增加代码难度而进行的杂糅，相信大家都会 $3$ 进制 $\rm fwt$ 吧。或者你可以点开这个东西的最后一节，里面有一些说法，看看和你的想法是不是一样的。

代码

由于是 $\rm Linux$ 环境测评， $r a n d ()$ 的范围高达 $2^{31}-1$ 。

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>
#include <cstdlib>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long int_;
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int Mod = 1e9+9;
const int gG = 13; // primitive
const int w1 = 115381398;
const int w2 = 884618610;
const int inv2 = (Mod+1)>>1;
const int inv3 = (Mod<<1|1)/3;
inline int qkpow(int_ b,int q){
	int a = 1;
	for(; q; q>>=1,b=b*b%Mod)
		if(q&1) a = a*b%Mod;
	return a;
}

int k; // stupid base
int bit[8]; // 3^x
void FWT(int a[],int n,int opt){
	if(k == 2){
		for(int w=1; w<(1<<n); w<<=1)
		for(int *p=a; p!=a+(1<<n); p+=(w<<1))
		for(int i=0,t; i<w; ++i){
			t = p[i], p[i] = (t+p[i+w])%Mod;
			p[i+w] = (t+Mod-p[i+w])%Mod;
			if(opt == 0){
				p[i] = 1ll*p[i]*inv2%Mod;
				p[i+w] = 1ll*p[i+w]*inv2%Mod;
			}
		}
	}
	else if(k == 3){
		for(int w=1; w<=n; ++w)
		for(int *p=a; p!=a+bit[n]; p+=bit[w])
		for(int i=0,x,y; i<bit[w-1]; ++i){
			int j = i+bit[w-1], k = j+bit[w-1];
			x = p[i]; p[i] = ((x+p[j])%Mod+p[k])%Mod;
			y = p[j]; p[j] = (x+1ll*y*w1+1ll*p[k]*w2)%Mod;
			p[k] = (x+1ll*y*w2+1ll*p[k]*w1)%Mod;
			if(opt == 0){
				swap(p[j],p[k]);
				p[i] = 1ll*p[i]*inv3%Mod;
				p[j] = 1ll*p[j]*inv3%Mod;
				p[k] = 1ll*p[k]*inv3%Mod;
			}
		}
	}
	else puts("You bastard!");
}

const int MaxN = 52;
const int MaxV = 2187;
int a[MaxN][MaxN][MaxV];
int tmp[MaxV], b[MaxN][MaxN];

int gauss(int n){
	int res = 1, f = 1;
	for(int i=1,id=0; i<=n; ++i){
		for(int j=i; j<=n; ++j)
			if(b[j][i]) id = j;
		if(!b[id][i]) return 0;
		if(id != i) res = -res,
			swap(b[id],b[i]);
		for(int j=i+1; j<=n; ++j){
			if(!b[j][i]) continue;
			rep(k,i+1,n) // designed for i
				b[j][k] = (1ll*b[j][k]*b[i][i]%Mod
					+Mod-1ll*b[i][k]*b[j][i]%Mod)%Mod;
			b[j][i] = 0; // of course
			f = 1ll*f*b[i][i]%Mod; // inversion
		}
		res = 1ll*res*b[i][i]%Mod;
	}
	res = 1ll*res*qkpow(f,Mod-2)%Mod;
	return (res+Mod)%Mod;
}

int main(){
	readint(); srand(5201314);
	int n = readint(); k = readint();
	rep(i,bit[0]=1,7)
		bit[i] = 3ll*bit[i-1]%Mod;
	rep(i,1,n) rep(j,1,n)
		b[i][j] = readint();
	rep(i,1,n) rep(j,1,n){
		a[i][j][b[i][j]] = rand()%Mod;
		FWT(a[i][j],7,1); // transform
	}
	drep(now,((k==2)?(1<<7):bit[7])-1,0){
		rep(i,1,n) rep(j,1,n)
			b[i][j] = a[i][j][now];
		tmp[now] = gauss(n);
	}
	FWT(tmp,7,0); // get back
	rep(i,0,MaxV-1) if(tmp[i])
		printf("%d ",i);
	return 0;
}