• 【学术篇】51nod 1238 最小公倍数之和


    这是一道杜教筛的入(du)门(liu)题目...


    题目大意

    [sum_{i=1}^nsum_{j=1}^nlcm(i,j) ]

    一看就是辣鸡反演一类的题目, 那就化式子呗..

    [sum_{i=1}^nsum_{j=1}^nlcm(i,j) \ =sum_{i=1}^nsum_{j=1}^nfrac{ij}{gcd(i,j)} \ =sum_{i=1}^nsum_{j=1}^nsum_{k=1}^nfrac{ij}k[gcd(i,j)=k] \ =sum_{k=1}^nsum_{i=1}^{leftlfloorfrac nk ight floor}sum_{j=1}^{leftlfloorfrac nk ight floor}ijk[gcd(i,j)=1] \ =sum_{k=1}^nksum_{i=1}^{leftlfloorfrac nk ight floor}isum_{j=1}^{leftlfloorfrac nk ight floor}j[gcd(i,j)=1] \ =sum_{k=1}^nksum_{i=1}^{leftlfloorfrac nk ight floor}isum_{j=1}^{leftlfloorfrac nk ight floor}jsum_{d=1}^{leftlfloorfrac nk ight floor}[d|i][d|j]mu(d) \ =sum_{k=1}^nksum_{d=1}^{leftlfloorfrac nk ight floor}mu(d) sum_{i=1}^{leftlfloorfrac nk ight floor}[d|i]i sum_{j=1}^{leftlfloorfrac nk ight floor}[d|j]j \ 设i=dq,j=dq, \ 原式=sum_{k=1}^nksum_{d=1}^{leftlfloorfrac nk ight floor}mu(d) sum_{p=1}^{leftlfloorfrac nk ight floor}dp sum_{q=1}^{leftlfloorfrac nk ight floor}dq \ =sum_{k=1}^nksum_{d=1}^{leftlfloorfrac nk ight floor}mu(d)*d^2sum_{p=1}^{leftlfloorfrac nk ight floor}p sum_{q=1}^{leftlfloorfrac nk ight floor}q \ ecausesum_{i=1}^ni=frac{n(n+1)}2, \ herefore原式=sum_{k=1}^nksum_{d=1}^{leftlfloorfrac nk ight floor}mu(d)*d^2(frac{{leftlfloorfrac nk ight floor}({leftlfloorfrac nk ight floor}+1)}{2})^2 \ 设t=kd, \ 原式=sum_{t=1}^n(frac{{leftlfloorfrac nt ight floor}({leftlfloorfrac nt ight floor}+1)}{2})^2sum_{d|t}mu(d)*td \ 设g(x)=sum_{d|x}mu(d)*xd, \ 则原式=sum_{t=1}^n(frac{{leftlfloorfrac nt ight floor}({leftlfloorfrac nt ight floor}+1)}{2})^2g(t) ]

    然后前面这一堆我们可以分块(O(1))搞出来, 所以问题的关键就变成了(g(x))的前缀和.
    因为题目中(nleq10^{10}), 所以这里要用(O(n^frac23))的杜教筛.
    但是杜教筛的框架我们是知道的,

    [设要求前缀和的函数为f(x), 其前缀和为s_f(x) \ 设有好求前缀和的积性函数g(x)和h(x), 使得(f*g)(x)=h(x), \ 设h(x)的前缀和为s_h(x), 则有 \ s_f(x)=frac{s_h(x)-sum_{d=2}^ns_f({leftlfloorfrac nd ight floor})g(d)}{g(1)} ]

    所以我们的任务又变成了找到合适的(g(x))(h(x)).
    在某篇歪果仁的blog上, 他说:

    With a little luck, we can notice that Id(l) = g * Id2(l),

    • 注: (id(x)=x)(以下称为(n))(id2(x)=x^2)(以下称为(n_2))

    But 这运气也太好了点吧? 我怎么就找不到这种函数呢?!
    不过我们还是要考虑一下的... (这波化式子的时候注意乘号(cdot)和卷积号(*))

    [g(x)=sum_{d|x}mu(d)xd \ =sum_{d|x}d^2mu(d)frac xd \ =(n^2cdotmu*n)(x) \ (g*n_2)(x) \ =((n^2cdotmu*n)*n_2)(x) cdots ① \ 积性函数有一个性质, 如果两边的乘积都含有某一项的时候, 可以把这一项提出来. \ 所以我们可以把x^2都提出来, (证明可以用定义式推咯~) \ ①=((n^2cdotmu*n_2)*n)(x) \ =(n^2cdot(mu*1)*n)(x) \ =(n^2cdotepsilon*n)(x) =n(x) ]

    所以我们就证明出来了((g*n_2)(x)=n(x))
    这样就可以直接杜教筛了= =
    然后就是注意细节了OvO
    比如非常坑人的(10^{10}*(1e9+7))会爆long long!!
    所以要开unsigned...
    就这样吧= =
    代码:

    #include <map>
    #include <cstring>
    #include <iostream>
    using namespace std;
    typedef unsigned long long LL;
    #define ri register int
    #define rll register unsigned long long 
    const int p=1e9+7;
    int qpow(int a,int b,int s=1){
    	for(;b;b>>=1,a=1LL*a*a%p)
    		if(b&1) s=1LL*s*a%p;
    	return s;
    }
    const int inv2=qpow(2,p-2),inv4=qpow(4,p-2),inv6=qpow(6,p-2);
    const int N=4800000,_N=N+10;
    const int SQ=233333;
    LL n;
    map<LL,int> mp;
    map<LL,int>::iterator it;
    int g[_N],pr[N/2],tot;
    bool notp[_N];
    void shai(){
    	g[1]=1; ri i,j; rll k;
    	for(i=2;i<=N;++i){
    		if(!notp[i])
    			pr[++tot]=i,g[i]=(-1LL*i*(i-1)%p+p)%p;
    		for(j=1;j<=tot&&(k=1LL*i*pr[j])<=N;++j){
    			notp[k]=1;
    			if(i%pr[j]==0){g[k]=1LL*g[i]*pr[j]%p;break;}
    			g[k]=1LL*g[i]*g[pr[j]]%p;
    		}
    	}
    	for(i=2;i<=N;++i) g[i]=(g[i]+g[i-1])%p;
    }
    inline int calcsqr(LL x){x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;}
    inline int calcsum(LL x){x%=p;return x*x%p*(x+1)%p*(x+1)%p*inv4%p;}
    int calc(LL x){
    	if(x<=N) return g[x];
    	it=mp.find(x);
    	if(it!=mp.end()) return it->second;
    	int ans=1LL*x%p*(x+1)%p*inv2%p;
    	LL last=0;
    	for(rll i=2;i<=x;i=last+1){
    		last=x/(x/i);
    		ans=(ans-1LL*calc(x/i)*(calcsqr(last)-calcsqr(i-1)+p)%p+p)%p;
    	}
    	return mp[x]=(ans%p+p)%p;
    }
    int main(){
    	shai();cin>>n;
    	LL last=0; int ans=0;
    	for(rll i=1;i<=n;i=last+1){
    		last=n/(n/i);
    		ans=(ans+1LL*calcsum(n/i)*(calc(last)-calc(i-1)+p)%p)%p;
    	}
    	cout<<ans;
    }
    

    这破题我连推带写加总结写了整整一天!!!
    窝还是太弱了..

  • 相关阅读:
    Agreementhasbeenupdated--EditPhoneNumber
    android 6.0之后动态获取权限
    appstore加速审核通道
    ListView多种item注意以及自己出现的莫名其妙的错误
    2020-10-15:mysql的双1设置是什么?
    2020-10-14:Redisson分布式锁超时自动释放,会有什么问题?
    2020-10-13:hash与B+tree的区别?
    2020-10-12:在做分布式集群时候一般会产生什么问题?
    2020-10-11:一条sql语句执行时间过长,应该如何优化?从哪些方面进行优化?
    2020-10-10:OOM都有哪些,说出几种?
  • 原文地址:https://www.cnblogs.com/enzymii/p/8505328.html
Copyright © 2020-2023  润新知