codeforces 1103E清华集训 石家庄的工人阶级队伍比较坚强

Posted 2021-02-10 yuhuger

希望复习高进制FWT的时候，能够快速回想起来。
FWT感觉就是每一维单独考虑，（虽然我不知道为什么这是对的）
分别对一个奇怪的东西做DFT，
那个奇怪的东西在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)<<'
';
}