该题果然是个好题啊!
题意来自上一题, ( 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;
}