• LOJ2085. 「NOI2016」循环之美


    定义一个(frac{x}{y})是“纯循环小数”,当且仅当:化成小数之后,可以写成(a.dot{c_1}c_2c_3c_4dotsdot{c_p})的形式。

    (xin[1,n])(yin[1,m]),有多少对((x,y),x perp y)满足(frac{x}{y})(k)进制下是纯循环小数。

    (n,mle 10^9,kle 2000)


    分析条件:如果是,当且仅当(exist p>0,xequiv k^pxpmod y)。也就是(kperp y)

    所以题目也就是要求(sum_{i=1}^msum_{j=1}^n [iperp j][iperp k])

    考虑拆式子。

    说在前面:作为数论白痴,下面估的时间复杂度往往是不准确的。如果有大佬路过希望指教一下。

    当拆([i perp j])

    经过推导得(sum_d [dperp k]mu(d)lfloorfrac{n}{d} floorsum_{i=1}^{lfloorfrac{m}{d} floor}[iperp k])

    (f(n)=sum_{i=1}^n[iperp k]=lfloorfrac{n}{k} floorphi(k)+f(nmod k))。可以快速求,不是复杂度瓶颈。

    对于(sum_d [dperp k]mu(d)lfloorfrac{n}{d} floor f(lfloorfrac{m}{d} floor))整除分块,现在要求前缀和,即(g(n,k)=sum_{d=1}^n mu(d)[dperp k])

    一种推法

    ([dperp k])拆开推一波式子可以得到(g(n,k)=sum_{t|k}mu(t)sum_{d=1}^{lfloorfrac{n}{t} floor}mu(dt)=sum_{t|k}mu(t)^2sum_{d=1}^{lfloorfrac{n}{t} floor}mu(d)[dperp t]=sum_{t|k}mu(t)^2g(lfloorfrac{n}{t} floor,t))。整除分块搞下去。

    边界是(g(n,1))(g(0,k))。其中(g(n,1))可以杜教筛求。

    时间复杂度(记忆化):(O(n^{2/3}+(k的因子个数)sqrt n))

    另一种推法

    (p)(k)的一个质因数,则(k)可以表示成(p^tq)的形式。

    (g(n,k)=sum_{i=1}^n[iperp q]mu(i)-sum_{p|i}[iperp q]mu(i)=g(n,q)+g(lfloorfrac{n}{p^t} floor,k))

    边界同上,时间复杂度(记忆化):(O(n^{2/3}+(k的质因子个数)sqrt n))

    当拆([iperp k])

    (f(m,n,k)=sum_{i=1}^msum_{j=1}^n [iperp j][iperp k])

    经过推导得(f(m,n,k)=sum_{d|k}mu(d)f(n,lfloorfrac{m}{d} floor,d))

    边界条件(k=1),化一下得(sum_d mu(d)lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor)。整除分块,杜教筛算(mu)的前缀和。

    这个东西不太会分析,但是实践表明递归的次数不是很多(没有对(f)进行记忆化的情况下)。杜教筛记忆化(lfloorfrac{n}{d} floor)(lfloorfrac{m}{d} floor)这样的前缀和,总时间(O(n^{2/3}+m^{2/3}))。整除分块的时间复杂度大概是(O(sqrt{max(n,m)})。总的时间复杂度不会算……

    代码是最后一种做法。


    using namespace std;
    #include <bits/stdc++.h>
    #define ll long long
    //#define ll __int128
    #define N 6000005
    int gcd(int a,int b){
    	int k;
    	while (b)
    		k=a%b,a=b,b=k;
    	return a;
    }
    int n,m,k;
    int p[N],np;
    bool inp[N];
    int mu[N],mnp[N];
    int smu[N];
    void init(int n=6000000){
    	mu[1]=1;
    	for (int i=2;i<=n;++i){
    		if (!inp[i])
    			p[++np]=i,mnp[i]=i,mu[i]=-1;
    		for (int j=1;j<=np && i*p[j]<=n;++j){
    			inp[i*p[j]]=1;
    			mnp[i*p[j]]=p[j];
    			if (i%p[j]==0)
    				break;
    			mu[i*p[j]]=-mu[i];
    		}
    	}
    	for (int i=1;i<=n;++i)
    		smu[i]=smu[i-1]+mu[i];
    }
    const int EMP=-1926081700;
    struct memory{
    	int n;
    	ll v0[N],v1[N];
    	void init(int _n){
    		n=_n;
    		for (int i=1;i*i<=n;++i)
    			v0[i]=v1[i]=EMP;
    	}
    	ll &operator[](int x){return (ll)x*x<=n?v0[x]:v1[n/x];}
    } bzn,bzm;
    ll S1(int n,memory &bzn){
    	if (n<=6000000)
    		return smu[n];
    	if (bzn[n]!=EMP)
    		return bzn[n];
    	ll sum=0;
    	for (int i=2;i<=n;){
    		int d=n/i,j=n/d;
    		sum+=(ll)(j-i+1)*S1(d,bzn);
    		i=j+1;
    	}
    	return bzn[n]=1-sum;
    }
    ll S2(int m,int n,int rev){
    	if (rev)
    		swap(m,n);
    	ll sum=0,pre=0;
    	for (int i=1;i<=m && i<=n;){
    		int dn=n/i,dm=m/i,j=min(n/dn,m/dm);
    		ll tmp=S1(j,n/dn<m/dm?bzn:bzm);
    		sum+=(ll)(tmp-pre)*dn*dm;
    		pre=tmp;
    		i=j+1;
    	}
    	return sum;
    }
    ll f(int m,int n,int k,int rev){
    	if (m==0 || n==0)
    		return 0;
    	ll s=0;
    	if (k==1){
    		s+=S2(m,n,rev);
    		return s;
    	}
    	int i=1;
    	for (;i*i<k;++i)
    		if (k%i==0){
    			if (mu[i]) 
    				s+=mu[i]*f(n,m/i,i,rev^1);
    			if (mu[k/i]) 
    				s+=mu[k/i]*f(n,m/(k/i),k/i,rev^1);
    		}
    	if (i*i==k && mu[i])
    		s+=mu[i]*f(n,m/i,i,rev^1);
    	return s;
    }
    int main(){
    	freopen("in.txt","r",stdin);
    //	freopen("out.txt","w",stdout);
    	init();
    	scanf("%d%d%d",&n,&m,&k);
    	bzn.init(n);
    	bzm.init(m);
    	printf("%lld
    ",(long long)f(m,n,k,0));
    //	printf("%lld
    ",cnt);
    	return 0;
    }
    
  • 相关阅读:
    容器操作--管理迭代器
    顺序容器--添加及访问元素
    日志记录-20151103
    顺序容器--容器库.迭代器
    使用-flat.vmdk恢复虚拟机
    H3C-交换机维护命令大全
    Centos6.5 安装zabbix-agent 3.0
    Linux系统调试工具之sysdig使用详解
    通过实例学习 tcpdump 命令
    系统之锹sysdig:Linux服务器监控和排障利器
  • 原文地址:https://www.cnblogs.com/jz-597/p/14310609.html
Copyright © 2020-2023  润新知