HDU2481 Toy (数论好题)

该题果然是个好题啊！

题意来自上一题， ( http://blog.csdn.net/jayye1994/article/details/37814965 )  BZOJ 1002: [FJOI2007]轮状病毒

上一题是旋转后相同视为不同情况，这题旋转后相同视为同一种情况。就这么一个小小的区别，

上一题用到了dp，这一题用到了dp、筛素数、二进制模拟乘法、矩阵、快速幂、欧拉函数、burnside引理。。

同样是人(题)。。做人(题)的差距怎么就这么大呢。。

数据范围：周围是n个点，答案取模m， n<=10^9, m <= 10^9

思路：

要做这题首先是肯定要会burnside引理的，不会的可以去做几题，刘汝佳的训练指南上有。

然后考虑有n种置换，分别是不动、旋转i个单位(i < n)。

首先解决不动的情况，实际上就是轮状病毒要求的结果，可是这里的n<=10^9，那题我计算的复杂度是O(n^2)的，

计算公式也不能快速幂，然后那题的dp公式可以推出dp[i][1] = 3*dp[i-1][1] - dp[i-2][1]，推结果ans的话我就推晕了，

那么打个表就可以发现规律。。ans[i] = 3*ans[i-1] - ans[i-2] + 2，于是不动的情况解决了。。。

接下来就是要求旋转的情况了，旋转i个单位，很容易可以得知循环节是n / gcd(n, i) ，循环的个数为gcd(i, n),所以每gcd(i,n)构成一块，每一块都是完全相同的，如下图所示

想了下还是把上一题的dp贴过来，这样子看起来比较清楚。

题目就是求最小生成树的种数，中心的点必须要和周围的一圈点连通，可以先不要管环，先考虑链的情况

dp[ i ][ 0 ]表示第i个点还没和中心点连通，并且前i-1个点和中心点或者第i个点是连通的

dp[ i ][ 1 ]表示前i个点全部都已经和中心点连通了

很容易可以推出状态转移方程 ：

dp[ i ][ 0 ] = dp[ i-1 ][ 0 ] + dp[ i-1 ][ 1 ]

dp[ i ][ 1 ] = dp[ i-1 ][ 0 ] + dp[ i-1 ][ 1 ]*2

那么对于n轮状病毒，可以枚举第一个点和多少个周围的点连通，其余的点就是一条链的情况了。

比如说i=2，n=6，那么gcd(2,6) = 2，循环长度是2，也就是说红色线段和绿色线段在重复出现，那么可以知道中间的所有的点都要是连通的，如果两侧的点也都和中心点连通了的话，绿色线段就要删掉，情况数就是dp[ gcd(i,n) ][ 1 ]。求循环长度是d的个数可以用欧拉函数来求，所以不需要枚举i(i<=n)，只需要枚举n的因子就可以了。

如果某一侧点没有和中心点连通的话，绿色线段是必要的，情况数就是2*dp[ gcd(i,n) ][0]，因为有两侧所以乘2，但是要注意如果这一段所有的点都没和中心点连通的话是不行了，所以情况数还需要减去 2 * 1。

然后问题差不多已经解决了，不过还有一个问题，最后算出的答案还需要除以n，但是n和m不一定互质，所以不能求逆元

所以是求(a/b) mod c，因为知道a/b是一个整数，所以ans = a/b + xc，然后ans * b = a + xbc = a mod (bc)，所以可以求出a mod (bc)然后除以b就可以了，中间乘法会爆long long，所以要二进制位上加法来模拟乘法。

code：

#include 
   
   
    
    
#include 
    
    
     
     
#include 
     
     
      
      
using namespace std;
typedef __int64 ll;

const int P = 100000+5;
int pri[P], pnum;
bool vis[P];
ll MOD;

void get_prime(int n) {
    for(int i = 2;i <= n; i++) {
        if(!vis[i]) {
            pri[pnum++] = i;
        }
        for(int j = 0;j < pnum; j++) {
            if(i*pri[j] > n)    break;
            vis[i*pri[j]] = 1;
            if(i%pri[j] == 0)   break;
        }
    }
}

ll Mul(ll a, ll b) {
    ll ret = 0;
    for(int i = 0;i < 63; i++) {
        if(b>>i&1)  ret += a;
        if(ret >= MOD)  ret -= MOD;
        a = a+a;
        if(a >= MOD)    a -= MOD;
    }
    return ret;
}

void Add(ll &x, ll y) {
    x += y;
    if(x >= MOD)    x -= MOD;
    if(x < 0)   x += MOD;
}

int euler(int n) {
    int ret = n;
    for(int i = 0;i < pnum && pri[i]*pri[i] <= n ;i++) if(n%pri[i] == 0) {
        while(n%pri[i] == 0)   n /= pri[i];
        ret -= ret/pri[i];
    }
    if(n > 1)   ret -= ret/n;
    return ret;
}

const int D = 3;
struct matrix {
    ll a[D][D];
    matrix() {
        memset(a, 0, sizeof(a));
    }
};

matrix mul(matrix &A, matrix &B) {
    matrix C;
    for(int i = 0;i < D; i++) {
        for(int j= 0;j < D; j++) {
            C.a[i][j] = 0;
            for(int k = 0;k < D; k++)
                Add(C.a[i][j], Mul(A.a[i][k],B.a[k][j]));
        }
    }
    return C;
}

ll pow_mod(int n, matrix &ret, matrix &x) {
    while(n) {
        if(n&1) ret = mul(ret, x);
        x = mul(x, x);
        n >>= 1;
    }
    return ret.a[0][1];
}

ll gao1(int n) {
    matrix ret, x;
    ret.a[0][0] = 1; ret.a[0][1] = 5; ret.a[0][2] = 2;
    x.a[0][1] = -1+MOD; x.a[1][0] = 1; x.a[1][1] = 3;
    x.a[2][1] = x.a[2][2] = 1;
    ll ans = pow_mod(n-2, ret, x);
    return ans;
}

ll gao2(int n, int i) {
    matrix ret, x;
    ret.a[0][0] = ret.a[0][1] = 1;

    x.a[0][0] = x.a[0][1] = x.a[1][0] = 1;
    x.a[1][1] = 2;
    pow_mod(i-1, ret, x);
    ll ans = ret.a[0][1];
    //printf("a0 = %d a1 = %d\n", ret.a[0][0], ret.a[0][1]);
    Add(ans, (ret.a[0][0]-1)*2%MOD);
    //printf("ans = %d\n", ans);
    return ans;
}

ll solve(int n) {
    ll ret = 0;
    Add(ret, euler(n));
    //printf("ret = %I64d\n", ret);
    Add(ret, gao1(n));
    //printf("ret = %I64d\n", ret);
    for(int i = 2;i*i <= n; i++) if(n%i == 0){
        int sum = euler(n/i);
        Add(ret, Mul(sum,gao2(n, i)));
      //  printf("ret = %d\n", ret);
        if(i*i != n) {
            sum = euler(i);
            Add(ret, Mul(sum,gao2(n, n/i)));
        //    printf("ret = %d\n", ret);
        }
    }
    return ret/n;
}

int main() {
    get_prime(100000);
    int n;
    while(scanf("%d%I64d", &n, &MOD) == 2) {
        MOD *= n;
        printf("%I64d\n", solve(n));
    }
    return 0;
}