BZOJ 4176 Lucas的数论

BZOJ 4176 Lucas的数论

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

题目是要求计算：

$$\sum_{i=1}^n\sum_{j=1}^n\sum_{d|ij}1$$

设：

$$nm=\prod_{i=1}^rP_i^{x_i+y_i}$$

其中：

$$n=\prod_{i=1}^rP_i^{x_i}$$

$$m=\prod_{i=1}^rP_i^{y_i}$$

则对于：

$$f(nm)=\sum_{d|nm}1\\=\sum_{k=1}^r\sum_{c\big|P_k^{x_k+y_k}}\sum_{d\big|\frac{nm}{P_{k}^{x_k+y_k}}}1\\=\sum_{k=1}^r\sum_{d\big|\frac{nm}{P_{k}^{x_k+y_k}}}\sum_{c\big|P_k^{x_k+y_k}}1\\=\sum_{k=1}^r\sum_{d\big|\frac{nm}{P_{k}^{x_k+y_k}}}\sum_{a\big|P_k^{x_k}}\sum_{b\big|P_k^{y_k}}[a\perp b]\\=\sum_{k=1}^r\sum_{a\big|P_k^{x_k}}\sum_{b\big|P_k^{y_k}}[a\perp b]\\=\sum_{a|n}\sum_{b|m}[a\perp b]$$

那么： $$answer=\sum_{i=1}^n\sum_{j=1}^n\sum_{a|i}\sum_{b|j}[a\perp b]\\=\sum_{a=1}^n\sum_{b=1}^n\sum_{a|i,i\leq n}\sum_{b|j,j\leq n}[a\perp b]$$

令：

$$F(k)=\sum_{a=1}^n\sum_{b=1}^n\sum_{a|i,i\leq n}\sum_{b|j,j\leq n}[k\big |\ gcd(a,b)]\\=\sum_{a=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{b=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{a|i,i\leq \lfloor\frac{n}{k}\rfloor}\sum_{b|j,j\leq \lfloor\frac{n}{k}\rfloor}1\\=\Big(\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\big\lfloor\frac{\lfloor\frac{n}{k}\rfloor}{i}\big\rfloor\Big)^2$$

记： $$S(n)=\Big(\sum_{i=1}^{n}\big\lfloor\frac{n}{i}\big\rfloor\Big)^2$$

则：

$$answer=\sum_{i=1}^n\mu(i)S(\big\lfloor\frac{n}{i}\big\rfloor)$$

令$m=\big\lfloor\sqrt{n}\big\rfloor$，$M$为梅藤斯函数

$$answer=\sum_{L=1}^mS(L)(M(\big\lfloor\frac{n}{L}\big\rfloor)-M(\big\lfloor\frac{n}{L+1}\big\rfloor))+\sum_{i=1}^{\lfloor\frac{n}{m+1}\rfloor}\mu(i)S(\big\lfloor\frac{n}{i}\big\rfloor)$$

显然计算$S(n)$复杂度是$O(\sqrt{n})$

<script type="math/tex; mode=display" id="MathJax-Element-15"></script>

而对于梅藤斯函数

因为$\mu*I=e$

直接套用公式有：

$$M(n)=1-\sum_{i=2}^nM(\big\lfloor\frac{n}{i}\big\rfloor)$$

注：对这一步有问题的可以参阅链接 :http://blog.csdn.net/ZLH_HHHH/article/details/77836062

进而有：

$$M(n)=1-\sum_{L=1}^mM(L)\big(\big\lfloor\frac{n}{L}\big\rfloor-\big\lfloor\frac{n}{L+1}\big\rfloor\big)+\sum_{i=2}^{\lfloor\frac{n}{m+1}\rfloor}M(\big\lfloor\frac{n}{i}\big\rfloor)$$

$m$不要取太大。BZOJ程序是逐个运行后算总时间的。。。

如果$m=\lfloor n^{\frac{2}{3}}\rfloor$,计算$M(n)$的时间为：

$$\sum_{i=1}^{n^{\frac{1}{3}}}O(\sqrt{\frac{n}{i}})=O(n^{\frac{2}{3}})$$

总时间复杂度：

$$T(n)=\sum_{i=1}^{\sqrt{n}}\Big(2*O(\sqrt{i})+O\big((\frac{n}{i})^{\frac{2}{3}}\big)\Big)=O(n^{\frac{5}{6}})\leq O\big(10^{7.5}\big)$$

实际上， m <script type="math/tex" id="MathJax-Element-24">m</script>略小反而更快。

下面是代码：

#include <algorithm>
#include <string.h>
#include <cmath>
#include <stdio.h>
#define MAXN 1000009
using namespace std;
typedef long long LL;
const int P=1e9+7;
const int P_=P-1;
int mu[100009];
int tmpm[MAXN]={0,-1};
int tmp[100009];
int ML[100009];
void slove(int n,int d)
{
    int ans=0,m=(int)sqrt(n)+1;
    for(int L=1;L<m;L++)
    {
        int u=n/L-n/(L+1);
        ans+=(LL)tmpm[L]*u%P;
        if(ans>P_)ans-=P;
    }
    for(int i=n/m;i>1;i--)
    {
        int u=n/i;
        if(u<MAXN)
            ans+=tmpm[u];
        else
            ans+=tmp[i*d];
        if(ans>P_)ans-=P;
    }
    tmp[d]=(P+1-ans)%P;
}

int M(int n)
{
    if(n<MAXN)return tmpm[n];
    for(int i=n/MAXN; i ; i--)    slove(n/i,i);
    return tmp[1];
}

int S(int n)
{
    int ans=0,m=(int)(sqrt((double)n)+0.0001)+1;
    for(int L=1;L<m;L++)
    {
        int u=n/L-n/(L+1);
        ans+=(LL)L*u%P;
        if(ans>P_)ans-=P;
    }
    for(int i=n/m;i;i--)
    {
        ans+=n/i;
        if(ans>P_)ans-=P;
    }
    ans=(LL)ans*ans%P;
    return ans;
}

void init()
{
    for(int i=1;i<MAXN;i++)
    {
        tmpm[i]=-tmpm[i];
        if(i<100000)mu[i]=tmpm[i];
        for(int j=i<<1;j<MAXN;j+=i) tmpm[j]+=tmpm[i];
        tmpm[i]+=tmpm[i-1];
        if(tmpm[i]<0)tmpm[i]+=P;
        if(tmpm[i]>P_)tmpm[i]-=P;
    }
}

int main ()
{
    init();
    int n;
    scanf("%d",&n);
    int m=(int)sqrt((double)n)+1,ans=0;
    ML[1]=M(n);
    for(int L=1;L<m;L++)
    {
        ML[L+1]=M(n/(L+1));
        ans+=(LL)S(L)*(ML[L]-ML[L+1]+P)%P;
        if(ans>P_)ans-=P;
    }
    for(int i=n/m;i;i--)
    {
        ans+=(LL)mu[i]*S(n/i);
        if(ans>P_)ans-=P;
        if(ans<0)ans+=P;
    }
    printf("%d\n",ans);
    return 0;
}