杂谈 2：论编程中的简单数学

杂谈 2：论编程中的简单数学

\(1^{st}\) 质数

因数只有 \(1\) 和数字本身的数称为质数。

素数计数函数：小于或等于 \(x\) 的素数的个数，用 \(\pi(x)\) 表示。随着 \(x\) 的增大，有这样的近似结果：\(\pi(x) \sim \dfrac{x}{\ln(x)}\)。

素性测试：素性测试（Primality test）是一类在不对给定数字进行素数分解（prime factorization）的情况下，测试其是否为素数的算法。

那我们怎么判断一个数是不是质数呢？我们可以枚举 \(2 \sim x\)，对于每个 \(i\)，如果 \(x \bmod i = 0\)，那么就不是质数。当然，\(x < 2\) 时 \(x\) 一定不是质数。

代码：

bool Prime(ll x) {
  if (x < 2) {
    return 0;
  }
  for (ll i = 2; i <= x; ++ i) {
    if (!(x % i)) {
      return 0;
    }
  }
  return 1;
}

可是，这样求解的时间复杂度是 \(O(x)\)，\(x\) 太大直接爆炸。我们发现如果 \(i\) 是 \(n\) 的约数，那么 \(\cfrac{n}{i}\) 也是 \(n\) 的约数（称 \(\Big(i,\ \cfrac{n}{i}\Big)\) 为 \(n\) 的约数对）。所以不必重复枚举 \(2\) 个，我们只要枚举最小的那个即可。由推理可得，所有 \(n\) 的约数对的最小数均 \(\le \sqrt{n}\)。所以只需枚举 \(2 \sim \sqrt{x}\) 即可。

代码：

bool Prime(ll x) {
  if (x < 2) {
    return 0;
  }
  for (ll i = 2; i * i <= x; ++ i) { // i * i <= x 等价于 i <= √x
    if (!(x % i)) {
      return 0;
    }
  }
  return 1;
}

那怎么一次性求出 \(1 \sim n\) 的所有数是不是质数呢？如果把每个数都检验一遍，那么时间复杂度为 \(O(n\sqrt{n})\)，大概能撑到 \(n = 2·10^5\)。下面介绍一种 \(O(n)\) 的算法：

找到一个质数后，就将它的倍数标记为合数，也就是把它的倍数筛掉，如果一个数没有被比它小的质数筛掉，那它就是质数。直到把 \(\le \sqrt{n}\) 的质数全部筛完，程序结束。时间复杂度为 \(O(n \log \log n)\)。

代码：

const int kMaxN = 1e6 + 10; // 筛掉 1 ~ 10^6+10 的质数
bool isp[kMaxN + 1]; // 多开 1 个保险
bool InitPrime(ll x) {
  for (int i = 2; i * i <= kMaxN; ++ i) {
    if (!isp[i]) { // 如果 i 是质数
      for (int j = i * i; j <= kMaxN; j += i) { // 标记
        isp[j] = 1;
      }
    }
  }
}

\(2^{\text{nd}}\) 最大公约数

最大公约数即为 Greatest Common Divisor，常缩写为 gcd。一组整数的最大公约数，是指所有公约数里面最大的一个。

\(2\) 个数的最大公约数：

考虑枚举 \(1 \sim \min(a,\ b)\)，找到最大的 \(i\) 有 \(i \bmod a = i \bmod b = 0\)，则 \(i\) 是 \(\gcd(a,\ b)\) 的值。

代码：

ll GCD(ll a, ll b) {
  ll ans = 0;
  for (ll i = 1; i <= min(a, b); ++ i) {
    if (!(i % a) && !(i % b)) {
      ans = i;
    }
  }
  return ans;
}

但是，这个算法的时间复杂度是 \(O(\min(a,\ b))\)，\(a\) 和 \(b\) 太大就会超时。设 \(a > b\)，我们发现如果 \(b\) 是 \(a\) 的约数，那么 \(b\) 就是二者的最大公约数。下面讨论不能整除的情况，即 \(a = bq + r\)（\(r < b\)）。我们可以得到 \(\gcd(a,\ b) = \gcd(b,\ a \bmod b)\)。

代码：

ll GCD(ll a, ll b) {
  if (!b) {
    return a;
  }
  return GCD(b, a % b);
}

但大整数（\(> 2^{64}\)）取模的时间复杂度较高，而加减法时间复杂度较低。针对大整数，我们可以用更相减损术。设 \(a \ge b\)，有：

\[\gcd(a, b) = \begin{cases}\text{if}\ a = b & a \\ \text{else}\ & \gcd(a - b, b) \end{cases} \]

因为当 \(a = b\) 时，\(\forall d\ |\ a, d\ |\ b\)，有 \(d\ |\ a - b\)，否则可以证明 \(a, b\) 的所有公因数都是 \(a - b,\ b\) 的公因数。

代码：

BigInt GCD(BigInt a, BigInt b) {
  int cnt1 = 0, cnt2 = 0;
  for (; !Mod(a, 2); ++ cnt1) { // Mod 为高精 a % b 函数
    RMove(a, 1); // a >>= 1 的高精版
  }
  for (; !Mod(b, 2); ++ cnt1) { // Mod 为高精 a % b 函数
    RMove(b, 1); // b >>= 1 的高精版
  }
  for (; ; ) {
    for (; !Mod(a, 2), RMove(a, 1));
    for (; !Mod(b, 2), RMove(b, 1));
    if (Same(a, b)) { // a == b 的高精版
      break;
    }
    if (a < b) {
      Swap(a, b); // 高精 swap(a, b)
    }
    Diff(a, b); // a -= b 高精版 
  }
  return LMove(a, min(cnt1, cnt2)); // LMove 为 a <<= b 的高精版
}

多个数的最大公约数：

显然答案一定是每个数的约数，那么也一定是每相邻两个数的约数。

代码：

ll GCD(ll a, ll b) {
  if (!b) {
    return a;
  }
  return GCD(b, a % b);
}

ll MGCD() {
  int _gcd = a[1];
  for (int i = 2; i <= n; ++ i) {
    _gcd = GCD(_gcd, a[i]);
  }
  return _gcd;
}

\(3^{\text{rd}}\) 快速幂

我们知道 \(x^a\) 乘上 \(x^b\)，等于 \(x^{a + b}\)，所以有 \(x^ax^b = x^{a + b}\)。因此，设 \(f(a, b)\) 表示 \(a\) 的 \(b\) 次方，我们可以这样算：

\[f(a, b) = \begin{cases}\text{if}\ b = 0 & 1 \\ \text{if}\ b = 1 & a \\ \text{if}\ b \bmod 2 = 1 & a·f(x, b - 1) \\ \text{if}\ b \bmod 2 = 0 & f\Big(a, \cfrac{b}{2}\Big)^2\end{cases} \]

代码：

ll Qkpow(ll a, ll b) {
  if (!b) {
    return 1;
  }
  if (b == 1) {
    return a;
  }
  int tmp = Qkpow(a, (b >> 1)); // b >> 1 相当于 b / 2
  return tmp * tmp * ((b & 1)? a : 1); // b & 1 相当于 b % 2
} 

但是，我们还可以使用二进制的方式进行运算。其实跟前面差不多，但是省去了递归步骤。

代码：

// 无模数版

ll Qkpow(ll b, ll p) {
  ll ans = 1;
  for (; p; p >>= 1) {  // p >>= 1 相当于 p /= 2
    if (p & 1) {
      ans = ans * b;
    }
    b = b * b;
  }
  return ans;
}   

// 有模数版

ll Qkpow(ll b, ll p, ll mod) {
  ll ans = 1;
  for (; p; p >>= 1) {  // p >>= 1 相当于 p /= 2
    if (p & 1) {
      ans = ans * b % mod;
    }
    b = b * b % mod;
  }
  return ans;
}