「笔记」莫比乌斯反演

Updated on 2020.8.6
巨幅更新,对积性函数和狄利克雷卷积部分进行重构。
新增对一类特殊求和式的讲解。

Updated on 2020.9.9
添加了几个例题。


前置知识

小碎骨

艾佛森括号 \([P] = \begin{cases} 1 &\text{If P is true}\\ 0 &\text{Otherwise} \end{cases}\)
此处 \(P\) 是一个可真可假的命题。

引理1

\[\forall a,b,c\in \mathbb{Z},\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c}}\right\rfloor \]

证明

\[\dfrac{a}{b} = \left\lfloor{\dfrac{a}{b}}\right\rfloor + r(0\le r < 1) \]

\[\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor\dfrac{a}{b}\times\dfrac{1}{c}\right\rfloor = \left\lfloor{\dfrac{1}{c}\times\left({\left\lfloor{\dfrac{a}{b} + r}\right\rfloor}\right)}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c} + \dfrac{r}{c}}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor{\dfrac{a}{b}}\right\rfloor}{c}}\right\rfloor \]

数论分块

内容独立了出来,详细内容见 数论分块 - Luckyblock

对于一类含有\(\left\lfloor\frac{n}{i}\right\rfloor\)的求和式 (\(n\) 为常数),由于\(\left\lfloor\frac{n}{i}\right\rfloor\)单调不增,故存在多个区间\([l,r]\), 使得\(\left\lfloor\frac{n}{i}\right\rfloor = \left\lfloor\frac{n}{j}\right\rfloor(i,j\in [l,r])\) 成立。

对于任意一个\(i\),最大的满足上式的 \(j=\left\lfloor{\dfrac{n}{\left\lfloor{\dfrac{n}{i}}\right\rfloor}}\right\rfloor\)


积性函数

定义

\(\gcd(x,y) = 1\)\(f(xy)=f(x)f(y)\),则 \(f(n)\) 为积性函数。

性质

\(f(x)\)\(g(x)\)均为积性函数,则以下函数也为积性函数:

\[\begin{aligned} & h(x) = f(x^p)\\ & h(x) = f^p(x)\\ & h(x) = f(x)g(x)\\ & h(x) = \sum_{d\mid x} f(d)g\left(\dfrac{x}{d}\right) \end{aligned}\]

常见积性函数

  • 单位函数 \(e(n) = [n = 1]\)
  • 幂函数 \(\operatorname{Id}_{k}(n) = n^k\)\(\operatorname{id}_1(n)\) 通常简记为\(\operatorname{id}(n)\)
  • 常数函数 \(1(n) = 1\)
  • 因数个数 \(\operatorname{d}(n) = \sum\limits_{d\mid n} 1\)
  • 除数函数 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)
    \(k=1\) 时为因数和函数,通常简记为 \(\sigma(n)\)
    \(k=0\) 时为因数个数函数 \(\sigma_{0}(n)\)
  • 欧拉函数 \(\varphi(n) = \sum\limits_{i=1}^{n} [\gcd(i,n) = 1]\)
  • 莫比乌斯函数 \(\mu(n) = \begin{cases}1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{为}\ n\ \text{的本质不同质因子个数} \end{cases}\)

不是很懂上面写的什么玩意?
不用深究,有个印象继续往下看就好。


莫比乌斯函数

定义

\(\mu\) 为莫比乌斯函数,定义为

\[\mu(n) = \begin{cases} 1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{为}\ n\ \text{的本质不同质因子个数} \end{cases}\]

解释

\(n = \prod\limits_{i=1}^{k} p_{i}^{c_i}\),其中\(p_i\)为质因子,\(c_i\ge 1\)

  1. \(n=1\)时,\(\mu (n) = 1\)
  2. \(n\not ={1}\)时 ,
    • \(\exist i\in [1,k], c_i > 1\) 时,\(\mu (n) = 0\)
      当某质因子出现次数大于\(1\)时,\(\mu (n) = 0\)
    • \(\forall i\in [1,k], c_i = 1\) 时,\(\mu (n) = (-1)^k\)
      当每个质因子只出现一次时,即\(n = \prod\limits_{i=1}^{k}p_i\)\(\{p_i\}\)中元素唯一。
      \(\mu (n) = (-1)^k\),此处\(k\)为质因子的种类数。

性质

莫比乌斯函数是积性函数,且具有以下性质

\[\large \sum_{d\mid n} \mu (d) = [n=1] \]

证明,设 \(n = \prod\limits_{i=1}^{k}{p_i^{c_i}}, n' = \prod\limits_{i=1}^{k}{p_i}\)

  • 根据莫比乌斯函数定义,则有:\(\sum\limits_{d\mid n}{\mu(d)} = \sum\limits_{d\mid n'}{\mu(d)}\)
  • \(n'\) 的某因子 \(d\),有 \(\mu (d) = (-1)^i\),则它由 \(i\) 个 本质不同的质因子组成。
    由于质因子总数为 \(k\),满足上式的因子数为 \(C_{k}^{i}\)
  • 对于原求和式,转为枚举 \(\mu(d)\) 的值。
    \(\sum\limits_{d\mid n'}{\mu(d)} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i\times 1^{k-i}}\)
    根据二项式定理,上式 \(= (1+(-1))^k\)
    易知该式在 \(k=0\),即 \(n=1\) 时为 \(1\),否则为 \(0\)

反演常用结论

一个反演常用结论:

\[[\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)} \]

证明 1:
\(n = \gcd(i,j)\),则右\(= \sum\limits_{d\mid n} {\mu (d)} = [n = 1] = [\gcd(i,j) = 1]=\) 左。

证明 2:
暴力展开:\([\gcd(i,j) = 1] = e(\gcd(i,j)) = \sum\limits_{d\mid \gcd(i,j)}\mu(d)\)

线性筛求莫比乌斯函数

\(\mu\) 为积性函数,因此可以线性筛莫比乌斯函数。

int pnum, mu[kMaxn], p[kMaxn];
bool vis[kMaxn];

void Euler(int lim_) {
  vis[1] = true, mu[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = - mu[i];
    }
  }
}

狄利克雷(Dirichlet)卷积

建议阅读 算法学习笔记(35): 狄利克雷卷积 By: Pecco

定义两个数论函数 \(f,g\) 的狄利克雷卷积为

\[\large(f\ast g) (n) = \sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right) \]

性质

  1. 显然满足 交换律,结合律,分配律。
    • \(f \ast g = g \ast f\)
    • \((f \ast g) \ast h = f\ast (g\ast h)\)
    • \(f\ast (g+h) = f\ast g + f\ast h\)
  2. \(e\) 为狄利克雷卷积的单位元,有\((f\ast e)(n) = f(n)\)
  3. \(f,g\) 为积性函数,则 \(f\ast g\) 为积性函数。

关于单位元 \(e\)

有:

\[e = \mu \ast 1=\sum\limits_{d\mid n} \mu (d) \]

证明

\[\begin{aligned} (f\ast e)(n) = \sum_{d\mid n} f(d)e(\dfrac{n}{d}) = \sum_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] \end{aligned}\]

  • 对于\([\dfrac{n}{d} = 1]\),当且仅当 \(\dfrac{n}{d}=1\),即 \(d=n\) 时为 \(1\),否则为\(0\)
  • 则当 \(d=n\) 时,\(f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\)
    \(d\not ={n}\) 时,\(f(d)\left[\dfrac{n}{d} = 1\right] = 0\)

综上,\((f\ast e)(n) = \sum\limits_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\),满足单位元的性质。
\(e = \mu \ast 1\) 成立。

除数函数与幂函数

幂函数 \(\operatorname{Id}_{k}(n) = n^k\)
除数函数 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)

显然有:

\[(\operatorname{Id}_k\ast 1)(n) = \sum_{d\mid n} \operatorname{Id_k}(d) = \sum_{d\mid n} d^k = \sigma_k(n) \]

\(k=0\) 时,\(\operatorname{Id_0}=1\)\(\sigma_0\) 为因数个数函数,有:

\[(1\ast 1)(n) = \sum_{d\mid n}1 = \sigma_0(n) \]

欧拉函数与恒等函数

\[\begin{aligned} \varphi \ast 1 =& \operatorname{Id}\\ \varphi =& \operatorname{Id}\ast \mu \end{aligned}\]

对于一式,当 \(n=p^m\) 时(\(p\) 为质数),有:

\[(\varphi \ast 1)(p^m) = \sum_{d\mid n}\varphi(d) = \varphi(1) +\sum_{i=1}^{m}\varphi(p^i) = 1 +\sum_{i=1}^{m}(p^i-p^{i-1}) = p^m \]

\(p^i\) 的因子有 \(p^{i-1}\) 个,为 \(1\sim p^{i-1}\),故 \(\varphi(p^i) = p^i-p^{i-1}\)

\((\varphi \ast 1)(n)\) 为积性函数,则对于任意正整数 \(n\),有:

\[(\varphi \ast 1)(n) = (\varphi \ast 1)\left(\prod p^m\right) = \prod\left(\varphi \ast 1\right)(p^m) = \prod p^m = n \]

得证。

对于 2 式,在 1 式基础上两侧同时 \(\ast \mu\) 即得。
左侧变为 \(\varphi \ast 1\ast \mu = \varphi \ast e = \varphi\)

计算

\[(f\ast g) (n) = \sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right) \]

考虑枚举 \(n\) 的因子,将贡献累加即可。
显然可以使用埃氏筛筛出所有前缀狄利克雷卷积,复杂度 \(O(nk\log n)\),其中 \(k\) 是计算函数值的复杂度。


莫比乌斯反演

反演是个啥?反演

公式

\(f(n),g(n)\)为两个数论函数。
如果有

\[\large f(n) = (g\ast 1)(n) = \sum\limits_{d\mid n}{g(d)} \]

那么有

\[\large g(n) = (\mu \ast f)(n)=\sum\limits_{d\mid n} {\mu(d)f\left(\dfrac{n}{d}\right)} \]

证明

方法一:运用卷积。

原问题为:已知 \(f = g\ast 1\),证明 \(g = f\ast \mu\)
易知如下转化:\(f\ast \mu = g\ast 1 \ast \mu \Longrightarrow f\ast \mu = g\ast e = g\)

方法二:对原式进行数论变换。

  1. \(\sum\limits_{d\mid n}g(d)\) 替换\(f\left(\dfrac{n}{d}\right)\)

    \[\sum_{d\mid n}{\mu(d)\sum_{k\mid \frac{n}{d}}g(k)} \]

  2. 变换求和顺序。

    \[\sum_{k\mid n}g(k)\sum_{d\mid \frac{n}{k}}{\mu(d)} \]

  3. 依据 \(\sum\limits_{d\mid n}{\mu(d)} = [n=1]\),仅当 \(\dfrac{n}{k} = 1\) 时,\(\sum\limits_{d\mid \frac{n}{k}}{\mu(d)} = 1\),否则为 \(0\)
    此时\(k=n\),故上式等价于 \(\sum\limits_{k\mid n} {[n=k]\cdot g(k)} = g(n)\)

举例

狄利克雷(Dirichlet)卷积 部分可以知道一些积性函数的关系。
尝试对它们进行反演:

\[e = \mu \ast 1\iff\mu = \mu\ast e = \mu \]

\[\sigma_k = \operatorname{Id}_k\ast 1\iff \operatorname{Id}_k = \mu \ast \sigma_k \]

\[\operatorname{Id}=\varphi \ast 1\iff \varphi = \operatorname{Id}\ast \mu \]


关于一类求和式

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

一般套路

考虑构造出函数 \(g\),满足下列形式:

\[f(n) = g\ast 1 = \sum_{d\mid n}g(d) \]

则求和式变为:

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}g(d) \]

考虑算术基本定理,发现若满足 \(d\mid \gcd (i,j)\),则 \(d\mid i\)\(d\mid j\) 成立。
考虑 \(g(d)\) 在何时做出贡献,调整枚举顺序:

\[\sum_{d=1}g(d)\sum_{i=1}^n[d\mid i]\sum_{j=1}^m[d\mid j] \]

\(\sum\limits_{i=1}^{n}[d\mid i]\) 等价于 \(1\sim n\)\(d\) 的倍数的个数,则上式等价于:

\[\sum_{d=1}g(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \]

数论分块求解即可。

例 1

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

发现此题的 \(f\) 等价于 \(\operatorname{Id}\),则上式等价于:

\[\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{Id}(\gcd(i,j))\\ =& \sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}\varphi(d)\\ =& \sum_{d=1}\varphi(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \end{aligned}\]

例 2

\[\begin{aligned} &\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}e(\gcd(i,j))\\ =& \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum_{d\mid \gcd(i,j)}\mu (d)\\ =& \sum_{d=1}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \end{aligned}\]

感悟

卷点什么东西,把 \(g\) 卷出来。
\(g\) 不一定是特殊意义的函数。


例题

[HAOI2011] Problem b

\(n\) 组询问,每次给定参数 \(a,b,c,d,k\),求:

\[\sum\limits_{x=a}^{b}\sum\limits_{y=c}^{d}[\gcd(x,y) = k] \]

\(1\le n,k,a,b,c,d\le 5\times 10^4\)\(a\le b,c\le d\)

\(f(n,m) = \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = k]\)
根据容斥原理,则原式等价于 \(f(b,d) - f(a-1,d) - f(b,d-1) + f(a-1,d-1)\)
\(f\) 变成了上述一类求和式的形式,考虑化简 \(f\)

易知原式等价于

\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\gcd(i,j) = 1] \]

代入反演常用结论 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),上式化为

\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum_{d\mid \gcd(i,j)}{\mu(d)} \]

变换求和顺序,先枚举\(d\mid \gcd(i,j)\),可得

\[\sum_{d=1}\mu(d)\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}[d\mid i]\sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[d\mid j] \]

对于上式后半的解释:当\(d\mid i\)\(d\mid j\) 时,\(d\mid \gcd(i,j)\)

易知\(1\sim \left\lfloor\dfrac{n}{k}\right\rfloor\)\(d\) 的倍数有 \(\left\lfloor\dfrac{\left\lfloor\dfrac{n}{k}\right\rfloor}{d}\right\rfloor = \left\lfloor\dfrac{n}{kd}\right\rfloor\) 个(由引理 1 可知),原式变为

\[\sum_{d=1}\mu(d)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \]

预处理 \(\mu\) 后,显然可以数论分块求解,复杂度\(O(n + T\sqrt{n})\)


代码

//知识点:莫比乌斯反演
/*
//By:Luckyblock
*/
#include <cstdio>
#include <cctype>
#include <algorithm>
#define ll long long
const int MARX = 6e4 + 10;
//=============================================================
int N, a, b, c, d, k;
int cnt, Mobius[MARX], Prime[MARX], Sum[MARX];
bool vis[MARX];
//=============================================================
inline int read()
{
    int f = 1, w = 0; char ch = getchar();
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
    for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
    return f * w;
}
void GetMobius()
{
    Mobius[1] = 1;
    int MAX = MARX - 10;
    for(int i = 2; i <= MAX; i ++)
    {
      if(! vis[i]) Mobius[i] = - 1, Prime[++ cnt] = i;
      for(int j = 1; j <= cnt && i * Prime[j] <= MAX; j ++)
      {
        vis[i * Prime[j]] = true;
        if(i % Prime[j] == 0) break;
        Mobius[i * Prime[j]] = - Mobius[i];
      }
    }
    for(int i = 1; i <= MAX; i ++) Sum[i] = Sum[i - 1] + Mobius[i];
}
ll Calc(int x, int y)
{
    ll ans = 0ll; int max = std :: min(x, y);
    for(int l = 1, r; l <= max; l = r + 1)
      r = std :: min(x / (x / l), y / (y / l)),
      ans += (1ll * x / (1ll * l * k)) * (1ll * y / (1ll * l * k)) * (Sum[r] - Sum[l - 1]);
    return ans;
}
ll Solve()
{
    a = read(), b = read(), c = read(), d = read(), k = read();
    return Calc(b, d) - Calc(b, c - 1) - Calc(a - 1, d) + Calc(a - 1, c - 1);
}
//=============================================================
int main()
{ 
    N = read(); GetMobius();
    while(N --) printf("%lld\n", Solve());
    return 0;
}

[国家集训队]Crash的数字表格

给定 \(n,m\) , 求:

\[\sum_{i=1}^n\sum_{j=1}^{m} \operatorname{lcm}(i,j)\bmod 20101009 \]

\(1\le n,m\le 10^7\)

易知原式等价于:

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

考虑枚举 \(\gcd(i,j)\),设枚举量为 \(d\)
\(d=\gcd(i,j)\) 的充要条件是满足 \(d|i, d|j\)\(\gcd(\dfrac{i}{d},\dfrac{j}{d}) = 1\),则原式等价于:

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

先枚举 \(d\),则原式等价于:

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

这个 \(d\) 很烦人,把 \(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\)\(\frac{j}{d}\)
消去 \(d\mid i\)\(d\mid j\) 的限定条件,则原式等价于:

\[\begin{aligned} &\sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]\dfrac{id\times jd}{d}\\ =& \sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]ijd\\ =& \sum_{d=1}d \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]ij \end{aligned}\]

单独考虑后半部分,设 \(f(x,y) = \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\)
发现 \(f(x,y)\) 的形式与 [HAOI2011] Problem b 中的式子类似,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\) 套路一波:

\[\begin{aligned} f(x,y) =& \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\\ =& \sum_{i=1}^{x} \sum_{j=1}^{y}\sum_{d\mid \gcd(i,j)}\mu(d) ij\\ =& \sum_{d=1}\mu(d)\sum_{i=1}^{x}[d\mid i] \sum_{j=1}^{y}[d\mid j] ij\\ =& \sum_{d=1}\mu(d) d^2\sum_{i=1}^{\left\lfloor \frac{x}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{y}{d}\right\rfloor}ij \end{aligned}\]

前半部分 \(\sum\limits_{d=1}\mu(d) d^2\),可以考虑筛出 \(\mu(d)\) 后求前缀和。
后半部分是等差数列乘等差数列的形式,设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算:

\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]

则对于 \(f(x,y)\),有:

\[f(x,y) = \sum_{d=1}\mu(d) d^2\cdot g(\left\lfloor \frac{x}{d}\right\rfloor, \left\lfloor\frac{y}{d}\right\rfloor) \]

数论分块求解即可。

再看回原式,原式等价于:

\[\sum_{d=1}d\cdot f(\left\lfloor\frac{n}{d}\right\rfloor, \left\lfloor\frac{m}{d}\right\rfloor) \]

又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块,复杂度 \(O(n + m)\),瓶颈是线性筛。


一些注意的点

处理时会出现求平方的运算,需要特别注意取模问题,ll 都会爆,被坑惨了。

在预处理前缀和的这个地方:

sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod

注意先对平方取模,在乘上 \(\mu\),否则就会爆掉。
以及可以仅令 \(mu + mod\)

以及这个地方:

int g(int n_, int m_) {
  return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}

平方计算,注意随时取模。

代码

//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void Getmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
  vis[1] = true, mu[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = - mu[i];
    }
  }
  sum[1] = 1ll;
  for (int i = 1; i <= lim_; ++ i) {
    sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
  }
}
int g(int n_, int m_) {
  return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
  int lim = std :: min(n_, m_), ret = 0;
  for (int l = 1, r; l <= lim; l = r + 1) {
    r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
    ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
  }
  return ret;
}
int Sum(ll l, ll r) {
  return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() { 
  int n = read(), m = read();
  int lim = std :: min(n, m), ans = 0;
  Euler(lim);
  for (int l = 1, r; l <= lim; l = r + 1) {
    r = std :: min(n / (n / l), m / (m / l));
    ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
  }
  printf("%d", ans);
  return 0;
}
/*
7718820 8445880
*/

SP5971 LCMSUM - LCM Sum

\(T\) 次询问,每次询问给定 \(n\),求:

\[\sum_{i=1}^{n}\operatorname{lcm}(i,n) \]

\(1<T\le 3\times 10^5\)\(1\le n\le 10^6\)

法一:无脑暴力

先拆 \(\operatorname{lcm}\),原式等价于:

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

套路的枚举 \(\gcd(i,n)\),调换枚举顺序,原式等价于:

\[\begin{aligned} &n\sum_{i=1}\sum_{d|i} [\gcd(i,n) = d]\dfrac{i}{d}\\ =& n\sum_{i=1}\sum_{d|i} [\gcd(\dfrac{i}{d},\dfrac{n}{d}) = 1]\dfrac{i}{d}\\ =& n\sum_{d|n}\sum_{i=1}^{n}[d|i][\gcd(\dfrac{i}{d},\dfrac{n}{d}) = 1]\dfrac{i}{d} \end{aligned}\]

\(i,n\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\),消去整除的条件,原式等价于:

\[n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\dfrac{n}{d}) = 1]i \]

代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:

\[n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} i \sum_{t|\gcd(i,\frac{n}{d})}\mu (t) \]

值得注意的是 \(t\) 的上界为 \(\frac{n}{d}\)\(dt\le n\)
调换枚举顺序,先枚举 \(t\),原式等价于:

\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t) \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} [t|i] i \]

套路地消去整除的条件,把 \(i\) 中的 \(t\) 提出来,原式等价于:

\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t)t \sum_{i=1}^{\left\lfloor\frac{n}{dt}\right\rfloor} i \]

对于最后的一个求和项,设 \(g(x) = \sum\limits_{i=1}^{x}i = \frac{x(x+1)}{2}\),显然可以 \(O(1)\) 求解,原式等价于:

\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t)t\cdot g(\left\lfloor\dfrac{n}{dt}\right\rfloor) \]

考虑枚举 \(T = dt\),显然 \(T\le n\)
\(\mu(t)t\)\(d\) 无关,可以直接考虑枚举 \(t|T\),原式等价于:

\[n\sum_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)\sum_{t|T}\mu(t)t \]

前半块是一个数论分块的形式,可以 \(O(\sqrt{n})\) 求解。
考虑后半块,设 \(f(n)=\sum\limits_{d|n}\mu(d)d\),发现它是一个积性函数,可以线性筛筛出,具体地:

\[f(n)= \begin{cases} 1-n &n\in \mathrm{primes} \\ f(\frac{x}{p}) &p^2\mid n\\ f(\frac{x}{p})f(p) &p^2\nmid n \end{cases} \]

其中 \(p\)\(n\) 的最小质因子。

此时已经可以线性筛 + 数论分块求解,复杂度 \(O(n+T\sqrt{n})\),比较菜鸡,时限 500ms 过不了。
考虑筛出 \(f\) 后再用埃氏筛预处理 \(\sum\limits_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)f(T)\),输出时乘上 \(n\),复杂度变为 \(O(n\log^2 n + n)\)


法二:

同样先拆 \(\operatorname{lcm}\),枚举 \(\gcd(i,n)\),调换枚举顺序,原式等价于:

\[\begin{aligned} &n\sum_{i=1}^{n}\dfrac{i}{\gcd(i,n)}\\ =& n\sum_{i=1}\sum_{d|i} [\gcd(i,n) = d]\dfrac{i}{d}\\ =& n\sum_{d|n}\sum_{i=1}^{n}[d|i][\gcd({i},{n}) = d]\dfrac{i}{d} \end{aligned}\]

\(i,n\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\),消去整除的条件,原式等价于:

\[\begin{aligned} &n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(id,n) = d]i\\ =& n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\dfrac{n}{d}) = 1]i \end{aligned}\]

调整枚举对象,上式等价于:

\[n\sum_{d|n}\sum_{i=1}^{d}[\gcd(i,d) = 1]i \]

考虑 \(\sum\limits_{i=1}^{d}[\gcd(i,d) = 1]i\) 的实际意义,表示 \([1,d]\) 中与 \(d\) 互质的数的和。
\(d>1\) 时,与 \(d\) 互质的数总是成对存在,即若 \(\gcd(i,d)=1\) 成立,则 \(\gcd(d-i,d)=1\) 成立。
每对这样的数的和为 \(d\),共有 \(\frac{\varphi(d)}{2}\) 对这样的数。
则原式等价于:

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

可以直接预处理答案。
预处理时先线性筛出 \(\varphi\),再埃氏筛枚举 \(i\) 的倍数,令它们的答案加上 \(\frac{\varphi(i)i}{2}\),最后输出时乘上 \(n\)
复杂度 \(O(n\log^2 n + T)\)


法二代码

//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 1e6;
//=============================================================
ll phi[kMaxn + 10], ans[kMaxn + 10];
int pnum, p[kMaxn + 10];
bool flag[kMaxn + 10];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void GetMax(int &fir, int sec) {
  if (sec > fir) fir = sec;
}
void GetMin(int &fir, int sec) {
  if (sec < fir) fir = sec;
}
void GetPrime() {
  phi[1] = 1, flag[1] = true; //注意初始化
  for (int i = 2; i <= kMaxn; ++ i) {
    if (! flag[i]) {
      p[++ pnum] = i;
      phi[i] = i - 1ll;
    }
    for (int j = 1; j <= pnum && i * p[j] <= kMaxn; ++ j) {
      flag[i * p[j]] = true;
      if (i % p[j]) {
        phi[i * p[j]] = phi[i] * phi[p[j]];
      } else {
        phi[i * p[j]] = phi[i] * p[j];
        break;
      }
    }
  }
  for (int i = 1; i <= kMaxn; ++ i) {
    for (int j = 1; i * j <= kMaxn; ++ j) {
      ans[i * j] += (i == 1 ? 1 : 1ll * phi[i] * i / 2);
    }
  }
}
//=============================================================
int main() { 
  GetPrime();
  int T = read();
  while (T --) {
    int n = read();
    printf("%lld\n", 1ll * ans[n] * n);
  }
  return 0; 
}

[SDOI2015]约数个数和

\(T\) 次询问,每次询问给定 \(n,m\)
定义 >\(\operatorname{d}(i)\)\(i\) 的约数个数,求:

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

\(1<T,n\le 5\times 10^4\)

一个结论:

\[\operatorname{d}(ij) = \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1] \]

证明

先考虑 \(i = p^a\)\(j=p^b(p\in \mathrm{primes})\) 的情况,有:

\[\operatorname{d}(p^{a+b})=\sum_{x|p^a}\sum_{y|p^b}[\gcd(x,y)=1] \]

对于等式左侧,\(p^{a+b}\) 的约数个数为 \(a+b+1\)
对于等式右侧,在保证 \(\gcd(x,y)=1\) 成立的情况下,有贡献的数对 \((x,y)\) 只能是下列三种形式:

  • \(x>0,y-0\)\(x\)\(a\) 种取值方法。
  • \(x=0,y>0\)\(y\)\(b\) 种取值方法。
  • \(x=0,y=0\)

则等式右侧贡献次数为 \(a+b+1\) 次,等于 \(p^{a+b}\) 的约数个数。
则当 \(i = p^a\)\(j=p^b(p\in \mathrm{primes})\) 时等式成立。

又不同质因子间相互独立,上述情况可拓展到一般的情况。


\(\operatorname{d}(i,j)\) 进行化简,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:

\[\begin{aligned} \operatorname{d}(ij) =& \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]\\ =& \sum_{x|i}\sum_{y|j}\sum\limits_{d\mid \gcd(x,y)} {\mu (d)} \end{aligned}\]

调换枚举顺序,先枚举 \(d\),原式等价于:

\[\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|i}[d|x]\sum_{y|j}[d|y] \]

把各项中的 \(d\) 提出来,消去整除的条件,原式等价于:

\[\begin{aligned} &\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|\frac{i}{d}}\sum_{y|\frac{j}{d}}1\\ =& \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \end{aligned}\]


\(\operatorname{d}(ij)\) 代回原式,原式等价于:

\[\sum_{i=1}^{n}\sum_{j=1}^m \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]

调换枚举顺序,先枚举 \(d\),原式等价于:

\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{n}[d|i]\sum_{j=1}^m [d|j]\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]

\(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}, \frac{j}{d}\),消去整除的条件,原式等价于:

\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\operatorname{d}(i)\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\operatorname{d}(j) \]

考虑预处理 \(S(x) = \sum\limits_{i=1}^{x}\operatorname{d}(i)\),则原式等价于:

\[\sum_{d=1}{\mu (d)}S\left({\left\lfloor\frac{n}{d}\right\rfloor}\right)S\left({\left\lfloor\frac{m}{d}\right\rfloor}\right) \]

线性筛预处理 \(\mu,\operatorname{d}\),数论分块求解即可,复杂度 \(O(n+T\sqrt{n})\)


代码

//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void Getmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
  vis[1] = true;
  mu[1] = d[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
      num[i] = 1;
      d[i] = 2;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        num[i * p[j]] = num[i] + 1;
        d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
        break;
      }
      mu[i * p[j]] = - mu[i];
      num[i * p[j]] = 1;
      d[i * p[j]] = 2ll * d[i]; //

    }
  }
  for (int i = 1; i <= lim_; ++ i) {
    summu[i] = summu[i - 1] + mu[i];
    sumd[i] = sumd[i - 1] + d[i];
  }
}

//=============================================================
int main() { 
  Euler(kMaxn - 10);
  int T = read();
  while (T --) {
    int n = read(), m = read(), lim = std :: min(n, m);
    ll ans = 0ll;
    for (int l = 1, r; l <= lim; l = r + 1) {
      r = std :: min(n / (n / l), m / (m / l));
      ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
    }
    printf("%lld\n", ans);
  }
  return 0;
}
/*
in
1
32 43
*/
/*
out
15420
*/

P3768 简单的数学题

给定参数 \(n\)\(p\),求:

\[\left(\sum_{i=1}^n\sum_{j=1}^ni\cdot j\cdot \gcd(i,j)\right)\bmod p \]

\(n\leq10^{10}\)\(5\times10^8\leq p\leq1.1\times10^9\)\(p\in \mathrm{primes}\)
时限 4s。

无脑套路暴力。

考虑先枚举 \(\gcd(i,j)\),原式等价于:

\[\begin{aligned} &\sum_{d=1}d\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)=d]ij\\ =& \sum_{d=1}d\sum_{i=1}^{n}[d|i]\sum_{j=1}^{n}[d|j][\gcd(\dfrac{i}{d},\dfrac{j}{d})=1]ij \end{aligned}\]

提出 \(d\),变为枚举 \(\frac{i}{d}\)\(\frac{j}{d}\),消去整除的条件,原式等价于:

\[\begin{aligned} &\sum_{d=1}d\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]id\cdot jd\\ =& \sum_{d=1} d^3\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]ij \end{aligned}\]

代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:

\[\sum_{d=1} d^3\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}ij\sum_{t|\gcd(i,j)}\mu (t) \]

值得注意的是 \(t\) 的上界为 \(\frac{n}{d}\)\(dt\le n\)
调换枚举顺序,先枚举 \(t\),原式等价于:

\[\sum_{d=1} d^3\sum_{t=1}\mu(t)\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[t|i] \sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[t|j]ij \]

和上面一样,提出 \(t\),套路地消去整除的条件,原式等价于:

\[\sum_{d=1} d^3\sum_{t=1}\mu(t)t^2\sum_{i=1}^{\left\lfloor\frac{n}{dt}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{dt}\right\rfloor}ij \]

发现后面两个求和是等差数列乘等差数列的形式。
\(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算:

\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]

代入原式,原式等价于:

\[\sum_{d=1} d^3\sum_{t=1}\mu(t)t^2\cdot g\left({\left\lfloor\frac{n}{dt}\right\rfloor},{\left\lfloor\frac{n}{dt}\right\rfloor}\right) \]

考虑枚举 \(T = dt\),显然 \(T\le n\)
再考虑枚举 \(d|T\),即可得到 \(t = \frac{T}{d}\),原式等价于:

\[\begin{aligned} &\sum_{T=1}^{n}g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right)\sum_{d|T}d^3\mu{\left(\dfrac{T}{d}\right)}\left(\dfrac{T}{d}\right)^2\\ =& \sum_{T=1}^{n}T^2 g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right)\sum_{d|T}d\cdot \mu{\left(\dfrac{T}{d}\right)} \end{aligned}\]

对于后面这一坨,用 \(\sum\limits_{d|T}d\cdot \mu{\left(\frac{T}{d}\right)} = \operatorname{Id} \ast \mu(T)= \varphi(T)\) 反演,则原式等价于:

\[\sum_{T=1}^{n}T^2 \varphi(T) \cdot g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right) \]

后半块可以数论分块,考虑前半块。
发现前半段即为 \(\operatorname{Id}^2(T)\varphi(T)\),又是前缀和形式,考虑杜教筛。

有:

\[f(n) = \operatorname{Id}^2\varphi(n) \]

考虑找到一个函数 \(g\),构造函数 \(h = f\ast g\) 使其便于求值,有:

\[h(n) = \sum_{d|n} d^2\varphi(d)\cdot g\left(\dfrac{n}{d}\right) \]

看到同时存在 \(d\)\(\frac{n}{d}\),考虑把 \(d^2\) 消去。
\(g = \operatorname{Id}^2\),有:

\[\begin{aligned} h(n) =& \sum_{d|n} d^2\varphi(d)\cdot \left(\dfrac{n}{d}\right)^2\\ =& n^2\sum_{d|n} \varphi(d)\\ =& n^2 \cdot \varphi \ast 1(n)\\ \end{aligned}\]

\(\varphi \ast 1 = \operatorname{Id}\),则有:

\[h(n) = n^3 \]

找到了合适的 \(g\),套杜教筛的公式。

\[\begin{aligned} g(1)S(n) &= \sum_{i=1}^{n}h(i) - \sum_{i=2}^{n}g(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\\ S(n) &= \sum_{i=1}^{n} i^3 - \sum_{i=2}^{n} i^2\cdot S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right) \end{aligned}\]

前一项是自然数的立方和,有 \(\sum\limits_{i=1}^{n} i^3 = (\frac{n(n+1)}{2})^2\)。证明详见:自然数前n项平方和、立方和公式及证明 - 百度文库
后一项直接等差数列求和 + 数论分块求解即可。


「SDOI2017」数字表格

\(f_{i}\) 表示斐波那契数列的第 \(i\) 项。
\(T\) 组数据,每次给定参数 \(n,m\),求:

\[\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)} \pmod {10^9 + 7} \]

\(1\le T\le 10^3\)\(1\le n,m\le 10^6\)
5S,256MB。

以下钦定 \(n\ge m\)
大力化式子,先套路地枚举 \(\gcd(i,j)\),用初中知识把两个 \(\prod\) 化到指数位置,原式等于:

\[\large\begin{aligned} &\prod_{d = 1}^{n}\prod_{i=1}^{n}\prod_{j=1}^{m}f_{d}[\gcd(i,j)=d]\\ =&\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=d]\right)} \end{aligned}\]

分母套路一波,有:

\[\begin{aligned} &\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \end{aligned}\]

代回原式,原式等于:

\[\large \begin{aligned} &\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor\right)}\\ =& \prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor} \end{aligned}\]

考虑再暴力拆一波,原式等于:

\[\large \begin{aligned} &\prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\\ =& \prod_{d = 1}^{n}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f_{d}^{\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\right) \end{aligned}\]

做不动了,但发现变量仅有 \(k,d,kd\),考虑更换枚举对象改为枚举 \(t = kd\)\(d\),则原式等于:

\[\large\prod_{t=1}^{n}\left(\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}}\right)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]

枚举对象变成了约数形式。从后面的式子推前面的式子是比较显然的,可以认为这种枚举 \(t=kd\) 的形式是一种逆向操作。

设:

\[\large g(t)=\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}} \]

\(g(t)\) 可以用类似埃氏筛的方法 \(O(n\log ^2 n)\) 地预处理出来。再把 \(g\) 代回原式,原式等于:

\[\large\prod_{t=1}^{n}g(t)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]

可以考虑预处理 \(g(t)\) 的前缀积,数论分块枚举指数求解即可。

总时间复杂度 \(O(n\log ^2 n + T\sqrt n)\),轻微卡常可以过。

//知识点:莫比乌斯反演 
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define LL long long
const LL mod = 1e9 + 7;
const int kN = 1e6;
//=============================================================
LL n, m, ans;
int p_num, p[kN + 10];
bool vis[kN + 10];
LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
LL invf[kN + 10], invp[kN];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) {
    w = (w << 3) + (w << 1) + (ch ^ '0');
  }
  return f * w;
}
void Chkmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
  x_ %= mod;
  LL ret = 1;
  for (; y_; y_ >>= 1ll) {
    if (y_ & 1) ret = ret * x_ % mod;
    x_ = x_ * x_ % mod;
  }
  return ret;
}
void Euler() {
  vis[1] = true, mu[1] = 1;
  for (int i = 2; i <= kN; ++ i) {
    if (! vis[i]) {
      p[++ p_num] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
}
void Prepare() {
  g[1] = g[2] = 1;
  f[1] = f[2] = 1;
  invf[1] = invf[2] = 1;
  for (int i = 3; i <= kN; ++ i) {
    g[i] = 1;
    f[i] = (f[i - 1] + f[i - 2]) % mod;
    invf[i] = QPow(f[i], mod - 2);
  }
  
  Euler();
  for (int d = 1; d <= kN; ++ d) {
    for (int j = 1; d * j <= kN; ++ j) {
      if (mu[j] == 1) {
        g[d * j] = g[d * j] * f[d] % mod;
      } else if (mu[j] == -1) {
        g[d * j] = g[d * j] * invf[d] % mod;
      }
    }
  }
  invp[0] = prod[0] = 1;
  for (int i = 1; i <= kN; ++ i) {
    prod[i] = prod[i - 1] * g[i] % mod;
    invp[i] = QPow(prod[i], mod - 2);
  }
}
//=============================================================
int main() {
  Prepare();
  int T = read();
  while (T -- ){
    n = read(), m = read(), ans = 1;
    if (n < m) std::swap(n, m);
    for (LL l = 1, r = 1; l <= m; l = r + 1) {
      r = std::min(n / (n / l), m / (m / l));
      ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
    }
    printf("%lld\n", ans);
  }
  return 0;
}

P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题

给定参数 \(p\),有 \(T\) 组数据,每次给定参数 \(A,B,C\),求:

\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,k)}\right)^{f(type)} \]

其中 \(f(type)\) 的取值如下:

\[f(type) = \begin{cases} 1 &type = 0\\ i\times j\times k &type = 1\\ \gcd(i,j,k) &type = 2 \end{cases}\]

\(1\le A,B,C\le 10^5\)\(10^7\le p\le 1.05\times 10^9\)\(p\in \mathbb{P}\)\(T=70\)
2.5S,128MB。

先化下原式,原式等于:

\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{i\times j }{\gcd(i,j)\times \gcd(i,k)}\right)^{f(type)} \]

发现每一项仅与两个变量有关,设:

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{f(type)}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{f(type)} \end{aligned}\]

发现 \(\prod\) 可以随意交换,则原式等价于:

\[\dfrac{f_1(a,b,c)\times f_1(b,a,c)}{f_2(a,b,c)\times f_2(a,c,b)} \]

考虑在 \(type\) 取值不同时,如何快速求得 \(f_1\)\(f_2\)
一共有 \(6\) 个需要推导的式子,不妨就叫它们 \(1\sim 6\) 面了(


type = 0

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j) \end{aligned}\]

对于 1 面,显然有:

\[\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i = \left(\prod_{i=1}^{a}i\right)^{b\times c} \]

预处理阶乘 + 快速幂即可,单次计算时间复杂度 \(O(\log n)\)


再考虑 2 面,套路地枚举 \(\gcd\),显然有:

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)\\ =&\left(\prod_{i=1}^{a}\prod_{j=1}^{b}\gcd(i,j)\right)^c\\ =& \left(\prod_{d=1} d^{\left(\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]\right)}\right)^c \end{aligned}\]

指数是个套路,可以看这里:P3455 [POI2007]ZAP-Queries。于是有:

\[\begin{aligned} &\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{a}{kd}\right\rfloor\left\lfloor\dfrac{b}{kd}\right\rfloor \end{aligned}\]

代回原式,略做处理,则原式等于:

\[\large \begin{aligned} &\left(\prod_{d=1} d^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)^c\\ =& \left(\prod_{d=1} \left(d^{\sum\limits_{k=1}\mu(k)}\right)^{\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor}\right)^c\\ =& \prod_{d=1} \left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{\left(\mu(k)\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)^c \end{aligned}\]

「SDOI2017」数字表格 一样,考虑枚举 \(t=kd\)\(d\),则原式等于:

\[\large \prod_{t=1}^{n}\left(\left(\prod_{d|t} d^{\mu{\left(\frac{t}{d}\right)}}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\right)^c \]

设:

\[\large g_0(t) = \prod_{d|t}d^{\mu\left(\frac{t}{d}\right)} \]

线性筛预处理 \(\mu\) 后,\(g_0(t)\) 可以用埃氏筛预处理,时间复杂度 \(O(n\log n)\)。再代回原式,原式等于:

\[\large \prod_{t=1}^{a}\left(g_0(t)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\right)^c \]

预处理 \(g_0(t)\) 的前缀积和前缀积的逆元,复杂度 \(O(n\log n)\)
数论分块 + 快速幂计算即可,单次时间复杂度 \(O(\sqrt n\log n)\)


type = 1

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{i\times j\times k}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{i\times j\times k} \end{aligned}\]

考虑 3 面,把 \(\prod k\) 扔到指数位置,有:

\[\large \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{i\times j\times k} = \prod_{i=1}^{a}\prod_{j=1}^{b}i^{\left(i\times j\times \sum\limits_{k = 1}^{c} k\right)} \]

再把 \(\prod j\) 也扔到指数位置,引入 \(\operatorname{sum}(n) = \sum_{i=1}^{n} i = \frac{n(n+1)}{2}\),原式等于:

\[\left(\prod_{i=1}^{a}i^i\right)^{\operatorname{sum}(b)\times \operatorname{sum}(c)} \]

预处理 \(i^i\) 的前缀积,复杂度 \(O(n\log n)\)
指数可以 \(O(1)\) 算出,然后快速幂,单次时间复杂度 \(O(\log n)\)

根据费马小定理,指数需要对 \(p - 1\) 取模。注意 \(p-1\) 不是质数,计算 \(\operatorname{sum}\) 时不能用逆元,但乘不爆 LL,直接算就行。


再考虑 4 面,发现 \(k\)\(\gcd\) 无关,则同样把 \(\prod k\) 扔到指数位置,则有:

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{i\times j\times k}\\ =& \left(\prod_{i=1}^a\prod_{j=1}^b\gcd(i,j)^{i\times j}\right)^{\operatorname{sum}(c)} \end{aligned}\]

套路地枚举 \(\gcd\),原式等于:

\[\large \left(\prod_{d=1}d^{\left(\sum\limits_{i=1}^a \sum\limits_{j=1}^b i\times j[\gcd(i,j)=d]\right)}\right)^{\operatorname{sum}(c)} \]

大力化简指数,有:

\[\large \begin{aligned} &\sum\limits_{i=1}^a \sum\limits_{j=1}^b i\times j[\gcd(i,j)=d]\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} i\times j[\gcd(i,j)=1\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j\sum\limits_{t|\gcd(i,j)}\mu(t)\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j\sum\limits_{k|\gcd(i,j)}\mu(k)\\ =& d^2 \sum\limits_{k=1}\mu(k)\sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i[k|i] \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j[k|j]\\ =& d^2 \sum\limits_{k=1}k^2\mu(k)\sum\limits_{i=1}^{\left\lfloor\frac{a}{kd}\right\rfloor} i\sum\limits_{j=1}^{\left\lfloor\frac{b}{kd}\right\rfloor} j\\ =& d^2\sum\limits_{k=1}k^2\mu(k)\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\\ \end{aligned}\]

指数化不动了,代回原式,原式等于:

\[\large \left(\prod_{d=1}d^{\left(d^2\sum\limits_{k=1}k^2\mu(k)\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)}\right)^{\operatorname{sum}(c)} \]

同 2 面的情况,先展开一下,再枚举 \(t=kd\)\(d\),原式等于:

\[\large \begin{aligned} &\left(\prod_{d=1}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{\left(d^2 k^2\mu(k)\right)}\right)^{\left(\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)}\right)^{\operatorname{sum}(c)}\\ =& \prod_{t=1}\left(\left(\prod_{d|t}d^{\left(d^2\left(\frac{t}{d}\right)^2\mu\left(\frac{t}{d}\right)\right)}\right)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)}\\ =& \prod_{t=1}\left(\left(\prod_{d|t}d^{\left(t^2\mu\left(\frac{t}{d}\right)\right)}\right)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)} \end{aligned}\]

与二面相同,设:

\[\large g_1(t) = \prod_{d|t}d^{\left(t^2\mu\left(\frac{t}{d}\right)\right)} \]

\(g_1(t)\) 可以用埃氏筛套快速幂筛出,时间复杂度 \(O(n\log^2 n)\)。再代回原式,原式等于:

\[\prod_{t=1}\left(g_1(t)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)} \]

同样预处理 \(g_1(t)\) 的前缀积及其逆元,时间复杂度 \(O(n\log n)\)
整除分块 + 快速幂即可,单次时间复杂度 \(O(\sqrt n\log n)\)

注意指数的取模。


type = 2

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{\gcd(i,j,k)}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{\gcd(i,j,k)} \end{aligned}\]

考虑 5 面,手段同上,大力反演化简一波,再调换枚举对象,则有:

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{\gcd(i,j,k)}\\ =&\prod_{d=1}\prod\limits_{i=1}^{a}i^{\left(\sum\limits_{j=1}^{b}\sum\limits_{k=1}^{c}[\gcd(i,j,k)=d]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{a}i^{\left(\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[\gcd(\frac{i}{d},j,k)=1]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[\gcd(i,j,k)=1]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}\sum\limits_{x|\gcd(i,j,k)}{\mu(x)}\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{x=1}\mu(x)[x|i]\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[x|j]\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[x|k]\right)}\\ =& \prod_{d=1}\prod_{x=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\times \mu(x)[x|i]\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[x|j]\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[x|k]\right)}\\ =& \prod_{d=1}\prod_{x=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{xd}\right\rfloor}(ixd)^{\left(d\times \mu(x){\left\lfloor\frac{b}{xd}\right\rfloor}{\left\lfloor\frac{c}{xd}\right\rfloor}\right)}\\ =& \prod_{t = 1}\prod_{d|T}\prod_{i=1}^{\left\lfloor\frac{a}{t}\right\rfloor}(it)^{\left(d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}\right)}\\ =& \prod_{t = 1}\prod_{d|T}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\prod_{i=1}^{\left\lfloor\frac{a}{t}\right\rfloor}i\right)^{d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ \end{aligned}\]

引入 \(\operatorname{fac}(n) = \prod_{i=1}^{n} i\),再根据枚举对象调整一下指数,原式等于:

\[\large \begin{aligned} &\prod_{t = 1}\prod_{d|t}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\left(d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}\right)}\\ =& \prod_{t = 1}\left(\prod_{d|t}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{d\times \mu\left(\frac{t}{d}\right)}\right)^{{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ =& \prod_{t = 1}\left(\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\sum\limits_{d|t}d\times \mu\left(\frac{t}{d}\right)}\right)^{{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}} \end{aligned}\]

指数中出现了一个经典的狄利克雷卷积的形式,对其进行反演。
\((\operatorname{Id}\ast \mu) (n)= \varphi (n)\) 代入原式,原式等于:

\[\large \begin{aligned} &\prod_{t = 1}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\varphi(t){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ =& \prod_{t = 1}\left(t^{\varphi(t)\left\lfloor\frac{a}{t}\right\rfloor{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)^{\varphi(t){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\right) \end{aligned}\]

预处理 \(t^{\varphi(t)}\) 的前缀积及逆元,阶乘的前缀积及阶乘逆元,\(\pmod {p-1}\) 下的 \(\varphi(t)\) 的前缀和(指数
),时间复杂度 \(O(n\log n)\)
同样整除分块 + 快速幂即可,单次时间复杂度 \(O(\sqrt n\log n)\)


然后是最掉 sans 的 6 面。有 \(\gcd(i,j,k) = \gcd(\gcd(i,j), k)\),考虑先枚举 \(\gcd(i,j)\),然后套路化式子,则有:

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{\gcd(i,j,k)}\\ =& \prod_{d=1}\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} [\gcd(i,j)=d] d^{\gcd(d,k)}\\ =& \prod_{d=1} \left(d^{\left(\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b} [\gcd(i,j)=d]\right)}\right)^{\sum\limits_{k=1}^{c}\gcd(d,k)} \end{aligned}\]

先考虑最外面的指数,这也是个套路,可以参考 一个例子。用 \(\operatorname{Id} = \varphi \ast 1\) 反演,显然有:

\[\large \begin{aligned} &\sum\limits_{k=1}^{c}\gcd(d,k)\\ =& \sum\limits_{k=1}^{c}\sum_{x|\gcd(d,k)}\varphi(x)\\ =& \sum_{x=1}\varphi(x)[x|d]\sum_{k=1}^{c}[x|k]\\ =& \sum_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor \end{aligned}\]

再考虑里面的指数,发现这式子在 2 面已经推了一遍了,于是直接拿过来用,有:

\[\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]=\sum_{y=1}\mu(y)\left\lfloor\dfrac{a}{yd}\right\rfloor\left\lfloor\dfrac{b}{yd}\right\rfloor \]

将化简后的两个指数代入原式,原式等于:

\[\large \begin{aligned} &\prod_{d=1} \left(d^{\left(\sum\limits_{y=1}\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor}\\ =& \prod_{d=1} \left(\prod\limits_{y=1}d^{\left(\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \end{aligned}\]

与 2、4 面同样套路地,考虑枚举 \(t=yd\)\(d\),再略作调整,原式等于:

\[\large \begin{aligned} &\prod_{d=1} \left(\prod\limits_{y=1}d^{\left(\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor}\\ =& \prod_{t=1}\prod_{d|t} d^{\left(\mu\left(\frac{t}{d}\right)\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\\ =& \prod_{t=1}\left(\prod_{d|t} d^{\left(\mu\left(\frac{t}{d}\right)\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{t=1}\left(\prod_{d|t} \prod_{x|d}d^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor} \end{aligned}\]

发现要同时枚举 \(d\)\(x\),化不动了。
从题解里学到一个比较神的技巧,考虑把 \(d\) 拆成 \(x\)\(\frac{d}{x}\) 分别计算贡献再相乘,即分别计算下两式:

\[\large \begin{aligned} &\prod_{t=1}\left(\prod_{d|t} \prod_{x|d}x^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ &\prod_{t=1}\left(\prod_{d|t} \prod_{x|d}{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor} \end{aligned}\]


先考虑 \(x\) 的情况,首先把枚举 \(x\) 调整到最外层。设 \(\operatorname{lim}=\max(a,b,c)\),则原式等于:

\[\large \begin{aligned} &\prod_{x=1} \prod_{t=1}^{\operatorname{lim}}[x|t]\left(\prod_{d|t} [x|d]{x}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|t} {x}^{\left(\mu\left(\frac{tx}{dx}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\prod_{d|t} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\mu\left(\frac{t}{d}\right)\right)} \end{aligned}\]

\(\prod {t}\) 挪到指数位置,原式等于:

\[\large \begin{aligned} &\prod_{x=1} {x}^{\left(\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\sum\limits_{d|t}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\mu\left(\frac{t}{d}\right)\right)}\\ =& \prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\sum\limits_{d|t}\mu\left(\frac{t}{d}\right)\right)} \end{aligned}\]

指数中又出现了一个经典的狄利克雷卷积的形式,对其进行反演。
\((\mu \ast 1) (n)= \epsilon (n)=[n=1]\) 代入原式,原式等于:

\[\large \begin{aligned} &\prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor[t=1]\right)}\\ =& \prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{a}{x}\right\rfloor\left\lfloor\frac{b}{x}\right\rfloor\left\lfloor\frac{c}{x}\right\rfloor\right)} \end{aligned}\]

得到了一个非常优美的式子,而且发现这个式子是 5 面最终答案的一部分。同 5 面的做法,直接整除分块即可。


再考虑 \(\frac{d}{x}\) 的情况,同上先把枚举 \(x\) 放到最外层,并调整一下指数,则原式等于:

\[\large \begin{aligned} &\prod_{x=1} \prod_{t=1}^{\operatorname{lim}}[x|t]\left(\prod_{d|t} [x|d]{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|tx} [x|d]{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{tx}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\\ =& \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|tx} [x|d]{\left(\frac{d}{x}\right)}^{\mu\left(\frac{tx}{d}\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \end{aligned}\]

考虑枚举 \(dx\),替换原来的 \(d\),注意一下这里的倍数关系。原式等于:

\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]

发现最内层的式子 \(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\),即为二面处理过的 \(g_0(t)\)。直接代入,原式等于:

\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}g_0(t)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]

一个小结论,证明可以看 这里

\[\forall a,b,c\in \mathbb{Z},\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c}}\right\rfloor \]

则原式等于:

\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}g_0(t)^{\left\lfloor\frac{\left\lfloor\frac{a}{x}\right\rfloor}{t}\right\rfloor\left\lfloor\frac{\left\lfloor\frac{b}{x}\right\rfloor}{t}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]

于是可以先对外层整除分块,再对内层整除分块。

然后就做完了,哈哈哈。


一些实现上的小技巧:

  • 逆元能预处理就预处理。
  • 注意对指数取模,模数为 \(p-1\)
//知识点:莫比乌斯反演 
/*
By:Luckyblock
用了比较清晰易懂的变量名,阅读体验应该会比较好。  
vsc 的自动补全真是太好啦!
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
using std::min;
using std::max;
#define LL long long
const int Lim = 1e5;
const int kN = 1e5 + 10;
//=============================================================
LL A, B, C, mod, ans;
int T, p_num, p[kN];
bool vis[kN];
LL mu[kN], phi[kN], fac[kN], g[2][kN];
LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) {
    w = (w << 3) + (w << 1) + (ch ^ '0');
  }
  return f * w;
}
void Chkmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
  x_ %= mod;
  y_ %= mod - 1;
  LL ret = 1;
  for (; y_; y_ >>= 1ll) {
    if (y_ & 1) ret = ret * x_ % mod;
    x_ = x_ * x_ % mod;
  }
  return ret;
}
LL Inv(LL x_) {
  return QPow(x_, mod - 2);
}
LL Sum(LL n_) {
  return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
}
void Euler() {
  vis[1] = true, mu[1] = phi[1] = 1; //初值
  for (int i = 2; i <= Lim; ++ i) {
    if (! vis[i]) {
      p[++ p_num] = i;
      mu[i] = -1;
      phi[i] = i - 1;
    }
    for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        phi[i * p[j]] = phi[i] * p[j];
        break;
      }
      mu[i * p[j]] = -mu[i];
      phi[i * p[j]] = phi[i] * (p[j] - 1);
    }
  }
}
void Prepare() {
  Euler();
  inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
  for (int i = 1; i <= Lim; ++ i) {
    g[0][i] = g[1][i] = 1;
    fac[i] = 1ll * fac[i - 1] * i % mod;
    sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
    prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
    if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;

    prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
    inv_prodt_phi[i] = Inv(prodt_phi[i]);
  }

  for (int d = 1; d <= Lim; ++ d) {
    for (int j = 1; d * j <= Lim; ++ j) {
      int t = d * j;
      if (mu[j] == 1) {
        g[0][t] = g[0][t] * d % mod;
        g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
      } else if (mu[j] == -1) {
        g[0][t] = g[0][t] * inv[d] % mod;
        g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
      }
    }
  }
  inv_prodg[0][0] = prodg[0][0] = 1;
  inv_prodg[1][0] = prodg[1][0] = 1;
  inv_prodt_phi[0] = 1;
  for (int i = 1; i <= Lim; ++ i) {
    for (int j = 0; j <= 1; ++ j) {
      prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
      inv_prodg[j][i] = Inv(prodg[j][i]);
    }
  }
}
LL f1(LL a_, LL b_, LL c_, int type) {
  if (! type) return QPow(fac[a_], b_ * c_);
  if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
  LL ret = 1, lim = min(min(a_, b_), c_);
  for (LL l = 1, r = 1; l <= lim; l = r + 1) {
    r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
    ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
    ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
  }
  return ret;
}
LL f2_2(LL a_, LL b_) { 
  LL ret = 1;
  for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
    r = min(a_ / (a_ / l), b_ / (b_ / l));
    ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
  }
  return ret;
}
LL f2(LL a_, LL b_, LL c_, int type) {
  LL ret = 1;
  if (! type) {
    for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
      r = min(a_ / (a_ / l), b_ / (b_ / l));
      LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
      ret = (ret * QPow(val, c_)) % mod;
    }
  } else if (type == 1) {
    for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
      r = min(a_ / (a_ / l), b_ / (b_ / l));
      LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
      ret = (ret * QPow(val, Sum(c_))) % mod;
    }
  } else {
    LL lim = min(min(a_, b_), c_);
    for (LL l = 1, r = 1; l <= lim; l = r + 1) {
      r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
      ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
      ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
    }
  }
  return ret;
}
//=============================================================
int main() {
  T = read(), mod = read();
  Prepare();
  while (T -- ) {
    A = read(), B = read(), C = read();
    for (int i = 0; i <= 2; ++ i) {
      ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
      ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
      printf("%lld ", ans);  
    }
    printf("\n");
  }
  return 0;
}

写在最后

参考资料:
Oi-Wiki-莫比乌斯反演
算法学习笔记(35): 狄利克雷卷积 By: Pecco
题解 SP5971 【LCMSUM - LCM Sum】 - BJpers2 的博客
题解 SP5971 【LCMSUM - LCM Sum】 - Venus 的博客
题解 P3327 【[SDOI2015]约数个数和】 - suncongbo 的博客

把 Oi-Wiki 上的内容进行了复制 整理扩充。
我是个没有感情的复读机(大雾)

posted @ 2020-04-07 17:41  Luckyblock  阅读(1653)  评论(15编辑  收藏  举报