题意：TTT 组数据，每组给定 n,an,an,a，求满足下列条件的项链数量：
有 nnn 个珠子。每个珠子上有三个 [1,a]∩Z[1,a]\cap \Z[1,a]∩Z 的数，且三个数 gcd\gcdgcd 为 111。相邻两个珠子不同。珠子旋转、翻转同构，项链旋转同构。对 109+710^9+7109+7 取模。
n≤1014,a≤107n\leq 10^{14},a\leq 10^7n≤1014,a≤107
有毒……
显然是道缝合怪题，先考虑求有多少个不同的珠子。
设 f(k)f(k)f(k) 表示 kkk 个 [1,a][1,a][1,a] 的数 gcd\gcdgcd 为 111 的方案数，用 Polya 定理数一下可知
ans=f(3)+3f(2)+2f(1)6ans=\frac{f(3)+3f(2)+2f(1)}{6}ans=6f(3)+3f(2)+2f(1)
然后就是个简单的反演
f(k)=∑i=1aμ(i)⌊ai⌋kf(k)=\sum_{i=1}^a\mu (i)\left\lfloor\frac ai\right\rfloor^kf(k)=i=1∑aμ(i)⌊ia⌋k
整除分块算就可以了（其实直接暴力也可以？）。
然后第二步，数项链。设第一步得到的珠子种数为 kkk，众所周知，F(n)F(n)F(n) 为 nnn 个珠子不循环同构的方案，答案就是 FFF 和 φ\varphiφ 的狄利克雷卷积除以 nnn。
考虑 FFF 怎么算。
考虑一个 nnn 个珠子的合法方案，我们删掉编号为 111 的点。如果两边的点不同，就对应了一种 F(n−1)F(n-1)F(n−1) 的方案，乘上 111 号珠子的方案 (k−2)(k-2)(k−2)。如果相同，任意删掉一个后就对应一种 F(n−2)F(n-2)F(n−2) 的方案，乘上 111 号珠子的方案 (k−1)(k-1)(k−1)。所以
F(n)=(k−2)F(n−1)+(k−1)F(n−2)F(n)=(k-2)F(n-1)+(k-1)F(n-2)F(n)=(k−2)F(n−1)+(k−1)F(n−2)
注意边界是 F(1)=0,F(2)=k(k−1)F(1)=0,F(2)=k(k-1)F(1)=0,F(2)=k(k−1),倒退出来 F(0)=mF(0)=mF(0)=m。
列出生成函数：
F=(k−2)xF+(k−1)x2F+k−k(k−2)xF=(k-2)xF+(k-1)x^2F+k-k(k-2)xF=(k−2)xF+(k−1)x2F+k−k(k−2)x
这里减 k(k−2)xk(k-2)xk(k−2)x 是为了平衡 (k−2)xF(k-2)xF(k−2)xF 对一次项的影响。
即
F=−m(m−2)x+m(1+x)(1−(m−1)x)F=\frac{-m(m-2)x+m}{(1+x)(1-(m-1)x)}F=(1+x)(1−(m−1)x)−m(m−2)x+m
直接待定系数拆掉
F=a1+x+b1−(k−1)xF=\frac{a}{1+x}+\frac{b}{1-(k-1)x}F=1+xa+1−(k−1)xb
列出方程：
a−a(k−1)x+b+bx=−k(k−2)x+ka-a(k-1)x+b+bx=-k(k-2)x+ka−a(k−1)x+b+bx=−k(k−2)x+k
解得 a=k−1,b=1a=k-1,b=1a=k−1,b=1
所以
F=k−11+x+11−(k−1)xF=\frac{k-1}{1+x}+\frac{1}{1-(k-1)x}F=1+xk−1+1−(k−1)x1
F(n)=(k−1)n+(−1)n(k−1)F(n)=(k-1)^n+(-1)^n(k-1)F(n)=(k−1)n+(−1)n(k−1)
然后就可以算了。
这里有点卡，就不用 Polya 模板那个垃圾做法，而是把因子预处理出来 dfs，可以做到 O(d(n))\Omicron(d(n))O(d(n)) 的优秀复杂度。
然后因为 nnn 很大，最后除 nnn 的时候又会有喜闻乐见的除 000 问题。所以要对 (109+7)2(10^9+7)^2(109+7)2 取模，最后再判一下。
总复杂度 O(n+T(a+d(n)log))\Omicron(n+T(\sqrt a+d(n)\log))O(n+T(a+d(n)log))
#include <iostream>#include <cstdio>#include <cstring>#include <cctype>#define int long longusing namespace std;const int mod=1e9+7,MOD=mod*mod;inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}inline int mul(int a,int b){return (a*b-(int)((long double)a/MOD*b+0.5)*MOD+MOD)%MOD;}inline int qpow(int a,int p){int ans=1;while (p){if (p&1) ans=mul(ans,a);a=mul(a,a),p>>=1;}return ans;}int INV6=833333345000000041ll;int n,a;const int N=1e7,MAXN=N+5;bool np[MAXN];signed pl[MAXN],cnt,mu[MAXN];void init(){np[1]=1,mu[1]=1;for (int i=2;i<=N;i++){if (!np[i]) mu[pl[++cnt]=i]=-1;for (int j=1,x;(x=i*pl[j])<=N;j++){np[x]=1;if (i%pl[j]==0) break;mu[x]=-mu[i];}}for (int i=2;i<=N;i++) mu[i]+=mu[i-1];}int calc(){int ans=0;for (int l=1,r;l<=a;l=r+1){r=a/(a/l);int t=mu[r]-mu[l-1];(t<0)&&(t+=MOD);ans=add(ans,mul(t,add(mul(a/l,mul(a/l,a/l)),3*mul(a/l,a/l)%MOD)));}return ans;}int k,ans;int lis[20005],siz[20005],tot;void dfs(int pos,int phi,int res){if (pos>tot){int t=qpow(k,res);if (res&1) t=dec(t,k);else t=add(t,k);ans=add(ans,mul(phi,t));return;}int t=qpow(lis[pos],siz[pos]);dfs(pos+1,phi,res*t);t/=lis[pos];dfs(pos+1,phi*=(lis[pos]-1),res*t);t/=lis[pos];if (!t) return;for (;t;dfs(pos+1,phi*=lis[pos],res*t),t/=lis[pos]);}signed main(){init();int T;cin>>T;while (T--){cin>>n>>a;k=dec(mul(calc()+2,INV6),1);tot=ans=0;int t=n;for (int i=1;i<=cnt;i++)if (t%pl[i]==0){lis[++tot]=pl[i],siz[tot]=0;while (t%pl[i]==0) t/=pl[i],++siz[tot];}if (t>1) lis[++tot]=t,siz[tot]=1;dfs(1,1,1);if (n%mod==0) n/=mod,ans/=mod;ans=mul(ans,qpow(n,mod-2));cout<<ans%mod<<'\n';}return 0;}
