分析
Attention！这道题的模数是(1e8+7)。
注意到循环同构会被认为是同一种方案，我们可以把顺时针旋转每个人的位置作为置换，容易发现这些置换一定会形成一个置换群，于是题目所求的所有合法的方案数便是这个置换群下等价类的数目，可以使用Burnside引理解决。
考虑如何求在置换“顺时针旋转(x)位置”下不动点的数目，可以发现这样的不动点的男女排列一定具有周期性，且(gcd(n,x))是它的一个周期，这里可以自行画图理解一下。这样我们又再次简化了题目，我们现在只需考虑这一个周期了。
令(d=gcd(n,x))，首先进行如下的分类讨论：
-
(n leq k)，这样不存在不合法的方案，返回(2^d)。
-
(n>k)且(d leq k)，这样只要这个周期内不全是女生即是合法方案，返回(2^d-1)。
-
(n>k)且(d>k)，DP解决。
怎么DP？
(f[i][j])表示不考虑循环的情况下，考虑到前(i)个人，最后(j)个人是女生的方案数。
(g[i][j])表示不考虑循环的情况下，考虑到前(i)个人，保证第一个人是男生，最后(j)个人是女生的方案数。
显然有转移：
(g[i][j])的转移与(f[i][j])类似。
于是对于长度为(i)的区间，其合法方案数为(sum_{j=1}^{min(i,k)}f[i][j]-sum_{j=k+1}^{min(i-1,2k)}g[i-j][0] imes (2k-j+1))，这个可以直接预处理出来。
总结一下上面所提到的，在置换“顺时针旋转(x)位置”下的不动点，男女排列一定具有周期性，(d=gcd(n,x))是其一个周期，我们只需要对这个区间进行讨论和计算，并统计到所有置换的不动点个数和中即可。
考虑到这样一个事实，(sum_{i=1}^{n}[gcd(i,n)=d]=varphi(frac{n}{d}))，所以可以直接枚举(n)的约数(d)，然后计算(d)对答案的贡献。（这一步可有可无）
最后根据Burnside引理乘上一个(n^{-1})。
代码
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,x) for(int i=head[(x)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
inline int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
return x;
}
const int MAXN=2005;
const int MOD=1e8+7;
int n,k;
LL f[MAXN][MAXN],g[MAXN][MAXN],cnt[MAXN];
int prm[MAXN],phi[MAXN],pcnt;
bool vis[MAXN];
inline LL qpow(LL x,LL y){
LL ret=1,tt=x%MOD;
while(y){
if(y&1) ret=ret*tt%MOD;
tt=tt*tt%MOD;
y>>=1;
}
return ret;
}
inline void pre_process1(){
phi[1]=1;
rin(i,2,n){
if(!vis[i]){
prm[++pcnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=pcnt&&i*prm[j]<=n;j++){
vis[i*prm[j]]=1;
if(i%prm[j]==0){
phi[i*prm[j]]=phi[i]*prm[j];
break;
}
phi[i*prm[j]]=phi[i]*(prm[j]-1);
}
}
}
inline void pre_process2(){
f[0][0]=f[1][0]=f[1][1]=g[1][0]=1;
rin(i,2,n){
f[i][0]=g[i][0]=0;
rin(j,0,std::min(i-1,k)){
f[i][0]+=f[i-1][j];
if(f[i][0]>=MOD) f[i][0]-=MOD;
}
rin(j,1,std::min(i,k)) f[i][j]=f[i-1][j-1];
rin(j,0,std::min(i-2,k)){
g[i][0]+=g[i-1][j];
if(g[i][0]>=MOD) g[i][0]-=MOD;
}
rin(j,1,std::min(i-1,k)) g[i][j]=g[i-1][j-1];
}
rin(i,1,n){
cnt[i]=0;
rin(j,0,std::min(i,k)){
cnt[i]+=f[i][j];
if(cnt[i]>=MOD) cnt[i]-=MOD;
}
rin(j,k+1,std::min(i-1,k<<1)){
cnt[i]-=g[i-j][0]*((k<<1)-j+1)%MOD;
if(cnt[i]<0) cnt[i]+=MOD;
}
}
}
inline LL solve(int x){
if(x<=k){
if(k<n) return (qpow(2,x)-1+MOD)%MOD;
else return qpow(2,x);
}
else return cnt[x];
}
int main(){
int T=read();
n=2000;
pre_process1();
while(T--){
n=read(),k=read();
pre_process2();
int lim=sqrt(n);
LL ans=0;
rin(i,1,lim){
if(n%i) continue;
ans=(ans+solve(i)*phi[n/i])%MOD;
if(i*i==n) continue;
ans=(ans+solve(n/i)*phi[i])%MOD;
}
printf("%lld
",ans*qpow(n,MOD-2)%MOD);
}
return 0;
}