cp-library

競技プログラミング用 C++ ライブラリ

View the Project on GitHub qandori/cp-library

:heavy_check_mark: 区間加算・最小値のノード演算

components/range_add_min_node.hpp · GitHub で見る

共通ページ

区間加算・区間最小値の 8 子ノードを処理する内部部品です。通常は data_structure/range_add_range_min.hpp を include します。

range_add_min_detail::Node<T> の T は 64 bit 以下の標準の符号付き整数型、U は同じ幅の符号なし型です。通常は、配列にノードの最小値からの差分を U で格納します。正規化済みの有効要素の最小値は 0 です。差分は U の全範囲を取り得るため、通常の経路は unsigned の順序を使います。build への入力は、差分ではなく元の符号付き値のビット表現です。

通常の values は、初期化済みの連続 8 要素の U 配列です。有効要素は先頭 count 個で、1 <= count <= 8 とします。更新には 8 要素とも読み書き可能な配列を渡します。build の count 0 と materialize の領域条件は後述します。

64 bit の SIMD 経路では、build・materialize を除き、配列の先頭アドレスを 32 バイト境界に揃えます。build・materialize と 32 bit 以下の操作には SIMD 用の追加のアラインメント条件はなく、T・U 本来のアラインメントを満たせば使えます。

操作 内容・返り値 計算量
build(values,count) U* の先頭 unsigned count 個を元の T のビット表現として解釈し、符号付き最小値を T で返す。有効要素を最小値からの unsigned 差分へ置き換え、残りは変更しない \(O(1)\)
minimum<Narrow>(values,l,r) const U* と unsigned l,r を渡し、0 <= l <= r <= count の最小差分を U で返す。空区間の返り値は本文参照 \(O(1)\)
normalize<Narrow>(values,count,negative,first,last,delta,x=0,dx=0,y=0,dy=0) U* の [first,last) に U delta、添字 x,y に U dx,dy を加え、有効要素を再び最小差分 0 に揃える。元の最小値に加える変化量を U で返す \(O(1)\)
normalize_boundaries<Narrow>(left,right,right_count,negative,first,last,delta,left_change,right_change) 左の添字 first に U left_change、[first+1,8) に U delta を加える。右の [0,last) に delta、添字 last に U right_change を加える。両ブロックを正規化し、左右の最小値の変化量を std::pair<U,U> で返す \(O(1)\)
change_one<Narrow>(values,count,index,delta,negative) U* の添字 index に U delta を加えて正規化する。元の最小値に加える変化量を U で返す。delta = 0 なら何もしない \(O(1)\)
full_minimum<Narrow>(values) const U* の 8 要素すべての最小差分を U で返す。パディングも参照する \(O(1)\)
subtract(values,count,change) U* の先頭 count 要素から U change を引く。パディングは変更しない。返り値は void \(O(1)\)
materialize(src,dst,count,base) const U* src の先頭 unsigned count 個に U base を巡回加算し、そのビット表現を T* dst へ出力。返り値は void \(O(1)\)
min_signed(x,y) AVX2 用の補助関数。__m256i の各 signed 64 bit レーンの最小値を __m256i で返す \(O(1)\)
broadcast_min(value) AVX2 用の補助関数。__m256i の 4 レーンの signed 64 bit 最小値を、全レーンに複製した __m256i で返す \(O(1)\)
reduce_min(x) AVX2 用の補助関数。__m256i の 4 レーンを signed 64 bit として比較し、最小値のビット表現を U で返す \(O(1)\)
broadcast_min32(value) AVX2 用の補助関数。__m256i の 8 レーンの unsigned 32 bit 最小値を、全レーンに複製した __m256i で返す \(O(1)\)
reduce_min32(value) AVX2 用の補助関数。__m256i の 8 レーンを unsigned 32 bit として比較し、最小値を U で返す \(O(1)\)
minimum_vector<Narrow>(values,l,r,neutral) AVX2 用の補助関数。const U* の [l,r) 以外を U neutral で置き換え、前半と後半の対応レーンを比較した 4 レーンの __m256i を返す。変換と中立値の条件は本文参照 \(O(1)\)

Narrow は省略できない bool のテンプレート引数で、true は 64 bit 整数用です。32 bit 以下では false を使います。添字と count,right_count,first,last,x,y,index は unsigned、negative は bool です。normalize は 0 <= first <= last <= count、change_one は 0 <= index < count とします。x,y は dx,dy が 0 でも常に 8 未満にし、非ゼロの差分を渡す添字は count 未満にします。normalize では範囲と一点の重複や x == y も指定できます。重なる位置には変化量を合計して加え、更新後の実際の値と Narrow の条件を満たすようにします。空の [first,last) も指定でき、dx,dy のみを適用できます。表に示したもの以外に既定値はありません。

build は 0 <= count <= 8 を指定します。先頭 count 個を std::bit_cast<T> で符号付き値として解釈し、その最小値を求めます。有効要素から最小値のビット表現を unsigned 演算で引き、差分へ正規化します。T の最小値と最大値を同時に含められ、Narrow の値域条件は不要です。

build は count が正なら、無効要素を含む初期化済みの読み書き可能な 8 要素領域を必要とします。無効要素は符号付き最小値の比較から除外し、元の値を保持します。通常は無効要素を std::numeric_limits<U>::max() で埋めておきます。これは raw の T では負の値ですが、最小値として選ばれません。count 0 は配列を参照せず std::numeric_limits<T>::max() を返し、values=nullptr も指定できます。

materialize も 0 <= count <= 8 を指定します。各出力は std::bit_cast<T>(U(src[i] + base)) です。src と base は任意の U の値を渡せ、加算は U の幅で巡回します。入力は変更せず、出力は先頭 count 個だけを書き換えます。count が正なら読み取り可能な src の count 要素と書き込み可能な dst の count 要素を用意し、両領域を重ねないでください。count 0 は両方のポインタに nullptr を指定できます。count が 8 未満なら、8 要素分の領域を確保する必要はありません。

normalize_boundaries の left,right は重ならない別の 8 要素ブロックを指す U* です。左の有効要素数は 8、右は 1 <= right_count <= 8 とし、0 <= first < 8、0 <= last < right_count を満たします。左右をそれぞれ正規化しておき、返された変化量をそれぞれの親へ反映します。左右の境界位置には共通の delta を重ねず、left_change,right_change だけを加えます。

差分への加減算は U の巡回演算で行います。負の変化量も U に変換して渡し、実際の増減方向を negative で別に指定します。1 回の呼び出しで渡す非ゼロの変化量はすべて同じ増減方向にし、減少なら negative=true、増加なら false とします。normalize_boundaries の左右の変化量にも同じ条件が必要です。差分の最上位ビットから増減方向を決めないでください。呼び出し元は、各子が表す更新後の値が T の範囲に収まることを保証します。更新前の有効要素は正規化しておき、返された変化量を親ノードの最小値に加えます。

Narrow=false では、パディングを std::numeric_limits<U>::max() にします。負更新で差分が wrap する場合も扱い、U 全域の差分を比較します。

32 bit の AVX2 経路は unsigned 32 bit の最小値を使い、符号付き整数の全範囲に対応します。negative=true では、更新後の unsigned 差分が元より大きくなるレーンを、0 を下回った差分として区別します。このレーンが存在すればその中の最小値を選び、なければ全有効レーンの最小値を選びます。normalize_boundaries では左右それぞれでこの判定を行ってから、両ブロックの最小値を同時に求めます。右のパディングには正規化量を引かず、U の最大値を維持します。

minimum の空区間は、32 bit 以下では std::numeric_limits<U>::max() を返します。64 bit も Narrow=false なら同じです。64 bit の Narrow=true では、SIMD 経路は std::numeric_limits<T>::max()、スカラー経路は std::numeric_limits<U>::max() を U で返します。空区間から通常の差分が得られるとみなして加算しないでください。

min_signed・broadcast_min・reduce_min・minimum_vector は AVX2 有効時だけ定義される、64 bit 整数用の補助関数です。broadcast_min32・reduce_min32 は同じ条件で定義される unsigned 32 bit 用です。

minimum_vector は 0 <= l <= r <= 8 とし、指定区間が有効要素だけを含むようにします。配列には前述の 64 bit 用のアラインメント条件が必要です。レーン j は添字 j と j+4 の小さい方を表します。l == r なら全レーンが neutral になります。レーン間の最小値はまだ取らず、reduce_min または broadcast_min と組み合わせます。

minimum_vector<true> は元の差分を signed 64 bit として比較します。minimum_vector<false> は読み込んだ各差分 value を key = value ^ sign_bit に変換し、key を signed 64 bit として比較します。これにより元の unsigned 順序と一致します。neutral はこの変換後の表現で渡し、signed の順序で全候補以上になる値を選びます。通常経路の中立値は U(std::numeric_limits<T>::max()) です。false の返り値に reduce_min を適用した後は、さらに sign_bit を XOR して元の差分に戻します。配列に保存する差分自体は変換しません。

Narrow=true は 64 bit 整数の内部用経路です。M = std::numeric_limits<T>::max()、初期列の最大値と最小値の差を \(S_0\)、全区間加算を除く更新の絶対変化量累計を \(A\) とし、\(\displaystyle L=\left\lfloor\frac{M}{3}\right\rfloor\) とします。一点代入では新旧の値の差の絶対値を累計に加えます。呼び出し元が \(S_0+A<L\) を保証し、構築後のパディングを \(2L\) にします。SIMD の normalize<true> はパディングも含めて signed の最小値を計算し、8 レーンすべてからその値を引きます。

有効要素の正規化後の差分は \([0,S_0+A]\) に収まります。各ノードの最小値の変化量の絶対値は、その更新で加えた変化量の絶対値以下です。したがって、パディングから差し引かれる正規化量の累計の絶対値も \(A\) 以下で、パディングの差分は \([2L-A,2L+A]\) に収まります。\(A<L\) なのでパディングは \(L\) より大きく、有効要素は \(L\) 未満です。一時的な負の差分も \(-L\) より大きく、パディングを含むすべての比較対象が signed 64 bit の範囲内です。これにより、パディングを含む signed の最小値計算と全レーンの正規化が成立します。

範囲条件を満たさなくなる更新の前に、各段の末尾のパディングを std::numeric_limits<U>::max() に戻し、通常の経路に切り替えます。この管理は RangeAddRangeMin が行い、利用者の値域制限は増えません。

子の数は 8 に固定されているため、各操作は \(O(1)\) 時間・追加空間です。

検証コード

コード

#ifndef RANGE_ADD_MIN_NODE_HPP
#define RANGE_ADD_MIN_NODE_HPP

#include <algorithm>
#include <array>
#include <bit>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <type_traits>
#include <utility>
#if (defined(__AVX2__) || defined(__AVX512F__)) && defined(__x86_64__)
#include <immintrin.h>
#endif

namespace range_add_min_detail {
template<class T> struct Node {
	static_assert(std::is_integral_v<T> && std::is_signed_v<T> && sizeof(T) <= 8);
	using U = std::make_unsigned_t<T>;
	static constexpr unsigned branch = 8;
	static constexpr U umax = std::numeric_limits<U>::max();
	static constexpr U sign_bit = U(1) << (sizeof(U) * 8 - 1);
#if defined(__AVX2__) && defined(__x86_64__)
	alignas(64) inline static constexpr auto masks_ = [] {
		std::array<std::array<U, branch>, branch + 1> masks{};
		for (unsigned i = 0; i <= branch; ++i)
			for (unsigned j = 0; j < branch; ++j) masks[i][j] = j < i ? umax : 0;
		return masks;
	}();
	alignas(64) inline static constexpr auto slots_ = [] {
		std::array<std::array<U, branch>, branch> slots{};
		for (unsigned i = 0; i < branch; ++i) slots[i][i] = umax;
		return slots;
	}();
	static __m256i min_signed(__m256i x, __m256i y) { return _mm256_blendv_epi8(x, y, _mm256_cmpgt_epi64(x, y)); }
	static __m256i broadcast_min(__m256i value) {
		value = min_signed(value, _mm256_permute2x128_si256(value, value, 1));
		return min_signed(value, _mm256_permute4x64_epi64(value, 0xb1));
	}
	static U reduce_min(__m256i value) { return U(_mm_cvtsi128_si64(_mm256_castsi256_si128(broadcast_min(value)))); }
	static __m256i broadcast_min32(__m256i value) {
		value = _mm256_min_epu32(value, _mm256_permute2x128_si256(value, value, 1));
		value = _mm256_min_epu32(value, _mm256_shuffle_epi32(value, 0x4e));
		return _mm256_min_epu32(value, _mm256_shuffle_epi32(value, 0xb1));
	}
	static U reduce_min32(__m256i value) { return U(std::uint32_t(_mm_cvtsi128_si32(_mm256_castsi256_si128(broadcast_min32(value))))); }
	template<bool Narrow>
	static __m256i minimum_vector(const U* values, unsigned l, unsigned r, U neutral) {
		auto identity = _mm256_set1_epi64x(std::bit_cast<T>(neutral));
		__m256i result;
		for (unsigned i = 0; i < branch; i += 4) {
			auto first = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[l].data() + i));
			auto last = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[r].data() + i));
			auto value = _mm256_load_si256(reinterpret_cast<const __m256i*>(values + i));
			if constexpr (!Narrow) value = _mm256_xor_si256(value, _mm256_set1_epi64x(std::bit_cast<T>(sign_bit)));
			value = _mm256_blendv_epi8(identity, value, _mm256_xor_si256(first, last));
			result = i ? min_signed(result, value) : value;
		}
		return result;
	}
#endif
	static T build(U* values, unsigned count) {
		if (!count) return std::numeric_limits<T>::max();
#if defined(__AVX512F__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto value = _mm512_loadu_si512(values);
			unsigned valid = (1u << count) - 1;
			T minimum = T(_mm512_mask_reduce_min_epi64(valid, value));
			value = _mm512_mask_sub_epi64(value, valid, value, _mm512_set1_epi64(minimum));
			_mm512_storeu_si512(values, value);
			return minimum;
		}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto first = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data()));
			auto last = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data() + 4));
			auto x = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values));
			auto y = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values + 4));
			auto identity = _mm256_set1_epi64x(std::numeric_limits<T>::max());
			auto a = _mm256_blendv_epi8(identity, x, first), b = _mm256_blendv_epi8(identity, y, last);
			auto minimum = broadcast_min(min_signed(a, b));
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(values), _mm256_sub_epi64(x, _mm256_and_si256(first, minimum)));
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(values + 4), _mm256_sub_epi64(y, _mm256_and_si256(last, minimum)));
			return std::bit_cast<T>(U(_mm_cvtsi128_si64(_mm256_castsi256_si128(minimum))));
		}
		if constexpr (sizeof(T) == 4) {
			auto valid = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data()));
			auto value = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values));
			auto sign = _mm256_set1_epi32(std::bit_cast<T>(sign_bit));
			auto key = _mm256_xor_si256(value, sign);
			key = _mm256_or_si256(key, _mm256_andnot_si256(valid, _mm256_set1_epi32(-1)));
			auto minimum = _mm256_xor_si256(broadcast_min32(key), sign);
			value = _mm256_sub_epi32(value, _mm256_and_si256(valid, minimum));
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(values), value);
			return std::bit_cast<T>(U(std::uint32_t(_mm_cvtsi128_si32(_mm256_castsi256_si128(minimum)))));
		}
#endif
		T minimum = std::bit_cast<T>(values[0]);
		for (unsigned i = 1; i < count; ++i) minimum = std::min(minimum, std::bit_cast<T>(values[i]));
		for (unsigned i = 0; i < count; ++i) values[i] = U(values[i] - U(minimum));
		return minimum;
	}
	template<bool Narrow>
	static U minimum(const U* values, unsigned l, unsigned r) {
#if defined(__AVX512F__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			unsigned mask = ((1u << r) - 1) ^ ((1u << l) - 1);
			auto value = _mm512_loadu_si512(values);
			if constexpr (Narrow) return U(_mm512_mask_reduce_min_epi64(mask, value));
			else return U(_mm512_mask_reduce_min_epu64(mask, value));
		}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			U result = reduce_min(minimum_vector<Narrow>(values, l, r, U(std::numeric_limits<T>::max())));
			if constexpr (Narrow) return result;
			else return result ^ sign_bit;
		}
		if constexpr (sizeof(T) == 4) {
			auto first = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[l].data()));
			auto last = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[r].data()));
			auto value = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values));
			auto excluded = _mm256_xor_si256(_mm256_xor_si256(first, last), _mm256_set1_epi32(-1));
			return reduce_min32(_mm256_or_si256(value, excluded));
		}
#endif
		U result = umax;
		for (unsigned i = l; i < r; ++i) result = std::min(result, values[i]);
		return result;
	}
	template<bool Narrow>
	static U normalize(U *values, unsigned count, bool negative, unsigned first, unsigned last, U delta, unsigned x = 0, U dx = 0, unsigned y = 0, U dy = 0) {
#if defined(__AVX512F__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto old = _mm512_loadu_si512(values);
			unsigned mask = ((1u << last) - 1) ^ ((1u << first) - 1), valid = (1u << count) - 1;
			auto next = _mm512_mask_add_epi64(old, mask, old, _mm512_set1_epi64(std::bit_cast<T>(delta)));
			next = _mm512_mask_add_epi64(next, 1u << x, next, _mm512_set1_epi64(std::bit_cast<T>(dx)));
			next = _mm512_mask_add_epi64(next, 1u << y, next, _mm512_set1_epi64(std::bit_cast<T>(dy)));
			U change;
			if constexpr (Narrow) {
				change = U(_mm512_reduce_min_epi64(next));
				next = _mm512_sub_epi64(next, _mm512_set1_epi64(std::bit_cast<T>(change)));
			} else {
				unsigned wrap = negative ? unsigned(_mm512_cmp_epu64_mask(next, old, _MM_CMPINT_GT)) & valid : 0;
				change = U(_mm512_mask_reduce_min_epu64(wrap ? wrap : valid, next));
				next = _mm512_mask_sub_epi64(next, valid, next, _mm512_set1_epi64(std::bit_cast<T>(change)));
			}
			_mm512_storeu_si512(values, next);
			return change;
		}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto identity = _mm256_set1_epi64x(std::numeric_limits<T>::max());
			__m256i next[2], result, wrapped = identity;
			auto underflow = _mm256_setzero_si256();
			for (unsigned i = 0; i < branch; i += 4) {
				auto a = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[first].data() + i));
				auto b = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[last].data() + i));
				auto old = _mm256_load_si256(reinterpret_cast<const __m256i*>(values + i));
				auto value = _mm256_add_epi64(old, _mm256_and_si256(_mm256_xor_si256(a, b), _mm256_set1_epi64x(std::bit_cast<T>(delta))));
				value = _mm256_add_epi64(value, _mm256_and_si256(_mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[x].data() + i)), _mm256_set1_epi64x(std::bit_cast<T>(dx))));
				value = _mm256_add_epi64(value, _mm256_and_si256(_mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[y].data() + i)), _mm256_set1_epi64x(std::bit_cast<T>(dy))));
				next[i / 4] = value;
				if constexpr (!Narrow) {
					auto biased = _mm256_xor_si256(value, _mm256_set1_epi64x(std::bit_cast<T>(sign_bit)));
					result = i ? min_signed(result, biased) : biased;
					if (negative) {
						auto wrap = _mm256_cmpgt_epi64(biased, _mm256_xor_si256(old, _mm256_set1_epi64x(std::bit_cast<T>(sign_bit))));
						underflow = _mm256_or_si256(underflow, wrap);
						wrapped = min_signed(wrapped, _mm256_blendv_epi8(identity, biased, wrap));
					}
				}
			}
			__m256i shift;
			U change;
			if constexpr (Narrow) {
				shift = broadcast_min(min_signed(next[0], next[1]));
				change = U(_mm_cvtsi128_si64(_mm256_castsi256_si128(shift)));
			} else {
				change = reduce_min(negative && _mm256_movemask_epi8(underflow) ? wrapped : result) ^ sign_bit;
				shift = _mm256_set1_epi64x(std::bit_cast<T>(change));
			}
			for (unsigned i = 0; i < branch; i += 4) {
				auto delta = shift;
				if constexpr (!Narrow) delta = _mm256_and_si256(delta, _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data() + i)));
				_mm256_store_si256(reinterpret_cast<__m256i*>(values + i), _mm256_sub_epi64(next[i / 4], delta));
			}
			return change;
		}
		if constexpr (sizeof(T) == 4) {
			auto first_mask = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[first].data()));
			auto last_mask = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[last].data()));
			auto old = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values));
			auto next = _mm256_add_epi32(old, _mm256_and_si256(_mm256_xor_si256(first_mask, last_mask), _mm256_set1_epi32(std::bit_cast<T>(delta))));
			next = _mm256_add_epi32(next, _mm256_and_si256(_mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[x].data())), _mm256_set1_epi32(std::bit_cast<T>(dx))));
			next = _mm256_add_epi32(next, _mm256_and_si256(_mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[y].data())), _mm256_set1_epi32(std::bit_cast<T>(dy))));
			auto candidates = next;
			if (negative) {
				auto sign = _mm256_set1_epi32(std::bit_cast<T>(sign_bit));
				auto wrapped = _mm256_cmpgt_epi32(_mm256_xor_si256(next, sign), _mm256_xor_si256(old, sign));
				if (_mm256_movemask_epi8(wrapped)) candidates = _mm256_or_si256(next, _mm256_andnot_si256(wrapped, _mm256_set1_epi32(-1)));
			}
			U change = reduce_min32(candidates);
			auto valid = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data()));
			next = _mm256_sub_epi32(next, _mm256_and_si256(valid, _mm256_set1_epi32(std::bit_cast<T>(change))));
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(values), next);
			return change;
		}
#endif
		U wrapped = umax, result = umax;
		bool underflow = false;
		for (unsigned i = 0; i < count; ++i) {
			U old = values[i], change = first <= i && i < last ? delta : 0;
			if (i == x) change += dx;
			if (i == y) change += dy;
			U next = U(old + change);
			values[i] = next;
			if (negative && next > old) { wrapped = std::min(wrapped, next); underflow = true; }
			result = std::min(result, next);
		}
		if (underflow) result = wrapped;
		for (unsigned i = 0; i < count; ++i) values[i] -= result;
		return result;
	}
	template<bool Narrow>
	static std::pair<U, U> normalize_boundaries(U* left, U* right, unsigned right_count, bool negative, unsigned first, unsigned last, U delta, U left_change, U right_change) {
#if defined(__AVX2__) && !defined(__AVX512F__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8 && Narrow) {
			__m256i x[2], y[2];
			auto common = _mm256_set1_epi64x(std::bit_cast<T>(delta));
			auto dx = _mm256_set1_epi64x(std::bit_cast<T>(left_change));
			auto dy = _mm256_set1_epi64x(std::bit_cast<T>(right_change));
			for (unsigned i = 0; i < branch; i += 4) {
				auto prefix = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[first].data() + i));
				auto suffix = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[last + 1].data() + i));
				auto lc = _mm256_andnot_si256(prefix, common);
				auto rc = _mm256_and_si256(suffix, common);
				lc = _mm256_blendv_epi8(lc, dx, _mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[first].data() + i)));
				rc = _mm256_blendv_epi8(rc, dy, _mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[last].data() + i)));
				x[i / 4] = _mm256_add_epi64(_mm256_load_si256(reinterpret_cast<const __m256i*>(left + i)), lc);
				y[i / 4] = _mm256_add_epi64(_mm256_load_si256(reinterpret_cast<const __m256i*>(right + i)), rc);
			}
			auto a = min_signed(x[0], x[1]), b = min_signed(y[0], y[1]);
			auto minimum = min_signed(_mm256_unpacklo_epi64(a, b), _mm256_unpackhi_epi64(a, b));
			minimum = min_signed(minimum, _mm256_permute2x128_si256(minimum, minimum, 1));
			auto lm = _mm256_permute4x64_epi64(minimum, 0x00), rm = _mm256_permute4x64_epi64(minimum, 0x55);
			_mm256_store_si256(reinterpret_cast<__m256i*>(left), _mm256_sub_epi64(x[0], lm));
			_mm256_store_si256(reinterpret_cast<__m256i*>(left + 4), _mm256_sub_epi64(x[1], lm));
			_mm256_store_si256(reinterpret_cast<__m256i*>(right), _mm256_sub_epi64(y[0], rm));
			_mm256_store_si256(reinterpret_cast<__m256i*>(right + 4), _mm256_sub_epi64(y[1], rm));
			return {U(_mm_cvtsi128_si64(_mm256_castsi256_si128(minimum))), U(_mm_extract_epi64(_mm256_castsi256_si128(minimum), 1))};
		}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 4) {
			auto common = _mm256_set1_epi32(std::bit_cast<T>(delta));
			auto prefix = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[first].data()));
			auto suffix = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[last + 1].data()));
			auto lc = _mm256_andnot_si256(prefix, common);
			auto rc = _mm256_and_si256(suffix, common);
			lc = _mm256_blendv_epi8(lc, _mm256_set1_epi32(std::bit_cast<T>(left_change)), _mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[first].data())));
			rc = _mm256_blendv_epi8(rc, _mm256_set1_epi32(std::bit_cast<T>(right_change)), _mm256_load_si256(reinterpret_cast<const __m256i*>(slots_[last].data())));
			auto old_left = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(left));
			auto old_right = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(right));
			auto x = _mm256_add_epi32(old_left, lc), y = _mm256_add_epi32(old_right, rc);
			auto a = x, b = y;
			if (negative) {
				auto sign = _mm256_set1_epi32(std::bit_cast<T>(sign_bit));
				auto lx = _mm256_cmpgt_epi32(_mm256_xor_si256(x, sign), _mm256_xor_si256(old_left, sign));
				auto ry = _mm256_cmpgt_epi32(_mm256_xor_si256(y, sign), _mm256_xor_si256(old_right, sign));
				if (_mm256_movemask_epi8(lx)) a = _mm256_or_si256(x, _mm256_andnot_si256(lx, _mm256_set1_epi32(-1)));
				if (_mm256_movemask_epi8(ry)) b = _mm256_or_si256(y, _mm256_andnot_si256(ry, _mm256_set1_epi32(-1)));
			}
			auto minimum = _mm256_min_epu32(_mm256_unpacklo_epi32(a, b), _mm256_unpackhi_epi32(a, b));
			minimum = _mm256_min_epu32(minimum, _mm256_shuffle_epi32(minimum, 0x4e));
			minimum = _mm256_min_epu32(minimum, _mm256_permute2x128_si256(minimum, minimum, 1));
			auto lm = _mm256_shuffle_epi32(minimum, 0x00), rm = _mm256_shuffle_epi32(minimum, 0x55);
			x = _mm256_sub_epi32(x, lm);
			y = _mm256_sub_epi32(y, _mm256_and_si256(rm, _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[right_count].data()))));
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(left), x);
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(right), y);
			auto low = _mm256_castsi256_si128(minimum);
			return {U(std::uint32_t(_mm_cvtsi128_si32(low))), U(std::uint32_t(_mm_extract_epi32(low, 1)))};
		}
#endif
		U a = first + 1 == branch ? change_one<Narrow>(left, branch, first, left_change, negative)
			: normalize<Narrow>(left, branch, negative, first + 1, branch, delta, first, left_change);
		U b = !last ? change_one<Narrow>(right, right_count, 0, right_change, negative)
			: normalize<Narrow>(right, right_count, negative, 0, last, delta, last, right_change);
		return {a, b};
	}
	template<bool Narrow>
	static U full_minimum(const U* values) {
#if defined(__AVX512F__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto value = _mm512_loadu_si512(values);
			if constexpr (Narrow) return U(_mm512_reduce_min_epi64(value));
			else return U(_mm512_reduce_min_epu64(value));
		}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto result = _mm256_load_si256(reinterpret_cast<const __m256i*>(values));
			if constexpr (!Narrow) result = _mm256_xor_si256(result, _mm256_set1_epi64x(std::bit_cast<T>(sign_bit)));
			for (unsigned i = 4; i < branch; i += 4) {
				auto value = _mm256_load_si256(reinterpret_cast<const __m256i*>(values + i));
				if constexpr (!Narrow) value = _mm256_xor_si256(value, _mm256_set1_epi64x(std::bit_cast<T>(sign_bit)));
				result = min_signed(result, value);
			}
			if constexpr (Narrow) return reduce_min(result);
			else return reduce_min(result) ^ sign_bit;
		}
		if constexpr (sizeof(T) == 4) return reduce_min32(_mm256_loadu_si256(reinterpret_cast<const __m256i*>(values)));
#endif
		U result = umax;
		for (unsigned i = 0; i < branch; ++i) result = std::min(result, values[i]);
		return result;
	}
	static void subtract(U* values, unsigned count, U change) {
#if defined(__AVX512F__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto value = _mm512_loadu_si512(values);
			value = _mm512_mask_sub_epi64(value, (1u << count) - 1, value, _mm512_set1_epi64(std::bit_cast<T>(change)));
			_mm512_storeu_si512(values, value);
			return;
		}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
		if constexpr (sizeof(T) == 8) {
			auto delta = _mm256_set1_epi64x(std::bit_cast<T>(change));
			for (unsigned i = 0; i < branch; i += 4) {
				auto valid = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data() + i));
				auto old = _mm256_load_si256(reinterpret_cast<const __m256i*>(values + i));
				_mm256_store_si256(reinterpret_cast<__m256i*>(values + i), _mm256_sub_epi64(old, _mm256_and_si256(valid, delta)));
			}
			return;
		}
		if constexpr (sizeof(T) == 4) {
			auto valid = _mm256_load_si256(reinterpret_cast<const __m256i*>(masks_[count].data()));
			auto value = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values));
			value = _mm256_sub_epi32(value, _mm256_and_si256(valid, _mm256_set1_epi32(std::bit_cast<T>(change))));
			_mm256_storeu_si256(reinterpret_cast<__m256i*>(values), value);
			return;
		}
#endif
		for (unsigned i = 0; i < count; ++i) values[i] -= change;
	}

	static void materialize(const U* values, T* result, unsigned count, U base) {
		if (count == branch) {
#if defined(__AVX512F__) && defined(__x86_64__)
			if constexpr (sizeof(T) == 8) {
				auto value = _mm512_loadu_si512(values);
				value = _mm512_add_epi64(value, _mm512_set1_epi64(std::bit_cast<T>(base)));
				_mm512_storeu_si512(result, value);
				return;
			}
#endif
#if defined(__AVX2__) && defined(__x86_64__)
			if constexpr (sizeof(T) == 8) {
				auto delta = _mm256_set1_epi64x(std::bit_cast<T>(base));
				for (unsigned i = 0; i < branch; i += 4) {
					auto value = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values + i));
					_mm256_storeu_si256(reinterpret_cast<__m256i*>(result + i), _mm256_add_epi64(value, delta));
				}
				return;
			}
			if constexpr (sizeof(T) == 4) {
				auto value = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(values));
				value = _mm256_add_epi32(value, _mm256_set1_epi32(std::bit_cast<T>(base)));
				_mm256_storeu_si256(reinterpret_cast<__m256i*>(result), value);
				return;
			}
			if constexpr (sizeof(T) == 2) {
				auto value = _mm_loadu_si128(reinterpret_cast<const __m128i*>(values));
				value = _mm_add_epi16(value, _mm_set1_epi16(std::bit_cast<T>(base)));
				_mm_storeu_si128(reinterpret_cast<__m128i*>(result), value);
				return;
			}
			if constexpr (sizeof(T) == 1) {
				auto value = _mm_loadl_epi64(reinterpret_cast<const __m128i*>(values));
				value = _mm_add_epi8(value, _mm_set1_epi8(std::bit_cast<T>(base)));
				_mm_storel_epi64(reinterpret_cast<__m128i*>(result), value);
				return;
			}
#endif
		}
		for (unsigned i = 0; i < count; ++i) result[i] = std::bit_cast<T>(U(values[i] + base));
	}
	template<bool Narrow>
	static U change_one(U *values, unsigned count, unsigned index, U delta, bool negative) {
		if (!delta) return 0;
		U old = values[index], next = U(old + delta);
		values[index] = next;
		if ((Narrow && sizeof(T) == 8) ? std::bit_cast<T>(next) < 0 : negative && next > old) {
			subtract(values, count, next);
			return next;
		}
		if (((Narrow && sizeof(T) == 8) || !negative) && !old) {
			U result = full_minimum<Narrow>(values);
			if (!result) return 0;
			subtract(values, count, result);
			return result;
		}
		return 0;
	}
};
}

#endif
Back to top page