进阶数学算法-破译密码

题目-破译密码

问题分析

给定$a,b,d$, 问$0\le x\le a$并且$0\le y\le b$使得$\gcd (x,y)=d$
问题等价于

$$0\le x\le \frac{a}{d},0\le y\le \frac{b}{d},\gcd (x,y)=1$$

满足上述条件$x,y$互质的数对个数, 可以用容斥原理

$M\ddot{o}bius$函数

假设对数字$x$分解质因数的结果是$p_1^{{\alpha }_1}{p}_2^{{\alpha }_2}...p_k^{{\alpha }_k}$

$$\mu (x)=\{\begin{array}{l}0{\alpha }_i>1\\ 1k\mathrm{mod}2=0,{\alpha }_i=1\\ -1k\mathrm{mod}2=1,{\alpha }_i=1\end{array}$$

假设计算与$N$互质的数的个数, $S_k$是$1\sim n$中$p$的倍数的个数
$$cnt=N-S_2-S_3-...-S_p+S_{2,3}+S_{2,5}+...-S_{2,3,5}-...$$

发现符号就是$\mu (x)$!, 可以在线性筛质数的情况下求出, 还是将情况分为三类

• 当前数字是质数, $\mu (i)=-1$
• $i$不是$p_j$的最小质因子, $\mu (i\times p_j)=-1\times \mu (i)$
• $i$是$p_j$的最小质因子, $\mu (i\times p_j)=0$

令$a^{\prime }=\frac{a}{d}$. $b^{\prime }=\frac{b}{d}$, 实际的数量就是
$$a^{\prime }{b}^{\prime }-S_2-S_3-...+S_{2,3}+S_{2,5}+...-S_{2,3,5}-...$$

结合$M\ddot{o}bius$函数, 得到
$$a^{\prime }{b}^{\prime }+\sum _{i=2}^{\min (a^{\prime },b^{\prime })}\lfloor \frac{a^{\prime }}{i}\rfloor \lfloor \frac{b^{\prime }}{i}\rfloor \times \mu (i)$$

现在令$a=a^{\prime }$, 展开得到
$$\frac{a}{1},\frac{a}{2},\frac{a}{3},....,\frac{a}{a}$$

引理: 上述数列不同数的数字个数不会超过$2\sqrt{n}$

证明:

将序列分为两部分, 第一部分有$\sqrt{a}$项
$$\frac{a}{1},\frac{a}{2},...,\frac{a}{\lfloor \sqrt{a}\rfloor }$$

第二部分是单调递减的, 取值最大不会超过$\sqrt{a}$
$$\frac{a}{\lfloor \sqrt{a}\rfloor +1}$$

因此第二部分不同的取值最多有$\sqrt{a}$个, 因此第二段最多$\sqrt{a}$项, 因此上述数列不同数的数字个数不会超过$2\sqrt{n}$, 证毕

在计算段内的时候, 因为数值$\frac{a^{\prime }}{i}$是相同的, 但是$\mu (x)$不是相同的, 可以预处理$\mu (x)$的前缀和, 然后再累计段内的值

分段的正确性已经证明, 现在问题变成了如何进行分段?

假设当前情况是
$$\lfloor \frac{a}{x}\rfloor$$
希望计算$x$最大到多少, 我们假设是$g(x)$, 使得下面的式子成立
$$\lfloor \frac{a}{x}\rfloor =\lfloor \frac{a}{g(x)}\rfloor ,\lfloor \frac{a}{g(x)+1}\rfloor <\lfloor \frac{a}{x}\rfloor$$

假设
$$g(x)=\lfloor \frac{a}{\lfloor \frac{a}{x}\rfloor }\rfloor$$
现在证明$g(x)$的正确性

首先证明$\lfloor \frac{a}{x}\rfloor =\lfloor \frac{a}{g(x)}\rfloor$
证明:

假设将$g(x)$内部的下取整符号去掉, $g(x{)}^{\prime }=x$, 因为去掉下取整$g(x)$会变大, 也就是$g(x)\ge g(x{)}^{\prime }$, 得到$g(x)\ge x$, 那么得到
$$\lfloor \frac{a}{x}\rfloor \ge \lfloor \frac{a}{g(x)}\rfloor \cdots (1)$$
现在需要证明的是
$$\lfloor \frac{a}{x}\rfloor \le \lfloor \frac{a}{g(x)}\rfloor$$
将式子展开
⌊ a x ⌋ ≤ ⌊ a ⌊ a ⌊ a x ⌋ ⌋ ⌋ \left \lfloor \frac{a}{x} \right \rfloor \le \left \lfloor \frac{a}{ \left \lfloor \frac{a}{\left \lfloor \frac{a}{x} \right \rfloor } \right \rfloor } \right \rfloor ⌊xa​⌋≤ ​⌊⌊xa​⌋a​⌋a​ ​
令
t = ⌊ a a ⌊ a x ⌋ ⌋ t = \left \lfloor \frac{a}{ \frac{a}{\left \lfloor \frac{a}{x} \right \rfloor } } \right \rfloor t= ​⌊xa​⌋a​a​ ​
那么$t=\lfloor \frac{a}{x}\rfloor$, 也就是说右边式子去掉一个下取整符号等于左边, 也就是原来的数变小等于左边, 那么右边大于等于左边
$$\lfloor \frac{a}{x}\rfloor \le \lfloor \frac{a}{g(x)}\rfloor \cdots (2)$$
由$(1)$和$(2)$, 证明出$\lfloor \frac{a}{x}\rfloor =\lfloor \frac{a}{g(x)}\rfloor$, 证毕

再证明$\lfloor \frac{a}{g(x)+1}\rfloor <\lfloor \frac{a}{x}\rfloor$
证明:

设$a=kx+r,(0\le r<x)$, 待入式子中得到
$$k>\lfloor \frac{a}{\lfloor \frac{a}{k}\rfloor +1}\rfloor$$
移项得到
$$k\cdot (\lfloor \frac{a}{k}\rfloor +1)>a$$
判断是否成立, 令$a=pk+q,(0\le q<k)$, 那么有
$$k\cdot (p+1)>pk+q$$
化简得到
$$k>q$$
成立, 因此$\lfloor \frac{a}{g(x)+1}\rfloor <\lfloor \frac{a}{x}\rfloor$, 证毕

$g(x)$就是$x$能达到的最大值, 通过该函数, 能计算出段的边界

算法步骤

• 线性筛法求$M\ddot{o}bius$函数
• 求$M\ddot{o}bius$函数的前缀和
• 计算答案

经过数论分块算法之后, 整体的算法时间复杂度降低到$O(n\sqrt{n})$

代码实现

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int N = 50010;

int T;
int a, b, d;
int primes[N], cnt;
bool st[N];
int mobius[N];
LL s[N];

void init(int n) {
   
    mobius[1] = 1;
    for (int i = 2; i <= n; ++i) {
   
        if (!st[i]) {
   
            primes[cnt++] = i;
            mobius[i] = -1;
        }
        for (int j = 0; primes[j] <= n / i; ++j) {
   
            st[i * primes[j]] = true;
            if (i % primes[j] == 0) {
   
                mobius[i * primes[j]] = 0;
                break;
            }
            mobius[i * primes[j]] = -1 * mobius[i];
        }
    }
    for (int i = 1; i <= n; ++i) s[i] = s[i - 1] + mobius[i];
}

int g(LL a, LL x) {
   
    return a / (a / x);
}

int main() {
   
    ios::sync_with_stdio(false);
    cin.tie(0);

    init(N - 1);

    cin >> T;
    while (T--) {
   
        cin >> a >> b >> d;
        a = a / d, b = b / d;
        
        LL ans = 0;
        int n = min(a, b);
        for (int l = 1, r; l <= n; l = r + 1) {
   
            r = min({
   n, g(a, l), g(b, l)});
            ans += (s[r] - s[l - 1]) * (LL) (a / l) * (b / l);
        }
        
        cout << ans << '\n';
    }

    return 0;
}