背景
一般に、除算命令の無いCPUで整数の除算を行うのは、重い処理になりがちである。
しかし、「10で割る」というのは、コンピュータ向けの二進数から人間向けの十進数に変換するときなど、よく行われる処理である。
2の正の整数乗での割り算は、ビットシフトで表現できるので、除算命令が無くても比較的簡単にできることが多い。
そこで、10で割った値を、愚直に除算を行うよりも簡単な処理で求めたい。
最初のアイデア
2の正の整数乗を十進数で表したとき、その1の位は必ず2・4・6・8のいずれかである。
偶数しか無いうえ、0 も無いため 10 の倍数にはならない。
しかし、1の位が「2」のときは、10で割ると「0.2」になる。
「0.2」を表現するのは難しいが、「0.25」を掛けるのであれば2ビット右シフトすることで表現できる。
ある程度大きい数を使ったほうが、誤差は少なくなるだろう。
という考察を踏まえ、819.25 / 8192 を計算してみると、0.100006103515625 となった。
0.1 のあとにある程度 0 が並んでいるので、これで10で割るのに近い値が得られるかな……？
最初のアイデアを検証
実際に「819.25を掛けて8192で割る」を行った値と、10で割った値を比較してみる。
大きい値を入力したときには誤差が出ても、小さい値の入力なら正しい値が出てくるのであれば、十分使えることが期待できる。
どのくらいの値まで耐えられるか、やってみよう。
#include <stdio.h>
#include <limits.h>
unsigned int div10(unsigned int a) {
return (a * 819ull + (a >> 2)) >> 13;
}
int main(void) {
for (unsigned int i = 0; ; i++) {
if (div10(i) != i / 10) {
printf("%u -> %u\n", i, div10(i));
break;
}
if (i == UINT_MAX) break;
}
return 0;
}
これを実行すると、
16389 -> 1639
と出力された。
すなわち、0～16388 を入力したときは正しい 10 で割った商が得られたが、16389 を入力すると正しい商 1638 より 1 大きい値が出てしまった。
せめて16ビット (65535) までは耐えてほしかったな……
大きすぎる値が出るなら減らしてみよう
もともと、819.25 / 8192 の値は 0.100006103515625 であり、0.1 より少し大きい。
入力を大きくすると出力も欲しい値より大きくなるのは、このためだろう。
出力が少し大きくなるのであれば、少し引いて補正してみたらどうだろうか？
やってみよう。
この章では、前章のコードから div10 関数のみを変えたコードを実行するので、div10 関数と出力結果のみを示す。
unsigned int div10(unsigned int a) {
return (a * 819ull + (a >> 2) - (a >> 3)) >> 13;
}
10 -> 0
引きすぎたようだ。
unsigned int div10(unsigned int a) {
return (a * 819ull + (a >> 2) - (a >> 4)) >> 13;
}
50 -> 4
まだ引きすぎのようだ。
unsigned int div10(unsigned int a) {
return (a * 819ull + (a >> 2) - (a >> 5)) >> 13;
}
43669 -> 4367
おっ、正しい値が得られる範囲が伸びたぞ！
unsigned int div10(unsigned int a) {
return (a * 819ull + (a >> 2) - (a >> 5) - (a >> 6)) >> 13;
}
261629 -> 26163
正しい値が得られる範囲がさらに伸び、目標の16ビットを突破したぞ！
unsigned int div10(unsigned int a) {
return (a * 819ull + (a >> 2) - (a >> 5) - (a >> 6) - (a >> 9) - (a >> 10)) >> 13;
}
4175869 -> 417587
さらに試すと、どうやらシフト幅を2連続で使い、その後2あけてまた2連続で使い……とすると、どんどん正しい値が得られる範囲が伸びるみたいだ！
unsigned int div10(unsigned int a) {
return (
a * 819ull + (a >> 2) - (a >> 5) - (a >> 6) - (a >> 9) - (a >> 10)
- (a >> 13) - (a >> 14)
) >> 13;
}
66650109 -> 6665011
いいぞいいぞ！
unsigned int div10(unsigned int a) {
return (
a * 819ull + (a >> 2) - (a >> 5) - (a >> 6) - (a >> 9) - (a >> 10)
- (a >> 13) - (a >> 14) - (a >> 17) - (a >> 18)
) >> 13;
}
1063780349 -> 106378035
もう少し！
unsigned int div10(unsigned int a) {
return (
a * 819ull + (a >> 2) - (a >> 5) - (a >> 6) - (a >> 9) - (a >> 10)
- (a >> 13) - (a >> 14) - (a >> 17) - (a >> 18) - (a >> 21) - (a >> 22)
) >> 13;
}
ついに、32ビット符号なし整数の全範囲で正しい値が得られるようになった！
考察と整理
この「2連続で使い、その後2あけてまた2連続で使い」という特徴的なパターンは、何を表しているのだろうか？
Interactive visualization of the floating-point format
で 0.1 を変換してみると、以下の結果が得られた。
| 要素 | 値 |
|---|---|
| Sign | 0 |
| Exponent | 0 1 1 1 1 0 1 1 |
| mantissa | 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 1 |
Sign や Exponent は大文字から始まるが、mantissa は小文字から始まる。
これは原文ママである。
これは、十進数の 0.1 は、二進数では 1.10011001100110011001101… × $2^{-4}$ と表せることを示している。
この「11001100…」というパターンが、「2連続で使い、その後2あけてまた2連続で使い」の正体であると考えられる。
さらに、十進数の 819 も、二進数で表すと 1100110011 となり、パターンに沿っていることがわかる。
ということは、今回やったのは二進数の 1100110011.01 - 0.00001100110011… を掛けて $2^{13}$ で割る、という操作だったということになる。
ところで、(a >> 2) のところだけ1ビットだけ単独で使っており、一見パターンに沿っていないように見える。
しかし、ここまでは足し算、この後は引き算である。
1ビットだけ単独で立っているところから、その下の2ビット連続で使うところを引くと、
11010000
- 11
----------
11001101
となり、「2ビット連続」が1セット増え、その下に新しい「1ビットだけ単独」が増えることがわかる。
2ビットの連続を引くと精度が上がったのは、このようにパターンを伸ばす効果があったからだろう。
さて、「2ビット連続」を引くことでパターンを伸ばせるなら、もっと上位からこれを用いてパターンを伸ばしてもよいだろう。
ついでに、「2ビット連続」をバラバラに書かず、まとめて処理しよう。
unsigned int div10(unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
return (
((unsigned long long)a << 10) - (a3 << 6) - (a3 << 2) - (a3 >> 2)
- (a3 >> 6) - (a3 >> 10) - (a3 >> 14) - (a3 >> 18) - (a3 >> 22)
) >> 13;
}
これでも、32ビット符号なし整数の全範囲で正しい値を得ることができた。
シフトの幅を減らす
最初は、十進数でみてよさそうな性質が見えた「819.25を掛けて8192で割る」を用いたが、二進数でのパターンが見えた今では「8192で割る」にこだわらなくて良さそうである。
もう少し小さい値で割るようにしても、精度を確保できるだろうか？
たとえば、シフト幅を1バイト捨てるだけで実現できる 8 にしてみた。
unsigned int div10(unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
return (
((unsigned long long)a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7)
- (a3 >> 11) - (a3 >> 15) - (a3 >> 19) - (a3 >> 23) - (a3 >> 27)
) >> 8;
}
これでも、32ビット符号なし整数の全範囲で正しい値を得ることができた。
シフト幅を 7 にすると、途中の値に対するシフトの幅が 4 の倍数になり、きれいである。
unsigned int div10(unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
return (
((unsigned long long)a << 4) - a3 - (a3 >> 4) - (a3 >> 8)
- (a3 >> 12) - (a3 >> 16) - (a3 >> 20) - (a3 >> 24) - (a3 >> 28)
) >> 7;
}
これでも、32ビット符号なし整数の全範囲で正しい値を得ることができた。
シフト幅をさらにもう1ビット削ってみる。
unsigned int div10(unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
return (
((unsigned long long)a << 3) - (a3 >> 1) - (a3 >> 5) - (a3 >> 9)
- (a3 >> 13) - (a3 >> 17) - (a3 >> 21) - (a3 >> 25) - (a3 >> 29)
) >> 6;
}
9786709 -> 978671
さすがに削りすぎたようである。
目標を16ビットに下げれば、これで十分である。
unsigned int div10(unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
return (
((unsigned long long)a << 3) - (a3 >> 1) - (a3 >> 5) - (a3 >> 9) - (a3 >> 13)
) >> 6;
}
120149 -> 12015
もう1ビット削ると、かなり役立たなくなってきた。
こうなると、引く項の数を増やしても、精度は伸びないようである。
unsigned int div10(unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
return (((unsigned long long)a << 2) - (a3 >> 2) - (a3 >> 6) - (a3 >> 10)) >> 5;
}
2389 -> 239
10で割った余りを求める
今回の方法で、32ビット符号なし整数を10で割った商を求めることができた。
しかし、二進数を十進数に変換する処理では、10で割った余りも欲しい。
そこで、最後の右シフトで捨てられる情報から、10で割った余りが得られないかを考えてみる。
シフト幅を8ビットにしたバージョンのプログラムを変更し、10で割った余りと「最後の右シフトで捨てられる情報」の関係を観察してみる。
#include <stdio.h>
#include <limits.h>
unsigned int div10(unsigned int *dropped, unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
unsigned long long value =
((unsigned long long)a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7)
- (a3 >> 11) - (a3 >> 15) - (a3 >> 19) - (a3 >> 23) - (a3 >> 27);
*dropped = value & 0xff;
return value >> 8;
}
unsigned int dropped_count[10][256];
int main(void) {
for (unsigned int i = 0; ; i++) {
unsigned int res, dropped;
res = div10(&dropped, i);
dropped_count[i % 10][dropped]++;
if (res != i / 10) {
printf("%u -> %u\n", i, res);
break;
}
if (i == UINT_MAX) break;
}
for (int a = 0; a < 256; a++) {
printf("%3d %02x", a, a);
for (int b = 0; b < 10; b++) {
printf("%10u", dropped_count[b][a]);
}
putchar('\n');
}
return 0;
}
以下の実行結果が得られた。
実行結果
0 00 1 0 0 0 0 0 0 0 0 0
1 01 68440 0 0 0 0 0 0 0 0 0
2 02 3079000 0 0 0 0 0 0 0 0 0
3 03 22056812 0 0 0 0 0 0 0 0 0
4 04 52269164 0 0 0 0 0 0 0 0 0
5 05 65706336 0 0 0 0 0 0 0 0 0
6 06 67094688 0 0 0 0 0 0 0 0 0
7 07 67103390 0 0 0 0 0 0 0 0 0
8 08 66210966 0 0 0 0 0 0 0 0 0
9 09 55320624 0 0 0 0 0 0 0 0 0
10 0a 26102512 0 0 0 0 0 0 0 0 0
11 0b 4351265 0 0 0 0 0 0 0 0 0
12 0c 133504 0 0 0 0 0 0 0 0 0
13 0d 28 0 0 0 0 0 0 0 0 0
14 0e 0 0 0 0 0 0 0 0 0 0
15 0f 0 0 0 0 0 0 0 0 0 0
16 10 0 0 0 0 0 0 0 0 0 0
17 11 0 0 0 0 0 0 0 0 0 0
18 12 0 0 0 0 0 0 0 0 0 0
19 13 0 0 0 0 0 0 0 0 0 0
20 14 0 0 0 0 0 0 0 0 0 0
21 15 0 0 0 0 0 0 0 0 0 0
22 16 0 0 0 0 0 0 0 0 0 0
23 17 0 0 0 0 0 0 0 0 0 0
24 18 0 0 0 0 0 0 0 0 0 0
25 19 0 0 0 0 0 0 0 0 0 0
26 1a 0 211 0 0 0 0 0 0 0 0
27 1b 0 244952 0 0 0 0 0 0 0 0
28 1c 0 5975796 0 0 0 0 0 0 0 0
29 1d 0 30327768 0 0 0 0 0 0 0 0
30 1e 0 57957875 0 0 0 0 0 0 0 0
31 1f 0 66556472 0 0 0 0 0 0 0 0
32 20 0 67107064 0 0 0 0 0 0 0 0
33 21 0 67076258 0 0 0 0 0 0 0 0
34 22 0 64996010 0 0 0 0 0 0 0 0
35 23 0 48829888 0 0 0 0 0 0 0 0
36 24 0 18278976 0 0 0 0 0 0 0 0
37 25 0 2112854 0 0 0 0 0 0 0 0
38 26 0 32606 0 0 0 0 0 0 0 0
39 27 0 0 0 0 0 0 0 0 0 0
40 28 0 0 0 0 0 0 0 0 0 0
41 29 0 0 0 0 0 0 0 0 0 0
42 2a 0 0 0 0 0 0 0 0 0 0
43 2b 0 0 0 0 0 0 0 0 0 0
44 2c 0 0 0 0 0 0 0 0 0 0
45 2d 0 0 0 0 0 0 0 0 0 0
46 2e 0 0 0 0 0 0 0 0 0 0
47 2f 0 0 0 0 0 0 0 0 0 0
48 30 0 0 0 0 0 0 0 0 0 0
49 31 0 0 0 0 0 0 0 0 0 0
50 32 0 0 0 0 0 0 0 0 0 0
51 33 0 0 0 0 0 0 0 0 0 0
52 34 0 0 21778 0 0 0 0 0 0 0
53 35 0 0 1728896 0 0 0 0 0 0 0
54 36 0 0 16513408 0 0 0 0 0 0 0
55 37 0 0 46979340 0 0 0 0 0 0 0
56 38 0 0 64548102 0 0 0 0 0 0 0
57 39 0 0 67061144 0 0 0 0 0 0 0
58 3a 0 0 67107912 0 0 0 0 0 0 0
59 3b 0 0 66682584 0 0 0 0 0 0 0
60 3c 0 0 59119013 0 0 0 0 0 0 0
61 3d 0 0 32476656 0 0 0 0 0 0 0
62 3e 0 0 6932284 0 0 0 0 0 0 0
63 3f 0 0 325144 0 0 0 0 0 0 0
64 40 0 0 469 0 0 0 0 0 0 0
65 41 0 0 0 0 0 0 0 0 0 0
66 42 0 0 0 0 0 0 0 0 0 0
67 43 0 0 0 0 0 0 0 0 0 0
68 44 0 0 0 0 0 0 0 0 0 0
69 45 0 0 0 0 0 0 0 0 0 0
70 46 0 0 0 0 0 0 0 0 0 0
71 47 0 0 0 0 0 0 0 0 0 0
72 48 0 0 0 0 0 0 0 0 0 0
73 49 0 0 0 0 0 0 0 0 0 0
74 4a 0 0 0 0 0 0 0 0 0 0
75 4b 0 0 0 0 0 0 0 0 0 0
76 4c 0 0 0 0 0 0 0 0 0 0
77 4d 0 0 0 7 0 0 0 0 0 0
78 4e 0 0 0 96384 0 0 0 0 0 0
79 4f 0 0 0 3673832 0 0 0 0 0 0
80 50 0 0 0 24051344 0 0 0 0 0 0
81 51 0 0 0 53845556 0 0 0 0 0 0
82 52 0 0 0 65981338 0 0 0 0 0 0
83 53 0 0 0 67099904 0 0 0 0 0 0
84 54 0 0 0 67099904 0 0 0 0 0 0
85 55 0 0 0 65981338 0 0 0 0 0 0
86 56 0 0 0 53845556 0 0 0 0 0 0
87 57 0 0 0 24051344 0 0 0 0 0 0
88 58 0 0 0 3673832 0 0 0 0 0 0
89 59 0 0 0 96384 0 0 0 0 0 0
90 5a 0 0 0 7 0 0 0 0 0 0
91 5b 0 0 0 0 0 0 0 0 0 0
92 5c 0 0 0 0 0 0 0 0 0 0
93 5d 0 0 0 0 0 0 0 0 0 0
94 5e 0 0 0 0 0 0 0 0 0 0
95 5f 0 0 0 0 0 0 0 0 0 0
96 60 0 0 0 0 0 0 0 0 0 0
97 61 0 0 0 0 0 0 0 0 0 0
98 62 0 0 0 0 0 0 0 0 0 0
99 63 0 0 0 0 0 0 0 0 0 0
100 64 0 0 0 0 0 0 0 0 0 0
101 65 0 0 0 0 0 0 0 0 0 0
102 66 0 0 0 0 0 0 0 0 0 0
103 67 0 0 0 0 5474 0 0 0 0 0
104 68 0 0 0 0 897898 0 0 0 0 0
105 69 0 0 0 0 11788240 0 0 0 0 0
106 6a 0 0 0 0 41006352 0 0 0 0 0
107 6b 0 0 0 0 62757599 0 0 0 0 0
108 6c 0 0 0 0 66975360 0 0 0 0 0
109 6d 0 0 0 0 67108752 0 0 0 0 0
110 6e 0 0 0 0 66926744 0 0 0 0 0
111 6f 0 0 0 0 61991961 0 0 0 0 0
112 70 0 0 0 0 38910104 0 0 0 0 0
113 71 0 0 0 0 10417072 0 0 0 0 0
114 72 0 0 0 0 707960 0 0 0 0 0
115 73 0 0 0 0 3214 0 0 0 0 0
116 74 0 0 0 0 0 0 0 0 0 0
117 75 0 0 0 0 0 0 0 0 0 0
118 76 0 0 0 0 0 0 0 0 0 0
119 77 0 0 0 0 0 0 0 0 0 0
120 78 0 0 0 0 0 0 0 0 0 0
121 79 0 0 0 0 0 0 0 0 0 0
122 7a 0 0 0 0 0 0 0 0 0 0
123 7b 0 0 0 0 0 0 0 0 0 0
124 7c 0 0 0 0 0 0 0 0 0 0
125 7d 0 0 0 0 0 0 0 0 0 0
126 7e 0 0 0 0 0 0 0 0 0 0
127 7f 0 0 0 0 0 0 0 0 0 0
128 80 0 0 0 0 0 0 0 0 0 0
129 81 0 0 0 0 0 32606 0 0 0 0
130 82 0 0 0 0 0 2112854 0 0 0 0
131 83 0 0 0 0 0 18278976 0 0 0 0
132 84 0 0 0 0 0 48829888 0 0 0 0
133 85 0 0 0 0 0 64996010 0 0 0 0
134 86 0 0 0 0 0 67076258 0 0 0 0
135 87 0 0 0 0 0 67107064 0 0 0 0
136 88 0 0 0 0 0 66556472 0 0 0 0
137 89 0 0 0 0 0 57957875 0 0 0 0
138 8a 0 0 0 0 0 30327768 0 0 0 0
139 8b 0 0 0 0 0 5975796 0 0 0 0
140 8c 0 0 0 0 0 244952 0 0 0 0
141 8d 0 0 0 0 0 211 0 0 0 0
142 8e 0 0 0 0 0 0 0 0 0 0
143 8f 0 0 0 0 0 0 0 0 0 0
144 90 0 0 0 0 0 0 0 0 0 0
145 91 0 0 0 0 0 0 0 0 0 0
146 92 0 0 0 0 0 0 0 0 0 0
147 93 0 0 0 0 0 0 0 0 0 0
148 94 0 0 0 0 0 0 0 0 0 0
149 95 0 0 0 0 0 0 0 0 0 0
150 96 0 0 0 0 0 0 0 0 0 0
151 97 0 0 0 0 0 0 0 0 0 0
152 98 0 0 0 0 0 0 0 0 0 0
153 99 0 0 0 0 0 0 0 0 0 0
154 9a 0 0 0 0 0 0 952 0 0 0
155 9b 0 0 0 0 0 0 426280 0 0 0
156 9c 0 0 0 0 0 0 7989851 0 0 0
157 9d 0 0 0 0 0 0 34632208 0 0 0
158 9e 0 0 0 0 0 0 60176580 0 0 0
159 9f 0 0 0 0 0 0 66783720 0 0 0
160 a0 0 0 0 0 0 0 67108394 0 0 0
161 a1 0 0 0 0 0 0 67040424 0 0 0
162 a2 0 0 0 0 0 0 64029864 0 0 0
163 a3 0 0 0 0 0 0 45052052 0 0 0
164 a4 0 0 0 0 0 0 14839700 0 0 0
165 a5 0 0 0 0 0 0 1402528 0 0 0
166 a6 0 0 0 0 0 0 14176 0 0 0
167 a7 0 0 0 0 0 0 0 0 0 0
168 a8 0 0 0 0 0 0 0 0 0 0
169 a9 0 0 0 0 0 0 0 0 0 0
170 aa 0 0 0 0 0 0 0 0 0 0
171 ab 0 0 0 0 0 0 0 0 0 0
172 ac 0 0 0 0 0 0 0 0 0 0
173 ad 0 0 0 0 0 0 0 0 0 0
174 ae 0 0 0 0 0 0 0 0 0 0
175 af 0 0 0 0 0 0 0 0 0 0
176 b0 0 0 0 0 0 0 0 0 0 0
177 b1 0 0 0 0 0 0 0 0 0 0
178 b2 0 0 0 0 0 0 0 0 0 0
179 b3 0 0 0 0 0 0 0 0 0 0
180 b4 0 0 0 0 0 0 0 8960 0 0
181 b5 0 0 0 0 0 0 0 1127526 0 0
182 b6 0 0 0 0 0 0 0 13263308 0 0
183 b7 0 0 0 0 0 0 0 43057520 0 0
184 b8 0 0 0 0 0 0 0 63435032 0 0
185 b9 0 0 0 0 0 0 0 67012480 0 0
186 ba 0 0 0 0 0 0 0 67108646 0 0
187 bb 0 0 0 0 0 0 0 66863912 0 0
188 bc 0 0 0 0 0 0 0 61133068 0 0
189 bd 0 0 0 0 0 0 0 36781096 0 0
190 be 0 0 0 0 0 0 0 9150989 0 0
191 bf 0 0 0 0 0 0 0 552392 0 0
192 c0 0 0 0 0 0 0 0 1800 0 0
193 c1 0 0 0 0 0 0 0 0 0 0
194 c2 0 0 0 0 0 0 0 0 0 0
195 c3 0 0 0 0 0 0 0 0 0 0
196 c4 0 0 0 0 0 0 0 0 0 0
197 c5 0 0 0 0 0 0 0 0 0 0
198 c6 0 0 0 0 0 0 0 0 0 0
199 c7 0 0 0 0 0 0 0 0 0 0
200 c8 0 0 0 0 0 0 0 0 0 0
201 c9 0 0 0 0 0 0 0 0 0 0
202 ca 0 0 0 0 0 0 0 0 0 0
203 cb 0 0 0 0 0 0 0 0 0 0
204 cc 0 0 0 0 0 0 0 0 0 0
205 cd 0 0 0 0 0 0 0 0 84 0
206 ce 0 0 0 0 0 0 0 0 182120 0
207 cf 0 0 0 0 0 0 0 0 5116903 0
208 d0 0 0 0 0 0 0 0 0 28198760 0
209 d1 0 0 0 0 0 0 0 0 56691792 0
210 d2 0 0 0 0 0 0 0 0 66400904 0
211 d3 0 0 0 0 0 0 0 0 67105650 0
212 d4 0 0 0 0 0 0 0 0 67087086 0
213 d5 0 0 0 0 0 0 0 0 65379968 0
214 d6 0 0 0 0 0 0 0 0 50595456 0
215 d7 0 0 0 0 0 0 0 0 20129524 0
216 d8 0 0 0 0 0 0 0 0 2560762 0
217 d9 0 0 0 0 0 0 0 0 47720 0
218 da 0 0 0 0 0 0 0 0 0 0
219 db 0 0 0 0 0 0 0 0 0 0
220 dc 0 0 0 0 0 0 0 0 0 0
221 dd 0 0 0 0 0 0 0 0 0 0
222 de 0 0 0 0 0 0 0 0 0 0
223 df 0 0 0 0 0 0 0 0 0 0
224 e0 0 0 0 0 0 0 0 0 0 0
225 e1 0 0 0 0 0 0 0 0 0 0
226 e2 0 0 0 0 0 0 0 0 0 0
227 e3 0 0 0 0 0 0 0 0 0 0
228 e4 0 0 0 0 0 0 0 0 0 0
229 e5 0 0 0 0 0 0 0 0 0 0
230 e6 0 0 0 0 0 0 0 0 0 0
231 e7 0 0 0 0 0 0 0 0 0 1800
232 e8 0 0 0 0 0 0 0 0 0 552392
233 e9 0 0 0 0 0 0 0 0 0 9150989
234 ea 0 0 0 0 0 0 0 0 0 36781096
235 eb 0 0 0 0 0 0 0 0 0 61133068
236 ec 0 0 0 0 0 0 0 0 0 66863912
237 ed 0 0 0 0 0 0 0 0 0 67108646
238 ee 0 0 0 0 0 0 0 0 0 67012480
239 ef 0 0 0 0 0 0 0 0 0 63435032
240 f0 0 0 0 0 0 0 0 0 0 43057520
241 f1 0 0 0 0 0 0 0 0 0 13263308
242 f2 0 0 0 0 0 0 0 0 0 1127526
243 f3 0 0 0 0 0 0 0 0 0 8960
244 f4 0 0 0 0 0 0 0 0 0 0
245 f5 0 0 0 0 0 0 0 0 0 0
246 f6 0 0 0 0 0 0 0 0 0 0
247 f7 0 0 0 0 0 0 0 0 0 0
248 f8 0 0 0 0 0 0 0 0 0 0
249 f9 0 0 0 0 0 0 0 0 0 0
250 fa 0 0 0 0 0 0 0 0 0 0
251 fb 0 0 0 0 0 0 0 0 0 0
252 fc 0 0 0 0 0 0 0 0 0 0
253 fd 0 0 0 0 0 0 0 0 0 0
254 fe 0 0 0 0 0 0 0 0 0 0
255 ff 0 0 0 0 0 0 0 0 0 0
この実行結果の行が「最後の右シフトで捨てられる情報」、列が10で割った余りに対応している。
この実行結果を見ると、入力が32ビット符号なし整数のとき、都合よく「最後の右シフトで捨てられる情報」の上位4ビットだけを見れば10で割った余りを求めることができそうなことがわかる。
というか、この部分は入力の値に 0.1 を掛けたときの小数部分に相当するので、10を掛けて同様に右シフトをすることで余りが得られそうだ。
この処理を組み込んでみよう。
#include <stdio.h>
#include <limits.h>
unsigned int div10(unsigned int *rem, unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
unsigned long long value =
((unsigned long long)a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7)
- (a3 >> 11) - (a3 >> 15) - (a3 >> 19) - (a3 >> 23) - (a3 >> 27);
*rem = ((value & 0xff) * 10) >> 8;
return value >> 8;
}
int main(void) {
for (unsigned int i = 0; ; i++) {
unsigned int res, rem;
res = div10(&rem, i);
if (res != i / 10 || rem != i % 10) {
printf("%u -> %u, %u\n", i, res, rem);
break;
}
if (i == UINT_MAX) break;
}
return 0;
}
結果は何も出力されず、任意の32ビット符号なし整数について正しい10で割った商と余りが求まることがわかった。
1個だけ掛け算を使っているのはカッコ悪いので、シフトと足し算に変換しておこう。
unsigned int div10(unsigned int *rem, unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
unsigned long long value =
((unsigned long long)a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7)
- (a3 >> 11) - (a3 >> 15) - (a3 >> 19) - (a3 >> 23) - (a3 >> 27);
unsigned int rem_work = value & 0xff;
*rem = (rem_work + (rem_work << 2)) >> 7;
return value >> 8;
}
結論
以下の関数で、任意の32ビット符号なし整数 (少なくとも UINT_MAX 以下) を10で割った商と余りを求めることができる。
unsigned int div10(unsigned int *rem, unsigned int a) {
unsigned long long a3 = ((unsigned long long)a << 1) + a;
unsigned long long value =
((unsigned long long)a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7)
- (a3 >> 11) - (a3 >> 15) - (a3 >> 19) - (a3 >> 23) - (a3 >> 27);
unsigned int rem_work = value & 0xff;
*rem = (rem_work + (rem_work << 2)) >> 7;
return value >> 8;
}
以下の関数で、任意の24ビット符号なし整数 (54,351,188 以下) を10で割った商と余りを求めることができる。
unsigned int div10(unsigned int *rem, unsigned int a) {
unsigned int a3 = (a << 1) + a;
unsigned int value =
(a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7)
- (a3 >> 11) - (a3 >> 15) - (a3 >> 19);
unsigned int rem_work = value & 0xff;
*rem = (rem_work + (rem_work << 2)) >> 7;
return value >> 8;
}
以下の関数で、任意の16ビット符号なし整数 (232,788 以下) を10で割った商と余りを求めることができる。
unsigned int div10(unsigned int *rem, unsigned int a) {
unsigned int a3 = (a << 1) + a;
unsigned int value = (a << 5) - (a3 << 1) - (a3 >> 3) - (a3 >> 7) - (a3 >> 11);
unsigned int rem_work = value & 0xff;
*rem = (rem_work + (rem_work << 2)) >> 7;
return value >> 8;
}
おまけ：GCC による「10で割る」の表現
Compiler Explorer を用いて、GCC が10 で割った商と余りをどう表現するかを見てみた。
以下のコードを入力とし、-O2 オプションをつけてコンパイルした。
unsigned int div10(unsigned int a) {
return a / 10;
}
unsigned int mod10(unsigned int a) {
return a % 10;
}
x86-64 gcc 14.2
div10:
movl %edi, %eax
movl $3435973837, %edx
imulq %rdx, %rax
shrq $35, %rax
ret
mod10:
movl %edi, %eax
movl $3435973837, %edx
imulq %rdx, %rax
shrq $35, %rax
leal (%rax,%rax,4), %edx
movl %edi, %eax
addl %edx, %edx
subl %edx, %eax
ret
商は入力にマジックナンバーを掛けて35ビット右シフトすることで、余りは入力から商の10倍を引くことで求めている。
マジックナンバーの 3435973837 を二進数で表すと 11001100110011001100110011001101 となり、パターンに沿っているらしいことがわかる。
ARM GCC 14.2.0
div10:
movw r3, #52429
movt r3, 52428
umull r3, r0, r3, r0
lsrs r0, r0, #3
bx lr
mod10:
movw r3, #52429
movt r3, 52428
movs r2, #10
umull r1, r3, r3, r0
lsrs r3, r3, #3
mls r0, r2, r3, r0
bx lr
「マジックナンバーを掛けて35ビット右シフトする」「余りは入力から商の10倍を引いて求める」方針のようである。
各命令の意味は、以下のものである。
| 命令 | 意味 |
|---|---|
movw a, b |
a ← b |
movt a, b |
a ← (a & 0xffff) | (b << 16) |
umull a, b, c, d |
b:a ← c * d |
lsrs a, b, c |
a ← b >> c (論理シフト) |
mls a, b, c, d |
a ← d - b * c |
参考：
- Assembly Programming on ARM Linux(06)
- ARM assembler for (dynabook AZ [AC100] and netwalker): armアセンブラ、掛け算命令MLS
- movw and movt in arm assembly - Stack Overflow
RISC-V (32-bits) gcc 14.2.0
div10:
li a5,-858992640
addi a5,a5,-819
mulhu a0,a0,a5
srli a0,a0,3
ret
mod10:
li a5,-858992640
addi a5,a5,-819
mulhu a5,a0,a5
srli a5,a5,3
slli a4,a5,2
add a5,a4,a5
slli a5,a5,1
sub a0,a0,a5
ret
「マジックナンバーを掛けて35ビット右シフトする」「余りは入力から商の10倍を引いて求める」方針のようである。
値 -858992640 を十六進数 (2の補数) で表すと 0xccccd000 であり、これに -819 を足すことでマジックナンバー 0xcccccccd (3435973837) を作っている。
mips gcc 14.2.0
div10:
li $2,10 # 0xa
bne $2,$0,1f
divu $0,$4,$2
break 7
mflo $2
jr $31
nop
mod10:
li $2,10 # 0xa
bne $2,$0,1f
divu $0,$4,$2
break 7
mfhi $2
jr $31
nop
除算命令を用いていることがわかる。
AVR GCC 14.2.0
__SP_H__ = 0x3e
__SP_L__ = 0x3d
__SREG__ = 0x3f
__tmp_reg__ = 0
__zero_reg__ = 1
div10:
.L__stack_usage = 0
ldi r22,lo8(10)
ldi r23,0
rcall __udivmodhi4
mov r24,r22
mov r25,r23
ret
mod10:
.L__stack_usage = 0
ldi r22,lo8(10)
ldi r23,0
rcall __udivmodhi4
ret
この環境の unsigned int は2バイトのようである。
肝心の処理内容は、ライブラリに隠蔽されている。