loj6271「长乐集训 2017 Day10」生成树求和 加强版

又是一个矩阵树套多项式的好题。

这里我们可以对每一位单独做矩阵树,但是矩阵树求的是边权积的和,而这里我们是要求加法,于是我们i将加法转化为多项式的乘法,其实这里相当于一个生成函数?之后如果我们暴力做的话,就是强行带入x插值,复杂度$O(8*2n*n^{3})$,还不够优秀,于是我们考虑用$dft$优化这个过程,这里我们需要找到一个三次单位根,于是我们考虑扩域的思想,我们把数表示为$(a+b*w_{3})$,这里$w_{3}$满足$w_{3}^{3}=1$且$w_{3}^{1}+w_{3}^{2}+w_{3}^{3}=0$,于是我们可以计算出这类数的计算法则,然后我们直接将矩阵dft之后做一遍矩阵树,之后再将答案逆变换回去,就是我们需要的答案了。

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <iostream>
  4 #include <algorithm>
  5 #include <cmath>
  6 #define N 111
  7 #define mod 1000000007
  8 #define inv3 333333336
  9 using namespace std;
 10 int qp(int a,int b){
 11     int c=1;
 12     for(;b;b>>=1,a=1ll*a*a%mod)
 13         if(b&1)c=1ll*c*a%mod;
 14     return c;
 15 }
 16 struct num{
 17     int a,b;
 18     num(){}
 19     num(int x,int y){a=x;b=y;}
 20     num operator + (num x){return num((a+x.a)%mod,(b+x.b)%mod);}
 21     num operator - (num x){return num((a-x.a+mod)%mod,(b-x.b+mod)%mod);}
 22     num operator * (num x){
 23         return num(
 24             ((1ll*a*x.a-1ll*b*x.b)%mod+mod)%mod,
 25             ((1ll*a*x.b+1ll*b*x.a-1ll*b*x.b)%mod+mod)%mod);
 26         }
 27     num inv(){
 28         int val=qp(((1ll*a*b-1ll*a*a-1ll*b*b)%mod+mod)%mod,mod-2);
 29         return num(1ll*(b-a+mod)*val%mod,1ll*b*val%mod);
 30     }
 31     bool operator ! (){return !a&&!b;}
 32 }A[N][N][3],w[3],det[3];
 33 int n,m,ans,du[N*N],dv[N*N],dw[N*N];
 34 int cnt;
 35 void work(){
 36     num t;
 37     for(int d=0;d<=2;d++){
 38         det[d]=num(1,0);
 39         for(int k=1;k<n;k++){
 40             if(!A[k][k][d]){
 41                 for(int i=k+1;i<=n;i++){
 42                     if(!(!A[i][k][d])){
 43                         det[d]=num(0,0)-det[d];
 44                         for(int j=k;j<=n;j++)swap(A[k][j][d],A[i][j][d]);
 45                         break;
 46                     }
 47                 }
 48                 if(!A[k][k][d]){det[d]=num(0,0);break;}
 49             }
 50             det[d]=det[d]*A[k][k][d];
 51             for(int i=k+1;i<n;i++){
 52                 t=A[i][k][d]*A[k][k][d].inv();
 53                 for(int j=k;j<=n;j++)
 54                     A[i][j][d]=A[i][j][d]-t*A[k][j][d];
 55             }
 56         }
 57     }
 58 }
 59 void dft(num *a){
 60     num b[3];
 61     b[0]=a[0]+a[1]+a[2];
 62     b[1]=a[0]+a[1]*w[1]+a[2]*w[2];
 63     b[2]=a[0]+a[1]*w[2]+a[2]*w[1];
 64     a[0]=b[0];a[1]=b[1];a[2]=b[2];
 65 }
 66 void add(int u,int v,int w){
 67     A[u][u][w]=A[u][u][w]+num(1,0);
 68     A[v][v][w]=A[v][v][w]+num(1,0);
 69     A[u][v][w]=A[u][v][w]-num(1,0);
 70     A[v][u][w]=A[v][u][w]-num(1,0);
 71 }
 72 void UPD(int &a,int b){
 73     a=(a+b>=mod)?(a+b-mod):(a+b);
 74 }
 75 int main(){
 76 freopen("sum.in","r",stdin);
 77 freopen("sum.out","w",stdout);
 78     w[0]=num(1,0);w[1]=num(0,1);
 79     w[2]=num(mod-1,mod-1);
 80     scanf("%d%d",&n,&m);
 81     for(int i=1;i<=m;i++)
 82         scanf("%d%d%d",&du[i],&dv[i],&dw[i]);
 83     for(int t=0,now=1;t<=8;t++,now=now*3){
 84         memset(A,0,sizeof A);
 85         for(int i=1;i<=m;i++){
 86             add(du[i],dv[i],dw[i]%3);
 87             dw[i]/=3;
 88         }
 89         for(int i=1;i<=n;i++)
 90             for(int j=1;j<=n;j++)
 91                 dft(A[i][j]);
 92         cnt=t;
 93         work();            
 94         swap(w[1],w[2]);
 95         dft(det);
 96         swap(w[1],w[2]);
 97         UPD(ans,1ll*det[1].a*inv3%mod*now%mod);
 98         UPD(ans,1ll*det[2].a*inv3%mod*2*now%mod);
 99     }
100     printf("%d\n",ans);
101     return 0;
102 }
View Code

 

posted @ 2018-04-21 07:46  Ren_Ivan  阅读(...)  评论(... 编辑 收藏