bzoj3202 SDOI2013项链

珍珠个数很好推不说了。本质不同项链个数用burnside得到\(\frac{1}{n}\sum\limits_{d|n}f(d,m)\varphi(\frac{n}{d})\)，\(f(n,m)\)表示长度为\(n\)，有\(m\)种本质不同的珍珠个数，相邻不同，首尾不同的项链个数。考虑容斥，用相邻不同，首尾任意的减去相邻不同，首尾相同的。任意一个相邻不同，首尾不同的相邻都可以考虑在开头加一个和末尾相同的珍珠得到一个相邻不同，首尾不同的项链，所以长度为\(n\)相邻不同首尾相同的就是\(f(n-1,m)\)。可得\(f(n,m)=m(m-1)^{n-1}-f(n-1,m)\)。边界为\(f(1,m)=0\)。那么有

\[f(n,m)=m\sum\limits_{i=1}^{n-1}(m-1)^i(-1)^{n-1-i}=m(-1)^{n-1}\sum\limits_{i=1}^{n-1}(1-m)^i=m(-1)^{n-1}\frac{(1-m)^n-(1-m)}{-m}=(m-1)^n+(m-1)(-1)^n \]

要求每个因数的\(\varphi(d)\)，这个东西可以考虑直接枚举\(n\)每个质因子的幂次求值。时间复杂度\(O(d(n)\log n)\)。

一个很恶心的东西就是\((1e9+7)|n\)。把模数变成\((1e9+7)^2\)去做，最后那个除以\(n\)先让他除以\(1e9+7\)，再除以\(\frac{n}{1e9+7}\)。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<queue>
#include<vector>
#include<map>
using namespace std;
typedef long long ll;
#define N 10000002
ll p;
int T,a;
int prim[N],cnt,mu[N],sm[N],ci[100002],tot;
ll b[100002],n,ans,m;
bool isp[N];
inline void upd(ll &x,ll y){x+=x+y>=p?y-p:y;}
inline ll gsc(ll d,ll k){ll ret=0;while(k){if(k&1)upd(ret,d);upd(d,d);k>>=1;}return ret;}
inline ll ksm(ll d,ll k){ll ret=1;while(k){if(k&1)ret=gsc(ret,d);d=gsc(d,d);k>>=1;}return ret;}
void dfs(int dep,ll fi,ll d)
{
    if(dep==tot+1)
    {ll te=ksm(m,n/d);
        if((n/d)&1)upd(te,(p-m)%p);
        else upd(te,m);
        te=gsc(te,fi);
        upd(ans,te);
        return;
    }
    dfs(dep+1,fi,d);
    fi*=(b[dep]-1);
    for(int i=1;i<=ci[dep];i++)
    {
        d*=b[dep];
        dfs(dep+1,fi,d);
        fi*=b[dep];
    }
}
int main()
{
    for(int i=2;i<=10000000;i++)isp[i]=1;mu[1]=1;
    for(int i=2;i<=10000000;i++)
    {
        if(isp[i])mu[i]=-1,prim[++cnt]=i;
        for(int j=1;j<=cnt&&prim[j]*i<=10000000;j++)
        {
            isp[prim[j]*i]=0;
            if(i%prim[j]==0){mu[prim[j]*i]=0;break;}
            mu[prim[j]*i]=-mu[i];
        }
    }
    for(int i=1;i<=10000000;i++)sm[i]=sm[i-1]+mu[i];
	scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%d",&n,&a);
        if(n%(1000000007)==0)p=1000000007ll*1000000007;
        else p=1e9+7;
        m=0;tot=0;ans=0;
        for(int i=1,j;i<=a;i=j+1)
        {
            j=a/(a/i);int d=a/i;
            ll te=(1ll*gsc(gsc(d,d),d)+3ll*gsc(d,d)%p+2ll*d%p)%p;
            upd(m,gsc((sm[j]-sm[i-1]+p)%p,te));
        }ll tn=n;
        if(n%1000000007==0)m=gsc(m,833333345000000041ll);
        else m=gsc(m,ksm(6,p-2));
        m=(m==0)?p-1:m-1;
        for(int i=1;i<=cnt&&prim[i]<=tn;i++)
        {
            if(tn%prim[i]==0)
            {
                b[++tot]=prim[i];ci[tot]=0;
                while(tn%prim[i]==0)ci[tot]++,tn/=prim[i];
            }
        }if(tn!=1)b[++tot]=tn,ci[tot]=1;
        dfs(1,1,1);
        if(n%1000000007==0)
        {
            ans/=1000000007;
            p=1000000007;
            ans=1ll*ans*ksm(n/p,p-2)%p;
        }
        else ans=1ll*ans*ksm(n%p,p-2)%p;
        printf("%lld\n",ans);
    }
}