[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;
}