模反元素 (Modular Multiplicative Inverse)

模板: Luogu P3811

求 \(1\) 到 \(n\) 每个数模 \(p\) 意义下最小正整数乘法逆元
\(n \leq 3*10^6\), \(p < 20000528\), \(500ms\)

概念

模反元素, 又称模逆元, 简称逆元, 其定义是在取模意义下的倒数

\[ax \% p = 1 \]

则称 \(x\) 是 \(a\) 模 \(p\) 意义下的逆元.

\(a\) 模 \(p\) 意义下的逆元存在当且仅当 \(a\), \(p\) 互质, 接下来给出证明.

Exgcd

根据定义, 整理式子

\[ax + py = 1 \]

根据最大公因数的定义, \(gcd(a, p)\) 是式子 \(ax + py\) 最小的正值. 如果 \(ax + py = 1\), 那么必须满足 \(gcd(a, b) = 1\), 即两数互质.

可以用 exgcd 快速算出使得 \(ax + py = 1\) 的 \(x\), \(y\) 的一组解, 通过使 \(x\) 对 \(p\) 的取模, 可以常数时间利用 \(x\) 求出最小的正整数解.

当然, 由于求 \(3 * 10^6\) 个数的逆元, 一次操作复杂度 \(O(logn)\), 最后是会超时的, 第一次提交总时间是 \(1.01s\), 在我在细节上压榨过常数之后, 下面这份代码的时间来到了 \(959ms\).

inline void Exgcd(int a, int b, int &x, int &y) {
  if(!b) {
    x = 1;
    y = 0;
  }
  else {
    Exgcd(b, a % b, y, x);
    y -= (a / b) * x;
  }
}
signed main() {
  n = RD();
  p = RD();
  for (register int i(1); i <= n; ++i) {
    Exgcd(p, i, A, B);
    B %= p;
    B += p;
    B %= p;
    printf("%d\n", B); 
  }
  return Wild_Donkey;
}

接下来, 分析代码, 发现递归是可以优化的, 于是改成循环, 优化到 \(920ms\)

inline void Exgcd(int a, int b, int &x, int &y) {
  int tmp;
  Cnt = 0;
  while(b) {
    a_b[++Cnt] = a / b;
    tmp = a;
    a = b;
    b = tmp % b;
  }
  x = 1;
  y = 0;
  for (register unsigned j(Cnt); j >= 1; --j) {
    tmp = x;
    x = y;
    y = tmp - a_b[j] * x;
  }
}

然后发现还有好多时间浪费到了输出上, 从网上嫖了个 fwrite, 并且进一步在细节上压榨 (如, 内层循环用 short 类型的控制器), 优化到 \(888ms\).

signed main() {
  n = RD();
  p = RD();
  for (register unsigned i(1); i <= n; ++i) {
    C = p;
    D = i;
    Cnt = 0;
    while(D) {
      a_b[++Cnt] = C / D;
      Tmp = C;
      C = D;
      D = Tmp % D;
    }
    A = 1;
    B = 0;
    for (register short j(Cnt); j >= 1; --j) {
      Tmp = A;
      A = B;
      B = Tmp - a_b[j] * A;
    }
    B %= p;
    B += p;
    B %= p;
    write(B);
  }
  fwrite(_d,1,_p-_d,stdout);
  return Wild_Donkey;
}

最后发现 Tmp 是 unsigned, 其它的变量是 int, 强转浪费时间了, 于是都改成 int 之后, 直接一波优化到 \(816ms\).

线性递推

用 \(a\) 表示 \(p\), 设 \(k = \left\lfloor\frac{a}{p}\right\rfloor\), \(r = p \% a\)

\[(ak + r) \% p = 0 \]

如果用 \(Inv_i\) 表示 \(i\) 模 \(p\) 意义下的逆元, 则式子除以 \(Inv_aInv_r\) 得

\[(Inv_rk + Inv_a) \% p = 0 \]

移项, 由于所求的 \(x = Inv_a\), 得到表达式

\[Inv_a = (-Inv_rk) \% p \]

可以递推求解

由于 \(r = p \% a\), 所以 \(r < a\), \(Inv_r\) 已经求出, 而 \(k = \left\lfloor\frac{a}{p}\right\rfloor\), 可以直接求, 所以可以线性时间求 \(Inv_1\) 到 \(Inv_n\), 下面这份代码只用了 \(253ms\)

signed main() {
  n = RD();
  p = RD();
  a[1] = 1;
  write(a[1]);
  for (register unsigned i(2); i <= n; ++i) {
    a[i] = ((long long)a[p % i] * (p - p / i)) % p;
    write(a[i]);
  }
  fwrite(_d,1,_p-_d,stdout);
  return Wild_Donkey;
}