HDU 3441 Rotation(Polya计数)

题意：将一个由A*A个方格组成的方块，分成k个B*B的方块和1个小方格后连乘一个由k+1个物体组成的环形物，如图，A=3时有两种分法。用C种颜色为每个方格染色，旋转得到的两种方案记为相同方案，问一共有多少种不同方案。

解法：如果将B*B的小正多边形看做一个整体，则B*B的方案数即外围“长度为k的项链”的颜色数，因此此题要嵌套两个polya分别对B*B和K+1个物体的方案数分别计算。

1.由于A巨大，因此不能提供枚举A*A-1的质因子确定B，首先将A*A-1化为(A+1)(A-1)然后分别求得A+1和A-1的质因子，然后将两个数的质因子分解结果合并，由于gcd(A+1,A-1)=1或2，因此对于A是奇数 的情况合并时需要合并2的幂数。由此得到若干对<B,K>；由于B*B*k=(A+1)(A-1),因此在枚举B是可以顺便处理k的质因子分解：初始化时与B*B的分解结果相同，每次dfs时B的质因子p增加delta倍，k的质因子p减少2delta倍。

2.对于一个B*B的正方形，共有四种置换(A,B,C,D)  (C,A,D,B) (D,C,B,A) (B,D,A,C)

          对于第1中置换，方案数为C^(B*B)

          对于第三种置换，如果B为偶数，根据前上面的B/2行便可得到下面的B/2行；若B为奇数，第B/2+1行任意染色，因此方案数C^(B*B+1)/2

          对于第二四种置换，如果B为偶数，左上角四分之一区域染色方案唯一确定整个B*B的染色方案；如果B为奇数，左上角(B+1)/2*(B-1)/2唯一确定其他三个角，同时中间格子随意染色，因此方案数为C^(B*B+3)/4

3.若B*B的不同方案是为c，则问题变为k个珠子的项链用c种颜色染色有多少种方案，通过k的质因子dfs求得所有可能的循环节个数得到方案数；同时中间的珠子有C中颜色，因此结果应乘C。

import java.util.Scanner;
public class Rotation3441 {
	long mod = 1000000007;
	long a,c,res,ans,k,color;
	long prime1[]=new long[100010],prime2[]=new long[100010];
	int count1[]=new int[100010],count2[]=new int[100010];
	void divide(long n,long[] prime,int[] count) {
		int cnt=count[0];//质因子数目
		long mx = (long) (Math.sqrt(n) + 3);
		if (n < mx)
			mx = n;
		for (int i = 2; i < mx; i++) {
			if (n == 1)
				break;
			if (n % i == 0) {
				prime[++cnt] = i;
				count[cnt] = 0;
				while (n % i == 0) {
					n /= i;
					count[cnt]++;
				}
			}
		}
		if (n != 1) {
			prime[++cnt] = n;
			count[cnt] = 1;
		}
		count[0]=cnt;
	}
	void dfs2(int now, long value, long euler) {//k
		if (now == count2[0] + 1) {		
			ans += (euler*pow(color, k / value))%mod;
			ans %= mod;
			return;
		}
		dfs2(now + 1, value, euler);
		long temp = prime2[now];
		for (int i = 1; i <= count2[now]; i++) {
			dfs2(now + 1, value * temp, euler * (prime2[now] - 1)
					* (temp / prime2[now])%mod);
			temp *= prime2[now];
			temp%=mod;
		}
	}

	void dfs1(int now, long value) {//b
		if (now == count1[0] + 1) {		
			res=(res+polya2(value))%mod;
			return;
		}
		dfs1(now + 1, value);
		long temp =prime1[now];
		for (int i = 1; i <= count1[now]; i++) {
			count2[now]-=2*i;
			dfs1(now + 1, value*temp);
			count2[now]+=2*i;
			temp *= prime1[now];
		}
	}
	long pow(long x, long p) {
		long ans = 1;
		while (p > 0) {
			if ((p & 1) == 1)
				ans = (ans * x) % mod;
			p >>= 1;
			x = (x * x) % mod;
		}
		return ans % mod;
	}
	long inv(long k) {// k在摸mod意义下的乘法逆元
		return pow(k%mod, mod - 2) % mod;
	}
	long solve() {
		if (a == 1)
			return c % mod;
		res=0;
		count1[0]=0;
		divide(a-1,prime1,count1);
		int temp=count1[0];
		divide(a+1,prime1,count1);
		if((a&1)==1)//合并a+1和a-1的质因子
		{
			 count1[1]+=count1[temp+1];
			 count1[temp+1]=count1[count1[0]];
			 prime1[temp+1]=prime1[count1[0]];
			 count1[0]--;
		}//初始时k的质因子等于b*b的质因子
		for(int i=0;i<=count1[0];i++){
			prime2[i]=prime1[i];
			count2[i]=count1[i];
		}//b*b的质因子只能取b的质因子的偶数倍次幂
		for(int i=1;i<=count1[0];i++)
			count1[i]/=2;
		dfs1(1,1);
		return res;
	}
	long polya1(long b) {// b*b
		long ans = 0;
		ans = pow(c, b * b);
		ans = (ans + pow(c, (b * b+1)/ 2 )) % mod;
		ans = (ans + 2 * pow(c, (b * b + 3) / 4)) % mod;
		return (ans * inv(4)) % mod;
	}
	long polya2(long b) {// k
		color = polya1(b);
		ans = 0;
		k=(a+1)*(a-1)/(b*b);
		dfs2(1, 1, 1);
		ans=(ans*c)%mod;
		return (ans * inv(k)) % mod;
	}
	void run() {
		Scanner scan = new Scanner(System.in);
		int cas = scan.nextInt();
		for (int i = 1; i <= cas; i++) {
			a = scan.nextLong();
			c = scan.nextLong();
			System.out.println("Case" + " " + i + ": " + solve());
		}
	}
	public static void main(String[] args) {
		new Rotation3441().run();
	}
}