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_易做图(LL a,LL b,LL &x,LL &y){
    if(b==0){x=1;y=0;return a;}
    LL 易做图=extend_易做图(b,a%b,x,y);
    LL t=x;x=y;
    y=t-a/b*y;
    return 易做图;
}
LL Get_inverse(LL num){
    LL x,y;
    LL 易做图=extend_易做图(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

补充：软件开发 , C++ ,