快速计算

快速计算

1. 快速幂

(1) 算法原理

求 \(a^b\bmod p\) 的结果。

我们可以将 \(b\) 进行二进制拆分，并构造如下算法：

\[a^b \bmod p=\begin{cases}(a^{\frac b 2})^2 \bmod p&\texttt{b is even}\\a(a^{\frac{b-1}2})^2 \bmod p&\texttt{b is odd}\end{cases} \]

---

为了更好地理解，我们举个栗子
以 \(5^{23}\) 为例：（先把取模省略掉）

• \(5^{23} = (5^{11})^2 ×5\)
• \(5^{11} = (5^5)^2 × 5\)
• \(5^5 = (5^2)^ 2 × 5\)
• \(5^2 = 5 × 5\)

再将式子从下往上递归即可

---

(2) 代码实现

递归版：

typedef long long ll;
ll ksm(ll x,ll k){
    if(k == 0) return 1;
	if(k == 1) return x;
	ll t = ksm(x * x % MOD,k / 2);
	if(k & 1) return (t * t % MOD) * x % MOD;
	else return t * t % MOD;
}

循环版（更常用）：

typedef long long ll;
ll ksm(ll x,ll k){
	ll res = 1;
	while(b){
		if(b&1) res = res * x % MOD;
		x = x * x % MOD;
        x >>= 1;
	}
	return res;
}

注：MOD为模数

(3) 时间复杂度：

\(O(\log b)\)

很显然每次 \(b\) 会减半

模板题：P1226 【模板】快速幂||取余运算

2. 光速幂

(1) 前置芝士

• 快速幂
• 欧拉函数 & 拓展欧拉定理

(2) 应用范围

在幂运算的底数和模数已经确定的情况下，可以用 \(O(\sqrt k)\)（\(k\) 为指数的值域）预处理，\(O(1)\) 查询（当然，你需要预处理 \(\varphi\)，或者说 \(p\) 为质数）

(3) 算法原理

根据拓展欧拉定理：

对于 \(a,b \in \mathbb{Z}\)，有 \(a^b \equiv \begin{cases} a^b & b < \varphi(p)\\a^{\varphi(p) + b ~mod~\varphi(p)} & b\geq\varphi(p)\end{cases}~\)

然后我们就可以将 \(b\) 降至 \(2\times\varphi(p)\)

接着，我们记 \(s = \lceil\sqrt p\rceil\)，则

\[a^b = a^{s\times{\lfloor\frac b s\rfloor}}\times a^{b~mod~s} \]

然后我们便只需要处理：\(a^{s\times i}\) 和 \(a^i\) 即可（\(i\) 从 \(1\sim s\)）

考虑本质上，就是把指数按根号下值域离线，空间换时间

(4) 代码实现

ll a,b,p,s;
ll p1[NN],p2[NN];//p1为a^i,p2位a^(i*s)
ll LSP(int k){
	k %= phi(p);
	return p1[k % s] * p2[k / s] % p;
}
int main(){
	scanf("%lld%lld%lld",&a,&b,&p);
	s = sqrt(p) + 1;
	p1[0] = p2[0] = 1;
	for(int i = 1; i <= s; ++i) p1[i] = p1[i-1] * a % p;
	for(int i = 1; i <= s; ++i) p2[i] = p2[i-1] * p1[s] % p;//预处理
	printf("%lld^%lld mod %lld=%lld\n",a,b,p,LSP(b));
}

3. 龟速乘

求 \(a\times b\bmod p\) 的结果。

\(a,b,p\) 都是 \(10^{18}\) 级别。

我们可以构造类似的算法（就是把乘法改成加法即可）

code:

long long ksm(long long a,long long b,long long p){
	long long res = 0,tmp = a;
	while(b){
		if(b&1) res = (res + tmp) % p;
		tmp = (tmp + tmp) % p;
	}
	return res;
}

时间复杂度 \(O(\log b)\)。

---

4. 更加快速的乘法

当然还有一种算法：

\[a×b\bmod p = a×b - \lfloor\frac{a×b}p\rfloor × p \]

我们不妨令 \(x = a × b,y = \lfloor\frac{a×b}p\rfloor × p\)：

\[a×b\bmod p = x - y \]

又因为 \(x - y = a×b\bmod p< p < 2^{64}\)

\[a×b\bmod p = x - y = (x-y)\bmod 2^{64} \]

我们先用 long double 对 \(\frac{a×b}p\) 进行计算，再下取整由于精度误差算出来的值可能会小 \(1\)，但是在模\(p\)意义下完全没有影响

我们再使用 unsigned long long 计算 \(x\) 和 \(y\) （unsigned long long 自然溢出即对 \(2^{64}\) 取模）

最后是代码：

typedef unsigned long long ull;
typedef long long ll;
ull ksm(ull a,ull b,ull p){
	ull x = a * b, y = (ull)((long double)a * b / p) * p;
	ll res = (ll) (x % p) - (ll)(y % p);
	return res + (res < 0 ? p : 0);
}

作者评价：不如直接 __int128