总的来说，这题要2次用到polya定理。
由题目条件A*A=B*B+1，变形为(A-1)*(A+1)=K*B*B;
分别分解A-1和A+1的质因数，在合并在一起。
第一步：搜索B，对B*B的正方形涂色，基本的polya定理搞定，即C^(B*B)+C^((B*B+1)/2)+2*C^((B*B+3)/4).
第二步：搜索K，在A-1和A+1的因子中搜索，这样不会超时，在用polya定理，最后在结果上乘C就可以了……
1 #include<iostream>
2 #include<stdio.h>
3 #include<algorithm>
4 #include<iomanip>
5 #include<cmath>
6 #include<string>
7 #include<vector>
8 using namespace std;
9 const long mod=1000000007;
10 int prime[40003],m,fac[1000][2],flen;
11 vector<int> v;
12 bool f[40003];
13 __int64 A,C,ans,ret_A,inverse_4,inverse_k;
14 void init()
15 {
16 __int64 i,j;
17 m=0;
18 for(i=2;i<=40000;i++)
19 if(f[i]==0)
20 {
21 prime[m++]=i;
22 for(j=i*i;j<=40000;j+=i)
23 f[j]=1;
24 }
25 }
26 __int64 extend_gcd(__int64 a,__int64 b,__int64 &x,__int64 &y)
27 {
28 __int64 d;
29 if(b==0)
30 {
31 x=1;y=0;
32 return a;
33 }
34 else
35 {
36 d=extend_gcd(b,a%b,x,y);
37 __int64 temp=x;
38 x=y;
39 y=temp-a/b*y;
40 }
41 return d;
42 }
43 __int64 pows(__int64 a,__int64 b)
44 {
45 __int64 ans=1;
46 a%=mod;
47 while(b)
48 {
49 if(b&1) ans=(ans*a)%mod;
50 b>>=1;
51 a=(a*a)%mod;
52 }
53 return ans%mod;
54 }
55 void factor(__int64 n)
56 {
57 int i;
58 for(i=0;i<m&&prime[i]*prime[i]<=n;i++)
59 while(n%prime[i]==0)
60 {
61 v.push_back(prime[i]);
62 n/=prime[i];
63 }
64 if(n>1) v.push_back(n);
65 }
66 void combine()
67 {
68 sort(v.begin(),v.end());
69 int i,j;
70 flen=0;
71 fac[flen][0]=v[0];
72 fac[flen++][1]=1;
73 for(i=1;i<v.size();i++)
74 {
75 if(v[i]==fac[flen-1][0])
76 fac[flen-1][1]++;
77 else
78 {
79 fac[flen][0]=v[i];
80 fac[flen++][1]=1;
81 }
82 }
83 }
84 __int64 inverse(__int64 a)
85 {
86 __int64 x,y;
87 extend_gcd(a,mod,x,y);
88 return (x%mod+mod)%mod;
89 }
90 __int64 euler(__int64 n)
91 {
92 int i;
93 __int64 ans=1;
94 for(i=0;i<flen&&fac[i][0]*fac[i][0]<=n;i++)
95 {
96 if(n%fac[i][0]==0)
97 {
98 ans*=fac[i][0]-1;
99 n/=fac[i][0];
100 while(n%fac[i][0]==0)
101 {
102 ans*=fac[i][0];
103 n/=fac[i][0];
104 }
105 }
106 }
107 if(n>1) ans*=n-1;
108 return ans%mod;
109 }
110 void dfs(__int64 d,__int64 num,__int64 cnt_B,__int64 k)
111 {
112 if(d>=flen)
113 {
114 ret_A=(ret_A+pows(cnt_B,k/num)*euler(num)%mod)%mod;
115 return ;
116 }
117 for(int i=0;i<=fac[d][1];i++)
118 {
119 dfs(d+1,num,cnt_B,k);
120 num*=fac[d][0];
121 }
122 }
123 __int64 get_A(__int64 k,__int64 cnt_B)
124 {
125 ret_A=0;
126 dfs(0,1,cnt_B,k);
127 return ((ret_A*inverse_k)%mod)*C%mod;
128 }
129 __int64 get_B(__int64 n)
130 {
131 __int64 ans=pows(C,n*n);
132 ans=(ans+pows(C,(n*n+1)/2)%mod)%mod;
133 ans=(ans+2*pows(C,(n*n+3)/4)%mod)%mod;
134 return (ans*inverse_4)%mod;
135 }
136 //枚举B
137 //d表示深度，num表示枚举因子的大小
138 void dfsB(__int64 d,__int64 num)
139 {
140 if(d>=flen)
141 {
142 __int64 cnt_B=get_B(num);
143 __int64 k=(A*A-1)/num/num;
144 inverse_k=inverse(k);
145 ans=(ans+get_A(k,cnt_B))%mod;
146 return ;
147 }
148 int temp=fac[d][1];
149 for(int i=0;i<=temp;i+=2,fac[d][1]-=2)
150 {
151 dfsB(d+1,num);
152 num*=fac[d][0];
153 }
154 fac[d][1]=temp;
155 }
156 __int64 solve()
157 {
158 v.clear();
159 factor(A-1);
160 factor(A+1);
161 combine();
162 ans=0;
163 dfsB(0,1);
164 return ans;
165 }
166 int main()
167 {
168 init();
169 int t,k=0;
170 inverse_4=inverse(4);
171 cin>>t;
172 while(t--)
173 {
174 scanf("%I64d%I64d",&A,&C);
175 if(A==1)
176 printf("Case %d: %I64d
",++k,C);
177 else printf("Case %d: %I64d
",++k,solve());
178 }
179 return 0;
180 }