[bzoj3202]项链 解题报告

这题wa了无数遍。。

先说做法，显然最终答案只与n和不同珠子数量有关系。
所以我们先求不同珠子数量。
如果两个珠子排完序后相同，那么这两个珠子就相同。
考虑珠子里最大的数是i，第二数与i的gcd是d，那么第三个数就必须与d互质。所以答案就是

$$\sum_{i=1}^a\sum_{d|i}\phi ({i \over d})\phi (d){i \over d}$$ 后面一坨显然是个积性函数，所以我们就可以直接筛出来了。但是考虑一下还是有点麻烦，所以简单起见我们把它化成这样： $$\sum_{i=1}^a\phi (i) \sum_{j=1}^{\lfloor {a \over i}\rfloor }\phi (j) j$$ ，这样我们就只需要筛 $\phi$，然后分块求就行了。
但是考虑一下上式，最大数+两个相同的与它互质的数这种情况会被考虑一遍，其余的都会被考虑两遍，所以再加上这种情况一起除以2就是答案了，即： ${\sum_{i=1}^a\phi (i) \sum_{j=1}^{\lfloor {a \over i}\rfloor }\phi (j) j+\sum_{i=1}^a\phi (i)}\over 2$

令$x={{\sum_{i=1}^a\phi (i) \sum_{j=1}^{\lfloor {a \over i}\rfloor }\phi (j) j+\sum_{i=1}^a\phi (i)}\over 2}$。
我们考虑求项链的数量，先考虑如果项链不考虑旋转同构的方案数，这可以dp。设$F_i$表示长度为i的项链不考虑旋转同构的方案数，则$F_i=(x-2)F_{i-1}+(x-1)F_{i-2}$（考虑倒数第二位与正数第一位是否相同），$F_1=0,F_2=x(x-1)$.这个就是一个比较常见的式子了，通常的做法是矩阵乘或者求通项。
矩乘的话就是乘

$$\begin{pmatrix} x-2 & x-1 \\ 1 & 0 \end{pmatrix}$$ ,时间复杂度 $O(2^3（logP+2logP））≈10^3$（设 $P=10^9+7$），考虑这题的数据范围还是有点虚的。
所以我们干脆把通项求出来，这种形如 $F_i=aF_{i-1}+bF_{i-2}(ab≠0)$的通项大致有3种求法。
（差比列）
设 $$F_n+xF_{n-1}$$ $$=(a+x)(F_{n-1}+xF_{n-2})$$ $$=(a+x)F_{n-1}+(a+x)xF_{n-2}$$ $$=(a+x)F_{n-1}+bF_{n-2}$$
$$∴(a+x)x=b$$ ,即 $$x^2+ax-b=0$$ 显然，若设方程解为 $x_1,x_2$,则 $x_1x_2≠0$, $a=-(x_1+x_2)$,将 $x=x_1$代入递推式，则有 $$F_n+x_1F_{n-1}$$ $$=(a+x_1)(F_{n-1}+x_1F_{n-2})$$ $$=(-x_2)(F_{n-1}+x_1F_{n-2}),$$ $$∴F_n+x_1F_{n-1}=(-x_2)^{n-1}(F_1+x_1F_0)$$ ,同除 $(-x_2)^n$，得 $${F_n \over (-x_2)^n}={x_1\over x_2}{F_{n-1}\over (-x_2)^{n-1}}-{1\over x_2}(F_1+x_1F_0)$$ ，令 $G_n={F_n \over (-x_2)^n}$，则 $G_n={x_1\over x_2}G_{n-1}-{1\over x_2}(F_1+x_1F_0)$,
若 ${x_1\over x_2}≠1$,即 $x_1≠x_2$，设 $$G_n+y$$ $$={x_1\over x_2}(G_{n-1}+y)$$ $$={x_1\over x_2}G_{n-1}+{x_1 \over x_2}y$$ $$={x_1 \over x_2}G_{n-1}-{1\over x_2}(F_1+x_1F_0)+y,$$ $$∴{x_1\over x_2}y=-{1\over x_2}(F_1+x_1F_0)+y$$ $$∴y={F_1+x_1F_0\over x_2-x_1},$$ $$∴G_n+{F_1+x_1F_0\over x_2-x_1}=({x_1\over x_2})^n(G_0+{F_1+x_1F_0\over x_2-x_1})$$ $$∴G_n=({x_1\over x_2})^n(G_0+{F_1+x_1F_0\over x_2-x_1})-{F_1+x_1F_0\over x_2-x_1},$$ $$∴F_n=(F_0+{F_1+x_1F_0\over x_2-x_1})(-x_1)^n+-{F_1+x_1F_0\over x_2-x_1}(-x_2)^n$$ $$={F_1+x_2F_0\over x_2-x_1}(-x_1)^n+-{F_1+x_1F_0 \over x_2-x_1}(-x_2)^n$$ 若 $x_1=x_2$，则 $G_n=G_{n-1}-{1\over x_2}(F_1+x_1F_0),$ $$∴G_n=-{F_1+x_1F_0\over x_2}n+G_0,$$ $$∴F_n=(-{F_1+x_1F_0\over x_2}n+G_0)(-x_2)^n$$
（母函数）
我并不是很会。。貌似需要很高端的数学知识。
（特征根法）
其实是上述做法的简化，考虑方程 $x^2=ax+b$的两根 $x_1,x_2$，是方程 $x^2+ax-b=0$的两根的相反数。注意到当 $x_1=x_2$时，通项公式为 $F_n=Ax_1^n+Bx_2^n$，当 $x_1≠x_2$时，通项公式为 $F_n=(An+B)x_1^n$.所以我们只需代入 $F_0,F_1$解出A,B即可。

所以考虑题目$F_n=(x-2)F_{n-1}+(x-1)F_{n-2}$,首项$F_1=0,F-2=x(x-1)$(注意$F_0$并非首项),由以上推导易得$F_n=(x-1)((x-1)^{n-1}+(-1)^n)$,注意到实际上，在$x-1≡-1 \ \mod p$的时候，这个式子依然是成立的。这样我们就能在$O(logn)$的时间内求出一个$F_n$了。

我们首先考虑一下关于循环节的问题。首先长度为n的项链的循环节长度必为n的约数。且若其有长为p的循环节和长为q的循环节，则其必然有长为(p,q)的循环节，因为若设p>q，该项链中第i个珠子为$x_i$，则有$x_i=x_{i+p-q}$，所以p-1也是该项链的循环节，所以(p,q)也是该项链的循环节。
那么显然，一条最短循环节为d的项链在不考虑循环的计算中会被计算d次，所以设最短循环节为d的项链有f(d)条，则$ans=\sum_{d|n}{1\over d}f(d)$，而考虑我们刚才求的F(n)，显然有$F(n)=\sum_{d|n}f(d)$,$f(n)=\sum_{d|n}F(d)\mu ({n\over d})$,$ans=\sum_{d|n}{1\over d}\sum_{i|d}F(i)\mu ({d \over i})=\sum_{d|n}F(d){1\over d}\sum_{i|{n \over d}}\mu i {1\over i}={\sum_{d|n}F(d)\phi ({n\over d})\over n}$

那么我们就可以在$O(sqrt(n)logn)≈5*10^8$的时间复杂度内回答单次询问了，注意到实际上这个$O(sqrt(n))$是很虚的，实际上应该只有$10^5~10^6$左右，所以实际上只会有$10^7$左右。

但是这道题就这么解决了么？不！注意到n非常非常大，有可能是p的倍数，所以如果我们要除n的话，就可能会面临没有逆元的问题。
那么这时候该怎么办呢？实际上有式子：${ab\ \mod \ ap\over a}=b\ \mod \ p$，这个还是比较显然的。但是这个n非常大啊，我要是直接乘n的话不就爆了么？但是注意到我们其实可以把除n分解成除p再除$n\over p$,注意到n中只会有一个p，而因为$n\over p$是有逆元的，所以我们可以对这两部分分别处理。
但是这样的话就需要在快速幂里再加一个快速乘，时间复杂度不就又爆了么？但是实际上，考虑到我们需要快速乘的时候是$p|n$之时，这是n的约数只会有$O(2\sqrt {n\over p})≈10^2$是非常少的。所以时间复杂度并不会增加，只是预处理更慢了而已。。

所以总的时间复杂度就是$O(a\log p+T(\sqrt n \log p+\sqrt {n \over p}\log ^2 p))≈10^8$.
然而似乎是预处理的时候太慢了跑了22s。。因为预处理的时候需要long long+快速乘，慢的飞起。。而如果把$\phi$换成$\mu$的话就会快很多(只跑了4s)。。这也是这题一个奇怪的地方，因为一般情况下用$\phi$都是会比用$\mu$更优的。

然后膜拜xyz大爷只跑了600ms。。应该是有什么更神的做法吧。。（难道是用杜教筛筛了$\mu$…如果用杜教筛筛$\mu$+数据较水，时间复杂度就成了$O(T(\sqrt a+\log p+\log n \log p+\log {n \over p}\log ^2 p))≈7*10^5$）
代码（$\phi$）：

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;
typedef long long LL;
const LL Mod[2]={(LL)(1e9+7)*(int)(1e9+7),1e9+7},phi[2]={(LL)(1e9+7)*(int)(1e9+6),1e9+6};
LL mul(LL a,LL b,bool o){
    if(a<=2e9&&b<=2e9)return a*b%Mod[o];
    else{
        LL ans=0,sum=a%Mod[o];
        for(b=(b+Mod[o])%Mod[o];b;b>>=1,sum=(sum<<1)%Mod[o])
            if(b&1)
                ans=(ans+sum)%Mod[o];
        return ans;
    }
}
LL pow(LL a,LL x,bool o){
    LL ans=1,prod=a%Mod[o];
    //cout<<"Pow("<<a<<"^"<<x<<")=";
    for(x%=Mod[o]-1;x;x>>=1,prod=mul(prod,prod,o))
        if(x&1)
            ans=mul(ans,prod,o);
    //cout<<ans<<" "<<prod<<endl;
    return ans;
}
const int A=1e7+5;
LL sphi[2][A],sphii[2][A];
int prime[1000005];
LL dvs[50];
int ind[50];
LL x;
LL ans;
void dfs(int i,LL f,LL phi,bool o){
    if(i<0){
        //cout<<"dfs("<<i<<","<<f<<","<<phi<<","<<o<<")\n";
        ans=(ans+mul(mul(x-1,pow(x-1,f-1,o)+(f&1?-1:1),o),phi,o))%Mod[o];
        return;
    }
    dfs(i-1,f,phi,o);
    dfs(i-1,f/=dvs[i],phi*=(dvs[i]-1),o);
    for(int j=ind[i];--j;)dfs(i-1,f/=dvs[i],phi*=dvs[i],o);
}
int main(){
    for(int i=2;i<=10000000;++i){
        if(!sphi[0][i]){
            sphi[0][i]=i-1;
            prime[++prime[0]]=i;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=10000000;++j)
            if(i%prime[j])sphi[0][i*prime[j]]=sphi[0][i]*(prime[j]-1);
            else sphi[0][i*prime[j]]=sphi[0][i]*prime[j];
    }
    sphi[0][1]=1;
    for(int i=1;i<=10000000;++i)
        for(int o=2;o--;){
            sphii[o][i]=(sphii[o][i-1]+sphi[0][i]*i)%Mod[o];
            sphi[o][i]=(sphi[o][i-1]+sphi[0][i])%Mod[o];
        }
    //cout<<sphi[0][10000000]<<" "<<sphii[0][10000000]<<endl;

    int t;
    cin>>t;
    while(t--){
        LL n;
        int a;
        cin>>n>>a;

        bool o=n%Mod[1];
        x=0;
        for(int i=1,j;i<=a;i=j+1){
            j=a/(a/i);
            //x=(x+mul(sphii[o][a/i],sphi[o][j]-sphi[o][i-1],o))%Mod[o];
            x=(x+sphii[o][a/i]*(sphi[o][j]-sphi[o][i-1]))%Mod[o];
            //printf("[%d,%d]=%I64d*%I64d=%I64d\n",i,j,sphii[o][a/i],sphi[o][j]-sphi[o][i-1],sphii[o][a/i]*(sphi[o][j]-sphi[o][i-1]));
        }
        x=mul((x+sphi[o][a])%Mod[o],pow(2,phi[o]-1,o),o);

        int dtot=0;
        for(int i=1;i<=prime[0]&&(LL)prime[i]*prime[i]<=n;++i)
            if(n%prime[i]==0){
                dvs[dtot]=prime[i];
                ind[dtot]=0;
                while(n%prime[i]==0)n/=prime[i],++ind[dtot];
                ++dtot;
            }
        if(n!=1){
            dvs[dtot]=n,ind[dtot]=1;
            ++dtot;
            n=1;
        }
        for(int i=dtot;i--;)
            for(int j=ind[i];j--;)
                n*=dvs[i];
        /*cout<<n<<"=";
        for(int i=dtot;i--;)cout<<dvs[i]<<"^"<<ind[i]<<"*";
        puts("");*/

        ans=0;
        dfs(dtot-1,n,1,o);
        if(!o)ans/=Mod[1],n/=Mod[1];
        ans=(ans*pow(n,Mod[1]-2,1)%Mod[1]+Mod[1])%Mod[1];
        cout<<ans<<endl;
    }
}

代码（$\mu$）：

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;
typedef long long LL;
const LL Mod[2]={(LL)(1e9+7)*(int)(1e9+7),1e9+7},phi[2]={(LL)(1e9+7)*(int)(1e9+6),1e9+6};
LL mul(LL a,LL b,bool o){
    if(a<=2e9&&b<=2e9)return a*b%Mod[o];
    else{
        LL ans=0,sum=a%Mod[o];
        for(b=(b+Mod[o])%Mod[o];b;b>>=1,sum=(sum<<1)%Mod[o])
            if(b&1)
                ans=(ans+sum)%Mod[o];
        return ans;
    }
}
LL pow(LL a,LL x,bool o){
    LL ans=1,prod=a%Mod[o];
    //cout<<"Pow("<<a<<"^"<<x<<")=";
    for(x%=Mod[o]-1;x;x>>=1,prod=mul(prod,prod,o))
        if(x&1)
            ans=mul(ans,prod,o);
    //cout<<ans<<" "<<prod<<endl;
    return ans;
}
const int A=1e7+5;
int smu[A];
bool p[A];
int prime[1000005];
LL dvs[50];
int ind[50];
LL x;
LL ans;
void dfs(int i,LL f,LL phi,bool o){
    if(i<0){
        //cout<<"dfs("<<i<<","<<f<<","<<phi<<","<<o<<")\n";
        ans=(ans+mul(mul(x-1,pow(x-1,f-1,o)+(f&1?-1:1),o),phi,o))%Mod[o];
        return;
    }
    dfs(i-1,f,phi,o);
    dfs(i-1,f/=dvs[i],phi*=(dvs[i]-1),o);
    for(int j=ind[i];--j;)dfs(i-1,f/=dvs[i],phi*=dvs[i],o);
}
int main(){
    smu[1]=1;
    for(int i=2;i<=10000000;++i){
        if(!p[i]){
            prime[++prime[0]]=i;
            smu[i]=-1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=10000000;++j){
            p[i*prime[j]]=1;
            if(i%prime[j])smu[i*prime[j]]=-smu[i];
            else break;
        }
    }
    for(int i=2;i<=10000000;++i)smu[i]+=smu[i-1];

    int t;
    cin>>t;
    while(t--){
        LL n;
        int a;
        cin>>n>>a;

        bool o=n%Mod[1];
        LL f2=0,f3=0,prod;
        x=0;
        for(int i=1,j;i<=a;i=j+1){
            j=a/(a/i);
            prod=(LL)(a/i)*(a/i)%Mod[o];
            f2=(f2+mul(prod,smu[j]-smu[i-1],o))%Mod[o];
            f3=(f3+mul(mul(prod,a/i,o),smu[j]-smu[i-1],o))%Mod[o];
        }
        x=mul((f3+mul(f2,3,o)+2)%Mod[o],pow(6,phi[o]-1,o),o);
        //cout<<mul(6,pow(6,phi[o]-1,o),o)<<endl;
        //cout<<f3+mul(f2,3,o)+2<<"*"<<pow(6,phi[o]-1,o)<<"="<<mul((f3+mul(f2,3,o)+2)%Mod[o],pow(6,phi[o]-1),o)<<endl;
        //cout<<f2<<" "<<f3<<"->"<<x<<endl; 

        int dtot=0;
        for(int i=1;i<=prime[0]&&(LL)prime[i]*prime[i]<=n;++i)
            if(n%prime[i]==0){
                dvs[dtot]=prime[i];
                ind[dtot]=0;
                while(n%prime[i]==0)n/=prime[i],++ind[dtot];
                ++dtot;
            }
        if(n!=1){
            dvs[dtot]=n,ind[dtot]=1;
            ++dtot;
            n=1;
        }
        for(int i=dtot;i--;)
            for(int j=ind[i];j--;)
                n*=dvs[i];
        /*cout<<n<<"=";
        for(int i=dtot;i--;)cout<<dvs[i]<<"^"<<ind[i]<<"*";
        puts("");*/

        ans=0;
        dfs(dtot-1,n,1,o);
        if(!o)ans/=Mod[1],n/=Mod[1];
        ans=(ans*pow(n,Mod[1]-2,1)%Mod[1]+Mod[1])%Mod[1];
        cout<<ans<<endl;
    }
}

这题因为不会写暴力，wa了无数遍，这里写下每一遍wa的原因：
①没有考虑n可能为p的倍数。
②快速乘的时候要保证分解的那个数是正数。
③指数应该模$\phi (p)$,而不是p。
④误以为珠子的翻转可以用旋转搞出来。。然而实际并不能。
⑤没有想清楚模$p^2$的时候该怎么搞。。一直在改一点wa一遍改一点wa一遍。。在一遍遍wa中渐渐发现要全改。。
⑥在开存质因子的数组时用了int。。忘了开long long。

总结：
①${1 \over {1-ax}}=\sum_{n=0}^{+∞}(ax)^n$
②一个字符串的所有循环节长度必然均为其最短循环节的倍数。
③在做模除法时，一定要考虑清楚是否会出现没有逆元的情况。
④除一个没有逆元的数时，可以对其没有逆元的部分和有逆元的部分分开处理。
⑤通常情况下用$\phi$会比用$\mu$更优，但是若预处理的时间复杂度成为主要矛盾，把$\phi$换成$\mu$通常会是更好的选择，因为筛$\mu$比筛$\phi$快很多。
⑥如果不知道怎么写暴力，就写两个标算对拍吧~
⑦快速乘的时候要保证分解的那个数是正数。
⑧指数应该模$\phi (p)$,而不是p。
⑨写代码前一定要仔细想清楚+冷静。
⑩在开变量之前一定要检查变量的范围！