投稿
  1. 高性能架构探索首页
  2. C/C++

杜教筛小记

admin • • C/C++ • 阅读 139

个人总结向博客。

杜教筛

用于求解积性函数前缀和。

设定我们当前要求的积性函数为 \(f(x)\) 。

我们需要求 \(\sum_{i=1}^n f(i)\) 设这东西为 \(S(n)\)。

这个时候的 \(n\) 一般很大范围达到 \(10^{10}\) 以上。

这样就不能暴力求了 qaq

杜教筛就是说帮助我们度过这个过程，其实就是一种构造，我们找到另外一个积性函数 \(g(x)\) 。

然后我们将他们执行狄利克雷卷积，也就是 \((g * f)(n) = \sum_{d|n} g(d) * f(\dfrac{n}{d})\) 。

然后我们考虑求一下这个卷积后的函数的前缀和。

\[\sum_{i=1}^n \sum_{d|i} g(d) * f(\dfrac{n}{d})
\]

考虑枚举变量。

\[\sum_{d=1}^n g(d) \sum_{i=1}^{n/d} f(i)=\sum_{d=1}^n g(d) * S(\dfrac{n}{d})
\]

我们发现这个卷积的前缀和居然和所求的 \(S\) 有关，那么考虑通过这个，我们能不能求出 \(S(n)\) 。

发现可以的！考虑容斥！

\[g(1) * S(n) = \sum_{d=1}^n g(d) * S(\dfrac{n}{d}) - \sum_{i=2}^{n} g(d) * S(\dfrac{n}{d})
\]

然后把 \(g(1)\) 除过去就好了/hanx 。

然后我们假设我们现在能够快速的求出狄利克雷卷积的前缀和。

就我们发现后面的东西其实是可以整除分块的，也就是我们只用求出 \(\sqrt{n}\) 个 \(S(\dfrac{n}{i})\) 的值就好了，然后就分治做就好了。

假设计算他的复杂度为 \(T(n)\) 。

那么 \(T(n) = O(\sqrt{n}) + \sum_{i=2}^{\sqrt{n}} T(i) + \sum_{i=2}^{\sqrt{n}} T(\dfrac{n}{i})\)

这个东西我们展开考虑只这一层，因为越往下越小，所以为。

\(T(n) = O(\sqrt{n}) + \sum_{i=2}^{\sqrt{n}} O(\sqrt{i}) + O(\sqrt{\dfrac{n}{i}})\)

通过均值不等式可得：

\(O(\sqrt{i}) + O(\sqrt{\dfrac{n}{i}}) \ge 2\sqrt{\sqrt{n}} = 2n^{\frac{1}{4}}\) 。

然后有 \(\sqrt{n}\) 项，于是为 \(\sqrt{n} \times \ 2n^{\frac{1}{4}} \ge 2n^{\frac{3}{4}}\)。

这就是直接做的复杂度惹，也就是 \(n^{\frac{3}{4}}\)。

考虑进一步分析复杂度，我们先可简化下式子。

将 \(\sqrt{i}\) 和 \(\sqrt{\dfrac{n}{i}}\) 视为同级。

那么复杂度可以看作为 \(\sum_{i=1}^{\sqrt{n}} \sqrt{\dfrac{n}{i}}\) 。

考虑我们还可以先用筛法预处理一个比较小范围的前缀和，假设这个范围为 \(k\) 。

那么这个时候的复杂度为 \(O(k)+O(\sqrt{n})+\sum_{i=2}^{\sqrt{n}} \left[ \sqrt{\dfrac{n}{i}} > k \right] O(\sqrt{\dfrac{n}{i}})\)

所以

\[\sqrt{\dfrac{n}{i}}>k
\]

\[\dfrac{n}{i} > k^2
\]

\[n > k^2 i
\]

\[i<\dfrac{n}{k^2}
\]

那么复杂度为 \(O(k)+O(\sqrt{n})+\sum_{i=2}^{\frac{n}{k^2}} O(\sqrt{\dfrac{n}{i}})\)

那也就是 \(O(k)+O(\sqrt{n})+O(\dfrac{n\sqrt{n}}{k^2})\)

在 \(k=n^{\frac{2}{3}}\) 取得最优复杂度为 \(n^{\frac{2}{3}}\) 。

那杜教筛就是这个东西了，这个东西主要是运用之类的，那么接下来刷几个数论板题。

不过在刷题前，需要牢记一件事狄利克雷卷积时点积没有结合律。

就是 \(\sum_{d|n} f(d)g(d)h(\dfrac{n}{d}) \ne \sum_{d|n} f(d)h(\dfrac{n}{d}) g(d)\)

这个东西可以通过自己卷看出（其实是问神仙 Bindir0 得出的orz）。

莫比乌斯函数之和

定义 \(S(a,b) = \sum_{i=a}^b \mu(i)\) 。

给出 \(a,b\) 求 \(S(a,b)\) 。

\(2 \leq a,b\leq 10^{10}\)

模板题，直接把杜教筛的式子拉过来。

\[g(1)S(n)=\sum_{i=1}^n(g*\mu)(i)-\sum_{i=2}^ng(i)S(\frac{n}{i})
\]

考虑找到一个积性函数 \(g\) 让 \(g\) 和 \(\mu\) 卷积很好算。

想起之前的学莫反时候的结论 \(\mu * I = e\)

那么卷 \(I\) 就是了。

那么那个 \(\sum_{i=1}^n(g*\mu)(i)\) 我们其实可以快速求出来，这东西就是求了一个 \(e\) 的前缀和。但是 \(e\) 的前缀和为 \(1\) /hanx

所以式子可以化简了。

\[S(n) = 1-\sum_{i=2}^n S(\dfrac{n}{i})
\]

然后做杜教筛就好了，为了保证最优复杂度我们提前预处理出一部分的前缀和，然后后面整除分块着做就好了，然后一定要加上记忆化，记忆化我们用一个 unordered_map 来进行，也可以手写 hash。

贴一个代码。

// 德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱
// Problem: P1244 莫比乌斯函数之和
// Contest: 51nod
// URL: https://www.51nod.com/Challenge/Problem.html#problemId=1244
// Memory Limit: 256 MB
// Time Limit: 3000 ms
// The Author : Pitiless0514
// 
// Powered by CP Editor (https://cpeditor.org)

#include <bits/stdc++.h>
#define int long long
using namespace std;
namespace IO {
	int len = 0;
	char ibuf[(1 << 20) + 1], *iS, *iT, out[(1 << 25) + 1];
	#define gh()                                                                   \
	  (iS == iT ? iT = (iS = ibuf) + fread(ibuf, 1, (1 << 20) + 1, stdin),         \
	   (iS == iT ? EOF : *iS++) : *iS++)
	inline int read() {
	  char ch = gh();
	  int x = 0;
	  char t = 0;
	  while (ch < '0' || ch > '9') t |= ch == '-', ch = gh();
	  while (ch >= '0' && ch <= '9') x = x * 10 + (ch ^ 48), ch = gh();
	  return t ? -x : x;
	}
	inline void putc(char ch) { out[len++] = ch; }
	template <class T> inline void write(T x) {
	  if (x < 0) putc('-'), x = -x;
	  if (x > 9) write(x / 10);
	  out[len++] = x % 10 + 48;
	}
	string getstr(void) {
	  string s = "";
	  char c = gh();
	  while (c == ' ' || c == '\n' || c == '\r' || c == '\t' || c == EOF) c = gh();
	  while (!(c == ' ' || c == '\n' || c == '\r' || c == '\t' || c == EOF))s.push_back(c), c = gh();
	  return s;
	}
	void putstr(string str, int begin = 0, int end = -1) {
	  if (end == -1)
	    end = str.size();
	  for (int i = begin; i < end; i++)
	    putc(str[i]);
	  return;
	}
	inline void flush() {
	  fwrite(out, 1, len, stdout);
	  len = 0;
	}
} // namespace IO by Macesuted
using IO::flush;
using IO::getstr;
using IO::putc;
using IO::putstr;
using IO::read;
using IO::write;
const int N = 3e6 + 5;
int a, b, n, B, lim, mu[N], pr[N], tot, vis[N];
unordered_map<int,int>sumu;
void pre(int n) {
  mu[1] = 1;
  int p;
  for(int i = 2; i <= n; i++) {
    if(!vis[i]) {pr[++tot] = i; mu[i] = -1;}
    for(int j = 1; j <= tot; j++) {
      p = i * pr[j];
      if(p > n) break;
      vis[p] = 1;
      mu[p] = -mu[i];
      if(i % pr[j] == 0) { mu[p] = 0; break; }
    }
  }
  for(int i = 2; i <= n; i++) mu[i] += mu[i - 1];
}//提前线性筛预处理
int calcsumu(int n) {
  if(n <= lim) return mu[n]; //预处理好了的一部分直接用
  if(sumu[n]) return sumu[n]; //记忆化
  int ret = 1;
  for(int l = 2, r; l <= n && n / l; l = r + 1) {
    r = n / (n / l) ;
    ret -= (r - l + 1) * calcsumu(n / l);
  }//整除分块
  return sumu[n] = ret; //记忆化
}
signed main () {
  a = read(), b = read();
  pre(lim = 3e6);
  cout << calcsumu(b) - calcsumu(a - 1) << "\n";
  return 0;
}

最小公倍数之和 V3

求 \(\sum_{i=1}^n \sum_{j=1}^n \rm{lcm}(i,j)\)

\(2 \le n \le 10^{10}\)

没有那么板子了，需要一定推式子能力了 /dk

\(\rm{lcm}\) 不好看，转化一下为 \(\dfrac{ij}{\gcd(i,j)}\)

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

考虑枚举 \(\gcd\) 。

\[\sum_{d=1}^n \sum_{i=1}^n \sum_{j=1}^n\dfrac{ij}{d} [(i,j)=d]
\]

将 \(i\) 变为 \(id\) ， 将 \(j\) 变为 \(jd\) 。

\[\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} ijd [(i,j)=1]
\]

\[\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} ij [(i,j)=1]
\]

\[\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} ij \sum_{k|(i,j)} \mu(k)
\]

\[\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} ij \sum_{k|i} \sum_{k|j} \mu(k)
\]

考虑交换枚举顺序，先枚举 \(d\)，后枚举 \(k\) ，最后枚举 \(i,j\) ，然后枚举 \(i,j\) 时直接枚举 \(k\) 的倍数。

\[\sum_{d=1}^n d\sum_{k=1}^{\lfloor \frac{n}{d} \rfloor} \mu(k) \sum_{i=1}^{\lfloor \frac{n}{dk} \rfloor} ik \sum_{j=1}^{\lfloor \frac{n}{dk} \rfloor} jk
\]

这样就变成了这个形式，进一步提取系数。

\[\sum_{d=1}^n d\sum_{k=1}^{\lfloor \frac{n}{d} \rfloor} \mu(k) k^2 \sum_{i=1}^{\lfloor \frac{n}{dk} \rfloor} i \sum_{j=1}^{\lfloor \frac{n}{dk} \rfloor} j
\]

令 \(g(n)=\sum_{i=1}^n\sum_{j=1}^nij=\dfrac{n^2(n+1)^2}{4}\) 继续化简。

考虑改变枚举变量，枚举 \(T=dk\) ，式子变为：

\[\sum_{T=1}^ng(\dfrac{n}{T})\sum_{x|T}x^2\mu(x)*\dfrac{T}{x}
\]

\[\sum_{T=1}^ng(\dfrac{n}{T})\sum_{x|T}xT\mu(x)
\]

发现前面的 \(g\) 可以整除分块，考虑快速求后面东西的前缀和。

令 \(f(n) = \sum_{d|n}dn\mu(d)\)

那么令 \(h(n)=n^2\mu(n)\) ，则 \(f=h*id\)

考虑快速计算 \(h\) ，把这个和 \(id^2\) 函数一卷可以消掉平方项。

观察 \((f*id^2)(n)\) ，也就是：

\[\sum_{d|n}d^2\mu(d)\dfrac{n^2}{d^2}
\]

\[n^2\sum_{d|n}\mu(d)
\]

\[n^2e=e
\]

那么 \(f*id^2=h * id * id^2=(h * id^2) * id = e * id=id\)

设 \(S\) 为 \(f\) 的前缀和，然后直接杜教筛就好了/se ， 拉一个式子过来。

\[id^2(1)S(n)=\sum_{i=1}^n(f*id^2)(i)-\sum_{i=2}^nid^2(i)S(\dfrac{n}{i})
\]

\[id^2(1)S(n)=\sum_{i=1}^ni-\sum_{i=2}^nid^2(i)S(\dfrac{n}{i})
\]

\[S(n)=\dfrac{n(n+1)}{2}-\sum_{i=2}^ni^2S(\dfrac{n}{i})
\]

直接杜教筛就好了。

然后代码应为我不知道我哪里漏取模了，于是我开了 _int128 ，然后手写了一个 unordered_map 的 hash ，然后就很丑，而且常数也很大，不过还是能过，毕竟实现 8s /hanx

// 德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱德丽莎你好可爱
// Problem: ?
// Contest: ?
// URL: ?
// Memory Limit: ?
// Time Limit: ?
// The Author : Pitiless0514
// 
// Powered by CP Editor (https://cpeditor.org)

#include <bits/stdc++.h>
#include<chrono>
#define int __int128
using namespace std;
struct custom_hash {
  static uint64_t splitmix64(uint64_t x) {
    x += 0x9e3779b97f4a7c15;
    x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
    x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
    return x ^ (x >> 31);
  }
  size_t operator()(uint64_t x) const {
    static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
    return splitmix64(x + FIXED_RANDOM);
  }
};
namespace IO {
  int len = 0;
  char ibuf[(1 << 20) + 1], *iS, *iT, out[(1 << 25) + 1];
  #define gh()                                                                   \
    (iS == iT ? iT = (iS = ibuf) + fread(ibuf, 1, (1 << 20) + 1, stdin),         \
    (iS == iT ? EOF : *iS++) : *iS++)
  inline int read() {
    char ch = gh();
    int x = 0;
    char t = 0;
    while (ch < '0' || ch > '9') t |= ch == '-', ch = gh();
    while (ch >= '0' && ch <= '9') x = x * 10 + (ch ^ 48), ch = gh();
    return t ? -x : x;
  }
  inline void putc(char ch) { out[len++] = ch; }
  template <class T> inline void write(T x) {
    if (x < 0) putc('-'), x = -x;
    if (x > 9) write(x / 10);
    out[len++] = x % 10 + 48;
  }
  string getstr(void) {
    string s = "";
    char c = gh();
    while (c == ' ' || c == '\n' || c == '\r' || c == '\t' || c == EOF) c = gh();
    while (!(c == ' ' || c == '\n' || c == '\r' || c == '\t' || c == EOF))s.push_back(c), c = gh();
    return s;
  }
  void putstr(string str, int begin = 0, int end = -1) {
    if (end == -1) end = str.size();
    for (int i = begin; i < end; i++) putc(str[i]);
    return;
  }
  inline void flush() {
    fwrite(out, 1, len, stdout);
    len = 0;
  }
} // namespace IO by Macesuted
using IO::flush;
using IO::getstr;
using IO::putc;
using IO::putstr;
using IO::read;
using IO::write;
const int N = 4e6, mod = 1000000007; 
int n, pr[N], vis[N], tot, lim, f[N], sto, ans;
unordered_map<int,int,custom_hash>s;
int power(int a,int b) {
  int ans = 1;
  while(b) {
    if(b & 1) ans = ans * a % mod;
    b >>= 1;
    a = a * a % mod;
  }
  return ans;
}
const int inv2 = power(2, mod - 2), inv6 = power(6, mod - 2);
void pre(int n) {
  f[1] = 1;
  int p;
  for(int i = 2; i <= n; i++) {
    if(!vis[i]) { pr[++tot] = i; f[i] = 1 - i + mod; }
    for(int j = 1; j <= tot; j++) {
      p = i * pr[j];
      if(p > n) break;
      vis[p] = 1;
      f[p] = f[i] * (1 - pr[j] + mod) % mod;
      if(i % pr[j] == 0) { f[p] = f[i]; break; }
    }
  }
  for(int i = 1; i <= n; i++) f[i] = f[i] * i + f[i - 1], f[i] %= mod;
}
int sum(int n) {
  n %= mod;
  return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}
int g(int n) {
  n %= mod;
  int s = (n + 1) * n % mod * inv2 % mod;
  return s * s % mod;
}
int calcS(int n) {
  if(n <= lim) return f[n];
  if(s[n]) return s[n];
  int ret = n * (n + 1) % mod * inv2 % mod;
  for(int l = 2, r; l <= n; l = r + 1) {
    r = n / (n / l);
    r = min(r, n);
    ret -= (sum(r) - sum(l - 1) + mod) % mod * calcS(n / l) % mod;
    ret %= mod; ret = (ret + mod) % mod;
  }
  return s[n] = ret;
}
signed main () {
  n = read();
  pre(lim = N - 1);
  for(int l = 1, r, val; l <= n && n / l; l = r + 1) {
    r = n / (n / l), val = n / l;
    r = min(r, n);
    sto = calcS(r) - calcS(l - 1); sto = (sto + mod) % mod;
    ans += g(val) * sto % mod; ans = (ans + mod) % mod;
  }
  write(ans), putc('\n'), flush();
  return 0;
}

原文链接: https://www.cnblogs.com/ptno/p/15759750.html

欢迎关注

微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍

    杜教筛小记

原创文章受到原创版权保护。转载请注明出处：https://www.ccppcoding.com/archives/183235

非原创文章文中已经注明原地址，如有侵权，联系删除

关注公众号【高性能架构探索】，第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处！

赞 (0)
adminadmin
0 0
上一篇 2023年2月12日 上午10:21
下一篇 2023年2月12日 上午10:21

相关推荐