■
VC++が定数による除算を乗算に変換する最適化をしてるのは、吐いたコードを見たことある人なら知ってることだと思うけど、なんとなく
このスレの400を見てて、確かに正確な変換方法は知らないので考えてみる。
mov eax, 40140141h ; r=esi/819とおく。esi=819r mul esi ;edx=205r sub esi, edx ;esi=614r shr esi, 1 ;esi=307r add esi, edx ;esi=512r shr esi, 9 ;esi=r これなんか、esi/=819; だぜ。人智を超えている。まあ、人智の生んだコンパイラだが。 40140141hは2^32*205/819で、205の大きさで精度を上げて9bitシフトで仕上げている。 一般的な式を出したいな。
まず、このコードについてるコメントは間違ってるので正確に書くと
mov eax, 40140141h ; r = esi/819とおく。esi = 819r + rem (0<=rem<819) mul esi ; edx = 205r + e (0<=e<=205) sub esi, edx ; esi = (819r + rem) - 205r - e = 614r + rem - e shr esi, 1 ; esi = 307r + (rem - e)/2 add esi, edx ; esi = (307r + (rem - e) / 2) + 205r + e= 512r + (rem + e) / 2 shr esi, 9 ; esi = r + ((rem + e) / 2)/512 = r
のようにちゃんと余りも書かないとだめ。esi = 2^32 - 1 - 2^32 % 819 のときに e = 205 になる。
e + rem < 1024 なのがポイント。
とりあえず符号なしとする。除数をMとして(3<=M<2^32)とする。
Mで割る代わりにTを掛けるようにしたいので、f(M) = F(M) * 2^32 を M で割った商+1を T とすることを考える
この f がみたすべき条件を考えると
f(M) = (T-1) * M + R, (0 <= R < f(M))
ここで、Tは32ビットに収まってないといけないので
f(M) < M * (2^32 - 1) F(M) <= M - 1
被除数をm(0<=m<2^32)として除算の結果を、m = k * M + r とする。
実際に T を掛けてみると、
m * T = m * ((f(M) - R) div M + 1)
= (m * (f(M) - R + M)) div M
= k * f(M) + floor((r * f(M) + m * (M - R)) / M)となる。
(M - R) >= 0 より
floor((r * f(M) + m * (M - R)) / M) >= 0
は常に成立している。さらに
floor((r * f(M) + m * (M - R)) / M) < f(M) + 2^32
も成立する。よって
k * F(M) <= m * T div 2^32 <= (k+1) * F(M)
である。
後ろの等号が成立しない場合に制限したければ
floor((r * f(M) + m * (M - R)) / M) < f(M) <=> f(M)/(M - f(M) mod M) > m / (M - m mod M)
をみたさなければならない。
右辺は m = 2^32-1 か m = 2^32-1-2^32%M で最大値をとる。
F(M)は2の冪になってるのが望ましいので、M-1以下の2の冪でこれを満たすものを見つければよい。
見つからない場合は、
a < M < 2*a となる2の冪 a に対して、F(M) = 2 * a - M とする。このときは後ろの等号が成立するかもしれないので
s := m * T div 2^32 x := s - k * F(M) y := m - s = k * (M - F(M)) + r - x = k * 2 * (M - a) + r - x z := y div 2 = k * (M - a) + (r-x) div 2 w := z + s = k * a + (r - x) div 2 + x = k * a + (r + x) div 2
と計算すると、ちょうど
0 <= x <= F(M) = 2 * a - M 0 <= r < M 0 <= r + x < 2 * a 0 <= (r + x) div 2 < a
となり、w / a = k を得る。
