unsigned の exact division などをうまくやる

背景としては、こちらの記事です。

rsk0315.hatenablog.com

1431655766, 1717986920, 1431655782, 1431655768, 1840700272, 1431655966, -1431655056 などさまざまなマジックナンバーが登場します。 「$x$ が $3$ の倍数であることがわかっているとき、$1431655766x \equiv \tfrac23 x\pmod{2^{32}}$ となる」のような性質を持ちます。

一般に除算は重い処理なので、結果が同じになるような別の命令（乗算など）に置き換える最適化が典型的に行われるのですが、その過程で出てきた値になります。

今回は、そうした値を自分で求められるようになってみようという回です（コンパイル時定数であれば、多くの場合はコンパイラが勝手に最適化してくれるとは思います）。

exact division というのは、あまりの出ない除算のことを指しています。

本題

まず、上記の記事でも書いていますが、$1431655766 = \tfrac13(2^{32}+2)$ です。$2^{32}+2\equiv 0\pmod{3}$ であることに注意しておきましょう。 $$ \begin{aligned} 3x\times 1431655766 &= 3x\times\tfrac13(2^{32}+2) \\ &= x\times (2^{32}+2) \\ &\equiv x\times(0+2) = 2x \pmod{2^{32}} \end{aligned} $$ となるわけです。

以降、$ax\times n\equiv bx \pmod{2^{32}}$ となることを $ax\xmapsto{n} bx$ と表すことにします。上記で言えば $3x\xmapsto{1431655766}2x$ です。

例

では、手始めに $3x\xmapsto{n}1x=x$ となるような $n$ について考えてみましょう。 合同式で $2^{32}$ が消えてくれることだけに目が行って、$\tfrac13(2^{32}+1)$ か？と思ってしまうと危ないです。$2^{32}+1\not\equiv 0\pmod{3}$ なので、$\tfrac13(2^{32}+1)$ が整数になってくれないんですね。

というわけで、$q\cdot 2^{32}+1 \equiv 0\pmod{3}$ なる $q$ があってくれるとうれしいわけです。 これは競プロ er ならもう全員わかっていて*1、 $$q \equiv (-1)\cdot (2^{32})^{-1} \pmod{3}$$ ですね。逆元を使いたくなったときは $\gcd(2^{32}, 3) = 1$ であることをちゃんと意識した上でやりましょう。

$3^{-1}\bmod 2^{32}$ であれば Newton 法を使ってうまくやる方法があるのですが、$(2^{32})^{-1} \bmod 3$ だとうまい方法はあるんでしょうか。わかりません。$x^{-1} \bmod y$ がわかっても $y^{-1}\bmod x$ を求めるのには役立たなさそうな気がします。

というわけで Euclid の互除法なりなんなりで求めると、$(2^{32})^{-1} \bmod 3 = 1$ とわかります（そもそも $2^{32}\equiv 1\pmod 3$ なので）。 よって、$q \equiv -1 \equiv 2$ とわかります。

すなわち、$3x\xmapsto{\tfrac13(2\cdot 2^{32}+1)}x$ というわけです。${\tfrac13(2\cdot 2^{32}+1)}=2863311531$ です。

一般化

$a=1$ の場合は単に $x\xmapsto{b} bx$ とできるので、以降は除外します。

奇数 $a\gt 1$ に対して、$2^k\cdot ax\xmapsto{n} bx$ なる $n$ を求めたいです。 そもそもシフト演算は高速ですから、$(2^k\cdot ax)\gg k = ax$ を用いて、$ax\xmapsto{n}bx$ に帰着させることができます。 $a$ は奇数なので、$\gcd(2^{32}, a) = 1$ が保証されて、うれしいです。

やりたいことは $\tfrac 1a(q\cdot 2^{32}+b)$ の形の整数を探すことで、それは $q\cdot 2^{32}+b \equiv 0 \pmod a$ を解くことによって得られます。 $$q \equiv (-b)\cdot (2^{32})^{-1} \pmod a$$ ですね。

気づいてしまったのですが、 $$ q \equiv (-b) \cdot (2^{-1})^{32} \pmod a $$ であり、$a$ が奇数であることから $2^{-1}\bmod a = \tfrac{a+1}2$ なので、 $$ q \equiv (-b) \cdot (\tfrac{a+1}2)^{32} \pmod a $$ ですね。32 乗の部分は繰り返し二乗法っぽくやりましょう。

warning: $a = 2^{32}-1$ のときは $\tfrac{a+1}2$ の計算過程でオーバーフローしうるので注意しましょう。$\floor{\tfrac a2}+1$ とするのがよいでしょうか。

不等式評価

より注意深くなるために、$\tfrac1a(q\cdot 2^{32}+b)$ の大きさについて見積もっておきましょうか。 $1\le q\le a-1$ であることと、$a$, $b$ が 32-bit 符号なし整数型であることから、 $$ \begin{aligned} \tfrac1a(q\cdot 2^{32} + b) &\le \tfrac1a( (a-1)\cdot 2^{32} + (2^{32}-1)) \\ &= \tfrac{a-1}a\cdot 2^{32} + \tfrac1a\cdot2^{32}-\tfrac1a \\ &\lt 2^{32} \end{aligned} $$ です。整数性から、$\tfrac1a(q\cdot 2^{32}+b)\le 2^{32}-1$ ですね。よって、予め前計算を 64-bit 整数などで行った後は、32-bit に収まる値になっていることがわかります*2。

なお、$a + b = 2^{32}$ となるケースでは $\tfrac1a(q\cdot 2^{32}+b) = 2^{32}-1$ となるようです。 $$ \begin{aligned} q &\equiv (-b) \cdot (2^{32})^{-1} \\ &= -(2^{32}-a)\cdot (2^{32})^{-1} \\ &\equiv -(2^{32})\cdot (2^{32})^{-1} \equiv -1 \pmod{a}. \end{aligned} $$ よって、 $$ \begin{aligned} \tfrac1a(q\cdot 2^{32}+b) &\equiv \tfrac1a(-2^{32}+b) \\ &= \tfrac1a(-a) \\ &= -1 \equiv 2^{32}-1 \pmod {2^{32}} \end{aligned} $$ となります。

具体例

というわけで、もう好きなようにできてしまうわけです。 試しに $271x\xmapsto{n} 314x$ となるような $n$ でも求めてみましょうか。 $$ \begin{aligned} q &\equiv (-314)\cdot (\tfrac{271+1}2)^{32} \\ &\equiv 228 \cdot 136^{32} \\ &\equiv 228\cdot 99 \equiv 79 \pmod {271} \end{aligned} $$ と出ました。よって、所望のマジックナンバー $n$ は $\tfrac1{271}(79\cdot 2^{32}+314) = 1252038438$ になりそうです。

>>> 123456760 % 271
0
>>> 123456760 * 1252038438 % 2**32
143045840
>>> 123456760 // 271 * 314
143045840

あっていそうですね。

おわり

こういう類のことを知っていると、コンパイラの出したアセンブリに不思議な数が入っていても驚かずに済んだりしそうです。

コンパイラの出したアセンブリを眺めて気になったところについて考えてみると、新しい発見があって面白いですね。

*1:わからなかった人はこれから競プロ er になりましょう。

*2:とはいえ、仮にオーバーフローしたとしても $2^{32}$ を法として乗算しているだけなので、$n\bmod 2^{32}$ を考えればいいだけのような気もしてきます？（たぶん）。