数论常用函数

数论常用函数 最大公因数

int gcd(int x,int y)
{
    return y == 0?x:gcd(y,x % y);
}

素数筛 暴力筛法

bool isPrime(a) {
  if (a < 2) return 0;
  for (int i = 2; i < a; ++i)
    if (a % i == 0) return 0;
  return 1;
}

时间复杂度 \(O(n^2)\)

埃氏筛

int Eratosthenes(int n) {
  int p = 0;
  for (int i = 0; i <= n; ++i) is_prime[i] = 1;
  is_prime[0] = is_prime[1] = 0;
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      prime[p++] = i;  // prime[p]是i,后置自增运算代表当前素数数量
      if ((long long)i * i <= n)
        for (int j = i * i; j <= n; j += i)
          // 因为从 2 到 i - 1 的倍数我们之前筛过了，这里直接从 i
          // 的倍数开始，提高了运行速度
          is_prime[j] = 0;  // 是i的倍数的均不是素数
    }
  }
  return p;
}

时间复杂度 \(O(n\log \log n)\)

线性筛（欧拉筛）

void init()
{
    int tot = 0;
    for(int i = 2;i <= n;i++)
    {
        if(!lst[i]) prime[++tot] = i;
        for(int j = 1;j <= tot && i * prime[j] <= n;j++)
        {
            lst[i * prime[j]] = 1;
            if(!i % prime[j]) break;
        }
    }
}

时间复杂度 \(O(n)\)

欧拉函数 求解

int euler_phi(int n) {
  int m = int(sqrt(n + 0.5));
  int ans = n;
  for (int i = 2; i <= m; i++)
    if (n % i == 0) {
      ans = ans / i * (i - 1);
      while (n % i == 0) n /= i;
    }
  if (n > 1) ans = ans / n * (n - 1);
  return ans;
}

优化后的代码：

int euler_phi(int n) {
  int ans = n;
  for (int i = 2; i * i <= n; i++)
    if (n % i == 0) {
      ans = ans / i * (i - 1);
      while (n % i == 0) n /= i;
    }
  if (n > 1) ans = ans / n * (n - 1);
  return ans;
}

筛法求解欧拉函数：

void pre() {
  memset(is_prime, 1, sizeof(is_prime));
  int cnt = 0;
  is_prime[1] = 0;
  phi[1] = 1;
  for (int i = 2; i <= 5000000; i++) {
    if (is_prime[i]) {
      prime[++cnt] = i;
      phi[i] = i - 1;
    }
    for (int j = 1; j <= cnt && i * prime[j] <= 5000000; j++) {
      is_prime[i * prime[j]] = 0;
      if (i % prime[j])
        phi[i * prime[j]] = phi[i] * phi[prime[j]];
      else {
        phi[i * prime[j]] = phi[i] * prime[j];
        break;
      }
    }
  }
}

欧拉定理(\(GCD\))

若 \(\gcd(a, m) = 1\)，则 \(a^{\varphi(m)} \equiv 1 \pmod{m}\)

扩展欧拉定理(\(exGCD\))

\[a^b\equiv \begin{cases} a^{b\bmod\varphi(p)},\,&\gcd(a,\,p)=1\\ a^b,&\gcd(a,\,p)\ne1,\,b<\varphi(p)\\ a^{b\bmod\varphi(p)+\varphi(p)},&\gcd(a,\,p)\ne1,\,b\ge\varphi(p) \end{cases} \pmod p \]

分解质因数 朴素算法

list<int> breakdown(int N) {
  list<int> result;
  for (int i = 2; i * i <= N; i++) {
    if (N % i == 0) {  // 如果 i 能够整除 N，说明 i 为 N 的一个质因子。
      while (N % i == 0) N /= i;
      result.push_back(i);
    }
  }
  if (N != 1) {  // 说明再经过操作之后 N 留下了一个素数
    result.push_back(N);
  }
  return result;
}

时间复杂度 \(O(\sqrt n)\)

裴蜀定理

设 \(a,b\) 是不全为零的整数，则存在整数 \(x,y\), 使得 \(ax+by=\gcd(a,b)\).

费马小定理

若 \(p\) 为素数，\(\gcd(a, p) = 1\)，则 \(a^{p - 1} \equiv 1 \pmod{p}\)

另一个形式：对于任意整数 \(a\)，有 \(a^p \equiv a \pmod{p}\)

求解乘法逆元

如果一个线性同余方程 \(ax \equiv 1 \pmod b\)，则 \(x\) 称为 \(a \bmod b\) 的逆元，记作 \(a^{-1}\)

朴素做法

void exgcd(int a, int b, int& x, int& y) {
      if (b == 0) {
        x = 1, y = 0;
        return;
      }
      exgcd(b, a % b, y, x);
      y -= a / b * x;
}

乘法逆元

因为 \(ax \equiv 1 \pmod b\)；

所以 \(ax \equiv a^{b-1} \pmod b\)

所以 \(x \equiv a^{b-2} \pmod b\)。

然后我们就可以用快速幂来求了。

inline int qpow(long long a, int b) {
      int ans = 1;
      a = (a % p + p) % p;
      for (; b; b >>= 1) {
        if (b & 1) ans = (a * ans) % p;
        a = (a * a) % p;
      }
      return ans;
}

线性求解

求出 \(1,2,...,n\) 中每个数关于 \(p\) 的逆元。

如果对于每个数进行单次求解，以上两种方法就显得慢了，很有可能超时，所以下面来讲一下如何线性（\(O(n)\)）求逆元。

首先，很显然的 \(1^{-1} \equiv 1 \pmod p\)；

证明:
对于 \(\forall p \in \mathbf{Z}\)，有 \(1 \times 1 \equiv 1 \pmod p\) 恒成立，故在 \(p\) 下 \(1\) 的逆元是 \(1\)，而这是推算出其他情况的基础。

其次对于递归情况 \(i^{-1}\)，我们令 \(k = \lfloor \frac{p}{i} \rfloor\)，\(j = p \bmod i\)，有 \(p = ki + j\)。再放到 \(\mod p\) 意义下就会得到：\(ki+j \equiv 0 \pmod p\)；

两边同时乘 \(i^{-1} \times j^{-1}\)：

\(kj^{-1}+i^{-1} \equiv 0 \pmod p\)

\(i^{-1} \equiv -kj^{-1} \pmod p\)

再带入 \(j = p \bmod i\)，有 \(p = ki + j\)，有：

\(i^{-1} \equiv -\lfloor\frac{p}{i}\rfloor (p \bmod i)^{-1} \pmod p\)

我们注意到 \(p \bmod i < i\)，而在迭代中我们完全可以假设我们已经知道了所有的模 \(p\) 下的逆元 \(j^{-1}, j < i\)。

故我们就可以推出逆元，利用递归的形式，而使用迭代实现：

\[i^{-1} \equiv \begin{cases} 1, & \text{if } i = 1, \\ -\lfloor\frac{p}{i}\rfloor (p \bmod i)^{-1}, & \text{otherwises}. \end{cases} \pmod p \]

代码实现：

// C++ Version
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}

中国剩余定理

可求解如下形式的一元线性同余方程组（其中 \(n_1, n_2, \cdots, n_k\) 两两互质）：

\[\begin{cases} x &\equiv a_1 \pmod {n_1} \\ x &\equiv a_2 \pmod {n_2} \\ &\vdots \\ x &\equiv a_k \pmod {n_k} \\ \end{cases} \]

算法流程

1. 计算所有模数的积 \(n\)；
2. 对于第 \(i\) 个方程：
   1. 计算 \(m_i=\frac{n}{n_i}\)；
   2. 计算 \(m_i\) 在模 \(n_i\) 意义下的 逆元 \(m_i^{-1}\)；
   3. 计算 \(c_i=m_im_i^{-1}\)（不要对 \(n_i\) 取模）。
3. 方程组的唯一解为：\(x=\sum_{i=1}^k a_ic_i \pmod n\)。

示例代码

// C++ Version
typedef long long LL;
LL CRT(int k, LL* a, LL* r) {
  LL n = 1, ans = 0;
  for (int i = 1; i <= k; i++) n = n * r[i];
  for (int i = 1; i <= k; i++) {
    LL m = n / r[i], b, y;
    exgcd(m, r[i], b, y);  // b * m mod r[i] = 1
    ans = (ans + a[i] * m * b % n) % n;
  }
  return (ans % n + n) % n;
}

同余定理 内容

1. \((a + b) \mod m = (a \mod m + b \mod m) \mod m\)

2. \((a*b) \mod m = ((a \mod m) \times (b \mod m)) \mod m\)

应用 高精度取余

int mod(string s,int md)
{
    int ans = 0;
    for(char c:s)
    {
        ans = 10 * ans + (c ^ '0');
        ans %= md;
    }
    return ans;
}

版权声明：

本文中有部分代码以及讲解来源于 oiwiki.com