[Polya计数] BZOJ 3202 [Sdoi2013]项链

大数论题
可以分成三部分解决
求合法的珠子的方案数 这个可以用反演解决
求本质不同的项链 这个直接Polya
要求相邻不能相同 这个可以递推出每个长度的环的方案数 但是不够快 可以直接求出通项
这里有篇超具体的题解 我就懒得一步步具体写了 数学公式太多绝对是大工程
最后有个小trick $n>P$ 那么$n$可能是$P$的倍数 那么就没有逆元 那么我们先模$P^2$ 如果$n$是$P$的倍数 那么当前答案一定能够整除$P$ 先除$P$然后在乘上$n\over P$的逆元

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const ll _P=1000000007LL,P=1000000007LL*1000000007LL;

inline ll mul(ll x,ll y){
  return (x*y-(ll)(((long double)x*y+0.5)/(long double)P)*P+P)%P;
}
inline ll Pow(ll a,ll b){
  ll ret=1; for (;b;b>>=1,a=mul(a,a)) if (b&1) ret=mul(ret,a); return ret;
}
inline ll Pow(ll a,ll b,ll P){
  ll ret=1; for (;b;b>>=1,a=a*a%P) if (b&1) ret=ret*a%P; return ret;
}
inline ll sqr(ll x) { return mul(x,x); }
inline ll cub(ll x) { return mul(x,mul(x,x)); }
inline void add(ll &x,ll y){ x+=y; if (x>=P) x-=P; }

const int mN=10000005,mP=mN/10;
bool vst[mN]; int miu[mN];
int prime[mP],num;
inline void Pre(int n){
  miu[1]=1;
  for (int i=2;i<=n;i++){
    if (!vst[i]) prime[++num]=i,miu[i]=-1;
    for (int j=1;j<=num && (ll)prime[j]*i<=n;j++){
      vst[prime[j]*i]=1;
      if (i%prime[j]==0){
    miu[i*prime[j]]=0; break;
      }else
    miu[i*prime[j]]=-miu[i];
    }
  }
  for (int i=1;i<=n;i++) miu[i]+=miu[i-1];
}
const ll Inv6=833333345000000041LL;
inline ll calc(ll n){
  ll A=0,B=0,C=1;
  for (int i=1,j;i<=n;i=j+1){
    j=n/(n/i);
    add(A,mul((miu[j]+P-miu[i-1])%P,cub(n/i)));
    add(B,mul((miu[j]+P-miu[i-1])%P,sqr(n/i)));
  }
  add(A,3*B%P); add(A,2); return mul(A,Inv6);
}

ll n,m,_m;

inline ll F(ll n){
  ll ret=Pow(_m,n);
  if (n&1) add(ret,P-_m); else add(ret,_m);
  return ret;
}

ll Ans=0;
ll cnt,p[mP],q[mP];

inline void dfs(int t,ll cur,ll phi){
  if (t==cnt+1){
    add(Ans,mul(phi,F(n/cur)));
    return;
  }
  dfs(t+1,cur,phi);
  cur=mul(cur,p[t]); phi=mul(phi,p[t]-1);
  dfs(t+1,cur,phi);
  for (int i=2;i<=q[t];i++){
    cur=mul(cur,p[t]); phi=mul(phi,p[t]);
    dfs(t+1,cur,phi);
  }
}

inline void Solve(){
  ll t=n; cnt=0;
  for (int i=1;i<=num && (ll)prime[i]*prime[i]<=n && t>1;i++)
    if (t%prime[i]==0){
      p[++cnt]=prime[i]; q[cnt]=0;
      while (t%p[cnt]==0) t/=p[cnt],q[cnt]++;
    }
  if (t!=1) p[++cnt]=t,q[cnt]=1;
  Ans=0; dfs(1,1,1);
  if (n%_P==0){
    Ans/=_P;
    Ans=Ans*Pow(n/_P,_P-2,_P)%_P;
  }else{
    Ans%=_P;
    Ans=Ans*Pow(n%_P,_P-2,_P)%_P;
  }
}

int main(){
  int Q,a;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  scanf("%d",&Q); Pre(1e7);
  while (Q--){
    scanf("%lld%d",&n,&a);
    m=calc(a); _m=(m+P-1)%P; 
    Solve();
    printf("%lld\n",Ans);
  }
  return 0;
}