高速な剰余演算
Mizar/みざー最適化の方法の例:
- コンパイラ等で用いられる除算最適化（商・剰余の計算、整除性の判定など）
- Montgomery Reduction
- Barrett Reduction
何故、汎用の除算命令・除算ルーチンを使う事が避けられるかというと、遅いからである。用途を絞れば、加減乗算や比較などの命令を上手く組み合わせた方が速いケースが多く存在する。
Mizar/みざーGarnerのアルゴリズム, 中国剰余定理, Montgomery Reduction
Garnerのアルゴリズム
(
関連: 中国剰余定理
x,y\in\text{整数}\mathbb{Z} M_x,M_y,\overline{M_x}\in\text{自然数}\mathbb{N} \gcd(M_x,M_y)=1 \overline{M_x}M_x\equiv 1\quad(\operatorname{mod}M_y) t=M_x(\overline{M_x}(x-y)\operatorname{mod}M_y) z=\begin{cases}x-t&(x \ge t)\\x-t+M_xM_y&(x \lt t)\\\end{cases}
とすると
0\le t\le M_x(M_y-1) x-M_x(M_y-1)\le x-t\le x+M_xM_y z\equiv x\quad(\operatorname{mod}M_x)\quad\because\quad t\equiv 0\quad(\operatorname{mod}M_x) z\equiv y\quad(\operatorname{mod}M_y)\quad\because\quad t\equiv x-y\quad(\operatorname{mod}M_y) - 特に、
0\le x\lt M_xM_y 0\le z\lt M_xM_y
Montgomery Reduction
(
上の Garner のアルゴリズムからさらに、
\overline{M_y}\in\text{自然数}\mathbb{N} \overline{M_y}M_y\equiv 1\quad(\operatorname{mod}M_x) M_x\lt M_y a,b\in\text{整数}\mathbb{Z} a'=M_ya\operatorname{mod}M_x=\overline{M_y}((M_y)^2\operatorname{mod}M_x)a\operatorname{mod}M_x b'=M_yb\operatorname{mod}M_x=\overline{M_y}((M_y)^2\operatorname{mod}M_x)b\operatorname{mod}M_x 0\le x=a'b'\le (M_x-1)^2 y=0
とすると
z\operatorname{mod}M_y=0 0\le z/M_y\lt M_x z/M_y=\overline{M_y}z\operatorname{mod}M_x=\overline{M_y}x\operatorname{mod}M_x z/M_y=\overline{M_y}a'b'\operatorname{mod}M_x=M_yab\operatorname{mod}M_x
この結果を使うと、
M_ya\times M_yb\times \overline{M_y} \equiv z/M_y\ (\operatorname{mod}M_x)
の形の計算を行うことができる。同じ
この式の前提を改めてまとめると以下の通り。
\gcd(M_x,M_y)=1 M_x\lt M_y \overline{M_x}M_x\equiv 1\quad(\operatorname{mod}M_y) \overline{M_y}M_y\equiv 1\quad(\operatorname{mod}M_x) a'=M_ya\operatorname{mod}M_x=\overline{M_y}((M_y)^2\operatorname{mod}M_x)a\operatorname{mod}M_x b'=M_yb\operatorname{mod}M_x=\overline{M_y}((M_y)^2\operatorname{mod}M_x)b\operatorname{mod}M_x 0\le x=a'b'\le (M_x-1)^2 y=0 t=M_x(\overline{M_x}(x-y)\operatorname{mod}M_y) z=\begin{cases}x-t&(x \ge t)\\x-t+M_xM_y&(x \lt t)\\\end{cases}
(C++, モンゴメリリダクション実装例)
(Rust, モンゴメリリダクション実装例)
Mizar/みざーBarrett Reduction
- AtCoder Library では汎用的な剰余算のライブラリとして使われている
- 同じ除数(AtCoder Library の場合は31bit整数までを想定)で剰余算を多く行う場合に有用、
\bar{m}' = \left\lfloor\frac{2^{64}-1}{m}\right\rfloor+1\mod2^{64}=\left\lceil\frac{2^{64}}{m}\right\rceil\operatorname{mod}2^{64} - 31bit整数 × 31bit整数
\operatorname{mod} -
ab-qm ab qm m \operatorname{mod}
-
-
\operatorname{mod} -
\lfloor x\rfloor \lceil x\rceil
(アルゴリズム)
- when
m=1 a=b=\bar{m}'=0 - when
2\le m\lt 2^{32} 2^{32}+2=\left\lceil\frac{2^{64}}{2^{32}-1}\right\rceil\le\bar{m}'=\left\lceil\frac{2^{64}}{m}\right\rceil\le \left\lceil\frac{2^{64}}{2}\right\rceil=2^{63} \bar{m}'\hspace{.1em}m=2^{64}+r\quad(\lbrace r\rbrace\in\text{整数(integer)}\mathbb{Z},\quad 0\le r\lt m) z = ab = cm + d\quad(\lbrace c,d\rbrace\in\text{整数(integer)}\mathbb{Z},\quad 0\le\lbrace c,d\rbrace\lt m) z\hspace{.1em}\bar{m}'=ab\hspace{.1em}\bar{m}'=(cm+d)\hspace{.1em}\bar{m}'=c(\bar{m}'\hspace{.1em}m)+d\hspace{.1em}\bar{m}'=2^{64}c+c\hspace{.1em}r+d\hspace{.1em}\bar{m}' -
2^{64}c\le ab\hspace{.1em}\bar{m}'\lt 2^{64}(c+2) ab\hspace{.1em}\bar{m}'=2^{64}c+c\hspace{.1em}r+d\hspace{.1em}\bar{m}' 0\le c\hspace{.1em}r\le (m-1)^2\le(2^{32}-2)^2=2^{64}-2^{34}+4\lt 2^{64} 0\le d\hspace{.1em}\bar{m}'\le\bar{m}'\hspace{.1em}(m-1)=2^{64}+r-\bar{m}'\le 2^{64}+(2^{32}-2)-(2^{32}+2)=2^{64}-4\lt 2^{64}
-
x=\left\lfloor\frac{ab\hspace{.1em}\bar{m}'}{2^{64}}\right\rfloor=\lbrace c (c+1)\rbrace -
ab-xm=ab-\left\lfloor\frac{ab\hspace{.1em}\bar{m}'}{2^{64}}\right\rfloor m=\lbrace d (d-m)\rbrace
(元コード)
(C++:
#ifdef _MSC_VER
#include <intrin.h>
#endif
// @param a `0 <= a < m`
// @param b `0 <= b < m`
// @return `a * b % m`
unsigned int barrett_mul_before(unsigned int a, unsigned int b, unsigned int _m, unsigned long long im) {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m)
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
unsigned long long z = a;
z *= b;
#ifdef _MSC_VER
unsigned long long x;
_umul128(z, im, &x);
#else
unsigned long long x =
(unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
unsigned int v = (unsigned int)(z - x * _m);
if (_m <= v) v += _m;
return v;
}
// @param a `0 <= a < m`
// @param b `0 <= b < m`
// @return `a * b % m`
unsigned int barrett_mul_after(unsigned int a, unsigned int b, unsigned int _m, unsigned long long im) {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m)
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
unsigned long long z = a;
z *= b;
#ifdef _MSC_VER
unsigned long long x;
_umul128(z, im, &x);
#else
unsigned long long x =
(unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
unsigned long long y = x * _m;
#ifdef __GNUC__
// https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html
unsigned long long v;
unsigned int w = __builtin_usubll_overflow(z, y, &v) ? _m : 0;
return (unsigned int)(v + w);
#elif defined(_MSC_VER) && defined(_M_AMD64)
// https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#text=_subborrow_u64&ig_expand=7252
unsigned long long v;
unsigned int w = _subborrow_u64(0, z, y, &v) ? _m : 0;
return (unsigned int)(v + w);
#else
return (unsigned int)((z - y) + (z < y ? _m : 0));
#endif
}
(Rust:
/// Calculates `a * b % m`.
///
/// * `a` `0 <= a < m`
/// * `b` `0 <= b < m`
/// * `m` `1 <= m <= 2^31`
/// * `im` = ceil(2^64 / `m`)
#[allow(clippy::many_single_char_names)]
pub fn mul_mod_before(a: u32, b: u32, m: u32, im: u64) -> u32 {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m)
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
let mut z = a as u64;
z *= b as u64;
let x = (((z as u128) * (im as u128)) >> 64) as u64;
let mut v = z.wrapping_sub(x.wrapping_mul(m as u64)) as u32;
if m <= v {
v = v.wrapping_add(m);
}
v
}
/// Calculates `a * b % m`.
///
/// * `a` `0 <= a < m`
/// * `b` `0 <= b < m`
/// * `m` `1 <= m < 2^32`
/// * `im` = ceil(2^64 / `m`) = floor((2^64 - 1) / `m`) + 1
#[allow(clippy::many_single_char_names)]
pub fn mul_mod_after(a: u32, b: u32, m: u32, im: u64) -> u32 {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m) = floor((2^64 - 1) / m) + 1
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
let z = (a as u64) * (b as u64);
let x = (((z as u128) * (im as u128)) >> 64) as u64;
match z.overflowing_sub(x.wrapping_mul(m as u64)) {
(v, true) => (v as u32).wrapping_add(m),
(v, false) => v as u32,
}
}
Mizar/みざー関数 function
fn barrett_reduction(z: u64, m: u64, im: u64) -> (u64, u64)
出力 output
return
入力の制約 input constraints
\left\{z\in\mathbb{Z}\;\middle|\; 0 \le z \lt\max((im-1)\times m,1)\right\} \left\{m\in\mathbb{Z}\;\middle|\; 2 \le m \lt 2^{32.5}\right\} \left\{im\in\mathbb{Z}\;\middle|\; im\equiv\left\lceil 2^{64}/m\right\rceil\mod 2^{64}\equiv\left\lfloor(2^{64}-1)/m\right\rfloor +1\mod 2^{64}\right\}
アルゴリズム algorithm
v=\lfloor(z\times im)/2^{64}\rfloor w=v\times m \lfloor z/m\rfloor=\begin{cases}v&(z\ge w)\\v-1&(z\lt w)\end{cases} z\mod m=\begin{cases}z-w&(z\ge w)\\z-w+m&(z\lt w)\end{cases}
実装例 implementation example
pub fn barrett_reduction(z: u64, m: u64, im: u64) -> (u64, u64) {
// return `(floor(z / m), z \mod m)`
// * `z` : `0 <= z < max((im - 1) * m, 1)`
// * `m` : `1 <= m < 2^32.5`
// * `im = ceil(2^64 / m) mod 2^64 = floor((2^64 - 1) / m) + 1 mod 2^64`
// [1] m = 1
// im = 0, so may fail (must be z = 0 to return correct result)
// [2] m >= 2
// im = ceil(2^64 / m) = floor((2^64 - 1) / m) + 1 >= 2
// -> im * m = 2^64 + r (0 <= r < m)
// let z = c*m + d (0 <= c < im, 0 <= d < m)
// z * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im = (c+1)*r - (m-d)*im + 2^64 <= im*(m-2) + 2^64 = 2 * (2^64 - im) < 2 * 2^64
// ((z * im) >> 64) == c or c + 1
debug_assert_eq!(im, ((!0u64) / m).wrapping_add(1));
debug_assert!(z < im.saturating_sub(1).saturating_mul(m).max(1));
let v = (((z as u128) * (im as u128)) >> 64) as u64;
let (t, f) = z.overflowing_sub(v.wrapping_mul(m));
(
v - (f as u64),
t.wrapping_add((f as u64).wrapping_neg() & m),
)
}
証明 proof
-
m=0 \quad\to\quad im: - ゼロ除算、
im - divide by zero, cannot satisfy the
im
- ゼロ除算、
-
m=1 \quad\to\quad im=0 \quad\to\quad z=0 - 制約を
0\le z\lt 2^{64} z\ne 0 z=0 - If the constraint is
0\le z<2^{64} z\ne 0 z=0
- 制約を
-
m\ge 2 \quad\to\quad im=\lceil 2^{64}/m\rceil=\lfloor(2^{64}-1)/m\rfloor+1 0\le z\lt(im-1)\times m=\lfloor(2^{64}-1)/m\rfloor\times m\le 2^{64}-1\lt 2^{64} 2^{31.5}\lt im\le 2^{63} -
m\times im=2^{64}+r\quad\left\{r\in\mathbb{Z}\;\middle|\; 0\le r\lt m\right\}
\to\quad r-m\times im+2^{64}=0 z=c\times m+d\quad\left\{c\in\mathbb{Z}\;\middle|\; 0\le c\lt im\right\},\left\{d\in\mathbb{Z}\;\middle|\; 0\le d\lt m\right\} \left\{r\in\mathbb{Z}\;\middle|\; 0\le r\lt m\right\}\quad\to\quad r\le m-1 \left\{c\in\mathbb{Z}\;\middle|\; 0\le c\lt im\right\}\quad\to\quad c+1\le im \left\{d\in\mathbb{Z}\;\middle|\; 0\le d\lt m\right\}\quad\to\quad d-m\le -1 z\times im=(c\times m+d)\times im=c\times (m\times im)+d\times im=2^{64}\times c+c\times r+d\times im -
c\times r+d\times im=(c\times r+d\times im)+(r-m\times im+2^{64})=(c+1)\times r+(d-m)\times im+2^{64}
\le im\times (m-1)+(-1)\times im+2^{64}=2\times 2^{64}-2\times im+m-1
\lt 2\times 2^{64}-2\times 2^{31.5}+2^{32.5}-1= 2\times 2^{64}-1\lt 2\times 2^{64} 2^{64}\times c\le z\times im\lt 2^{64}\times(c+2) v=\lfloor(z\times im)/2^{64}\rfloor=c\text{ or }c+1 -
w=\lfloor(z\times im)/2^{64}\rfloor\times m\le\lfloor z/m\rfloor\times m+m
\le (im-1)\times m=\lfloor(2^{64}-1)/m\rfloor\times m\le 2^{64}-1\lt 2^{64}
Mizar/みざー(下書き中)
a. 丁度の値を求める (
b. 丁度か1小さい値を求める (
c. 丁度か1大きい値を求める (
Mizar/みざー(下書き中)
a. 丁度の値を求める (
mulh(m, x)
-
d=1 \displaystyle\left\lfloor\frac{x}{d}\right\rfloor=x - 擬似コード:
x
-
\log_2d=\lceil\log_2d\rceil d 2 \displaystyle t=\log_2d,\quad\left\lfloor\frac{x}{d}\right\rfloor=\left\lfloor\frac{x}{2^{\log_2d}}\right\rfloor - 前計算:
s=\log_2d - 擬似コード:
x >> s
-
\displaystyle\log_2d\ne\lceil\log_2d\rceil,\quad t=\lfloor\log_2((2^W-1)d)\rfloor,\quad\left(-2^t\right)\operatorname{mod}d\lt 2^{(t-W)} \displaystyle t=\lfloor\log_2((2^W-1)d)\rfloor,\quad\left\lfloor\frac{x}{d}\right\rfloor=\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^t}\right\rfloor=\left\lfloor\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^W}\right\rfloor\frac{1}{2^{(t-W)}}\right\rfloor - 前計算:
\displaystyle t=\lfloor\log_2((2^W-1)d)\rfloor,\quad s=t-W,\quad m=\left\lceil\frac{2^t}{d}\right\rceil - 擬似コード:
mulh(m, x) >> s
- 上記以外の場合:
-
\displaystyle t=W+\lceil\log_2d\rceil,\quad y=\left\lfloor\left(\left\lceil\frac{2^t}{d}\right\rceil-2^W\right)\frac{x}{2^W}\right\rfloor,
\displaystyle\left\lfloor\frac{x}{d}\right\rfloor=\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^t}\right\rfloor=\left\lfloor\left(\left\lfloor\frac{x-y}{2}\right\rfloor+y\right)\frac{1}{2^{(t-W-1)}}\right\rfloor - 前計算:
\displaystyle t=W+\lceil\log_2d\rceil,\quad s=t-W-1,\quad m=\left\lceil\frac{2^t}{d}\right\rceil-2^W - 擬似コード:
y = mulh(m, x); (((x - y) >> 1) + y) >> s
-
b. 丁度か
前計算:
擬似コード: ((x >> s) * m) >> t
(x * m) >> t
c. 丁度か
前計算:
擬似コード: (((x + u) >> s) * m) >> t
(x * m) >> t
(a * b) >> 64, x / 998244353, x mod 998244353, x / 7, x mod 7 を計算する関数とそのコンパイル結果の例 (C++): https://godbolt.org/z/6fbWGsEhr
#include <cstdint>
#ifdef _MSC_VER
#include <intrin.h>
#endif
// u64mulh_shim(a, b) == floor(a * b / 2**64)
uint64_t u64mulh_shim(uint64_t a, uint64_t b) {
#ifdef _MSC_VER // msvc https://learn.microsoft.com/ja-jp/cpp/intrinsics/umulh?view=msvc-170
return __umulh(a, b);
#else // gcc/clang
return (unsigned __int128)a * (unsigned __int128)b >> 64;
#endif
}
// div998244353_shim(x) == floor(x / 998244353)
uint64_t div998244353_shim(uint64_t x) {
// 0 <= x < (2**93 / ((-2**93) mod 998244353)) --> floor(ceil(2**93 / 998244353) * x / 2**93) == floor(x / 998244353)
// ceil(2**93 / 998244353) == floor((2**93 - 1) / 998244353) + 1 == 9920937979283557439
// floor(x / 998244353) == floor((x * 9920937979283557439) / 2**93) == floor(floor((x * 9920937979283557439) / 2**64) / 2**29)
return u64mulh_shim(x, 9920937979283557439ULL) >> 29;
}
uint64_t div998244353(uint64_t x) {
return x / 998244353ULL;
}
// rem998244353_shim(x) == (x mod 998244353)
uint64_t rem998244353_shim(uint64_t x) {
return x - div998244353_shim(x) * 998244353ULL;
}
uint64_t rem998244353(uint64_t x) {
return x % 998244353ULL;
}
// div7_shim(x) == floor(x / 7)
uint64_t div7_shim(uint64_t x) {
// 0 <= x < (2**67 / ((-2**67) mod 7)) --> floor(ceil(2**67 / 7) * x / 2**67) == floor(x / 7)
// y := floor((ceil(2**67 / 7) - 2**64) * x / 2**64)
// floor(x / 7) == floor((x + y) / 2**3) == floor((floor((x - y) / 2) + y) / 2**2)
// ceil(2**67 / 7) - 2**64 == 2635249153387078803
uint64_t y = u64mulh_shim(x, 2635249153387078803ULL);
return (((x - y) >> 1) + y) >> 2;
}
uint64_t div7(uint64_t x) {
return x / 7;
}
// rem7_shim(x) == (x mod 7)
uint64_t rem7_shim(uint64_t x) {
return x - div7_shim(x) * 7ULL;
}
uint64_t rem7(uint64_t x) {
return x % 7ULL;
}
# x86-64 gcc
u64mulh_shim(unsigned long, unsigned long):
mov rax, rdi
mul rsi
mov rax, rdx
ret
div998244353_shim(unsigned long):
movabs rax, -8525806094425994177
mul rdi
mov rax, rdx
shr rax, 29
ret
div998244353(unsigned long):
movabs rax, -8525806094425994177
mul rdi
mov rax, rdx
shr rax, 29
ret
rem998244353_shim(unsigned long):
movabs rax, -8525806094425994177
mul rdi
mov rax, rdi
shr rdx, 29
imul rdx, rdx, 998244353
sub rax, rdx
ret
rem998244353(unsigned long):
movabs rax, -8525806094425994177
mul rdi
mov rax, rdx
shr rax, 29
imul rdx, rax, 998244353
mov rax, rdi
sub rax, rdx
ret
div7_shim(unsigned long):
movabs rax, 2635249153387078803
mul rdi
sub rdi, rdx
shr rdi
lea rax, [rdi+rdx]
shr rax, 2
ret
div7(unsigned long):
movabs rax, 2635249153387078803
mul rdi
sub rdi, rdx
shr rdi
lea rax, [rdx+rdi]
shr rax, 2
ret
rem7_shim(unsigned long):
movabs rax, 2635249153387078803
mov rcx, rdi
mul rdi
mov rax, rdi
sub rcx, rdx
shr rcx
add rdx, rcx
shr rdx, 2
lea rcx, [0+rdx*8]
sub rcx, rdx
sub rax, rcx
ret
rem7(unsigned long):
movabs rax, 2635249153387078803
mul rdi
mov rax, rdi
sub rax, rdx
shr rax
add rax, rdx
shr rax, 2
lea rdx, [0+rax*8]
sub rdx, rax
mov rax, rdi
sub rax, rdx
ret
# ARM64 gcc
u64mulh_shim(unsigned long, unsigned long):
umulh x0, x0, x1
ret
div998244353_shim(unsigned long):
mov x1, 52287
movk x1, 0x5de0, lsl 16
movk x1, 0x4087, lsl 32
movk x1, 0x89ae, lsl 48
umulh x0, x0, x1
lsr x0, x0, 29
ret
div998244353(unsigned long):
mov x1, 52287
movk x1, 0x5de0, lsl 16
movk x1, 0x4087, lsl 32
movk x1, 0x89ae, lsl 48
umulh x0, x0, x1
lsr x0, x0, 29
ret
rem998244353_shim(unsigned long):
mov x2, 52287
movk x2, 0x5de0, lsl 16
movk x2, 0x4087, lsl 32
movk x2, 0x89ae, lsl 48
umulh x2, x0, x2
lsr x2, x2, 29
lsl x1, x2, 4
sub x1, x1, x2
lsl x1, x1, 3
sub x1, x1, x2
add x1, x2, x1, lsl 23
sub x0, x0, x1
ret
rem998244353(unsigned long):
mov x2, 52287
movk x2, 0x5de0, lsl 16
movk x2, 0x4087, lsl 32
movk x2, 0x89ae, lsl 48
umulh x2, x0, x2
lsr x2, x2, 29
lsl x1, x2, 4
sub x1, x1, x2
lsl x1, x1, 3
sub x1, x1, x2
add x1, x2, x1, lsl 23
sub x0, x0, x1
ret
div7_shim(unsigned long):
mov x1, 9363
movk x1, 0x9249, lsl 16
movk x1, 0x4924, lsl 32
movk x1, 0x2492, lsl 48
umulh x1, x0, x1
sub x0, x0, x1
add x0, x1, x0, lsr 1
lsr x0, x0, 2
ret
div7(unsigned long):
mov x1, 9363
movk x1, 0x9249, lsl 16
movk x1, 0x4924, lsl 32
movk x1, 0x2492, lsl 48
umulh x1, x0, x1
sub x0, x0, x1
add x0, x1, x0, lsr 1
lsr x0, x0, 2
ret
rem7_shim(unsigned long):
mov x2, 9363
movk x2, 0x9249, lsl 16
movk x2, 0x4924, lsl 32
movk x2, 0x2492, lsl 48
umulh x2, x0, x2
sub x1, x0, x2
add x1, x2, x1, lsr 1
lsr x1, x1, 2
lsl x2, x1, 3
sub x1, x2, x1
sub x0, x0, x1
ret
rem7(unsigned long):
mov x2, 9363
movk x2, 0x9249, lsl 16
movk x2, 0x4924, lsl 32
movk x2, 0x2492, lsl 48
umulh x2, x0, x2
sub x1, x0, x2
add x1, x2, x1, lsr 1
lsr x1, x1, 2
lsl x2, x1, 3
sub x1, x2, x1
sub x0, x0, x1
ret
(下書き中)
除算の商を求める3つのアプローチ
この節では、除算の商
- a. 丁度の値を求める
\lfloor x/d\rfloor=f_d(x)
- b. 丁度か
1 \lfloor x/d\rfloor-1\leq f_d(x)\leq\lfloor x/d\rfloor
- c. 丁度か
1 \lfloor x/d\rfloor\leq f_d(x)\leq\lfloor x/d\rfloor+1
コンパイラによる最適化で使われるパターンは主に 「a. 丁度の値を求める」 ですが、
後者の2つのアプローチは、求めた除算の商の近似値から、仮に除算の剰余を計算してみて、商が大きすぎたり小さすぎたりしたと分かれば商の近似値から1だけ調整して、正しい除算の商とするものです。後ほど、 「b. 丁度か
また、以下のような略記を所々で使うので定義しておきます。
-
mulh(m, x)\displaystyle=\left\lfloor\frac{mx}{2^W}\right\rfloor\quad\left(m,x\in\mathbb{Z},\quad 0\le m\lt 2^W,\quad 0\le x\lt 2^W\right) -
W W W=64
-
-
\displaystyle\left(-2^t\right)\operatorname{mod}d=(d-1)-\left(\left(2^t-1\right)\operatorname{mod}d\right)\quad\left(d,t\in\mathbb{Z},\quad 1\leq d\lt 2^W,\quad 0\leq t\right) - 負の被除数に対する剰余計算。ここでは主に、「
((2^t\text{ 以上で }d\text{ の倍数となる最小の整数}) - 2^t) - 逆に、「
(2^t-(2^t\text{ 以下で }d\text{ の倍数となる最大の整数})) 2^t\operatorname{mod}d
- 負の被除数に対する剰余計算。ここでは主に、「
a. 丁度の値を求める (d: x: W:
2^W
-
d=1 \displaystyle\left\lfloor\frac{x}{d}\right\rfloor=x - 擬似コード:
x
-
\log_2d=\lceil\log_2d\rceil d 2 \displaystyle t=\log_2d,\quad\left\lfloor\frac{x}{d}\right\rfloor=\left\lfloor\frac{x}{2^{\log_2d}}\right\rfloor - 前計算:
s=\log_2d - 擬似コード:
x >> s
-
\displaystyle 2\leq d\lt 2^W,\quad t=\lfloor\log_2(d-1)\rfloor+W,\quad\left((-2^t)\operatorname{mod}d\right)x\lt 2^t 0\leq x\lt 2^{W-1}\quad\longrightarrow\quad\left((-2^t)\operatorname{mod}d\right)x\lt 2^t 0\leq x\lt 2^W,\quad (-2^t)\operatorname{mod}d\lt 2^{(t-W)}\quad\longrightarrow\quad \left((-2^t)\operatorname{mod}d\right)x\lt 2^t \displaystyle 2^{(t-W)}\leq(d-1)\lt d\leq 2^{(t-W+1)},\quad 2^{(W-1)}\leq\lceil 2^t/d\rceil\lt 2^W,\quad 0\leq t-W=\lfloor\log_2(d-1)\rfloor\lt W, \displaystyle \left\lfloor\frac{x}{d}\right\rfloor=\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^t}\right\rfloor=\left\lfloor\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^W}\right\rfloor\frac{1}{2^{(t-W)}}\right\rfloor - 前計算:
\displaystyle s=\lfloor\log_2(d-1)\rfloor,\quad t=s+W,\quad m=\left\lceil 2^t/d\right\rceil - 擬似コード:
mulh(m, x) >> s
-
\left(2\leq d\lt 2^W,\quad 0\leq x\lt 2^W\right) \displaystyle t=\lfloor\log_2(d-1)\rfloor+W+1,\quad \displaystyle 2^{(t-W-1)}\leq(d-1)\lt d\leq 2^{(t-W)},\quad 2^W\leq\lceil 2^t/d\rceil\lt 2^{(W+1)},\quad \displaystyle y=\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^W}\right\rfloor-x=\left\lfloor\left(\left\lceil\frac{2^t}{d}\right\rceil-2^W\right)\frac{x}{2^W}\right\rfloor,\quad 0\leq y\leq x\lt 2^W,\quad \displaystyle\left\lfloor\frac{x}{d}\right\rfloor=\left\lfloor\left\lceil\frac{2^t}{d}\right\rceil\frac{x}{2^t}\right\rfloor=\left\lfloor\frac{x+y}{2^{(t-W)}}\right\rfloor=\left\lfloor\left(\left\lfloor\frac{x-y}{2}\right\rfloor+y\right)\frac{1}{2^{(t-W-1)}}\right\rfloor - 前計算:
\displaystyle s=\lfloor\log_2(d-1)\rfloor,\quad t=s+W+1,\quad m=\left\lceil 2^t/d\right\rceil-2^W - 擬似コード:
y = mulh(m, x); (((x - y) >> 1) + y) >> s
例
b. 丁度か 1 d: x: s,t: x
- 前計算:
\displaystyle m=\left\lfloor\frac{2^{(s+t)}}{d}\right\rfloor - 擬似コード:
((x >> s) * m) >> t -
s=0 (m * x) >> t -
s=0,t=W mulh(m, x)※d=1 m=\lfloor 2^W/d\rfloor=2^W
前半の導出
例
c. 丁度か 1 d: x: s,t: x
- 前計算:
\displaystyle u=2^s-1,\quad m=\left\lceil\frac{2^{(s+t)}}{d}\right\rceil - 擬似コード:
(((x + u) >> s) * m) >> t -
s=0 (x * m) >> t -
s=0, t=W mulh(m, x)※d=1 m=\lceil 2^W/d\rceil=2^W
前半の導出