[Sdoi2013]项链

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=3202

解法：读完题目，很明显我们有两个任务：（1）计算符合提交的颜色数目；（2）计算有m种颜色，长度为n相邻两个颜色不相同的本质不同项链数目（旋转相同被认为是本质相同。）

Part 1：

我们设c2,c3分别为gcd=1的有序整数对以及有序三元组的数目。

那么我们可以得出答案即为m=(c3-(c2-1)*3-1)/6+c2。

计算c2,c3需要莫比乌斯反演的知识，这里直接给出公式好了（不会的同学可以参考《具体数学》）

c2=simga(miu(i)*floor(a/i)^2),c3=simga(miu(i)*floor(a/i)^3),1<=i<=a。

Part 2：

我们设上面所求答案m。

接下来需要用到polya计数相关知识。

我们设f(i)表示长度为i的项链，相邻颜色两两不同的方案数。

那么答案就是sigma(f(gcd(i,n)),1<=i<=n)/n=sigma(phi(d)*f(n/d),d|n)/n。

如何求f(i)?

若第i-1个和第1个的颜色不同，情况数有f(i-1)，那么第i个可以放m-2种颜色，方案数为f(i-1)*(m-2)

若第i-1个和第1个的颜色相同，那么i-2肯定和第1个不同，情况数有f(i-2),第i个可以放m-1种颜色。

所以f(i)=f(i-1)*(m-2)+f(i-2)*(m-1)。于是就可以在log(n)的时间求出f(n)了（我的程序里是求过通项的）。

于是这题差不多解决了= =。我们可以通过官方数据了。。

但是。。。vfk大神在bzoj上加了一组n%mod!=0的数据。。这样我们到最后就不能直接乘逆元了。。

于是对于这种情况我们先对答案%mod^2，然后这个结果必然被mod整除，于是我们直接除掉，再乘上(n/mod)的逆元。（n<mod^2，所以这下n/mod的逆元存在了。。）

下面是代码（在bzoj上目前rank1）

/**************************************************************
    Problem: 3202
    User: hta
    Language: C++
    Result: Accepted
    Time:4636 ms
    Memory:93072 kb
****************************************************************/
 
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <string>
using namespace std;
typedef long long LL;
const int maxn=10000003;
const int mod=1000000007;
int A[11],miu[maxn],T,a,sum[maxn];
bool check[maxn];
int totp,prime[maxn/10];
LL N[11],p;
inline int get()
{
    int f=0,v=0;char ch;
    while(!isdigit(ch=getchar()))if(ch=='-')break;
    if(ch=='-')f=1;else v=ch-48;
    while(isdigit(ch=getchar()))v=v*10+ch-48;
    if(f==1)return -v;else return v;
}
 
void init()
{
    T=get();int n=0;
    for(int i=1;i<=T;i++)cin>>N[i],A[i]=get(),n=max(n,A[i]);
    miu[1]=sum[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!check[i])prime[++totp]=i,miu[i]=-1;
        for(int j=1,tp;(tp=prime[j]*i)<=n&&j<=totp;j++)
        {
            check[tp]=1;
            if(i%prime[j]==0){miu[tp]=0;break;}else miu[tp]=-miu[i];
        }
    }
    for(int i=2;i<=n;i++)sum[i]=sum[i-1]+miu[i];
}
/*LL Mul(LL a,LL b)
{
    LL res=0;a%=p,b%=p;
    if(b<0)b+=p;
    for(;b>0;b/=2,a+=a,a-=(a>=p)*p)
        if(b%2==1)res+=a,res-=(res>=p)*p;
    return res;
}*/
 
LL Mul(LL a,LL b)
{
  return (a*b-(LL)((long double)a/p*b)*p+p)%p;
}
void extgcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0){x=1,y=0;return;}
    extgcd(b,a%b,x,y);
    LL tp=x;x=y;y=tp-a/b*x;
}
 
LL Pow(LL a,LL b)
{
    LL res=1;
    for(a%=mod;b>0;b/=2,a=a*a%mod)
        if(b&1)res=res*a%mod;
    return res;
}
LL calc(int a)
{
    int i=1,j,tp;
    LL c2=0,c3=0;
    while(i<=a)
    {
        j=a/(a/i),tp=a/i;
        LL tmp=LL(tp)*tp%mod*(sum[j]-sum[i-1])%mod;
        c2=(c2+tmp)%mod,c3+=tmp*tp,c3%=mod;
        i=j+1;
    }
    LL res=(c3-(c2-1)*3-1)%mod*Pow(6,mod-2)+c2;
    res%=mod;if(res<0)res+=mod;
    return res;
}
LL phi(LL n)
{
    LL res=n;
    for(LL i=2;i*i<=n;i++)
    {
        if(n%i!=0)continue;
        res-=res/i;
        while(n%i==0)n/=i;
    }
    if(n>1)res-=res/n;
    return res;
}
LL f(LL n,LL m)
{
    LL res;
    if(n%2==1)res=(Pow(m-1,n)-m+1+mod)%mod;
    else res=(Pow(m-1,n)+m-1+mod)%mod;
    return res;
}
 
 
LL Pow2(LL a,LL b)
{
    LL res=1;
    for(a%=p;b>0;b/=2,a=Mul(a,a))
        if(b&1)res=Mul(res,a);
    return res;
}
LL calc2(int a)
{
    int i=1,j;
    LL c2=0,c3=0,tp;
    while(i<=a)
    {
        j=a/(a/i),tp=a/i;
        LL tmp=Mul(tp*tp%p,sum[j]-sum[i-1]);
        c2=(c2+tmp)%p,c3+=Mul(tmp,tp),c3%=p;
        i=j+1;
    }
    LL x,y;
    extgcd(6,p,x,y);if(x<0)x+=p;
    LL res=Mul((c3-(c2-1)*3-1)%p,x)+c2;
    res%=p;if(res<0)res+=p;
    return res;
}
LL f2(LL n,LL m)
{
    LL res;
    if(n%2==1)res=(Pow2(m-1,n)-m+1+p)%p;
    else res=(Pow2(m-1,n)+m-1+p)%p;
    return res;
}
 
void work(LL n,int a)
{
    if(n%mod==0)
    {
        p=LL(mod)*mod;
        LL col=calc2(a);
        LL ans=0;
        for(LL i=1;i*i<=n;i++)
        {
            if(n%i!=0)continue;
            ans+=Mul(phi(i)%p,f2(n/i,col));
            if(i*i!=n)ans+=Mul(phi(n/i)%p,f2(i,col));
            ans%=p;
        }
        ans/=mod;ans%=mod;
        LL x,y;extgcd(n/mod,p,x,y);
        x%=mod;if(x<0)x+=mod;
        printf("%d\n",int(ans*x%mod));
        return;
    }
    LL col=calc(a);
    LL ans=0;
    for(LL i=1;i*i<=n;i++)
    {
        if(n%i!=0)continue;
        ans+=phi(i)%mod*f(n/i,col);
        if(i*i!=n)ans+=phi(n/i)%mod*f(i,col);
        ans%=mod;
    }
    printf("%d\n",int(ans*Pow(n,mod-2)%mod));
}
 
int main()
{
    init();
    for(int i=1;i<=T;i++)work(N[i],A[i]);
    return 0;
}