• [20191003机房测试] 太阳神


    太阳神拉很喜欢最小公倍数,有一天他想到了一个关于最小公倍数的题目
    求满足如下条件的数对(a,b)对数:
    a,b 均为正整数且 a,b<=n 而lcm(a,b)>n
    其中的 lcm 当然表示最小公倍数
    答案对 1,000,000,007取模
    

    数据范围是1e10的,打表找了半天规律发现没用……
    那就莫比乌斯反演呗

    题目求:

    [sum_{i=1}^{n}sum_{j=1}^{n}[lcm(i,j)> n] ]

    但是大于的太多了,那我们就反过来求,最后用总的来减

    也就是求:

    [sum_{i=1}^{n}sum_{j=1}^{n}[lcm(i,j)leq n] ]

    先把(lcm)拆开

    [sum_{i=1}^{n}sum_{j=1}^{n}[dfrac{i·j}{gcd(i,j)}leq n] ]

    看到(gcd),枚举 (i,j) 的约数 (d)

    [sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{n}[dfrac{i·j}{d}leq n] ]

    拉出去反演:

    [sum_{d=1}^{n}sum_{i'=1}^{lfloor{frac{n}{d}} floor}sum_{j'=1}^{lfloor{frac{n}{i'·d}} floor}[gcd(i',j')=1] ]

    引入莫比乌斯函数:

    [sum_{d=1}^{n}sum_{i'=1}^{lfloor{frac{n}{d}} floor}sum_{j'=1}^{lfloor{frac{n}{i'·d}} floor}mu(gcd(i',j')) ]

    再枚举 (i',j') 的约数 (d'),继续反演

    [sum_{d=1}^{n}sum_{d'=1}^{sqrt{frac{n}{d}}}mu(d')sum_{i''=1}sum_{j''=1}lfloor{frac{n}{d·d'^2i''}} floor ]

    (d')拿出来:

    [sum_{d'=1}^{sqrt{n}}mu(d')[d·i''·j''leq frac{n}{d'^2} ext{的组数}] ]

    然后就是求满足(d·i''·j''leq frac{n}{d'^2})的三元组(d,i'',j'')的个数

    这个可以(Theta(n^{frac{2}{3}}))算出来

    就可以了

    代码:

    #include<bits/stdc++.h>
    #define ll long long
    #define N 100005
    #define M 100000
    #define mod 1000000007
    using namespace std;
    
    ll n,ans,siz,sum,maxn,x;
    ll mu[N],prime[N],primenum;
    bool isprime[N];
    
    template<class T>inline void read(T &res)
    {
    	char c;T flag=1;
    	while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;res=c-'0';
    	while((c=getchar())>='0'&&c<='9')res=res*10+c-'0';res*=flag;
    }
    
    void init()
    {
    	mu[1]=1;
    	for(register int i=2;i<=M;++i)
    	{
    		if(!isprime[i])
    		{
    			prime[++primenum]=i;
    			mu[i]=-1;
    		}
    		for(register int j=1;j<=primenum;++j)
    		{
    			if(i*prime[j]>M) break;
    			isprime[i*prime[j]]=1;
    			if(!(i%prime[j]))
    			{
    				mu[i*prime[j]]=0;
    				break;
    			}
    			mu[i*prime[j]]=-mu[i];
    		}
    	}
    }
    
    int main()
    {
    	freopen("ra.in","r",stdin);
    	freopen("ra.out","w",stdout);
    	read(n);
    	init();
    	siz=(ll)(sqrt(n)+0.5);
    	for(register ll i=1;i<=siz;++i)
    	{
    		sum=0;
    		maxn=n/i/i;
    		for(register ll j=1;j*j*j<=maxn;++j)
    		{ 
    			for(register ll k=j;k*k<=maxn/j;++k)
    			{
    				x=(maxn/j/k-k+1);
    				if(j==k) sum=(sum+1+(x-1)*3)%mod;
    				else sum=(sum+3+(x-1)*6)%mod;
    			}
    		}
    		sum%=mod;
    		ans=(ans+mu[i]*sum)%mod;
    	}
    	ans=((n%mod)*(n%mod)-ans)%mod;
    	while(ans<0) ans+=mod;
    	printf("%lld
    ",ans);
    	return 0;
    }
    /*
    10000000000
    210705255
    */
    
  • 相关阅读:
    Eclipse背景颜色修改
    Android动画效果translate、scale、alpha、rotate详解
    代理上网的方法
    ubuntu系统使用SSH免密码登陆
    Git的思想和基本工作原理
    GitHub详细教程
    Ubuntu和Redhat(Debian)的差别
    T2: 一种能累积计算积分的EC2实例类型
    win server 2008 r2 iis+php 500错误内部服务器错误。
    从OTF字体文件里查找字体名称
  • 原文地址:https://www.cnblogs.com/tqr06/p/11619827.html
Copyright © 2020-2023  润新知