zoukankan      html  css  js  c++  java
  • 51Nod 1238

    题目

    戳这里

    推导

       i=1nj=1nlcm(i,j)~~~sum_{i=1}^{n}sum_{j=1}^{n}lcm(i,j)

    =i=1nj=1nijgcd(i,j)=sum_{i=1}^{n}sum_{j=1}^{n}frac{ij}{gcd(i,j)}

    =i=1nd1i=1nj=1nij[gcd(i,j)==d]=sum_{i=1}^{n}d^{-1}sum_{i=1}^{n}sum_{j=1}^{n}ij[gcd(i,j)==d]

    =i=1ndi=1ndj=1ndij[gcd(i,j)==1]=sum_{i=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor}ij[gcd(i,j)==1]

    =i=1ndi=1ndij=1ndj[gcd(i,j)==1]=sum_{i=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}isum_{j=1}^{lfloorfrac{n}{d} floor}j[gcd(i,j)==1]

    =i=1nd(2i=1ndij=1ij[gcd(i,j)==1]1)=sum_{i=1}^{n}d(2sum_{i=1}^{lfloorfrac{n}{d} floor}isum_{j=1}^{i}j[gcd(i,j)==1]-1)

    =i=1nd(2i=1ndiiφ(i)+[i==1]21)=sum_{i=1}^{n}d(2sum_{i=1}^{lfloorfrac{n}{d} floor}ifrac{ivarphi(i)+[i==1]}{2}-1)

    =i=1ndi=1ndi2φ(i)=sum_{i=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}i^2varphi(i)

    子问题:

    i=1ni2φ(i)sum_{i=1}^{n}i^2varphi(i)
    f(i)=i2φ(i)f(i)=i^2varphi(i)
    使用狄利克雷卷积,卷一个g(i)=i2g(i)=i^2
    那么:

        i=1n(fg)(i)~~~~sum_{i=1}^{n}(f*g)(i)

    =i=1ndif(d)g(id)=sum_{i=1}^{n}sum_{d|i}^{}f(d)g(frac{i}{d})

    =i=1ndid2φ(d)(id)2=sum_{i=1}^{n}sum_{d|i}^{}d^2varphi(d)(frac{i}{d})^2

    =i=1ni2diφ(d)=sum_{i=1}^{n}i^2sum_{d|i}^{}varphi(d)

    =i=1ni3=sum_{i=1}^{n}i^3

    =n2(n+1)24=frac{n^2(n+1)^2}{4}

    又因为:

        i=1ndid2φ(d)(id)2~~~~sum_{i=1}^{n}sum_{d|i}^{}d^2varphi(d)(frac{i}{d})^2

    =i=1ni2d=1nid2φ(d)=sum_{i=1}^{n}i^2sum_{d=1}^{lfloorfrac{n}{i} floor}d^2varphi(d)

    =i=2ni2d=1nid2φ(d)+i=1ni2φ(i)=sum_{i=2}^{n}i^2sum_{d=1}^{lfloorfrac{n}{i} floor}d^2varphi(d)+sum_{i=1}^{n}i^2varphi(i)

    =n2(n+1)24=frac{n^2(n+1)^2}{4}

    所以:

    i=1ni2φ(i)=n2(n+1)24i=2ni2d=1nid2φ(d)sum_{i=1}^{n}i^2varphi(i)=frac{n^2(n+1)^2}{4}-sum_{i=2}^{n}i^2sum_{d=1}^{lfloorfrac{n}{i} floor}d^2varphi(d)

    使用杜教筛将时间复杂度降到O(n23)O(n^{frac{2}{3}})

    数学太难了QAQ

    代码:

    #include<cstdio>
    #include<cmath>
    #include<cstring>
    #include<map>
    #include<algorithm>
    
    #define maxn 5000000
    #define INF 0x3f3f3f3f
    #define MOD 1000000007
    #define two 500000004
    #define six 166666668
    
    using namespace std;
    
    inline long long getint()
    {
    	long long num=0,flag=1;char c;
    	while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;
    	while(c>='0'&&c<='9')num=num*10+c-48,c=getchar();
    	return num*flag;
    }
    
    long long n;
    bool not_prime[maxn+5];
    int prime[maxn+5],cnt;
    long long phi[maxn+5];
    map<long long,long long>M;
    
    inline void init()
    {
    	phi[1]=1;
    	for(int i=2;i<=maxn;i++)
    	{
    		if(!not_prime[i])prime[++cnt]=i,phi[i]=i-1;
    		for(int j=1;j<=cnt&&i*prime[j]<=maxn;j++)
    		{
    			not_prime[i*prime[j]]=1;
    			if(i%prime[j])phi[i*prime[j]]=phi[i]*phi[prime[j]];
    			else{phi[i*prime[j]]=phi[i]*prime[j];break;}
    		}
    	}
    	for(int i=1;i<=maxn;i++)(phi[i]*=1ll*i*i%MOD)%=MOD;
    	for(int i=1;i<=maxn;i++)(phi[i]+=phi[i-1])%=MOD;
    }
    
    inline long long getsqr(long long x)
    {return x%MOD*((x+1)%MOD)%MOD*((2*x+1)%MOD)%MOD*six%MOD;}
    
    inline long long solve(long long x)
    {
    	if(x<=maxn)return phi[x];
    	if(M.count(x))return M[x];
    	long long sum=x%MOD*((x+1)%MOD)%MOD*two%MOD;
    	(sum*=sum)%=MOD;
    	for(long long i=2,j;i<=x;i=j+1)
    	{
    		j=x/(x/i);
    		(sum-=(getsqr(j)-getsqr(i-1))%MOD*solve(x/i)%MOD)%=MOD;
    		(sum+=MOD)%=MOD;
    	}
    	return M[x]=sum;
    }
    
    int main()
    {
    	init();
    	n=getint();
    	long long sum=0;
    	for(long long i=1,j;i<=n;i=j+1)
    	{
    		j=n/(n/i);
    		(sum+=1ll*(j+i)%MOD*(j-i+1)%MOD*two%MOD*solve(n/i)%MOD)%=MOD;
    	}
    	printf("%lld
    ",sum);
    }
    
  • 相关阅读:
    SoapUI 使用笔记
    git 使用笔记(二)
    git 使用笔记(一)
    jquery 拓展
    hdu 1024 Max Sum Plus Plus (DP)
    hdu 2602 Bone Collector (01背包)
    hdu 1688 Sightseeing (最短路径)
    hdu 3191 How Many Paths Are There (次短路径数)
    hdu 2722 Here We Go(relians) Again (最短路径)
    hdu 1596 find the safest road (最短路径)
  • 原文地址:https://www.cnblogs.com/Darknesses/p/12002540.html
Copyright © 2011-2022 走看看