初等数论之同余总结

上周很颓废的就完成了一个数论的同余专题……总结一下吧，也是备忘。

扩展欧几里德算法

可以求解如Ax+By=C的二元一次不定方程的整数解。

考虑不定方程Ax+By=C，当(A,B)|C时一定有整数解。考虑如何解Ax+By=(A, B)

考虑欧几里德GCD算法GCD(A,B)的最后一步：
当GCD(A, 0)时，(x, y) = (1, 0)为满足Ax+0y=(A, 0)的一组解。考虑倒着推回去：

设
   A1x1 + B1y1 = (A1, B1)
   A2x2 + B2y2 = (A2, B2) = (A1, B1)
由
   (A, B) = (B, A % B)  (B, A - [A / B] * B)
可得
   A1x1 + B1y1 = A2x2 + B2y2 = B1x2 + (A1 - [A1 / B1] * B1)y2
移项得
   (x1-y2)A1 + (y1-x2+[A1/B1]y2)B1 = 0
因此
   x1 = y2, y1 = x2 - [A1 / B1] * y2
如此递推就可以得到Ax+By=(A, B)的解了

可以得到如下算法：

typedef long long Int64;
struct Triple
{
  Int64 x, y, z;
  Triple(const Int64 a, const Int64 b, const Int64 c): x(a), y(b), z(c) { }
};
Triple ExtendedGCD(const Int64 a, const Int64 b)
{
  if (b == 0) return Triple(1, 0, a);
  Triple last(ExtendedGCD(b, a%b));
  return Triple(last.y, last.x - a / b * last.y, last.z);
}

其中Triple是一个三元组，(x, y, z)分别保存(解x, 解y, GCD)。

在我们解出Ax+By=(A, B)的解后，要求Ax+By=(A, B)的解，只需把(x, y)乘上C/(A, B)即可！

Ax+By=C的x的最小非负整数解

其实二元一次不定方程有一个定理：

若(x0, y0)为一组解，那么(x0+B/(A,B)*t, y0-A/(A,B)*t), t∈Z，也是一组解

由此可以知道，x的最小非负整数解就是把x0不断加减B/(A,B)使得它是最小非负整数解

简单地来说就是：x = x0 % (B / (A,B))

注意这里的模是数学模，用C++的写法a%b就是：{cpp}(a % b + b) % b{/cpp}

乘法逆元

Multiplicative inverse

大家都知道a / b = a * b^(-1)，大家也都知道在模上不能直接计算除法。因此模上的除法一般转换为乘上其乘法逆元。

所谓a在模p下的乘法逆元x也就是：ax ≡ 1(mod p)

很简单，设ax = py + 1，扩展欧几里德套上解之

注意，这里有解的条件，(a, p)|1，也就是(a, p) = 1。因此，乘法逆元不一定存在。

线性模方程组

Linear congruence theorem

考虑如下方程组：

A1x ≡ B1(mod M1)
A2x ≡ B2(mod M2)
A3x ≡ B3(mod M3)
...
Anx ≡ Bn(mod Mn)

先把系数A除掉（乘上乘法逆元），变成如下形式：

x ≡ C1(mod M1)
x ≡ C2(mod M2)
x ≡ C3(mod M3)
...
x ≡ Cn(mod Mn)

求解的方法是两两合并：

先考虑
   x ≡ C1(mod M1)
   x ≡ C2(mod M2)
设x = M1y + C1 = M2z + C2
即M1y + M2z = C2 - C1
扩展欧几里德套上解之，得(y, z)
那么x = M1y + C1
此时两个方程合并为：
   x ≡ M1y + C1(mod LCM(M1, M2))
再与
   x ≡ C3(mod M3)
合并，如此反复，最后的得到我们要的解：
   x ≡ C(mod M)

其中只要任意一步(y, z)无解，方程组无解。

程序如下：

Int64 RemainderEquations(const Int64* const a, const Int64* const m, int n)
{
  Int64 lasta = a[0], lastm = m[0];
  for (int i = 1; i < n; ++i)
  {
    Triple gcd(ExtendedGCD(lastm, m[i]));
    if ((a[i]-lasta)%gcd.z) return -1;
    const Int64 times = (a[i]-lasta)/gcd.z;
    const Int64 tmp = m[i] / gcd.z;
    gcd.x = (times*gcd.x%tmp+tmp)%tmp;
    lasta = gcd.x*lastm+lasta;
    lastm = lastm/(gcd.z)*m[i];
  }
  return (lasta%lastm+lastm)%lastm;
}

中国剩余定理

Chinese remainder theorem

如果方程组：

x ≡ C1(mod M1)
x ≡ C2(mod M2)
x ≡ C3(mod M3)
...
x ≡ Cn(mod Mn)

满足Mi两两互质，可以使用中国剩余定理来解决。

中国剩余定理的核心思想是对于每一个模方程求出一个Ei，符合：Ei ≡ 0(mod Mj) j != i 以及 Ei ≡ 1(mod Mj) j == i

如此一来，x ≡ Sigma(EiCi)(mod M)就是方程组的解了

下面看如何求解Ei：

令
   M = M1*M2*...*Mn
   mi = M/Mi
因为Ei满足：
  Ei ≡ 0(mod mi)
  Ei ≡ 1(mod Mi)

扩展欧几里德套上解之即可得到Ei

注意这里必定有解

因为，设Ei = mi*y又设Ei = Mi*z+1，那么mi*y - Mi*z = 1

由于(mi, Mi) = 1，因此必定有解

long long M = 1, ans = 0;
for (int i = 0; i < n; ++i) M *= m[i];
for (int i = 0; i < n; ++i)
{
  const long long mi = M/m[i];
  ans = (ans + a[i] * mi * Extended(mi, m[i]).x) % M;
}

模上除法

计算A/B%M，如果有逆元，乘上逆元就好了。没有的话比较麻烦：

把M拆成P1^C1 + P2^C2 + ... + Pk^Ck的形式，分别计算A / B % (Pi^Ci)的结果，再合并方程组

由于B在Pi^Ci仍然不一定存在逆元，所以我们把A和B都拆成如下形式：Pi^m * n

这样一来A/B%Pi变成了：(Pi^m * r) / (Pi^n * s) % Pi

显然，m < n时等于0，否则就变成了：(Pi^(m-n) * r) / s % Pi

此时有(s, Pi) = 1，因此乘上s的逆元就好了

到最后把方程组合并就好了

阶乘取模分解

n! % P假设模P的是素数（不是的话分解质因数拆掉，分别处理，再中国剩余定理合并），还是得把n!转换成P^a * b的形式

令
   k = [n/p]
   m = n%p
n! mod p =  1      * 2     * ... * (p-1)  * p
          * (p+1)  * (p+2) * ... * (2p-1) * 2p
          * ...
          * (kp+1) * (kp+2) * ... * (kp+m)      mod p

         =  1 * 2 * ... * (p-1) * p
          * 1 * 2 * ... * (p-1) * 2p
          * ...
          * 1 * 2 * ... * m          mod p

         =  [(n%p)! * (p-1)!^(n/p)] * [(n/p)! * p^(n/p)]    mod p

其中(n%p)!和(p-1)!都可以预处理，(p-1)!^(n/p)和p^(n/p)都可以快速幂，剩下的只需要递归处理(n/p)!即可

可以参考[BZOJ2313]分组，以及[BZOJ2142]礼物。

组合数取模

Lucas’ theorem

C(n, k) % p仍然假设p为质数。这时候就有神奇的Lucas定理：

C(n, k) % p = Lucas(n, k, p) = Lucas(n/p, k/p, p) * C(n%p, k%p)

原理就是阶乘取模

另外注意费马小定理：

p为质数那么，a^(p-1) ≡ 1(mod p)

也就是说a^(p-2)在mod p下为a的乘法逆元！如此一来，配合快速幂，计算C(n%p, k%p)就非常快了，p的范围能上10^4

程序如下：

//fac[i] = i! % p
long long Power(long long n, long long k)
{
  long long ans = 1;
  for (; k; k >>= 1)
  {
    if (k & 1) ans = ans*n%p;
    n = n*n%p;
  }
  return ans;
}
inline long long C(const long long n, const long long k)
  { return n >= k ? fac[n]*Power(fac[k]*fac[n-k]%p, p-2)%p : 0; }
long long Lucas(long long n, long long k, const long long p)
{
  long long ans = 1;
  while (n && k && ans)
  {
    ans = ans*C(n%p, k%p)%p;
    n /= p;
    k /= p;
  }
  return ans;
}

欧拉函数

[Euler’s totient function](http://en.wikipedia.org/wiki/Euler’s_totient_function)
线性时间内筛素数和欧拉函数

phi(x)或者φ(x)，表示1~x中和x互质的有多少个数。定义/性质如下：

phi(x) = n * product{p|n & p is prime, 1 - 1/p}
phi(p^k) = (p-1) * p^(k-1)
phi(n^m) = phi(n) * n^(m-1)
phi(mn) = phi(m) * phi(n), (gcd(m, n) = 1)
phi(mn) = phi(m) * phi(n) * gcd(m, n)/phi(gcd(m, n))
n = sum{d|n, phi(d)}
若gcd(a, p) = 1，则a^phi(p) ≡ 1(mod p)

阶

Multiplicative order

若a^n ≡ 1 (mod p)，那么称最小的符合要求的那个n为a关于mod p的阶。因为a^phi(p) ≡ 1 (mod p)，所以n一定是phi(p)的约数。算的时候就是从小到大枚举phi(p)约数然后快速幂验证。

inline Int64 order(Int64 a, const Int64 p)
{
  static int b[100];
  const Int64 phi = p-1;
  int cnt = 0;
  for (Int64 i = 1; i*i < p; ++i) if (phi % i == 0) b[cnt++] = i;
  for (int i = cnt-1; i >= 0; --i) b[cnt++] = phi / b[i];
  for (int i = 0; i < cnt; ++i) if (power(a, b[i], p) == 1) return b[i];
  return -1;
}

原根

Primitive root modulo n

如果g的阶是phi(p)的话，那么称g是原根。只有**1、2、4、奇质数的幂、2*奇质数的幂**才有原根。因为原根都很小，所以一般就暴力枚举了。

Int64 PrimitiveRoot(const int p, const int pc)
{
  static int d[SQRTN*2];
  const int phi = pc/p * (p-1);
  int cnt = 0;
  for (int i = 1; i*i <= phi; ++i) if (phi % i == 0) d[cnt++] = i;
  for (int i = cnt-1; i >= 0; --i) d[cnt++] = phi/d[i];
  for (int g = 1, ord; ; ++g)
  {
    for (ord = 0; ord < cnt && Power(g, d[ord], pc) != 1; ++ord);
    if (d[ord] == phi) return g;
  }
}

指标（离散对数）

Discrete logarithm

若(b, p) = 1, 0 <= a < phi(p)，满足g^a ≡ b (mod p)，那么a称为b关于p的指标，记为ind(g)b，简写为ind b。有下列性质：

a ≡ b (mod p) <=> ind a ≡ ind b (mod phi(p))
ind (a*b) ≡ ind a + ind b (mod phi(p))
ind (a^n) ≡ n ind a (mod phi(p))

很明显，指标要用BSGS求了，因为p是质数，所以BSGS还是非常好写的。参考：Baby Step Giant Step 及其扩展算法

QuarterGeek

Lequn Chen (abcdabcd987) 's Blog