ulong（64bit符号なし整数型）同士の掛け算を128bitで計算して上位64bitを取得する必要に迫られました。C#のサンプルが見当たらなかったので、作成したものをパブリックドメインに置いておきます。下位64bitは普通に掛け算すれば得られます。
public ulong umulh(ulong x, ulong y) { ulong xh = x >> 32, xl = x & 0xffffffff; ulong yh = y >> 32, yl = y & 0xffffffff; ulong a = xh * yl, ah = a >> 32, al = a & 0xffffffff; ulong b = xl * yh, bh = b >> 32, bl = b & 0xffffffff; return ((((xl * yl) >> 32) + al + bl) >> 32) + ah + bh + xh * yh; }
参考までに導出過程を書いておきます。
10進数2桁×10進数2桁
取っ掛かりとして小学校で習う10進数での2桁×2桁の計算を考えます。
1文字1桁として変数を用意します。
ab * cd = wxyz 例: 67 * 89 = 5963 → a = 6, b = 7, c = 8, d = 9
int a = ab / 10, b = ab % 10, c = cd / 10, d = cd % 10;
筆算で必要となる組み合わせを調べます。
a b (例: 67)
* c d (例: 89)
---------------
[a*d][b*d]
[a*c][b*c]列挙した積を1文字1桁の形式にします。
b * d = ef (例: 7 * 9 = 63) a * d = gh (例: 6 * 9 = 54) b * c = ij (例: 7 * 8 = 56) a * c = kl (例: 6 * 8 = 48)
int ef = b * d, gh = a * d, ij = b * c, kl = a * c; int e = ef / 10, f = ef % 10; int g = gh / 10, h = gh % 10; int i = ij / 10, j = ij % 10; int k = kl / 10, l = kl % 10;
筆算で足す桁を調べます。
ef (例: 63) gh (例: 54 ) ij (例: 56 ) + kl (例: 48 ) ------ wxyz (例: 5963)
2桁ごとに区切ってwxとyzを求めます。wxはyの繰り上がりから影響を受けるため、最初にyを求めます。
e + h + j = my (例: 6 + 4 + 6 = 16) m + g + i + kl = wx (例: 1 + 5 + 5 + 48 = 59) f + y * 10 = yz (例: 3 + 6 * 10 = 63)
int my = e + h + j; int m = my / 10, y = my % 10; int wx = m + g + i + kl; int yz = f + y * 10;
以上で10進数2桁の積が求められました。まとめると以下の通りです。
public static int[] umul(int ab, int cd) { int a = ab / 10, b = ab % 10; int c = cd / 10, d = cd % 10; int ef = b * d, gh = a * d, ij = b * c, kl = a * c; int e = ef / 10, f = ef % 10; int g = gh / 10, h = gh % 10; int i = ij / 10, j = ij % 10; int my = e + h + j; int m = my / 10, y = my % 10; return new[] { m + g + i + kl, f + y * 10 }; }
32bit2桁×32bit2桁
1文字を10進数1桁から32bit1桁に拡張します。桁ごとに分解する部分を変更します。
xy → x, y 10進数: x = xy / 10, y = xy % 10 32 bit: x = xy >> 32, y = xy & 0xffffffff
これを適用すると以下のようになります。
public static ulong[] umul(ulong ab, ulong cd) { ulong a = ab >> 32, b = ab & 0xffffffff; ulong c = cd >> 32, d = cd & 0xffffffff; ulong ef = b * d, gh = a * d, ij = b * c, kl = a * c; ulong e = ef >> 32, f = ef & 0xffffffff; ulong g = gh >> 32, h = gh & 0xffffffff; ulong i = ij >> 32, j = ij & 0xffffffff; ulong my = e + h + j; ulong m = my >> 32, y = my & 0xffffffff; return new[] { m + g + i + kl, f + (y << 32) }; }
下位は普通の掛け算で得られるため省略して、単純な計算はインライン展開して変数名を整理したものが冒頭のコードです。
