ベクトル化されたCLZおよびCTZ命令

ベクトル化されたCLZおよびCTZ 2026年10月9日 clzとctzは、固定サイズの整数における先頭（末尾）のゼロビット数を計算する命令です。これらは最新のCPUでネイティブにサポートされていますが、常に高速であるとは限りません。例えば、tzcntはArrow Lakeでレイテンシが3です。 私は現在開発中のFPUエミュレータでctzを使用していましたが、浮動小数点演算のトリックで回避できることに気づき、それが間接的にベクトル化可能なctzのポリフィルに一般化できることを発見しました。完全性のためにclzも実装しました。 clz まず、2つのうち簡単な方であるclzから始めましょう。 fn clz(x: u32) -> u32 { let a = 2.0f64.powi(-970); 32 - ((f64::from_bits(a.to_bits() | x as u64) - a).to_bits() >> 52) as u32 } 一般的な考え方は以下の通りです。 浮動小数点数の指数は、その値のバイアスされた対数です。32ビットの数値xを値2kのマントッサに代入すると、2k(1+2−52x)を表すdoubleが得られます。次に、2kをdoubleとして減算すると、2k−52xが得られます。指数を抽出すると、k−52+31−clz(x)が得られ、ここからビット単位の減算でclzを計算できます。そこから、x=0でclzが正しく動作するようにkを選択できます。 32ビットの入力値を扱うには64ビットのdoubleが必要です。残念ながら、これはこのトリックが任意の64ビット入力値には適用できず、52ビットまでしか機能しないことを意味します。 入力と出力がu64x4に格納されていると仮定すると、これは以下のようにコンパイルされます。 vpbroadcastq ymm1, [rip + bias] vorpd ymm0, ymm0, [rip + a] vsubpd ymm0, ymm0, [rip + a] vpsrlq ymm0, ymm0, 52 vpsubq ymm0, ymm1, ymm0 a: .quad 0x350000000000000, 0x350000000000000, 0x350000000000000, 0x350000000000000 bias: .quad 32 私のHaswellでは、これは1イテレーションあたり0.45ナノ秒で実行されます。スカラーバージョンでは1ナノ秒です。レイテンシバウンドの場合、数値は2ナノ秒対1ナノ秒に上昇します（ただし、ベクトル化されたctzでレイテンシバウンドになっている場合は、おそらく何か間違っています）。 Ian Qvist氏がAlder Lakeでこれをテストしてくれましたが（ありがとう！）、1イテレーションあたり0.29ナノ秒、スカラーバージョンでは0.85ナノ秒、レイテンシバウンドでは1.3ナノ秒対0.85ナノ秒でした。最新のIntel CPUでは、数値は同じかそれ以上になるはずです。 AMD CPUではlzcntが非常に安価であるため、スカラーバージョンが有利になる可能性が高いです。ただし、Zen CPUはAVX-512をサポートしており、vplzcntdがあるため、それも選択肢となります。 ctz 次にctzについてです。 fn ctz(x: u32) -> u32 { let a = f64::from_bits((0x340000100000001 ^ x as u64) ^ (x as u64 + u32::MAX as u64)) - f64::from_bits(0x340000000000000); (a.to_bits() >> 52) as u32 } まずx⊕︎(x−1)から始めて、最も低いセットビットを分離します。ctzは、その値の対数に等しく、ビット単位で2kを加え、doubleとして2kを減算し、指数を調べます。適切に選択されたkにより、バイアスされていないctzが含まれます。x=0を正しく処理するために、マントッサに232を事前に混合し、x−1の代わりにx+232−1を使用します。また、奇数xが遅いサブノーマルではなくa=0を生成するように、マントッサに1を事前に混合します。（どれだけの時間を費やしたか想像できますか？） この関数は以下のようにコンパイルされます。 vpxor ymm1, ymm0, [rip + c1] vpaddq ymm0, ymm0, [rip + u32_max] vpxor ymm0, ymm1, ymm0 vsubpd ymm0, ymm0, [rip + c2] vpsrlq ymm0, ymm0, 52 c1: .qword 0x340000100000001, 0x340000100000001, 0x340000100000001, 0x340000100000001 u32_max: .qword 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff c2: .qword 0x340000000000000, 0x340000000000000, 0x340000000000000, 0x340000000000000 Haswellでは、これは1イテレーションあたり0.49ナノ秒、レイテンシバウンドでは2.3ナノ秒で実行されます。Alder Lakeでは、1イテレーションあたり0.35ナノ秒、レイテンシバウンドでは1.3ナノ秒です。clzと比較して遅延が増加するのは、命令が1つ多く使用されているためです。AVX-512が存在する場合は、vpternlogqを使用することで回避できますが、その場合は(x - 1) & !xに対してvpopcntdを実行するのと同じくらい簡単です。スカラーバージョンはclzと違いはありません。 10月10日追加： Nikolay Malkovsky氏が、デ・ブルイン系列が別のベクトル化可能なアプローチを提供すると指摘しました。いくつかのテストの後、私は以下のコードにたどり着きました。 const char table[32] = { 0, 4, 5, 6, 11, 9, 7, 12, 15, 3, 10, 8, 14, 2, 13, 1, 0, 4, 5, 6, 11, 9, 7, 12, 15, 3, 10, 8, 14, 2, 13, 1, }; __m256i bit = _mm256_andnot_si256(x, _mm256_sub_epi32(x, _mm256_set1_epi32(1))); __m256i high = _mm256_madd_epi16( _mm256_cmpeq_epi16(bit, _mm256_set1_epi16(-1)), _mm256_set1_epi16(-16) ); __m256i index = _mm256_srli_epi32(_mm256_mullo_epi32(bit, _mm256_set1_epi32(0xf0a6f0a7)), 28); __m256i low = _mm256_shuffle_epi8(_mm256_loadu_si256((__m256i*)table), index); return _mm256_add_epi32(low, high); vpshufbが16バイトのレーンを越えることができないため、真の32バイトLUTを使用することはできません。代わりに私が使用したアプローチは説明が難しいですが、基本的に16ビットのデ・ブルイン系列を2回繰り返してctzのビット0-3を計算し、ゼロの半分に応じて16または32を加算します。0xf0a6f0a7の67は、これが機能する4つのマジック定数のうちの1つです。 これはHaswellで1ナノ秒（Alder Lakeでは0.7ナノ秒）かかりますが、スループットが2倍になるため、シャッフルを回避するのに役立つ場合はFPベースのアプローチよりもわずかに高速になる可能性があります。 x=0を扱う必要がない場合（またはctz(0)が32ではなく0になることを望む場合）、 __m256i high = _mm256_and_si256( _mm256_cmpgt_epi32(bit, _mm256_set1_epi32(0x7fff)), _mm256_set1_epi32(16) ); を使用すると、時間を0.82ナノ秒に短縮できます。 10月11日追加： arya dradjica氏と私は、アルゴリズムの64ビットバージョンについてブレインストーミングを行い、私が最適だと考えるアルゴリズムを見つけることができました。 const char table[32] = { // x -> 4 - popcnt(x) 4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0, 4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0, }; __m256i lut = _mm256_loadu_si256((__m256i*)table); __m256i bit = _mm256_andnot_si256(x, _mm256_sub_epi64(x, _mm256_set1_epi64x(1))); __m256i low = _mm256_shuffle_epi8(lut, bit); __m256i high = _mm256_shuffle_epi8(lut, _mm256_srli_epi16(bit, 4)); return _mm256_sad_epu8(_mm256_add_epi8(low, high), _mm256_set1_epi8(8)); これはctz(0)=64を正しく返し、私のHaswellでは0.63ナノ秒（レイテンシバウンドでは2.8ナノ秒）かかります。これは、(x - 1) & !xの各バイトの2つのニブルの反転されたポップカウントを個別に計算し（pshufbのいくつかの特殊性のために反転されています）、次にvpsadbwを使用して水平合計を実行し、一度に反転を修正します。 vpsadbwは正確に8バイトの水平合計しかできないため、このアルゴリズムは64ビットに特化しており、32ビット入力では機能しません。デ・ブルイン系列に基づいたアプローチは、そのシナリオでは依然として最速です。