希望复习高进制FWT的时候，能够快速回想起来。
FWT感觉就是每一维单独考虑，（虽然我不知道为什么这是对的）
分别对一个奇怪的东西做卷积，
那个奇怪的东西在k进制下就是关于k次单位根的范德蒙特矩阵。
范德蒙特矩阵的逆矩阵大致就是每行除了第一个数之外翻转一下，然后除以矩阵的阶。
也可理解为原矩阵把k次单位根取反，反正每一维上就是个DFT，按FFT的做法搞就行了。
然后这两道题毒瘤的地方就在于不能用实数来表示一个复数，否则会爆精。
而且你也不一定可以算出单位根的值，十分毒瘤。
于是你就需要搞事，
就是把一个数用向量表示，表示单位根的某些次幂前的系数。
单位根的性质 https://blog.csdn.net/DT_Kang/article/details/79944113
分圆多项式 https://yhx-12243.github.io/OI-transit/memos/17.html#pr-1-1
可以证明在模分圆多项式的意义下对答案没有影响，虽然我也不会证。
然后考虑清华集训这题。
当(k=3)时，分圆多项式为(x^2+x+1)
并且把(w^1_3)带进去由分圆多项式的定义可得(w^1_3+w^2_3+1=0)
然后你就只需要存两个系数了。
有一个优化是DFT时乘上的是单位根，所以可以不用乘法，直接移系数。
#pragma GCC optimize(3)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace io {
const int SIZE = 1.1e6;
char buff[SIZE];
char *l = buff, *r = buff;
void init() {
l = buff;
r = l + fread(buff, 1, SIZE, stdin);
}
char gc() {
if (l == r) init();
if(l==r)return EOF;
return *(l++);
}
void read(int &x) {
x = 0;
char ch = gc();
while(!isdigit(ch)) ch = gc();
while(isdigit(ch)) x = x * 10 + ch - '0', ch = gc();
}
}
using io::read;
namespace{
int M;
void FIX(int &x){
x>=M?x-=M:(x<0?x+=M:0);
}
int ADD(int x,int y){
return (x+=y)>=M?x-M:x;
}
int SUB(int x,int y){
return (x-=y)<0?x+M:x;
}
void ADDT(int &x,int y){
(x+=y)>=M?x-=M:0;
}
void SUBT(int &x,int y){
(x-=y)<0?x+=M:0;
}
int MUL(int x,int y){
return (ll)x*y%M;//(ll)x*y-(ll)x*y/M*M
}
int fp(int x,int y){
int ret=1;
for (; y; y>>=1,x=MUL(x,x))
if (y&1) ret=MUL(ret,x);
return ret;
}
void exgcd(int a,int b,int &x,int &y){
if (b==0) return (void)(x=1,y=0);
exgcd(b,a%b,y,x);
y-=a/b*x;
}
}
template<int N>
struct Number{
int a[N-1];
Number(){
memset(a,0,sizeof(a));
}
friend ostream & operator <<(ostream &out,const Number<N> &b){
out<<'{';
for (int i=0; i<N-2; ++i) out<<b.a[i]<<',';
return out<<b.a[N-2]<<'}';
}
Number operator +(const Number &b) const{
Number ret;
for (int i=0; i<N-1; ++i) ret.a[i]=ADD(a[i],b.a[i]);
return ret;
}
Number operator -(const Number &b) const{
Number ret;
for (int i=0; i<N-1; ++i) ret.a[i]=SUB(a[i],b.a[i]);
return ret;
}
Number operator *(const Number &b) const{
Number ret;
int tmp=0;
for (int i=0; i<N-1; ++i)
for (int j=0; j<N-1; ++j)
ADDT(i+j==N-1?tmp:ret.a[i+j>=N?i+j-N:i+j],MUL(a[i],b.a[j]));
for (int i=0; i<N-1; ++i) SUBT(ret.a[i],tmp);
return ret;
}
bool operator <(const Number &b) const{
for (int i=0; i<N-1; ++i) if (a[i]!=b.a[i]) return a[i]<b.a[i];
return 0;
}
Number operator <<(const int &c) const{
//less than N
Number ret;
int tmp=0;
for (int i=0; i<N-1; ++i) ADDT(i+c==N-1?tmp:ret.a[i+c>=N?i+c-N:i+c],a[i]);
for (int i=0; i<N-1; ++i) SUBT(ret.a[i],tmp);
return ret;
}
};
namespace std{
template<int N>
struct hash<Number<N> >{
int operator()(const Number<N> &a) const{
int ret=0;
for (int i=0; i<N-1; ++i) ret^=std::hash<int>{}(a.a[i]);
return ret;
}
};
template<int N>
struct equal_to<Number<N> >{
bool operator()(const Number<N> &a,const Number<N> &b) const{
for (int i=0; i<N-1; ++i) if (a.a[i]!=b.a[i]) return 0;
return 1;
}
};
}
unordered_map<Number<3>,Number<3> > mp;
template<int N>
Number<N> fp(Number<N> x,int y){
Number<N> ret;
ret.a[0]=1;
for (; y; y>>=1,x=x*x) if (y&1) ret=ret*x;
return ret;
}
template<int N,int LEN>
struct FWT_helper{
typedef Number<N> Num;
Num w[N][N];
int pwd[LEN+1];
FWT_helper(){
assert(N>1);
assert(M);
for (int i=0; i<N; ++i) w[i][0].a[0]=1;
w[0][1].a[0]=1;
w[1][1].a[1]=1;
for (int i=2; i<N; ++i) w[i][1]=w[i-1][1]*w[1][1];
for (int i=0; i<N; ++i)
for (int j=2; j<N; ++j)
w[i][j]=w[i][j-1]*w[i][1];
pwd[0]=1;
for (int i=1; i<=LEN; ++i) pwd[i]=pwd[i-1]*N;
}
void DFT(Num *a){
Num b[N];
for (int i=0; i<N; ++i){
//for (int j=0; j<N; ++j) b[i]=b[i]+a[j]*w[i][j];???
for (int j=0,d=0; j<N; ++j,(d+=i)>=N?d-=N:0)
b[i]=b[i]+(a[j]<<d);
}
for (int i=0; i<N; ++i) a[i]=b[i];
}
void IDFT(Num *a){
reverse(a+1,a+N);
DFT(a);
//Don't forget to div
}
void FWT(Num *a,int m,bool rev){
for (int i=0; i<m; ++i)
for (int cj=0; cj<pwd[m]; cj+=pwd[i+1])
for (int j=cj; j<cj+pwd[i]; ++j){
Num tmp[N];
for (int k=0; k<N; ++k) tmp[k]=a[j+pwd[i]*k];
if (rev) IDFT(tmp);
else DFT(tmp);
for (int k=0; k<N; ++k) a[j+pwd[i]*k]=tmp[k];
}
if (rev){
int inv,tmp;
exgcd(pwd[m],M,inv,tmp);
FIX(inv);
for (int i=0; i<pwd[m]; ++i) a[i].a[0]=MUL(a[i].a[0],inv);
}
}
};
void test(){
M=998244353;
FWT_helper<3,15> a;
Number<3> done[9];
for (int i=0; i<a.pwd[2]; ++i) done[i].a[0]=i;
a.FWT(done,2,0);
for (int i=0; i<a.pwd[2]; ++i) done[i]=done[i]*done[i];
a.FWT(done,2,1);
for (int i=0; i<a.pwd[2]; ++i) cerr<<done[i]<<" "; cerr<<endl;
}
int m,t;
int cnt1[531441],cnt2[531441];
Number<3> f[531441],g[531441];
int b[13][13];
int main(){
//test();
read(m); read(t); read(M);
FWT_helper<3,12> ups;
for (int i=0; i<ups.pwd[m]; ++i){
int x; read(x);
f[i].a[0]=x;
}
for (int i=0; i<=m; ++i)
for (int j=0; i+j<=m; ++j)
read(b[i][j]);
for (int i=0; i<ups.pwd[m]; ++i){
cnt1[i]=cnt1[i/3]+(i%3==1);
cnt2[i]=cnt2[i/3]+(i%3==2);
}
for (int i=0; i<ups.pwd[m]; ++i)
g[i].a[0]=b[cnt1[i]][cnt2[i]];
ups.FWT(f,m,0);
ups.FWT(g,m,0);
for (int i=0; i<ups.pwd[m]; ++i)
f[i]=f[i]*(mp.count(g[i])?mp[g[i]]:(mp[g[i]]=fp(g[i],t)));
ups.FWT(f,m,1);
for (int i=0; i<ups.pwd[m]; ++i) cout<<f[i].a[0]<<'
';
}
然后考虑cf题，朴素想法是直接模(x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1)，最后再模10的分圆多项式，这样也可以过。
#pragma GCC optimize(3)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
namespace io {
const int SIZE = 1.1e6;
char buff[SIZE];
char *l = buff, *r = buff;
void init() {
l = buff;
r = l + fread(buff, 1, SIZE, stdin);
}
char gc() {
if (l == r) init();
if(l==r)return EOF;
return *(l++);
}
void read(int &x) {
x = 0;
char ch = gc();
while(!isdigit(ch)) ch = gc();
while(isdigit(ch)) x = x * 10 + ch - '0', ch = gc();
}
}
using io::read;
namespace{
ull ADD(ull x,ull y){
return x+y;
}
void ADDT(ull &x,ull y){
x+=y;
}
ull SUB(ull x,ull y){
return x-y;
}
void SUBT(ull &x,ull y){
x-=y;
}
ull MUL(ull x,ull y){
return x*y;
}
}
template<int N>
struct Number{
ull a[N-1];
Number(){
memset(a,0,sizeof(a));
}
friend ostream & operator <<(ostream &out,const Number<N> &b){
out<<'{';
for (int i=0; i<N-2; ++i) out<<b.a[i]<<',';
return out<<b.a[N-2]<<'}';
}
Number operator +(const Number &b) const{
Number ret;
for (int i=0; i<N-1; ++i) ret.a[i]=ADD(a[i],b.a[i]);
return ret;
}
Number operator -(const Number &b) const{
Number ret;
for (int i=0; i<N-1; ++i) ret.a[i]=SUB(a[i],b.a[i]);
return ret;
}
Number operator *(const Number &b) const{
Number ret;
ull tmp=0;
for (int i=0; i<N-1; ++i)
for (int j=0; j<N-1; ++j)
ADDT(i+j==N-1?tmp:ret.a[i+j>=N?i+j-N:i+j],MUL(a[i],b.a[j]));
for (int i=0; i<N-1; ++i) SUBT(ret.a[i],tmp);
return ret;
}
bool operator <(const Number &b) const{
for (int i=0; i<N-1; ++i) if (a[i]!=b.a[i]) return a[i]<b.a[i];
return 0;
}
Number operator <<(const int &c) const{
//less than N
Number ret;
ull tmp=0;
for (int i=0; i<N-1; ++i) ADDT(i+c==N-1?tmp:ret.a[i+c>=N?i+c-N:i+c],a[i]);
for (int i=0; i<N-1; ++i) SUBT(ret.a[i],tmp);
return ret;
}
};
template<int N>
Number<N> fp(Number<N> x,int y){
Number<N> ret;
ret.a[0]=1;
for (; y; y>>=1,x=x*x) if (y&1) ret=ret*x;
return ret;
}
const ull inv5=14757395258967641293ull;
template<int N,int LEN>
struct FWT_helper{
typedef Number<N> Num;
Num w[N][N];
int pwd[LEN+1];
FWT_helper(){
assert(N>1);
for (int i=0; i<N; ++i) w[i][0].a[0]=1;
w[0][1].a[0]=1;
w[1][1].a[1]=1;
for (int i=2; i<N; ++i) w[i][1]=w[i-1][1]*w[1][1];
for (int i=0; i<N; ++i)
for (int j=2; j<N; ++j)
w[i][j]=w[i][j-1]*w[i][1];
pwd[0]=1;
for (int i=1; i<=LEN; ++i) pwd[i]=pwd[i-1]*N;
}
void DFT(Num *a){
Num b[N];
for (int i=0; i<N; ++i){
//for (int j=0; j<N; ++j) b[i]=b[i]+a[j]*w[i][j];???
for (int j=0,d=0; j<N; ++j,(d+=i)>=N?d-=N:0)
b[i]=b[i]+(a[j]<<d);
}
for (int i=0; i<N; ++i) a[i]=b[i];
}
void IDFT(Num *a){
reverse(a+1,a+N);
DFT(a);
//Don't forget to div
}
void FWT(Num *a,int m,bool rev){
for (int i=0; i<m; ++i)
for (int cj=0; cj<pwd[m]; cj+=pwd[i+1])
for (int j=cj; j<cj+pwd[i]; ++j){
Num tmp[N];
for (int k=0; k<N; ++k) tmp[k]=a[j+pwd[i]*k];
if (rev) IDFT(tmp);
else DFT(tmp);
for (int k=0; k<N; ++k) a[j+pwd[i]*k]=tmp[k];
}
}
};
int n;
Number<10> a[100000];
int main(){
FWT_helper<10,5> ups;
read(n);
for (int i=0; i<n; ++i){
int x;
read(x);
++a[x].a[0];
}
ups.FWT(a,5,0);
for (int i=0; i<ups.pwd[5]; ++i) a[i]=fp(a[i],n);
ups.FWT(a,5,1);
for (int i=0; i<n; ++i){
//cerr<<a[0]<<" "<<(~0ull)<<endl;
for (int j=8; j>=7; --j){
a[i].a[j-5]-=a[i].a[j];
}
//a[i].a[3]-=a[i].a[6];
a[i].a[1]-=a[i].a[6];
a[i].a[0]-=a[i].a[5];
a[i].a[3]+=a[i].a[4];
a[i].a[2]-=a[i].a[4];
a[i].a[1]+=a[i].a[4];
a[i].a[0]-=a[i].a[4];
//cerr<<a[i].a[0]<<endl;
//exit(0);
}
ull inv=inv5*inv5*inv5*inv5*inv5;
for (int i=0; i<ups.pwd[5]; ++i) a[i].a[0]=MUL(a[i].a[0],inv);
for (int i=0; i<n; ++i){
//for (int j=1; j<9; ++j) assert(a[i].a[j]==0);
cout<<(a[i].a[0]>>5)%(1ll<<58)<<'
';
}
}
另一种做法是直接模10的分圆多项式，但是比较难写，因为乘法高次项的贡献比较恶心。
然后考虑将10次单位根转化为5次单位根，有(w^5_{10}=-1)
转化了之后直接对长度为4的搞，也可以。
#pragma GCC optimize(3)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace io {
const int SIZE = 7e5+10;
char buff[SIZE];
char *l = buff, *r = buff;
void init() {
l = buff;
r = l + fread(buff, 1, SIZE, stdin);
}
char gc() {
if (l == r) init();
if(l==r)return EOF;
return *(l++);
}
void read(int &x) {
x = 0;
char ch = gc();
while(!isdigit(ch)) ch = gc();
while(isdigit(ch)) x = x * 10 + ch - '0', ch = gc();
}
}
using io::read;
typedef unsigned long long ull;
template<int N>
struct Number{
ull a[N-1];
Number(){
memset(a,0,sizeof(a));
}
friend ostream & operator <<(ostream &out,const Number<N> &b){
out<<'{';
for (int i=0; i<N-2; ++i) out<<b.a[i]<<',';
return out<<b.a[N-2]<<'}';
}
Number operator +(const Number &b) const{
Number ret;
for (int i=0; i<N-1; ++i) ret.a[i]=a[i]+b.a[i];
return ret;
}
Number operator -(const Number &b) const{
Number ret;
for (int i=0; i<N-1; ++i) ret.a[i]=a[i]-b.a[i];
return ret;
}
Number operator *(const Number &b) const{
Number ret;
ull tmp=0;
for (int i=0; i<N-1; ++i)
for (int j=0; j<N-1; ++j)
(i+j==N-1?tmp:ret.a[i+j>=N?i+j-N:i+j])+=a[i]*b.a[j];
for (int i=0; i<N-1; ++i) ret.a[i]-=tmp;
return ret;
}
bool operator <(const Number &b) const{
for (int i=0; i<N-1; ++i) if (a[i]!=b.a[i]) return a[i]<b.a[i];
return 0;
}
Number operator <<(const int &c) const{
//less than N
Number ret;
ull tmp=0;
for (int i=0; i<N-1; ++i) (i+c==N-1?tmp:ret.a[i+c>=N?i+c-N:i+c])+=a[i];
for (int i=0; i<N-1; ++i) ret.a[i]-=tmp;
return ret;
}
};
typedef Number<5> Num;
template<int N>
Number<N> fp(Number<N> x,int y){
Number<N> ret;
ret.a[0]=1;
for (; y; y>>=1,x=x*x) if (y&1) ret=ret*x;
return ret;
}
void DFT(Num *a){
Num b[10];
for (int i=0; i<10; ++i)
for (int j=0,t=0; j<10; ++j,(t+=i)>=10?t-=10:0){
if (~t&1) b[i]=b[i]+(a[j]<<(t>>1));
else b[i]=b[i]-(a[j]<<((t+5)>>1)%5);
}
for (int i=0; i<10; ++i) a[i]=b[i];
}
const int pwd[6]={1,10,100,1000,10000,100000};
const ull inv5=14757395258967641293ull;
void FWT(Num *a,bool rev){
for (int i=0; i<5; ++i){
for (int cj=0; cj<100000; cj+=pwd[i+1])
for (int j=cj; j<cj+pwd[i]; ++j){
Num t[10];
for (int k=0; k<10; ++k) t[k]=a[j+k*pwd[i]];
if (rev) reverse(t+1,t+10);
//for (int k=0; k<10; ++k) cout<<t[k]<<endl;
DFT(t);
//for (int k=0; k<10; ++k) cout<<t[k]<<endl;
//exit(0);
for (int k=0; k<10; ++k) a[j+k*pwd[i]]=t[k];
}
}
if (rev){
ull c=inv5*inv5*inv5*inv5*inv5;
for (int i=0; i<100000; ++i) a[i].a[0]*=c;
}
}
const int N=100000;
int n;
Num a[N];
int main(){
read(n);
for (int i=0; i<n; ++i){
int x; read(x);
++a[x].a[0];
}
FWT(a,0);
for (int i=0; i<100000; ++i) a[i]=fp(a[i],n);
FWT(a,1);
for (int i=0; i<n; ++i)
cout<<(a[i].a[0]>>5)%(1ll<<58)<<'
';
}