• Luogu6788 「EZEC-3」四月樱花


    2020.8.29更新

    被命题人@了,原来的做法已经过不去了……

    这题是一道很不错的数论题,对整除分块的考察及其时间复杂度分析与优化较为深入。

    算法考察

    整除分块,常见数论函数性质。

    算法分析

    我们来推一推柿子。

    显然(d(y)=sum_{z|y}1),因此(y^{d(y)}=prod_{z|y}y)

    我们把(y)拆开:

    [prod_{z|y}y=prod_{z|y}z imesprod_{z|y}frac{y}{z}=left(prod_{z|y}z ight)^2=prod_{z|y}z^2 ]

    因此原式可推导如下:

    [prod_{x=1}^tprod_{y|x}prod_{z|y}frac{z^2}{(z+1)^2} ]

    改为枚举(y)

    [left(prod_{y=1}^tprod_{z|y}frac{z^2}{(z+1)^2} ight)^{leftlfloorfrac{t}{y} ight floor} ]

    改为枚举(z)

    [left(prod_{z=1}^tfrac{z^2}{(z+1)^2} ight)^{sum_{y=1}^{leftlfloorfrac{t}{z} ight floor}leftlfloorfrac{t}{yz} ight floor} ]

    我们设(f(n)=sum_{i=1}^n leftlfloorfrac{n}{i} ight floor)

    [left(prod_{z=1}^tfrac{z^2}{(z+1)^2} ight)^{fleft(leftlfloorfrac{t}{z} ight floor ight)} ]


    整除分块基础——求解(f(n))

    考虑(f(n))如何求解。(有整除分块基础的同学可以跳过)

    通过观察,我们可以发现(leftlfloorfrac{n}{i} ight floor)最多有(O(sqrt n))种不同的取值,且随(i)增大单调不减。

    单调性显然,取值个数简单证明一下:

    • (ile sqrt n)时,有(sqrt n)个不同的(i),对应至多(sqrt n)个不同的(leftlfloorfrac{n}{i} ight floor)
    • (i>sqrt n)时,发现(leftlfloorfrac{n}{i} ight floor<sqrt n),因此也至多有(sqrt n)个不同的取值。

    因此我们考虑枚举这(sqrt n)种不同的取值即可求解,单次求解的时间复杂度为(O(sqrt n))

    求解(f(n))的代码如下:

    int get_f(int n){
    	int ret=0;
    	for(int l=1,r;l<=n;l=r+1){
    		r=n/(n/l);
    		// l,r表示区间[l,r]内的整数i,floor(n/i)相等(均为 n/l)
    		// l,r 的值的推导本处不展开
    		// 可以参考下面相关题目的题解 
    		ret=(ret+1ll*(r-l+1)*(n/l))%mod;
    		// 区间[l,r]内有(r-l+1)个整数,它们的整除值都是 n/l 
    	}
    	return ret;
    }
    

    相关题目:[CQOI2007]余数求和


    回到原问题,现在我们已经可以用(O(sqrt n))的时间复杂度求出(f(n))的值。

    整除分块进阶——整除分块嵌套

    我们再观察式子:

    [left(prod_{z=1}^tfrac{z^2}{(z+1)^2} ight)^{fleft(color{red}{leftlfloorfrac{t}{z} ight floor} ight)} ]

    我们发现(f(n))的自变量也出现整除,也就意味着会出现一段连续的(z),使得(leftlfloorfrac{t}{z} ight floor)的值相同,即(fleft({leftlfloorfrac{t}{z} ight floor} ight))也一定相同

    因此对于指数相同的部分,我们只需要知道底数的部分积即可。

    观察底数部分的部分积。

    [prod_{z=l}^rfrac{z^2}{(z+1)^2}=frac{l^2}{(r+1)^2} ]

    发现这个值可以配合费马小定理求逆元快速求出。

    类比整除分块的过程,我们枚举(O(sqrt t))个不同的(leftlfloorfrac{t}{z} ight floor)的值,指数部分可以用整除分块求解。

    这就是整除分块套整除分块。

    类比杜教筛,时间复杂度为(O(t^{frac{3}{4}}log t)),前面是整除分块套整除分块的复杂度,即(sum_{i=1}^{sqrt t}O(sqrt i)+O(sqrt{frac{n}{i}})=O(t^{frac{3}{4}})),后面(log t)是过程中求逆元。


    算法优化

    本算法还有进一步优化的空间。

    继续类比杜教筛。在杜教筛中,我们预处理出部分函数值及其前缀和,从而将时间复杂度由(O(n^frac{3}{4}))优化到(O(n^frac{2}{3}))。对于本题,我们能不能也预处理出一部分(f(n)),从而优化时间复杂度呢?

    因数与倍数——(f(n))求法优化

    我们考虑(leftlfloorfrac{a}{b} ight floor)的实际含义,其可以表示不超过(a)的数中有多少个是(b)的倍数,即(sum_{i=1}^a [b|i])

    因此我们把(f(n))改写如下:

    [f(n)=sum_{i=1}^n leftlfloorfrac{n}{i} ight floor ]

    [=sum_{i=1}^n sum_{j=1}^n [i|j] ]

    [=sum_{j=1}^n sum_{i=1}^n [i|j] ]

    [=sum_{j=1}^n sum_{i|j}1 ]

    [=sum_{j=1}^n sigma(j) ]

    其中(sigma(n))表示(n)有多少个因数。

    根据(sigma(n))的定义,我们不难写出单次询问(O(sqrt n))的算法求出单个(sigma(n)),但这个显然不能满足我们的需求。

    现在我们考虑如何批量求出(sigma(n))

    注意到对每一个数(i)都会对(i)的倍数产生(1)的贡献,因此我们再次枚举倍数,即可批量求出(sigma(n)),最后求一遍前缀和即可批量求出(f(n))

    具体实现如下:

    for(int i=1;i<=N;i++){
    	for(int j=1;i*j<=N;j++){
    		sig[i*j]++;
    	}
    }
    

    时间复杂度为(sum_{i=1}^nfrac{n}{i}=O(nlog n)),可用调和级数证明。

    上述过程也可以用狄利克雷卷积解释,即(1 * 1),而求狄利克雷卷积的通用方法正是枚举倍数。

    我们也不难看出(sigma(n))是一个积性函数,并且(sigma(p^k)=k+1),这意味着我们可以线性筛批量求出(sigma(n)),时间复杂度为(O(n))

    相关题目:[AHOI2005]约数研究


    综上,我们可以用(O(nlog n))的时间复杂度求出前(n)(f(n))。假设我们预处理前(M)(f(n))的值,那么总时间复杂度为

    [Mlog M+sqrt tlog t+sum_{i=1}^{t/M}sqrt{frac{t}{i}}approx(2sqrt{frac{t^2}{M}}+M)log t ]

    由不等式相关知识得(M=t^frac{1}{3}),原式取最小值为(O(t^{frac{2}{3}}log t)),实测(M=500000)较优。

    代码实现

    1. 注意(2.5 imes10^9)超出了 int 范围,要使用unsigned int
    2. 注意对指数上的整除分块求和模数应为(varphi(p)=p-1)
    3. 注意常数优化。
    #include<bits/stdc++.h>
    using namespace std;
    #define maxn 1000005
    #define ui unsigned
    template <typename Tp> void read(Tp &x){
    	int fh=1;char c=getchar();x=0;
    	while(c>'9'||c<'0'){if(c=='-'){fh=-1;}c=getchar();}
    	while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+(c&15);c=getchar();}x*=fh;
    }
    ui n,mod;
    ui ans=1;
    ui ksm(ui B,ui P){ui ret=1;while(P){if(P&1)ret=1llu*ret*B%mod;B=1llu*B*B%mod;P>>=1;}return ret;}
    ui f[maxn];
    ui calc(ui x){//整除分块求f(n) 
    	if(x<=1000000)return f[x];
    	ui ret=0;
    	for(ui l=1,r;l<=x;l=r+1){
    		r=x/(x/l);
    		ret=(ret+1ll*(r-l+1)*(x/l))%(mod-1);
    	}
    	return ret;
    }
    void preprocess(int N){//预处理求f(n) 
    	for(int i=1;i<=N;i++){
    		for(int j=1;i*j<=N;j++){
    			f[i*j]++;
    		}
    	}
    	for(int i=1;i<=N;i++)f[i]=(f[i-1]+f[i])%mod;
    }
    signed main(){
    	read(n);read(mod);
    	preprocess(1000000);
    	for(ui l=1,r;l<=n;l=r+1){
    		r=n/(n/l);//整除分块嵌套 
    		ui invr=1llu*l*ksm((r+1)%mod,mod-2)%mod;
    		invr=1llu*invr*invr%mod;
    		ans=1llu*ans*ksm(invr,calc(n/l))%mod;
    	}
    	printf("%u
    ",ans);
    	return 0;
    }
    

    后记

    本题解中附有一些相关题目,其本身可能难度不是很大,但其思想与本题相同。

    对于求解不同的数论问题,枚举因数与枚举倍数各有所长,有时可能要反复转化。

  • 相关阅读:
    sharepoint无法到搜索人员
    sharepoint 自定义操作 CustomAction
    sharepoint 自定列表总结
    向数据表中插入字符序号,例如"001""002"
    通过sharepoint 客户脚本,验证列表添加页面上的时间
    WebService入门(1)
    Java虚拟机1:什么是Java
    SFTP服务搭建
    Shell学习——特殊变量
    js知识点~扁平数据结构转Tree
  • 原文地址:https://www.cnblogs.com/ZigZagKmp/p/13581931.html
Copyright © 2020-2023  润新知