[SDOI2013]项链

description

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

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

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

solution

先求本质不同的珠子数(m),发现要用(Burnside)。
对于一个置换考虑计算其循环长度为(k)是对应的不动点个数。
可以知道其为(sum_{x_1}^a...sum_{x_k}^a[gcd_{i=1}^kx_i=1])。
莫比乌斯反演即得(sum_{i=1}^amu(i)lfloorfrac{a}{i} floor^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(0le r<p^2)),
如果(p|Ans)那么显然有(p|r),直接除掉,发现(frac{Ans}{p}=kp+frac{r}{p}equivfrac{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
",mul(ans,inv(n))%p);
    else ans/=p,printf("%lld
",mul(ans,inv(n/p))%p);
  }
  return 0;
}