[SDOI2013]项链

description

luogu
最近,铭铭迷恋上了一种项链。与其他珍珠项链基本上相同,不过这种项链的珠子却与众不同,是正三菱柱的泰山石雕刻而成的。
三菱柱的侧面是正方形构成的,上面刻有数字。 能够让铭铭满意的项链必须满足下面的条件:

• 这串项链由\(n\)颗珠子构成的。
• 每一个珠子上面的数字\(x\),必须满足\(0<x\le a\),且珠子上面的数字的最大公约数要恰 好为\(1\)。两个珠子被认为是相同的,当且仅当他们经过旋转,或者翻转后能够变成一样的。
• 相邻的两个珠子必须不同。
• 两串项链如果能够经过旋转变成一样的,那么这两串项链就是相同的。

铭铭很好奇如果给定\(n\)和\(a\),能够找到多少不同串项链。
由于答案可能很大,所以对输出的答案\(mod\ 10^9+7\)。
多组数据,\(T\le 10,n\le 10^{14},a\le 10^7\),时限\(3s\)。

solution

先求本质不同的珠子数\(m\),发现要用\(Burnside\)。
对于一个置换考虑计算其循环长度为\(k\)是对应的不动点个数。
可以知道其为\(\sum_{x_1}^a...\sum_{x_k}^a[gcd_{i=1}^kx_i=1]\)。
莫比乌斯反演即得\(\sum_{i=1}^a\mu(i)\lfloor\frac{a}{i}\rfloor^k\)。
直接\(O(n)\)应该也能满足要求,可以用数论分块做到\(O(\sqrt n)\)卡常。
再求项链本质不同的方案数。
根据[POJ2888]Magic Bracelet中的老套路转化为\(\frac{1}{n}\sum_{i=1}^nf(gcd(i,n))=\frac{1}{n}\sum_{d|n}f(d)\varphi(\frac{n}{d})\),
其中\(f(n)\)表示长度为\(n\)的手环不考虑旋转时的合法方案数。
发现\(f(1)=0,f(2)=1,f(n)=(m-2)f(n-1)+(m-1)f(n-2)\)。
使用生成函数递推,设\(G(x)=\sum_{i=1}^{\infty}f(i)x^i\),
那么根据递推式有\(G(x)-f(1)x-f(2)x^2=(m-2)x(G(x)-f(1)x)+(m-1)x^2G(x)\)
化简/使用待定系数法解得\(G(x)=\frac{m-1}{1-(m-1)x}-\frac{m-1}{1+x}=\sum_{i=1}^{\infty}(m-1)^i-(m-1)(-1)^i\)。
于是\(O(\sqrt n)\)地\(dfs\)所有约数并递推算\(\varphi\)。
最后需要注意\(n\)可能是\(p\)的倍数,然后就gg了。
解决方法是将原来的\(mod\ p=10^9+7\)变成\(mod\ p^2\),这样求出来的\(Ans=kp^2+r(0\le r<p^2)\),
如果\(p|Ans\)那么显然有\(p|r\),直接除掉,发现\(\frac{Ans}{p}=kp+\frac{r}{p}\equiv\frac{r}{p}(mod\ p)\),仍然是对的。
因为\(n<p^2\),所以\(n\)里面最多只含一个\(p\)为因子,最后乘一下\(\frac{n}{p}\)对于\(p^2\)的逆元即可。

code

#include<bits/stdc++.h>
#define FL "a"
using namespace std;
typedef long long ll;
typedef long double dd;
const int N=1e7+10;
const int p=1e9+7;
const ll mod=1ll*p*p;
const int inv6=166666668;
inline ll read(){
  ll data=0,w=1;char ch=getchar();
  while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
  if(ch=='-')w=-1,ch=getchar();
  while(ch<='9'&&ch>='0')data=data*10+ch-48,ch=getchar();
  return data*w;
}
inline void file(){
  freopen(FL".in","r",stdin);
  freopen(FL".out","w",stdout);
}
inline void upd(ll &a,ll b){a+=b;if(a>=mod)a-=mod;}
inline void dec(ll &a,ll b){a-=b;if(a<0)a+=mod;}
inline ll mul(ll a,ll b){ll k=(dd)a/mod*b;return a*b-k*mod;}
inline ll poww(ll a,ll b){
  ll res=1;
  for(;b;b>>=1,a=mul(a,a))
    if(b&1)res=mul(res,a);
  return res;
}
ll n,a,m,m1,m2,mu[N],ans;
inline void init(){
  static int pri[N],cnt;static bool vis[N];
  memset(vis,0,sizeof(vis));cnt=0;vis[1]=mu[1]=1;
  for(int i=2;i<N;i++){
    if(!vis[i])pri[++cnt]=i,mu[i]=mod-1;
    for(int j=1;j<=cnt&&1ll*i*pri[j]<N;j++){
      vis[i*pri[j]]=1;if(mu[i])mu[i*pri[j]]=mod-mu[i];
      if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
    }
  }
  for(int i=2;i<N;i++)upd(mu[i],mu[i-1]);  
}
inline ll calc(ll x){return (poww(m1,x)-mul(m1,poww(mod-1,x-1))+mod)%mod;}
void exgcd(ll a,ll b,ll &x,ll &y,ll &d){
  if(!b){d=a;x=1;y=0;return;}exgcd(b,a%b,y,x,d);y-=a/b*x;
}
inline ll inv(ll a){
  ll x,y,d;exgcd(a,mod,x,y,d);x=(x%mod+mod)%mod;return x;
}
ll pri[52];int sum[52],cnt;
inline void fact(ll n){
  cnt=0;
  for(int i=2;1ll*i*i<=n;i++)
    if(n%i==0){
      cnt++;pri[cnt]=i;sum[cnt]=0;
      while(n%i==0)n/=i,sum[cnt]++;
    }
  if(n!=1)cnt++,pri[cnt]=n,sum[cnt]=1;
}
void dfs(int x,ll res,ll phi){
  if(x==cnt+1){upd(ans,mul(calc(n/res),phi));return;}
  for(int i=0;i<=sum[x];i++){
    if(i==1)phi=mul(pri[x]-1,phi);
    else if(i)phi=mul(pri[x],phi);
    if(i)res=mul(res,pri[x]);dfs(x+1,res,phi);
  }
}
int main()
{
  int T=read();init();
  while(T--){
    n=read();a=read();m=2;ans=0;
    for(ll i=1,j;i<=a;i=j+1){
      j=a/(a/i);
      upd(m,mul(mul(3,(mu[j]-mu[i-1]+mod)%mod),poww((a/i),2)));
      upd(m,mul((mu[j]-mu[i-1]+mod)%mod,poww((a/i),3)));
    }
    m=m1=m2=mul(m,inv(6));dec(m1,1);dec(m2,2);
    fact(n);dfs(1,1,1);
    if(n%p)printf("%lld\n",mul(ans,inv(n))%p);
    else ans/=p,printf("%lld\n",mul(ans,inv(n/p))%p);
  }
  return 0;
}