zoukankan      html  css  js  c++  java
  • 【POJ2888】Magic Bracelet-Burnside引理+数论+DP矩阵优化

    测试地址:Magic Bracelet
    题目大意:要用N(109)颗珠子连成环形的手镯,共有M(10)种不同的珠子,规定K个条件,每个条件规定某两种珠子不能相邻,旋转后相同的方案视作相同,问有多少种本质不同的方案,对9973取模(gcd(N,9973)=1)。
    做法:这题非常经典,思路很有代表性,是进阶Burnside引理和Polya定理题的一块敲门砖,需要用到:Burnside引理,欧拉函数,DP+矩阵快速幂,乘法逆元。
    首先一看题目是计数,就自然而然地联想到Burnside引理和Polya定理,然而这题有不能相邻的条件,所以不能直接用Polya定理,那么我们考虑使用Burnside引理来解决。
    我们知道环形的旋转置换有N种,第i种置换(这里定义为旋转i颗珠子的置换)的循环数为gcd(N,i),观察这样的置换,我们发现它很有规律:第1个元素属于第1个循环,第2个元素属于第2个循环,…,第gcd(N,i)个元素属于第gcd(N,i)个循环,第gcd(N,i)+1个元素属于第1个循环,以此类推。由于我们计算|C(f)|时每个循环内元素着色应该相同,那么我们实际上只需确定前gcd(N,i)个元素的填色方案数即可。设对前i个元素进行着色,且最后一个元素颜色为j的方案数为f(i,j),并设二元函数m(i,j),表示如果i,j两种颜色的珠子可以相邻,那么m(i,j)=1,否则m(i,j)=0,那么我们可以得到状态转移方程:
    f(i,j)=Mk=1m(j,k)×f(i1,k)
    但是最后一个元素的填色除了和它前一个元素有关,还和第一个元素有关,那么这个状态转移方程是不是错的呢?实际上,我们只需一开始枚举第一个元素的颜色k,然后这一种情况的答案就是第gcd(N,i)+1个元素与第一个元素同色的方案数,即:f(gcd(N,i)+1,k),累加每种情况即可得出|C(f)|。带进Burnside公式计算之后,最外面还要乘一个1/N,由于gcd(N,9973)=1,所以我们直接求N关于模9973的乘法逆元,然后再乘即可。
    上面的方法的时间复杂度是O(N2M),显然不能通过此题,需要考虑优化。此题需要从两个方面入手,一是计算Burnside公式的复杂度,而是动态规划的复杂度。
    首先优化计算Burnside公式的时间复杂度。我们知道gcd(N,i)|N,而且当gcd(N,i)相同时,|C(f)|相同,那么我们完全可以枚举N的因子d,计算gcd(N,i)=d时的|C(f)|,再乘上使得gcd(N,i)=di的个数,累加起来得到答案。后面那个个数显然就是φ(N/d)了(不懂的可以去看看我写的POJ2154的题解),那么枚举的时间复杂度就大大降低了。
    接下来优化DP的时间复杂度,注意到M很小,那么二元函数m显然可以构造成一个矩阵M,如果我们把f(i,1),f(i,2),...,f(i,M)从上到下排成一个列向量Fi,那么可以看出:Fi=MFi1,所以Fi=Mi1F1,于是我们可以先用矩阵快速幂算出Md,然后就可以算Fd+1了,但其实我们没有必要算出全部的Fd+1,因为我们只要求f(d+1,k),这是第k行第1列的元素,那么就用Md的第k行与F1的第1列做运算,根据定义,F1中只有第k行为1,其余都为0,所以f(d+1,k)=Mdkk,那么|C(f)|就应该等于Md对角线上的元素和。
    经过了以上两个优化,外层枚举的复杂度降到差不多O(N),DP的时间复杂度从O(NM)降到O(M3logN),因此总的时间复杂度为O(M3NlogN),可以通过此题。
    以下是本人代码:

    #include <cstdio>
    #include <cstdlib>
    #include <cstring>
    #include <iostream>
    #include <algorithm>
    #define ll long long
    using namespace std;
    ll n,ans;
    int T,m,k;
    struct matrix {ll s[11][11];} M[40];
    
    void exgcd(ll a,ll b,ll &x,ll &y)
    {
      ll x0=1,y0=0,x1=0,y1=1;
      while(b)
      {
        ll tmp,q;
        q=a/b;
        tmp=x0,x0=x1,x1=tmp-q*x1;
        tmp=y0,y0=y1,y1=tmp-q*y1;
        tmp=a,a=b,b=tmp%b;
      }
      x=x0,y=y0;
    }
    
    matrix mult(matrix A,matrix B)
    {
      matrix S;
      memset(S.s,0,sizeof(S.s));
      for(int i=1;i<=m;i++)
        for(int j=1;j<=m;j++)
          for(int k=1;k<=m;k++)
            S.s[i][j]=(S.s[i][j]+A.s[i][k]*B.s[k][j])%9973;
      return S;
    }
    
    matrix power(ll x)
    {
      int i=0;
      matrix S;
      memset(S.s,0,sizeof(S.s));
      for(int j=1;j<=m;j++) S.s[j][j]=1;
      while(x)
      {
        if (x&1) S=mult(S,M[i]);
        i++;x>>=1;
      }
      return S;
    }
    
    ll phi(ll x)
    {
      ll s=x;
      for(ll i=2;i*i<=x;i++)
        if (!(x%i))
        {
          s=s/i*(i-1);
          while(!(x%i)) x/=i;
        }
      if (x>1) s=s/x*(x-1);
      return s;
    }
    
    void solve(ll x)
    {
      matrix S=power(x);
      ll sum=0;
      for(int i=1;i<=m;i++)
        sum=(sum+S.s[i][i])%9973;
      sum=(sum*phi(n/x))%9973;
      ans=(ans+sum)%9973;
    }
    
    int main()
    {
      scanf("%d",&T);
      while(T--)
      {
        scanf("%lld%d%d",&n,&m,&k);
        ans=0;
        memset(M[0].s,0,sizeof(M[0].s));
        for(int i=1;i<=m;i++)
          for(int j=1;j<=m;j++)
            M[0].s[i][j]=1;
        for(int i=1,a,b;i<=k;i++)
        {
          scanf("%d%d",&a,&b);
          M[0].s[a][b]=M[0].s[b][a]=0;
        }
    
        for(int i=1;i<=35;i++) M[i]=mult(M[i-1],M[i-1]);
        for(ll i=1;i*i<=n;i++)
          if (n%i==0)
          {
            solve(i);
            if (i!=n/i) solve(n/i);
          }
    
        ll x0,y0;
        exgcd(n,9973,x0,y0);
        x0=(x0%9973+9973)%9973;
        printf("%lld
    ",(x0*ans)%9973);
      }
    
      return 0;
    }
    
  • 相关阅读:
    leetcode 268. Missing Number
    DBSCAN
    python二维数组初始化
    leetcode 661. Image Smoother
    leetcode 599. Minimum Index Sum of Two Lists
    Python中的sort() key含义
    leetcode 447. Number of Boomerangs
    leetcode 697. Degree of an Array
    滴滴快车奖励政策,高峰奖励,翻倍奖励,按成交率,指派单数分级(1月3日)
    北京Uber优步司机奖励政策(1月2日)
  • 原文地址:https://www.cnblogs.com/Maxwei-wzj/p/9793641.html
Copyright © 2011-2022 走看看