P4213 【模板】杜教筛（Sum）

【模板】杜教筛（Sum）

rqy 姐姐！

Orz _rqy。

Orz Apia Du。

修改了原来不太规范的公式写法。

不想改 int。。。所以我选择调整预处理的大小。

效果拔群。

这里简单把杜教筛的式子推一遍。

已知函数 \(f\) 和 \(g\)。我们要求的是 \(S_f\)。这里 \(S_f(n)=\sum_{i=1}^{n}f(i)\)。

假设我们可以比较轻松地求出 \(S_g\) 和 \(S_{(f\ast g)}\) 及它们的前缀和。

这里 \((f\ast g)(i)=\sum_{xy=i}f(x)g(y)\)（好像看杜教筛的人不会有不知道的。。。）。

先从 \(f*g\) 的前缀和开始推导：

\[\sum_{i=1}^n(f\ast g)(i)=\sum_{i=1}^n\sum_{xy=i}f(x)g(y) \]

然后把 \(y\) 提前枚举：

\[\sum_{i=1}^n(f\ast g)(i)=\sum_{y=1}^ng(y)\sum_{xy\le n}f(x) \]

实际上就是：

\[\sum_{i=1}^n(f\ast g)(i)=\sum_{y=1}^ng(y)\sum_{x=1}^{\lfloor n/y \rfloor}f(x) \]

即：

\[\sum_{i=1}^n(f\ast g)(i)=\sum_{y=1}^ng(y)S_f\left (\left \lfloor\dfrac{n}{y}\right \rfloor\right ) \]

然后我们把右边除 \(g(1)S_f(n)\) 的部分往另一边移：

\[g(1)S_f(n)=\sum_{i=1}^n(f\ast g)(i)-\sum_{y=2}^ng(y)S_f\left (\left \lfloor\dfrac{n}{y}\right \rfloor\right ) \]

一般情况下我们找的 \(g\) 是积性函数并且 \(g(1)=1\)，所以我们可以直接得到 \(S_f(n)\) 的一个表达式。如果不是积性函数把 \(g(1)\) 除过去就完了。

然后我们显然得到了一个可以递推解决的东西，这个东西很像 DP，而且后面那一部分显然可以数论分块。

所以我们可以采取类似 DP 的方式记忆化瞎搞。

主定理计算这个复杂度似乎有一点点玄学。

我们发现函数调用的 \(\lfloor n/y \rfloor\)（设为 \(x\)），分别是 \(1,2,3,\cdots,\lfloor\sqrt{n}\rfloor\) 和 \(\lfloor n / 1\rfloor,\lfloor n / 2\rfloor,\cdots,\lfloor n / \sqrt{n}\rfloor\)。于是我们的复杂度就是求个和：

\[O \left(\sum_{i=1}^{\sqrt{n}}\sqrt{i}+\sum_{i=1}^{\sqrt{n}}\sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor} \right) = O\left(\int_{1}^{\sqrt{n}}x^{1/2}\, \mathrm{d}x+\int_{1}^{\sqrt{n}}\sqrt{\dfrac{n}{x}}\, \mathrm{d}x\right) \]

都是初等函数，直接积分：

\[O\left(\dfrac{2}{3}n^{3/4}-\dfrac{2}{3}+2n^{3/4}-2n^{1/2}\right) = O\left(n^{3/4}\right) \]

也就是说我们成功找到了低于线性的一种方法。

但我们不能满足于此。

发现这个东西我们可以先预处理出 \(S_f\) 在 \(1,2,\cdots,n^c\)（\(c> 1 / 2\)，其实可以尽量大，这里只是图证明省事）的值，于是函数调用的 \(x\) 就是 \(\lfloor n/1\rfloor,\lfloor n/2\rfloor,\cdots,\lfloor n/n^{1-c} \rfloor\)。

我们再求一次和：

\[O\left(n^c+\sum_{i=1}^{n^{1-c}}\sqrt{\left\lfloor\dfrac{n}{i}\right\rfloor}\right) \]

积分近似，就是：

\[O\left(n^c+\int_{1}^{n^{1-c}}\sqrt{\dfrac{n}{x}}\,\mathrm{d}x\right) = O(n^c+n^{1-c/2}) \]

现在我们要求复杂度最小，所以 \(c=2/3\)。此时时间复杂度 \(O(n^{2/3})\)。

现在我们开始着手考虑实现这个想法。

首先这是个递归调用，还需要加记忆化。但是一个很烦人的问题就是这个记忆化的数组可能需要很大。

一个比较 naive 的方法是用 STL map，开了 O2 跑得巨快。但可惜的是理论复杂度会多一个 \(\log\)。

我们发现，其实真正调用的 \(x\)，都是一些形如 \(\lfloor n/ y \rfloor\) 的东西，我们已经预处理了 \(y>n^{1-c}=n^{1/3}\) 的情况，现在调用的全部都是 \(y\le n^{1/3}\) 的情况，也就是说，理论上需要记忆化的数组量不会超过 \(n^{1/3}\)。

让杜教筛复杂度能过的题目，其 \(n\) 基本能满足 \(n^{1/3}\) 个数组够开。

所以我们把记忆化数组的下标设为 \(y\) 即可。每次调用函数的时候我们找到满足 \(x=\lfloor n/y\rfloor\) 的最大的 \(y\)。显然 \(y=\lfloor n/x\rfloor\)（数论分块的一个小结论）。返回函数值的时候记忆化在这个数组中即可。

好了，讨论完杜教筛的一般套路，我们来看一看这道模板题的做法。

首先先线性筛肯定是不用说了。我们直接看这个递归怎么做：

对于求 \(S_{\mu}\)：
- 我们显然知道一个 \(1*\mu=\varepsilon\)。
- 那么设 \(f=\mu\)，\(g=1\)，就会有 \(f*g=\varepsilon\)。
- \(1\) 和 \(\varepsilon\) 的值及其前缀和实在是太好求了。
- 于是很轻松的套上了杜教筛。
给出这一部分的代码：
```
inline ll SMu(ll x) {
  if(x<=M) return SMu1[x];ll tmp=n/x;
  if(vis[tmp]) return SMu2[tmp];
  vis[tmp]=1;SMu2[tmp]=1;
  for(ll i=2,j;i<=x;i=j+1) {
      j=x/(x/i);SMu2[tmp]-=(j-i+1)*SMu(x/i);
  }
  return SMu2[tmp];
}
```
对于求 \(S_{\varphi}\)：
- 首先我们知道 \(\varphi*1=\operatorname{id}\)。
- 我们设 \(f=\varphi\)，\(g=1\)，则 \(f*g=\operatorname{id}\)。
- \(1\) 很好求，\(\operatorname{id}\) 的值很好求，前缀和就是等差数列求和，也非常好求。
- 于是也很轻松地套上了杜教筛。
给出这一部分的代码：
```
inline ll SPhi(ll x) {
  if(x<=M) return SPhi1[x];ll tmp=n/x;
  if(vis[tmp]) return SPhi2[tmp];
  vis[tmp]=1;SPhi2[tmp]=(x+1)*x/2;
  for(ll i=2,j;i<=x;i=j+1) {
      j=x/(x/i);SPhi2[tmp]-=(j-i+1)*SPhi(x/i);
  }
  return SPhi2[tmp];
}
```

好了我们做完了。

那么肯定有人问如果题目没那么模板，\(g\) 和 \(f*g\) 不太能 \(O(1)\) 求怎么办？

对于 \(g\) 来说我实在是没辙，换成其它函数试试吧（因为 \(g\) 要再套杜教筛一类的东西很容易复杂度爆炸）。

对于 \(f\ast g\)，大胆杜教筛没有问题，理论复杂度仍然是 \(O(n^{2/3})\)。

（算是我写得相当用心的 blog 了。。。）

代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;

const ll N=(1ll<<31)-1,M=2e6,L=1.3e3;

ll T,n,cnt;

ll phi[M+5],mu[M+5],SMu1[M+5],SPhi1[M+5],prime[M+5];
ll SMu2[L+5],SPhi2[L+5];
bool vis[L+5],f[M+5];

inline void Init() {
  mu[1]=1;phi[1]=1;f[1]=1;
  for(ll i=2;i<=M;i++) {
    if(!f[i]) {prime[++cnt]=i;mu[i]=-1;phi[i]=i-1;}
    for(ll j=1;j<=cnt&&prime[j]*i<=M;j++) {
      f[i*prime[j]]=1;
      if(i%prime[j]==0) {
        mu[i*prime[j]]=0;phi[i*prime[j]]=phi[i]*prime[j];
        break;
      }
      mu[i*prime[j]]=-mu[i];phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
  }
  for(ll i=1;i<=M;i++) {
    SMu1[i]=SMu1[i-1]+mu[i];SPhi1[i]=SPhi1[i-1]+phi[i];
  }
}

inline ll SMu(ll x) {
  if(x<=M) return SMu1[x];ll tmp=n/x;
  if(vis[tmp]) return SMu2[tmp];
  vis[tmp]=1;SMu2[tmp]=1;
  for(ll i=2,j;i<=x;i=j+1) {
    j=x/(x/i);SMu2[tmp]-=(j-i+1)*SMu(x/i);
  }
  return SMu2[tmp];
}

inline ll SPhi(ll x) {
  if(x<=M) return SPhi1[x];ll tmp=n/x;
  if(vis[tmp]) return SPhi2[tmp];
  vis[tmp]=1;SPhi2[tmp]=(x+1)*x/2;
  for(ll i=2,j;i<=x;i=j+1) {
    j=x/(x/i);SPhi2[tmp]-=(j-i+1)*SPhi(x/i);
  }
  return SPhi2[tmp];
}

inline ll read() {
  ll ret=0,f=1;char ch=getchar();
  while(ch<'0'||ch>'9') {if(ch=='-') f=-f;ch=getchar();}
  while(ch>='0'&&ch<='9') {ret=(ret<<3)+(ret<<1)+ch-'0';ch=getchar();}
  return ret*f;
}

inline void write(ll x) {
  static char buf[22];static ll len=-1;
  if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
  else {putchar('-');do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
  while(len>=0) putchar(buf[len--]);
}

int main() {

  Init();

  T=read();

  while(T--) {
    n=read();
    memset(vis,0,sizeof(vis));
    write(SPhi(n));putchar(' ');
    memset(vis,0,sizeof(vis));
    write(SMu(n));putchar('\n');
  }

  return 0;
}

posted @ 2023-02-06 22:56 Aryper 阅读(25) 评论(0) 编辑收藏举报

刷新页面返回顶部

Aryper.

P4213 【模板】杜教筛（Sum）

公告