【hdu3441】Rotation，Polya套Polya

传送门
思路：
首先，
这题我TM做了一天半
最后发现问题是没有特判A=1

说一下题意
给定$A,C$,在$A×A$个格子中，每个格子都可以染成$C$种颜色中任意一种，然后把这些格子重新组合成$n$个$B×B$的矩形和一个剩下的格子（即$n×B^2+1=A^2$），把这k个矩形连在那个小格子上，求方案数（旋转后重合视为同一种方案）
一开始我直接套Polya，发现有$4^n×n$种置换，然后就写成了这个样子

$$\sum_{n^2|A^2-1}[B=\frac {A^2-1}{n^2}]\sum^n_{i=1}\sum^3_{a_1=0}\sum^3_{a_2=0}..\sum^3_{a_{(i,n)}=0}C^{\sum^{(i,n)}_{j=1}f(a_j,\frac n {(i,n)})}$$
$f(x,y)=4^{y-1}×\lceil \frac {B^2} {p[x]} \rceil$
$p[]=\{1,4,2,4\}$
按理来说我们不可以愉快地直接求了
因为这个复杂度连 $A=4$都跑不过去

其实标题早已暗示了一切
我们先对 $B×B$的矩形套Polya
四种置换的循环节数分别是 $B^2,\lceil \frac {B^2} 4 \rceil,\lceil \frac {B^2} 2 \rceil,\lceil \frac {B^2} 4 \rceil$
颜色数就是 $C$了
我们定义一个新的$C'$是这个Polya算出来的答案
就是这个东西
$$C'=\frac {C^{B^2}+C^{\lceil \frac {B^2} 2\rceil}+C^{\lceil \frac {B^2} 2\rceil}+C^{\lceil \frac {B^2} 4\rceil}} 4$$
因而这 $k$个$B×B$的矩形又变成了新的染色对象，颜色数为 $C'$，置换方式是旋转1到k格，最终答案就是这个
$$\sum_{B^2|A^2-1}[n=\frac {A^2-1} {B^2}]\sum^n_{i=1}[i|n]\varphi(\frac n i)C'^i$$
按理来说我们就可以愉快地直接求了吧
但是 $A\leq 10^9$，总不能枚举 $\sqrt {A^2-1}$以内的因数吧？

我想的是分解质因数
$A^2-1=(A-1)(A+1)$
所以只要考虑 $A-1,A+1$就可以了
而且 $B^2$的限制使得符合条件的数中质因子次数都是偶数
小于 $10^9$的数不同的质因子最多有九个，所以考虑暴搜来求符合条件的数
如果 $A^2-1=\prod a_i^{p_i}$
复杂度是 $O(\prod p_i)$
因而只要把 $A-1,A+1$分别分解质因数，然后合到一起去就可以了
按理来说
但是n的数量级也是 $10^{18}$，显然 $F(n)=\sum^n_{i=1}[i|n]\varphi(\frac n i)C'^i$这个非积性函数连杜教筛都做不了

我只能说：继续暴搜吧
反正因数又不大，想不出什么好方法就直接dfs呗
第一遍时算出 $n$的质因子及次数，第二遍dfs时记录$i$，算出 $\frac n i$的质因子及次数，然后就可以求 $\varphi$了
$\varphi$的数量级也是 $10^{18}$,我没有算它的逆元啥的，所以乘的过程中好几次炸了long long
因为要算 $A-1$的质因子，所以要特判 $A=1$的情况，不然会T死(各位也可能写的做法比较优美，不像我这么思博)
我还在那里一直压常数，找网上的程序下来拍，发现随机 $10^9$数据下我跑的比它还快……

至此这个问题终于算是解决了
最后贴张图感受一下

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#define LL long long
#define mo 1000000007
using namespace std;
LL A,B,col;
int tot,T,C,ans,cnt;
int prime[3505],num[3500][2];
bool vis[100005];
struct node{
    int num[3500],pr;
    void clr()
    {
        memset(num,0,sizeof(num));
        pr=0;
    }
    void cal(LL x)
    {
        clr();
        for (int i=1;i<=prime[0];++i)
            while (x%prime[i]==0)
                x/=prime[i],
                ++num[i];
        if (x!=1) pr=x;
    }
}X,Y;
LL qr(LL x,LL y)
{
    LL t=1;
    for (x%=mo;y;y>>=1,x=x*x%mo)
        if (y&1) t=t*x%mo;
    return t;
}
LL ph()
{
    LL t=1;
    for (int i=1;i<=cnt;++i)
        if (num[i][1])
        {
            t=t*(num[i][0]-1);
            for (int j=1;j<num[i][1];++j)
                t=t*num[i][0];
        }
    return t;
}
void dfs2(int x,LL sum)
{
    if (x>cnt)
    {
        ans=(ans+qr(col,sum)*(ph()%mo)%mo)%mo;
        return;
    }
    LL t=1;
    int k=num[x][1];
    for (;num[x][1]>=0;--num[x][1],t=t*num[x][0])
        dfs2(x+1,sum*t);
    num[x][1]=k;
}
void dfs(int x,LL n)
{
    if (x>cnt)
    {
        B=sqrt(A/n);
        ans=0;
        bool flag=B&1;
        col=((qr(C,1LL*B*B)+qr(C,1LL*B*B/2+flag)+2*qr(C,1LL*B*B/4+flag))%mo)*qr(4,mo-2)%mo;
        dfs2(1,1);
        tot=(tot+ans*qr(n,mo-2)%mo)%mo;
        return;
    }
    LL t=1;
    int k=num[x][1];
    for (;num[x][1]>=0;num[x][1]-=2,t=t*num[x][0]*num[x][0])
        dfs(x+1,n/t);
    num[x][1]=k;
}
void work()
{
    scanf("%lld%d",&A,&C);
    if (A==1) {printf("%d\n",C);return;}
    X.cal(A-1);Y.cal(A+1);
    tot=cnt=0;
    for (int i=1;i<=prime[0];++i)
        if (X.num[i]+Y.num[i])
            num[++cnt][0]=prime[i],
            num[cnt][1]=X.num[i]+Y.num[i];
    if (X.pr) num[++cnt][0]=X.pr,num[cnt][1]=1;
    if (Y.pr) num[++cnt][0]=Y.pr,num[cnt][1]=1;
    if (X.pr==Y.pr&&X.pr) --cnt,num[cnt][1]=2;
    A=1LL*A*A-1;
    dfs(1,A);
    printf("%d\n",1LL*tot*C%mo);
}
main()
{
    for (int i=2;i<=31625;++i)
    {
        if (!vis[i])
            prime[++prime[0]]=i;
        for (int j=1;prime[j]*i<=31625&&j<=prime[0];++j)
        {
            vis[prime[j]*i]=1;
            if (i%prime[j]);
            else
                break;
        }
    }
    scanf("%d",&T);
    for (int i=1;i<=T;++i)
        printf("Case %d: ",i),
        work();
}