BZOJ 4176 Lucas的数论

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

题目是要求计算：

∑i=1n∑j=1n∑d|ij1

设：

nm=∏i=1rPxi+yii

其中：

n=∏i=1rPxii

m=∏i=1rPyii

则对于：

f(nm)=∑d|nm1=∑k=1r∑c∣∣∣Pxk+ykk∑d∣∣∣nmPxk+ykk1=∑k=1r∑d∣∣∣nmPxk+ykk∑c∣∣∣Pxk+ykk1=∑k=1r∑d∣∣∣nmPxk+ykk∑a∣∣∣Pxkk∑b∣∣∣Pykk[a⊥b]=∑k=1r∑a∣∣∣Pxkk∑b∣∣∣Pykk[a⊥b]=∑a|n∑b|m[a⊥b]

那么：

answer=∑i=1n∑j=1n∑a|i∑b|j[a⊥b]=∑a=1n∑b=1n∑a|i,i≤n∑b|j,j≤n[a⊥b]

令：

F(k)=∑a=1n∑b=1n∑a|i,i≤n∑b|j,j≤n[k∣∣ gcd(a,b)]=∑a=1⌊nk⌋∑b=1⌊nk⌋∑a|i,i≤⌊nk⌋∑b|j,j≤⌊nk⌋1=(∑i=1⌊nk⌋⌊⌊nk⌋i⌋)2

记：

S(n)=(∑i=1n⌊ni⌋)2

则：

answer=∑i=1nμ(i)S(⌊ni⌋)

令m=⌊n‾‾√⌋，M为梅藤斯函数

answer=∑L=1mS(L)(M(⌊nL⌋)−M(⌊nL+1⌋))+∑i=1⌊nm+1⌋μ(i)S(⌊ni⌋)

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

而对于梅藤斯函数

因为μ∗I=e

直接套用公式有：

M(n)=1−∑i=2nM(⌊ni⌋)

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

进而有：

M(n)=1−∑L=1mM(L)(⌊nL⌋−⌊nL+1⌋)+∑i=2⌊nm+1⌋M(⌊ni⌋)

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

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

∑i=1n13O(ni‾‾√)=O(n23)

总时间复杂度：

T(n)=∑i=1n√(2∗O(i√)+O((ni)23))=O(n56)≤O(107.5)

实际上，m略小反而更快。

下面是代码：

#include <algorithm>#include <string.h>#include <cmath>#include <stdio.h>#define MAXN 1000009using 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;}