题解 P6788 【「EZEC-3」四月樱花】

· · 题解

本题解同步发表于我的cnblog

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\times\prod_{z|y}\frac{y}{z}=\left(\prod_{z|y}z\right)^2=\prod_{z|y}z^2

因此原式可推导如下：

\prod_{x=1}^t\prod_{y|x}\prod_{z|y}\frac{z^2}{(z+1)^2}

改为枚举y

\left(\prod_{y=1}^t\prod_{z|y}\frac{z^2}{(z+1)^2}\right)^{\left\lfloor\frac{t}{y}\right\rfloor}

改为枚举z

\left(\prod_{z=1}^t\frac{z^2}{(z+1)^2}\right)^{\sum_{y=1}^{\left\lfloor\frac{t}{z}\right\rfloor}\left\lfloor\frac{t}{yz}\right\rfloor}

我们设f(n)=\sum_{i=1}^n \left\lfloor\frac{n}{i}\right\rfloor：

\left(\prod_{z=1}^t\frac{z^2}{(z+1)^2}\right)^{f\left(\left\lfloor\frac{t}{z}\right\rfloor\right)}

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

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

通过观察，我们可以发现\left\lfloor\frac{n}{i}\right\rfloor最多有O(\sqrt n)种不同的取值，且随i增大单调不减。

单调性显然，取值个数简单证明一下：

因此我们考虑枚举这\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}^t\frac{z^2}{(z+1)^2}\right)^{f\left(\color{red}{\left\lfloor\frac{t}{z}\right\rfloor}\right)}

我们发现f(n)的自变量也出现整除，也就意味着会出现一段连续的z，使得\left\lfloor\frac{t}{z}\right\rfloor的值相同，即f\left({\left\lfloor\frac{t}{z}\right\rfloor}\right)也一定相同。

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

观察底数部分的部分积。

\prod_{z=l}^r\frac{z^2}{(z+1)^2}=\frac{l^2}{(r+1)^2}

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

类比整除分块的过程，我们枚举O(\sqrt t)个不同的\left\lfloor\frac{t}{z}\right\rfloor的值，指数部分可以用整除分块求解。

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

类比杜教筛，时间复杂度为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)求法优化

我们考虑\left\lfloor\frac{a}{b}\right\rfloor的实际含义，其可以表示不超过a的数中有多少个是b的倍数，即\sum_{i=1}^a [b|i]。

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

f(n)=\sum_{i=1}^n \left\lfloor\frac{n}{i}\right\rfloor =\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}^n\frac{n}{i}=O(n\log n)，可用调和级数证明。

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

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

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

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

M\log M+\sqrt t\log t+\sum_{i=1}^{t/M}\sqrt{\frac{t}{i}}\approx(2\sqrt{\frac{t^2}{M}}+M)\log t

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

代码实现

  1. 注意2.5\times10^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\n",ans);
    return 0;
}

后记

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

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