BZOJ 4176 Lucas的数论

BZOJ 4176 Lucas的数论

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

题目是要求计算:

i=1nj=1nd|ij1

设:

nm=i=1rPxi+yii

其中:

n=i=1rPxii

m=i=1rPyii

则对于:

f(nm)=d|nm1=k=1rcPxk+ykkdnmPxk+ykk1=k=1rdnmPxk+ykkcPxk+ykk1=k=1rdnmPxk+ykkaPxkkbPykk[ab]=k=1raPxkkbPykk[ab]=a|nb|m[ab]

那么:
answer=i=1nj=1na|ib|j[ab]=a=1nb=1na|i,inb|j,jn[ab]

令:

F(k)=a=1nb=1na|i,inb|j,jn[k gcd(a,b)]=a=1nkb=1nka|i,inkb|j,jnk1=(i=1nknki)2

记:
S(n)=(i=1nni)2

则:

answer=i=1nμ(i)S(ni)

m=n M 为梅藤斯函数

answer=L=1mS(L)(M(nL)M(nL+1))+i=1nm+1μ(i)S(ni)

显然计算 S(n) 复杂度是 O(n)

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

而对于梅藤斯函数

因为 μI=e

直接套用公式有:

M(n)=1i=2nM(ni)

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

进而有:

M(n)=1L=1mM(L)(nLnL+1)+i=2nm+1M(ni)

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

如果m=n23,计算 M(n) 的时间为:

i=1n13O(ni)=O(n23)

总时间复杂度:

T(n)=i=1n(2O(i)+O((ni)23))=O(n56)O(107.5)

实际上, 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;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值