HDU 3441 Rotation （两次Polya!!!!）

题目：有一个A*A的正方形，拆成A*A个1*1的小正方形，然后组成k个B*B的正方形，而且剩下一个小正方形，也就是A*A=K*B*B+1。中间小小正方形连到K个B*B正方形的形状有多少种，有C种颜色，而且旋转视为等价。
 拿下AC的身体，开心吖，不过跪舔 了好久。
思路： 找到B，然后对B*B的正方形进行染色，然后将每一种染色方案视为一种颜色。将K个B*B的正方形看成K个物品的环，用之前得到的颜色进行染色。中间的小方块有C种颜色可选 。
那么第一步：找到B，直接分解肯定不行，A*A-1达到10^18左右，可以分解成(A-1)(A+1),分别找到因子，然后再合并在一起，然后搜索所有的因子，得到可能的B。
第二步：对于得到的B，进行B*B的正方形的染色，这个很基础了，4种旋转，C^(B*B)+2*C^((B*B+3)/4)+C^((B*B+1)/2);
然后再乘以4的逆元，若这步的方案是Cnt_B
第三步：那么剩下的相当于对K个物品的环进行染色，颜色数量为Cnt_B,也是基础的Polya,但是有一点，当B比较小，那么K就会非常大，可能达到10^18的级别，那么即使用欧拉函数优化，sqrt(K)也接受不了，之前我们已经找过A*A-1的因子，一部分给了B*B，剩下的为K的，那么直接对剩下的质因子搜索即可。所以在第一步搜索的时候把B*B的因子去掉，不过记得还原现场。这一步的实现就是进行第二次搜索，第二次Polya。
注意：纠结了好久，由于范围很大，时刻注意溢出问题，级别大部分都是10^18，记得取模后再乘
在欧拉函数的时候，用素数表还是会TLE，同样是利用之前分解因子得到的因子表。
[cpp]
#include<iostream> 
#include<cstdio> 
#include<vector> 
#include<cstring> 
#include<cmath> 
#include<algorithm> 
#define inf 1<<29 
#define LL long long 
#define N 1000000000 
#define MOD 1000000007 
#define pb(a) push_back(a) 
using namespace std; 
int tot; 
LL A,c; 
LL inverse_4,inverse_k; 
int prime[40000],cnt=0; 
LL fac[1000][2]; 
bool flag[40000]={0}; 
vector<int>v; 
void Prime(){ 
    for(int i=2;i<=sqrt(N+1.0);i++){ 
        if(flag[i])  continue; 
        prime[cnt++]=i; 
        for(int j=2;j*i<=sqrt(N+1.0);j++) 
            flag[i*j]=true; 
    } 
} 
//以上素数表 
LL extend_gcd(LL a,LL b,LL &x,LL &y){ 
    if(b==0){x=1;y=0;return a;} 
    LL gcd=extend_gcd(b,a%b,x,y); 
    LL t=x;x=y; 
    y=t-a/b*y; 
    return gcd; 
} 
LL Get_inverse(LL num){ 
    LL x,y; 
    LL gcd=extend_gcd(num,MOD,x,y); 
    return (x%MOD+MOD)%MOD; 
} 
//以上求逆元 
LL Eular(LL n){ 
    LL ret=1; 
   for(int i=0;i<tot&&fac[i][0]*fac[i][0]<=n;i++){ 
        if(n%fac[i][0]==0){ 
            n/=fac[i][0];ret*=fac[i][0]-1; 
            while(n%fac[i][0]==0){n/=fac[i][0];ret*=fac[i][0];} 
        } 
    } 
    if(n>1) ret*=n-1; 
    return (ret%MOD); 
} 
//以上求欧拉函数，注意要直接用之前分解的因子，不然会TLE 
LL PowMod(LL a,LL b){ 
    LL ret=1; 
    a%=MOD; 
    while(b){ 
        if(b&1) ret=((LL)ret*a)%MOD; 
        a=((LL)a*a)%MOD; 
        b>>=1; 
    } 
    return ret; 
} 
//快速幂 
void get_fact(LL t){ 
    for(int i=0;i<cnt&&prime[i]*prime[i]<=t;i++){ 
        while(t%prime[i]==0){ 
            t/=prime[i]; 
            v.pb(prime[i]); 
        } 
    } 
    if(t>1) v.pb(t); 
} 
//分解因子 
void get_union(){ 
    sort(v.begin(),v.end()); 
    tot=0; 
    fac[tot][0]=v[0];fac[tot++][1]=1; 
    for(int i=1;i<v.size();i++){ 
        if(v[i]==fac[tot-1][0]) 
            fac[tot-1][1]++; 
        else{ 
            fac[tot][0]=v[i]; 
            fac[tot++][1]=1; 
        } 
    } 
} 
//将A-1和A+1的因子整合在一起 
LL ret_A; 
void dfs(int idx,LL num,LL cnt_B,LL K){ 
    if(idx>=tot){ 
        ret_A=(ret_A+PowMod(cnt_B,K/num)*Eular(num)%MOD)%MOD; 
        return ; 
    } 
    for(int i=0;i<=fac[idx][1];i++){ 
        dfs(idx+1,num,cnt_B,K); 
        num*=fac[idx][0]; 
    } 
} 
//搜索K的因子，欧拉函数优化，第二个Polya 
LL get_A(LL K,LL cnt_B){ 
    ret_A=0; 
    dfs(0,1,cnt_B,K); 
    return (((ret_A*inverse_k)%MOD)*c)%MOD; 
} 
//用B*B的数量给K个环染色 
LL get_B(LL B){ 
    LL ans=PowMod(c,(LL)B*B); 
    ans=(ans+2*PowMod(c,((LL)B*B+3)/4))%MOD; 
    ans=(ans+PowMod(c,((LL)B*B+1)/2))%MOD; 
    return (ans*inverse_4)%MOD; 
} 
//B*B的正方形染色,4种旋转 
LL ans; 
void dfsB(int idx ,LL nowB){ 
    if(idx>=tot){ 
        LL cnt_B=get_B(nowB); 
        LL K=(A*A-1)/nowB/nowB; 
        inverse_k=Get_inverse(K); 
        ans=(ans+get_A(K,cnt_B))%MOD; 
        return ; 
    } 
    LL temp=fac[idx][1]; 
    //因子每次减少2，因为是B*B，而且剩余的用作搜索K的因子 
    for(int i=0;i<=temp;i+=2,fac[idx][1]-=2){ 
        dfsB(idx+1,nowB); 
        nowB*=fac[idx][0]; 
    } 
    fac[idx][1]=temp; 
} 
//以上搜索B 
LL slove(){ 
    v.clear(); 
    get_fact(A-1); 
    get_fact(A+1); 
    get_union(); 
    ans=0; 
    dfsB(0,1); 
    return ans; 
} 
int main(){ 
    int t,cas=0; 
    Prime(); 
    inverse_4=Get_inverse(4); 
    scanf("%d",&t); 
    while(t--){ 
        scanf("%I64d%I64d",&A,&c); 
        printf("Case %d: ",++cas); 
        if(A==1)  printf("%I64d/n",c); 
        else 
            printf("%I64d/n",slove()); 
    } 
    return 0; 
}