[组合数学][数论][矩阵][优化][Polya][状压DP][x] XDOJ 1296 敬老师的手环II

去年校赛的时候感觉此题难的一批，，，昨天做polya练习题的时候做到了poj 2888去翻了下matrix67大神的十个利用矩阵乘法解决的经典题目中的第9题才恍然大悟。。。Oooooor2 做完这俩题再来看XDOJ的这个题已经是自然而然就可以搞出结果的了。

这个题的思想就是：把每一列3*1的小方格当成一个“珠子”，把他可能的8种情况当成“颜色”，然后这些颜色同时也代表状态，类似于poj2888的颜色相邻限制，这里的状态之间“能不能相邻”其实就是在状压dp里“能不能转移”的限制，不同状态之间的转移图如下（来自matrix67的那篇文章

根据转移图可以做出来一个转移矩阵（用矩阵可以快速幂加速啊这个题n那么大）。然后这个矩阵的作用就是对于循环节个数为k的置换求出其方案数，求的原理就是polya的思想：循环节内部涂上同一种颜色，说白了就是把循环节内部缩成一个点，然后求有多少可以在指定的限定条件下正好走k步走到自身的方法，也就是把矩阵^k的对角线上的数字求和（这一步如果不懂的话移步poj2888的题解，各路大佬已经讲的很清楚了）

此外，polya中常用的枚举循环节长度然后用欧拉函数求的方法自然也不必说肯定是要用的，但是这个题有个非常坑爹的问题！时限实在是太紧了！！！！ 1s时限处理1000个1e9的数字我优化到吐血好吗？？？？出题人真的是丧心病狂。

所以，第一版的算法就正确的情况下tle了8次。。。各种优化都用上了。。。真心给跪了，最后参考了同校tdl的题解中的几个优化之处终于给过掉了。。。主要是

1.欧拉函数的求法和遍历因子的处理互相考虑减少运算，只利用一次质因数分解

2.矩阵的乘法中有些位置的元素恒为0从而可以不用计算

3.矩阵快速幂可以把中间过程的矩阵存下来，因为每次的底数矩阵都是同一个，可以减少大量重复计算

最后终于用242ms的速度通过了此题………

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <map>
using namespace std;
 
typedef long long LL;
const int N = 8;
const int mod = 1e9 + 7;
const int NUM = 40004;
int Num_Prime = 0;
bool visit[NUM];
LL n;
bool check[NUM + 10];
LL phi[NUM + 10];
LL prime[NUM + 10];
int tot;
void phi_and_prime_table()
{
    int N = NUM - 1;
    //phi[1] = 1;
    tot = 0;
    for (int i = 2; i <= N; i++)
    {
        if (!check[i])
        {
            prime[tot++] = i;
            //phi[i] = i - 1;
        }
        for (int j = 0; j < tot; j++)
        {
            if (i * prime[j] > N)
                break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                //phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            // else
            // {
            //     phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            // }
        }
    }
}
 
LL mul(LL a, LL b)
{
    if (a >= mod)
        a %= mod;
    if (b >= mod)
        b %= mod;
    a *= b;
    return a > mod ? a % mod : a;
}
 
LL add(LL a, LL b)
{
    return a + b > mod ? a + b - mod : a + b;
}
 
struct Matrix
{
    LL a[8][8];
    void init(int x)
    {
        memset(a, 0, sizeof(a));
        if (x)
            for (int i = 0; i < 8; i++)
                a[i][i] = 1;
    }
} a, b, c,d[32];
 
 
 
 
/*
1 0 0 1 0 0 1 0
0 1 0 0 1 0 0 1
0 0 1 0 0 0 0 0
1 0 0 1 0 0 1 0
0 1 0 0 1 0 0 1
0 0 0 0 0 1 0 0
1 0 0 1 0 0 1 0
0 1 0 0 1 0 0 1
*/
Matrix fuck(Matrix a, Matrix b)
{
    Matrix p;
    p.init(0);
    // for (int i = 0; i < 8; i++)
    //     for (int j = 0; j < 8; j++)
    //     {
    //         for (int k = 0; k < 8; k++)
    //             p.a[i][j] = add(p.a[i][j], mul(a.a[i][k], b.a[k][j]));
    //     }
    for(int k=0;k<3;k++)
        for(int i=3*k;i<3*k+2;i++)
            for(int j=i%3;j<8;j+=3)
                for(int k=0;k<8;k++)
                {
                    p.a[i][j] = add(p.a[i][j], mul(a.a[i][k], b.a[k][j]));
                }
    for(int k=0;k<8;k++)
    {
            p.a[2][2] = add(p.a[2][2], mul(a.a[2][k], b.a[k][2]));
            p.a[5][5] = add(p.a[5][5], mul(a.a[5][k], b.a[k][5]));
    }
    return p;
}
 
void power(int t)
{
    int tp=0;
    b.init(1);
    while (t)
    {
        if (t & 1)
            b = fuck(b, d[tp]);
        t >>= 1;
        tp++;
    }
}
 
LL modPow(LL s, LL index, LL mod)
{
    LL ans = 1;
    s %= mod;
    while (index)
    {
        if ((index & 1))
            ans = mul(ans, s);
        index >>= 1;
        s = mul(s, s);
    }
    return ans;
}
 
 
 
map<int, LL> PP;
 
LL Cal(int L)
{
    if (PP[L])
    {
        return PP[L];
    }
    power(L);
    LL res = 0;
    for (int i = 0; i < 8; i++)
        res = add(res, b.a[i][i]);
    return PP[L] = res % mod;
}
 
map<int, LL> mp;
 
LL factor[100];
LL divs[5000],dtot;
int fatCnt,pt;
void getFactors(LL x)
{
    dtot=fatCnt = 0;
    LL tmp = x;
    divs[dtot++]=1;
    for (int i = 0; prime[i]*prime[i] <= tmp ; i++)
    {
        //factor[fatCnt][1] = 0;
        if (tmp % prime[i] == 0)
        {
            factor[fatCnt] = prime[i];
            pt=dtot;
            while (tmp % prime[i] == 0)
            {
                //factor[fatCnt][1]++;
                tmp /= prime[i];
                for(int j=0;j<pt;j++)
                {
                    divs[dtot]=divs[dtot-pt]*prime[i];dtot++;
                }
            }
            fatCnt++;
        }
    }
    if (tmp != 1)
    {
        factor[fatCnt] = tmp;
        fatCnt++;
        //factor[fatCnt++][1] = 1;
        pt=dtot;
        for(int j=0;j<pt;j++)
        {
            divs[dtot]=divs[dtot-pt]*tmp;dtot++;
        }
    }
    //return fatCnt;
}
 
 
 
int main()
{
    phi_and_prime_table();
    a.a[0][0] = 0, a.a[0][1] = 0, a.a[0][2] = 0, a.a[0][3] = 0, a.a[0][4] = 0, a.a[0][5] = 0, a.a[0][6] = 0, a.a[0][7] = 1;
    a.a[1][0] = 0, a.a[1][1] = 0, a.a[1][2] = 0, a.a[1][3] = 0, a.a[1][4] = 0, a.a[1][5] = 0, a.a[1][6] = 1, a.a[1][7] = 0;
    a.a[2][0] = 0, a.a[2][1] = 0, a.a[2][2] = 0, a.a[2][3] = 0, a.a[2][4] = 0, a.a[2][5] = 1, a.a[2][6] = 0, a.a[2][7] = 0;
    a.a[3][0] = 0, a.a[3][1] = 0, a.a[3][2] = 0, a.a[3][3] = 0, a.a[3][4] = 1, a.a[3][5] = 0, a.a[3][6] = 0, a.a[3][7] = 1;
    a.a[4][0] = 0, a.a[4][1] = 0, a.a[4][2] = 0, a.a[4][3] = 1, a.a[4][4] = 0, a.a[4][5] = 0, a.a[4][6] = 0, a.a[4][7] = 0;
    a.a[5][0] = 0, a.a[5][1] = 0, a.a[5][2] = 1, a.a[5][3] = 0, a.a[5][4] = 0, a.a[5][5] = 0, a.a[5][6] = 0, a.a[5][7] = 0;
    a.a[6][0] = 0, a.a[6][1] = 1, a.a[6][2] = 0, a.a[6][3] = 0, a.a[6][4] = 0, a.a[6][5] = 0, a.a[6][6] = 0, a.a[6][7] = 1;
    a.a[7][0] = 1, a.a[7][1] = 0, a.a[7][2] = 0, a.a[7][3] = 1, a.a[7][4] = 0, a.a[7][5] = 0, a.a[7][6] = 1, a.a[7][7] = 0;
    d[0]=a;
    for(int i=1;i<=30;i++)
    {
        d[i]=fuck(d[i-1],d[i-1]);
    }
    while (~scanf("%lld", &n))
    {
        if (n & 1)
        {
            printf("0\n");
            continue;
        }
        if (mp[n])
        {
            printf("%lld\n", mp[n]);
            continue;
        }
        LL ans=0,eular;
        getFactors(n);
        for(int i=0;i<dtot;i++)
        {
            eular=divs[i];
            for(int j=0;j<fatCnt;j++)
            {
                if(divs[i]%factor[j]==0)
                {
                    eular=eular*(factor[j]-1)/factor[j];
                }
            }
            ans=add(ans,mul(eular,Cal(n/divs[i])));
        }
        ans = mul(ans, modPow(n % mod, mod - 2, mod));
        printf("%lld\n", ans);
        mp[n] = ans;
    }
    return 0;
}