快速计算

一、快速幂

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 【模板】快速幂||取余运算

二、光速幂

1. 前置芝士

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

2. 应用范围

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

3. 算法原理

根据拓展欧拉定理：

对于\(~a,b \in \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));
}

三、龟速乘

求 \(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)\)。

---

四、更加快速的乘法

当然还有一种算法：

\[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} \]

我们先用longdouble对\(\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