__declspec(align(16)) struct MultiCrossMsg final {...

__declspec(align(16)) struct MultiCrossMsg final {
float intervals[15]; uint8_t count;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& intervals& count; }
};

__declspec(align(16)) struct BestSolutionMsg final {
float bestF, bestX, bestY, bestQ[32]; uint8_t dim;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& bestF& bestX& bestY& bestQ& dim; }
};

__declspec(align(16)) struct CtrlMsgND final {
uint8_t kind; CrossMsg xchg; MultiCrossMsg multiXchg; BestSolutionMsg bestSol;
template<typename Archive> __declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& kind; if (kind == 1) ar& xchg; else if (kind == 2) ar& multiXchg; else if (kind == 3) ar& bestSol;
}
};

static __forceinline int generate_lhs_seeds_lite(const MortonND& map, const int dim, float* __restrict S, int stride, uint32_t seed) {
int temp_dim = dim; const int ns = --temp_dim * temp_dim; uint32_t st = seed;
alignas(32) int permutations[32][256];
for (int d = 0; d < dim; ++d) {
for (int s = 0; s < ns; ++s) permutations[d][s] = s;
for (int s = ns - 1; s > 0; --s) { st ^= st << 13; st ^= st >> 17; st ^= st << 5; int j = st % (s + 1); std::swap(permutations[d][s], permutations[d][j]); }
}
for (int s = 0; s < ns; ++s) {
for (int d = 0; d < dim; ++d) {
st ^= st << 13; st ^= st >> 17; st ^= st << 5;
float u = (st & 0xFFFFFF) * 5.9604645e-8f;
int stratum = permutations[d][s];
float pos = ((float)stratum + u) / (float)ns;
int pd = map.perm[d];
float lo = map.low[pd], hi = map.high[pd];
S[s * stride + d] = fmaf(pos, (hi - lo), lo);
}
}
return ns;
}

static __forceinline int generate_heuristic_seeds(const ManipCost& cost, const MortonND& map, int dim, float* __restrict S, int stride, uint32_t seed) {
const int n = cost.n; const bool VL = cost.variableLen;
const float tx = cost.targetX, ty = cost.targetY;
int total_seeds = 0;

text
{
	float* s0 = S + total_seeds * stride;
	float phi = atan2f(ty, tx);
	float rho = sqrtf(fmaf(tx, tx, ty * ty));
	float len = fminf(fmaxf(rho / (float)n, 0.5f), 2.0f);
	for (int i = 0; i < n; ++i) s0[i] = phi / (float)n;
	if (VL) for (int i = 0; i < n; ++i) s0[n + i] = len;
	total_seeds++;
}

{
	float* s1 = S + total_seeds * stride;
	float phi = atan2f(ty, tx);
	for (int i = 0; i < n; ++i) s1[i] = 0.5f * phi * ((i & 1) ? -1.0f : 1.0f);
	if (VL) for (int i = 0; i < n; ++i) s1[n + i] = 1.0f * (0.8f + 0.4f * (float)i / (float)n);
	total_seeds++;
}

{
	float* s2 = S + total_seeds * stride;
	const float inv = (n > 1) ? 1.0f / (float)(n - 1) : 0.0f;
	float phi = atan2f(ty, tx);
	for (int i = 0; i < n; ++i) { float pr = (float)i * inv; s2[i] = phi * (1.0f - 0.3f * pr); }
	if (VL) { for (int i = 0; i < n; ++i) { float si = sinf(1.5f * (float)i); s2[n + i] = (1.0f + 0.2f * si); } }
	total_seeds++;
}

int lhs_count = generate_lhs_seeds_lite(map, dim, S + total_seeds * stride, stride, seed);
total_seeds += lhs_count;
return total_seeds;

}

static __forceinline void agp_run_branch_mpi(
const MortonND& map, const ManipCost& cost, int maxIter, float r, bool adaptive, float eps, unsigned seed,
std::vector<IntervalND*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY, float M_prior = 1e-3f)
{
const int n = cost.n;
const int dim = n + (cost.variableLen ? n : 0);
uint32_t exchange_counter_500 = 0;
uint32_t exchange_counter_T = 0;

text
alignas(32) float M_by_span[12]; for (int i = 0; i < 12; ++i) M_by_span[i] = M_prior;
float Mmax = M_prior;

alignas(64) float q_local[32], phi[32], s_arr[32], c_arr[32], sum_s[32], sum_c[32], q_try[32];
bestQ.reserve(dim);
float x = 0.0f, y = 0.0f;
int no_improve = 0;

auto t_to_idx = [&](float t) -> uint64_t {
	float tt = (t <= 0.0f) ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t);
	uint64_t idx = (uint64_t)((double)tt * (double)map.scale);
	if (idx >= map.scale) idx = map.scale - 1ull;
	return idx;
	};

auto update_pockets_and_Mmax = [&](IntervalND* I) {
	const int k = I->span_level; if (I->M > M_by_span[k]) M_by_span[k] = I->M; if (M_by_span[k] > Mmax) Mmax = M_by_span[k];
	};

float a = 0.0f, b = 1.0f;

auto evalAt = [&](float t) -> float {
	map.map01ToPoint(t, q_local);
	float f = cost(q_local, x, y);

	if (f < bestF * 1.25f) {
		float acc = 0.0f; for (int i = 0; i < n; ++i) { acc += q_local[i]; phi[i] = acc; }
		FABE13_SINCOS(phi, s_arr, c_arr, n);
		float as = 0.0f, ac = 0.0f;
		for (int k = n - 1; k >= 0; --k) {
			const float Lk = cost.variableLen ? q_local[n + k] : 1.0f;
			as += Lk * s_arr[k]; ac += Lk * c_arr[k]; sum_s[k] = as; sum_c[k] = ac;
		}
		const float dx = x - cost.targetX, dy = y - cost.targetY;
		float dist = sqrtf(fmaf(dx, dx, dy * dy)) + 1e-8f;

		float eta = 0.125f;
		for (int stepI = 0; stepI < 3; ++stepI) {
			for (int i = 0; i < n; ++i) {
				float gpen = 0.0f;
				{
					float ai = fabsf(q_local[i]);
					float v = cost.minTheta - ai;
					if (v > 0.0f) {
						float scale = 2.0f / (cost.minTheta + 1e-6f);
						float e = exp2f(scale * v);
						float dpen_dtheta = cost.sharpW * (e * 0.69314718055994530941723212145818f * scale) * (-copysignf(1.0f, q_local[i]));
						gpen += dpen_dtheta;
					}
				}
				{
					float tsg = -q_local[i] * cost.archBiasK;
					float sig = 1.0f / (1.0f + expf(-tsg));
					gpen += -(cost.archBiasW * cost.archBiasK) * sig;
				}

				float g = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen;
				q_try[i] = q_local[i] - eta * g;

				const float deg2rad = 3.14159265358979323846f / 180.0f;
				const float lo0 = -60.0f * deg2rad, hi0 = 150.0f * deg2rad;
				const float lo = -150.0f * deg2rad, hi = 150.0f * deg2rad;
				const float Lb = (i == 0) ? lo0 : lo;
				const float Hb = (i == 0) ? hi0 : hi;
				if (q_try[i] < Lb) q_try[i] = Lb; else if (q_try[i] > Hb) q_try[i] = Hb;
			}
			if (cost.variableLen) {
				for (int i = 0; i < n; ++i) {
					float g = (dx * c_arr[i] + dy * s_arr[i]) / dist;
					float v = q_local[n + i] - eta * g;
					if (v < 0.5f) v = 0.5f; else if (v > 2.0f) v = 2.0f;
					q_try[n + i] = v;
				}
			}
			float x2, y2; float f2 = cost(q_try, x2, y2);
			if (f2 < f) { memcpy(q_local, q_try, dim * sizeof(float)); f = f2; x = x2; y = y2; break; }
			eta *= 0.5f;
		}

		const int last = n - 1;
		const float deg2rad = 3.14159265358979323846f / 180.0f;
		const float lo = (last == 0) ? (-60.0f * deg2rad) : (-150.0f * deg2rad);
		const float hi = 150.0f * deg2rad;
		float bestLocF = f;
		float saved = q_local[last];
		for (float delta = 0.05f; delta >= 0.00625f; delta *= 0.5f) {
			for (int sgn = -1; sgn <= 1; sgn += 2) {
				float cand = saved + sgn * delta;
				if (cand < lo) cand = lo; else if (cand > hi) cand = hi;
				float backup = q_local[last]; q_local[last] = cand;
				float x2, y2; float f2 = cost(q_local, x2, y2);
				if (f2 < bestLocF) { bestLocF = f2; x = x2; y = y2; saved = cand; }
				q_local[last] = backup;
			}
		}
		if (bestLocF < f) { q_local[last] = saved; f = bestLocF; }
	}

	if (f < bestF) { bestF = f; bestQ.assign(q_local, q_local + dim); bestX = x; bestY = y; no_improve = 0; }
	else { ++no_improve; }
	return f;
	};

float f_a = evalAt(a), f_b = evalAt(b);
const int K = (std::min)((std::max)(2 * dim, 8), 128);

H.reserve((size_t)maxIter + K + 16);
const int rank = g_world->rank();
const int world = g_world->size();

alignas(64) float seeds[256 * 32];
const int seedCnt = generate_heuristic_seeds(cost, map, dim, seeds, 32, seed + rank * 7919u);

for (int i = 0; i < seedCnt; ++i) {
	const float* s = seeds + i * 32;
	float t_seed = map.pointToT(s);
	float interval_size = (i < 3) ? (0.0004f * (float)dim)
		: (0.00031f * (float)dim) * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(0.00025f / 0.00031f) * (float)(i - 3));
	float t1 = fmaxf(a, t_seed - interval_size), t2 = fminf(b, t_seed + interval_size);
	if (t2 <= t1) continue;
	alignas(64) float q1[32], q2[32]; float x1, y1, x2, y2;
	map.map01ToPoint(t1, q1); float f1 = cost(q1, x1, y1);
	map.map01ToPoint(t2, q2); float f2 = cost(q2, x2, y2);
	IntervalND* I = new IntervalND(t1, t2, f1, f2);
	I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2);
	I->diam = map.block_diameter(I->i1, I->i2);
	I->compute_span_level(map);
	I->set_metric(I->diam);
	update_pockets_and_Mmax(I);
	I->ChangeCharacteristic(r * Mmax);
	if (i < 3) I->R *= fmaf(0.01f, (float)dim, 0.85f);
	else {
		float start_mult = 0.214f * (float)dim;
		float end_mult = 0.174f * (float)dim;
		float mult = start_mult * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_mult / start_mult) * (float)(i - 3));
		I->R *= mult;
	}
	H.emplace_back(I);
	std::push_heap(H.begin(), H.end(), ComparePtrND);
	if (f1 < bestF) { bestF = f1; bestQ.assign(q1, q1 + dim); bestX = x1; bestY = y1; }
	if (f2 < bestF) { bestF = f2; bestQ.assign(q2, q2 + dim); bestX = x2; bestY = y2; }
}

float prev_t = a, prev_f = f_a;
for (int k = 1; k <= K; ++k) {
	float t = a + (b - a) * ((float)k / (K + 1)) + (float)rank / (float)(world * (K + 1));
	float f = evalAt(t);
	IntervalND* I = new IntervalND(prev_t, t, prev_f, f);
	I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t);
	I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); I->set_metric(I->diam);
	update_pockets_and_Mmax(I);
	I->ChangeCharacteristic(r * Mmax);
	H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND);
	prev_t = t; prev_f = f;
}
IntervalND* tail = new IntervalND(prev_t, b, prev_f, f_b);
tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b);
tail->diam = map.block_diameter(tail->i1, tail->i2); tail->compute_span_level(map); tail->set_metric(tail->diam);
update_pockets_and_Mmax(tail);
tail->ChangeCharacteristic(r * Mmax);
H.emplace_back(tail); std::push_heap(H.begin(), H.end(), ComparePtrND);

float dmax = b - a, initial_len = dmax, thr03 = 0.3f * initial_len, inv_thr03 = 1.0f / thr03;
int it = 0;

auto kickEveryByDim = [&](int dim) -> int {
	float z = 120.0f * exp2f(-0.05f * (float)dim);
	if (z < 60.0f) z = 60.0f;
	return (int)z;
	};

auto noImproveThrByDim = [&](int dim) -> int {
	float z = 80.0f * exp2f(-0.08f * (float)dim);
	if (z < 30.0f) z = 30.0f;
	return (int)z;
	};

while (it < maxIter) {
	if ((it % kickEveryByDim(dim)) == 0 && no_improve > noImproveThrByDim(dim)) {
		float t_best = map.pointToT(bestQ.data());
		for (int i = 0; i < 2; ++i) {
			float off = (i == 0) ? 0.01f : -0.01f;
			float t_seed = fminf(b, fmaxf(a, t_best + off));
			float f_seed = evalAt(t_seed);
			IntervalND* J = new IntervalND(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed);
			J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f);
			J->diam = map.block_diameter(J->i1, J->i2); J->compute_span_level(map); J->set_metric(J->diam);
			update_pockets_and_Mmax(J);
			J->ChangeCharacteristic(r * Mmax); J->R *= 0.9f;
			H.emplace_back(J); std::push_heap(H.begin(), H.end(), ComparePtrND);
		}
		no_improve = 0;
	}

	const float p = fmaf(-1.0f / initial_len, dmax, 1.0f);
	bool   stagnation = (no_improve > 100) && (it > 270);

	float A = 200.0f + 64.0f * exp2f(-0.06f * (float)dim);
	float B = 210.0f + 67.0f * exp2f(-0.06f * (float)dim);
	const int T = (int)fmaf(-expm1f(p), A, B);

	float  r_eff = fmaxf(1.0f, r * (0.7f + 0.3f * (1.0f - p)));

	std::pop_heap(H.begin(), H.end(), ComparePtrND);
	IntervalND* cur = H.back(); H.pop_back();

	const float x1 = cur->x1, x2 = cur->x2, y1 = cur->y1, y2 = cur->y2;
	float m = r_eff * Mmax;
	float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r);
	tNew = fminf(fmaxf(tNew, a), b);
	float fNew = evalAt(tNew);

	IntervalND* L = new IntervalND(x1, tNew, y1, fNew);
	IntervalND* Rv = new IntervalND(tNew, x2, fNew, y2);

	L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew);
	Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2);
	L->diam = map.block_diameter(L->i1, L->i2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2);
	L->compute_span_level(map); Rv->compute_span_level(map);
	L->set_metric(L->diam); Rv->set_metric(Rv->diam);

	float Mloc = (std::max)(L->M, Rv->M);
	update_pockets_and_Mmax(L); update_pockets_and_Mmax(Rv);

	const float prevMmax = Mmax; if (Mloc > Mmax) Mmax = Mloc;
	m = r_eff * Mmax;

	if (adaptive) {
		float len1 = tNew - x1, len2 = x2 - tNew;
		if (len1 + len2 == dmax) {
			dmax = (std::max)(len1, len2);
			for (auto pI : H) { float Ls = pI->x2 - pI->x1; if (Ls > dmax) dmax = Ls; }
		}
		if ((thr03 > dmax && !(it % 3)) || (10.0f * dmax < initial_len)) {
			const float progress = fmaf(-dmax, inv_thr03, 1.0f);
			const float alpha = progress * progress;
			const float beta = fmaf(-alpha, 1.0f, 2.0f);
			const float MULT = (1.0f / dmax) * Mmax;
			const float global_coeff = fmaf(MULT, r_eff, -MULT);
			const float GF = fmaf(beta, global_coeff, 0.0f);
			L->ChangeCharacteristic(fmaf(GF, len1, L->M * alpha));
			Rv->ChangeCharacteristic(fmaf(GF, len2, Rv->M * alpha));
			size_t sz = H.size();
			RecomputeR_AffineM_AVX2_ND(H.data(), sz, GF, alpha);
			std::make_heap(H.begin(), H.end(), ComparePtrND);
		}
		else {
			if (Mloc > prevMmax) {
				L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m);
				if (Mloc > 1.15f * prevMmax) {
					size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND);
				}
			}
			else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); }
		}
	}
	else {
		if (Mloc > prevMmax) {
			L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m);
			if (Mloc > 1.15f * prevMmax) { size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); }
		}
		else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); }
	}

	H.push_back(L);  std::push_heap(H.begin(), H.end(), ComparePtrND);
	H.push_back(Rv); std::push_heap(H.begin(), H.end(), ComparePtrND);
	if (!H.empty()) {
		_mm_prefetch((const char*)H[0], _MM_HINT_T0);
		if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0);
	}

	IntervalND* top = H.front();
	float interval_len = top->x2 - top->x1;

	bool want_term = (exp2f(log2f(interval_len) / (float)dim) < eps) || (it == maxIter - 1);
	if (!(it % T) || want_term) {
		CtrlMsgND out; out.kind = want_term ? 0 : 2;
		if (!want_term) {
			uint8_t cnt = (uint8_t)((H.size() >= 3) ? 3 : H.size());
			out.multiXchg.count = cnt;
			float* dest = out.multiXchg.intervals;
			IntervalND* t1 = H[0];
			IntervalND* t2 = (H.size() > 1 ? H[1] : H[0]);
			IntervalND* t3 = (H.size() > 2 ? H[2] : H[H.size() - 1]);
			IntervalND* tops[3] = { t1, t2, t3 };
			for (uint8_t i2 = 0; i2 < cnt; ++i2) {
				IntervalND* Tt = tops[i2];
				dest[0] = Tt->x1; dest[1] = 0.0f; dest[2] = Tt->x2; dest[3] = 0.0f; dest[4] = Tt->R;
				dest += 5;
			}
			const size_t iterations = std::bit_width(static_cast<size_t>(world - 1));
			bool active = true;
			bool invert_T = ++exchange_counter_T & 1;

			for (size_t i = 0; i < iterations && active; ++i) {
				const size_t step = 1ULL << i;
				int partner = rank ^ step;

				if (partner >= world) continue;

				bool am_sender = !!(rank & step) ^ invert_T;

				if (am_sender) {
					g_world->isend(partner, 0, out);
					active = false;
				}
			}
		}
		else {
			for (int i2 = 0; i2 < world; ++i2) if (i2 != rank) g_world->isend(i2, 0, out);
			return;
		}
	}

	if (!(it % 500)) {
		CtrlMsgND out; out.kind = 3; out.bestSol.bestF = bestF; out.bestSol.bestX = bestX; out.bestSol.bestY = bestY;
		out.bestSol.dim = (uint8_t)bestQ.size(); memcpy(out.bestSol.bestQ, bestQ.data(), bestQ.size() * sizeof(float));
		const size_t iterations = std::bit_width(static_cast<size_t>(world - 1));
		bool active = true;
		bool invert_T = ++exchange_counter_500 & 1;

		for (size_t i = 0; i < iterations && active; ++i) {
			const size_t step = 1ULL << i;
			int partner = rank ^ step;

			if (partner >= world) continue;

			bool am_sender = !!(rank & step) ^ invert_T;

			if (am_sender) {
				g_world->isend(partner, 0, out);
				active = false;
			}
		}
	}

	while (g_world->iprobe(boost::mpi::any_source, 0)) {
		CtrlMsgND in; g_world->recv(boost::mpi::any_source, 0, in);
		if (in.kind == 0) return;
		else if (in.kind == 1) {
			float sx = in.xchg.s_x1, ex = in.xchg.e_x1;
			if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f;
			if (ex > sx) {
				alignas(64) float tmp[32]; float tx, ty;
				map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty);
				map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty);
				IntervalND* inj = new IntervalND(sx, ex, y1i, y2i);
				inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex);
				inj->diam = map.block_diameter(inj->i1, inj->i2);
				inj->compute_span_level(map); inj->set_metric(inj->diam);
				update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax);
				if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); }
				IntervalND* topH = H.front();
				if (inj->R > 1.15f * topH->R) {
					float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f);
					float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f)
						: fmaf(0.3491860241215959f, expm1f(p2), 0.6f);
					inj->R = in.xchg.Rtop * k;
					H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND);
				}
			}
		}
		else if (in.kind == 2) {
			const MultiCrossMsg& mX = in.multiXchg;
			for (uint8_t ii = 0; ii < mX.count; ++ii) {
				const float* d = &mX.intervals[ii * 5];
				float sx = d[0], ex = d[2];
				if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f;
				if (ex > sx) {
					alignas(64) float tmp[32]; float tx, ty;
					map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty);
					map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty);
					IntervalND* inj = new IntervalND(sx, ex, y1i, y2i);
					inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex);
					inj->diam = map.block_diameter(inj->i1, inj->i2);
					inj->compute_span_level(map); inj->set_metric(inj->diam);
					update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax);
					if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); }
					IntervalND* topH = H.front();
					if (inj->R > 1.15f * topH->R) {
						float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f);
						float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f)
							: fmaf(0.3491860241215959f, expm1f(p2), 0.6f);
						inj->R = d[4] * k;
						H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND);
					}
				}
			}
		}
		else if (in.kind == 3) {
			const BestSolutionMsg& bm = in.bestSol;
			if (bm.bestF < bestF * 1.15f) {
				alignas(64) float tmp_q[32]; memcpy(tmp_q, bm.bestQ, bm.dim * sizeof(float));
				float t_best = map.pointToT(tmp_q);
				float t1 = fmaxf(a, t_best - 0.001f), t2 = fminf(b, t_best + 0.001f);
				if (t2 > t1) {
					alignas(64) float tq1[32], tq2[32]; float xx1, yy1, xx2, yy2;
					map.map01ToPoint(t1, tq1); float f1 = cost(tq1, xx1, yy1);
					map.map01ToPoint(t2, tq2); float f2 = cost(tq2, xx2, yy2);
					IntervalND* I = new IntervalND(t1, t2, f1, f2);
					I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2);
					I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map);
					I->set_metric(I->diam); update_pockets_and_Mmax(I);
					I->ChangeCharacteristic(r * Mmax); I->R *= 0.90f;
					H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND);
				}
				if (bm.bestF < bestF) { bestF = bm.bestF; bestX = bm.bestX; bestY = bm.bestY; bestQ.assign(bm.bestQ, bm.bestQ + bm.dim); }
			}
		}
	}
	++it;
}

}

extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_Manip2D(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon,
unsigned int seed, float** out_bestQ, size_t* out_bestQLen, float* out_bestX,
float* out_bestY, float* out_bestF)
{
Slab* const __restrict slab = tls.local(); slab->current = slab->base;
while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
const int dim = nSegments + (variableLengths ? nSegments : 0);

text
g_mc.permCache.resize(dim); for (int i = 0; i < dim; ++i) g_mc.permCache[i] = i;
uint32_t s = g_mc.baseSeed;
for (int i = dim - 1; i > 0; --i) { s ^= s << 13; s ^= s >> 17; s ^= s << 5; uint32_t j = s % (uint32_t)(i + 1); std::swap(g_mc.permCache[i], g_mc.permCache[j]); }
g_mc.invMaskCache.resize(dim); for (int k = 0; k < dim; ++k) { s ^= s << 13; s ^= s >> 17; s ^= s << 5; g_mc.invMaskCache[k] = (uint64_t)s; }

const float deg2rad = 3.14159265358979323846f / 180.0f;
const float theta0Min = -60.0f * deg2rad, theta0Max = 150.0f * deg2rad;
const float thetaMin = -150.0f * deg2rad, thetaMax = 150.0f * deg2rad;
const float lenMin = 0.5f, lenMax = 2.0f;

std::vector<float> low;  low.reserve(dim);
std::vector<float> high; high.reserve(dim);
for (int i = 0; i < nSegments; ++i) { low.push_back(i == 0 ? theta0Min : thetaMin); high.push_back(i == 0 ? theta0Max : thetaMax); }
if (variableLengths) { for (int i = 0; i < nSegments; ++i) { low.push_back(lenMin); high.push_back(lenMax); } }

ManipCost cost(nSegments, variableLengths, targetX, targetY, minTheta);

const int rank = g_world->rank(), world = g_world->size();
std::vector<float> bestQ; float bestF = FLT_MAX, bx = 0.0f, by = 0.0f;

const int levels0 = (std::min)(peanoLevels, 8);
const int maxIter0 = (int)(maxIterPerBranch * 0.2f);
MortonND map0(dim, levels0, low.data(), high.data(), g_mc);

std::vector<IntervalND*> H_coarse; std::vector<float> bestQ_coarse;
float bestF_coarse = FLT_MAX, bx_coarse = 0.0f, by_coarse = 0.0f;

float M_prior = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << levels0));
if (variableLengths) M_prior += 1.41421356237f * (1.0f / (float)(1u << levels0));

agp_run_branch_mpi(map0, cost, maxIter0, r, adaptiveMode, epsilon, seed,
	H_coarse, bestQ_coarse, bestF_coarse, bx_coarse, by_coarse, M_prior);

if (bestF_coarse < bestF) { bestF = bestF_coarse; bestQ = bestQ_coarse; bx = bx_coarse; by = by_coarse; }

if (levels0 < peanoLevels) {
	while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
	MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc);
	std::vector<IntervalND*> H_fine;
	std::vector<float> bestQ_fine = bestQ;
	float bestF_fine = bestF, bx_fine = bx, by_fine = by;

	float M_prior_fine = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << peanoLevels));
	if (variableLengths) M_prior_fine += 1.41421356237f * (1.0f / (float)(1u << peanoLevels));

	std::sort(H_coarse.begin(), H_coarse.end(), [](const IntervalND* a, const IntervalND* b) { return a->R < b->R; });
	const size_t topCount = (size_t)(H_coarse.size() * 0.3f);
	auto t_to_idx_fine = [&](float t) -> uint64_t {
		float tt = (t <= 0.0f) ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t);
		uint64_t idx = (uint64_t)((double)tt * (double)map1.scale);
		if (idx >= map1.scale) idx = map1.scale - 1ull;
		return idx;
		};
	for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) {
		const IntervalND* C = H_coarse[i];
		alignas(64) float q1[32], q2[32]; float x1, y1, x2, y2;
		map1.map01ToPoint(C->x1, q1); float f1 = cost(q1, x1, y1);
		map1.map01ToPoint(C->x2, q2); float f2 = cost(q2, x2, y2);
		IntervalND* I = new IntervalND(C->x1, C->x2, f1, f2);
		I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2);
		I->diam = map1.block_diameter(I->i1, I->i2);
		I->set_metric(I->diam);
		H_fine.push_back(I);
		if (f1 < bestF_fine) { bestF_fine = f1; bestQ_fine.assign(q1, q1 + dim); bx_fine = x1; by_fine = y1; }
		if (f2 < bestF_fine) { bestF_fine = f2; bestQ_fine.assign(q2, q2 + dim); bx_fine = x2; by_fine = y2; }
	}
	std::make_heap(H_fine.begin(), H_fine.end(), ComparePtrND);
	agp_run_branch_mpi(map1, cost, maxIterPerBranch - maxIter0, r, adaptiveMode, epsilon, seed,
		H_fine, bestQ_fine, bestF_fine, bx_fine, by_fine, M_prior_fine);

	if (bestF_fine < bestF) { bestF = bestF_fine; bestQ = bestQ_fine; bx = bx_fine; by = by_fine; }
}

BestSolutionMsg best;
best.bestF = bestF;
best.bestX = bx;
best.bestY = by;
best.dim = (uint8_t)bestQ.size();
memcpy(best.bestQ, bestQ.data(), best.dim * sizeof(float));

const size_t iterations = std::bit_width(static_cast<size_t>(world - 1));
bool active = true;

for (size_t i = 0; i < iterations && active; ++i) {
	const size_t step = 1ULL << i;
	int partner = rank ^ step;

	if (partner >= world) continue;

	bool am_sender = (rank & step) != 0U;

	if (am_sender) {
		g_world->isend(partner, 2, best);
		active = false;
	}
	else {
		BestSolutionMsg in;
		g_world->recv(partner, 2, in);
		if (in.bestF < best.bestF) best = in;
	}
}

if (rank == 0) {
	*out_bestQLen = (size_t)best.dim;
	*out_bestQ = (float*)CoTaskMemAlloc(sizeof(float) * (*out_bestQLen));
	memcpy(*out_bestQ, best.bestQ, sizeof(float) * (*out_bestQLen));
	*out_bestX = best.bestX;
	*out_bestY = best.bestY;
	*out_bestF = best.bestF;
}

}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline int AgpInit(int peanoLevel, float a, float b, float c, float d) {
g_env = new boost::mpi::environment();
g_world = new boost::mpi::communicator();
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
const int rank = g_world->rank();
const int world_size = g_world->size();
if (world_size == 4) { new (&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, rank & 3); }
g_mc.baseSeed = fmaf(0x9E3779B9u, rank, 0x9E3779B9u);
return rank;
}

__forceinline float ShekelFunc(float x, float seed) {
int i = 0; float st = seed, r1, r2, res = 0.0f;
while (i < 10) {
XOR_RAND(st, r1); float xp = fmaf(-r1, 10.0f, x); XOR_RAND(st, r1); XOR_RAND(st, r2);
float d = fmaf(fmaf(r1, 20.0f, 5.0f), xp * xp, fmaf(r2, 0.2f, 1.0f)); d = copysignf(fmaxf(fabsf(d), FLT_MIN), d); res -= 1.0f / d; ++i;
}
return res;
}

__forceinline float RastriginFunc(float x1, float x2) {
const float t = fmaf(x1, x1, x2 * x2); float c1, c2;
FABE13_COS(6.28318530717958647692f * x1, c1); FABE13_COS(6.28318530717958647692f * x2, c2);
return (t - fmaf(c1 + c2, 10.0f, -14.6f)) * fmaf(-t, 0.25f, 18.42f);
}

__forceinline float HillFunc(float x, float seed) {
int j = 0; __declspec(align(32)) float ang[14u];
float st = 6.28318530717958647692f * x; while (j < 14) { ang[j] = st * (float)(j + 1); ++j; }
__declspec(align(32)) float sv[14u], cv[14u]; FABE13_SINCOS(ang, sv, cv, 14u);
float state = seed, r1, r2; XOR_RAND(state, r1); float res = fmaf(r1, 2.0f, -1.1f); --j;
while (j >= 0) { XOR_RAND(state, r1); XOR_RAND(state, r2); res += fmaf(fmaf(r1, 2.0f, -1.1f), sv[j], fmaf(r2, 2.0f, -1.1f) * cv[j]); --j; }
return res;
}

__forceinline float GrishaginFunc(float x1, float x2, float seed) {
int j = 0; __declspec(align(32)) float aj[8u], ak[8u];
while (j < 8) { float pj = 3.14159265358979323846f * (float)(j + 1); aj[j] = pj * x1; ak[j] = pj * x2; ++j; }
__declspec(align(32)) float sj[8u], cj[8u], sk[8u], ck[8u]; FABE13_SINCOS(aj, sj, cj, 8u); FABE13_SINCOS(ak, sk, ck, 8u);
--j; float p1 = 0.0f, p2 = 0.0f; float st = seed, r1, r2;
while (j >= 0) {
size_t k = 0u; while (k < 8u) {
float s = sj[j] * sj[j]; float c = ck[k] * ck[k];
XOR_RAND_GRSH(st, r1); XOR_RAND_GRSH(st, r2); p1 = fmaf(r1, s, fmaf(r2, c, p1));
XOR_RAND_GRSH(st, r1); XOR_RAND_GRSH(st, r2); p2 = fmaf(-r1, c, fmaf(r2, s, p2)); ++k;
}
--j;
}
return -sqrtf(fmaf(p1, p1, p2 * p2));
}

extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_1D(float global_iterations, float a, float b, float r, bool mode, float epsilon, float seed,
float** out_data, size_t* out_len) {
Slab* const __restrict slab = tls.local(); slab->current = slab->base; int schetchick = 0;
const float initial_length = b - a; float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length, inv_threshold_03 = 1.0f / threshold_03;
const float start_val = ShekelFunc(a, seed); float best_f = ShekelFunc(b, seed);
float x_Rmax_1 = a, x_Rmax_2 = b; float y_Rmax_1 = start_val, y_Rmax_2 = best_f;
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval1D*, boost::alignment::aligned_allocator<Interval1D*, 64u>> R;
Extr.reserve((size_t)global_iterations << 2u); R.reserve((size_t)global_iterations << 1u);
R.emplace_back(new Interval1D(a, b, start_val, best_f, 1.0f)); float Mmax = R.front()->M; float m = r * Mmax;

text
while (true) {
	const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 1.0f, r);
	const float new_value = ShekelFunc(new_point, seed);
	if (new_value < best_f) { best_f = new_value; Extr.emplace_back(best_f); Extr.emplace_back(new_point); }
	std::pop_heap(R.begin(), R.end(), ComparePtr1D);
	const Interval1D* pro = R.back();
	const float new_x1 = pro->x1, new_x2 = pro->x2;
	const float len2 = new_x2 - new_point, len1 = new_point - new_x1;
	const float interval_len = (len1 < len2 ? len1 : len2);
	if (++schetchick == (int)global_iterations || interval_len < epsilon) {
		Extr.emplace_back((float)schetchick); Extr.emplace_back(interval_len);
		*out_len = Extr.size(); *out_data = (float*)CoTaskMemAlloc(sizeof(float) * (*out_len));
		memcpy(*out_data, Extr.data(), sizeof(float) * (*out_len)); return;
	}
	Interval1D* curr = new Interval1D(new_x1, new_point, pro->y1, new_value, 1.0f);
	Interval1D* curr1 = new Interval1D(new_point, new_x2, new_value, pro->y2, 1.0f);
	const float currM = curr->M > curr1->M ? curr->M : curr1->M;
	const size_t r_size = R.size();
	if (mode) {
		if (len2 + len1 == dmax) { dmax = len2 > len1 ? len2 : len1; for (auto p : R) { float L = p->x2 - p->x1; if (L > dmax) dmax = L; } }
		if (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) {
			if (currM > Mmax) { Mmax = currM; m = r * Mmax; }
			const float progress = fmaf(-inv_threshold_03, dmax, 1.0f);
			const float alpha = progress * progress; const float betta = 2.0f - alpha;
			const float MULT = (1.0f / dmax) * Mmax;
			const float global_coeff = fmaf(MULT, r, -MULT);
			const float GF = betta * global_coeff;
			curr->ChangeCharacteristic(fmaf(GF, len1, curr->M * alpha));
			curr1->ChangeCharacteristic(fmaf(GF, len2, curr1->M * alpha));
			RecomputeR_AffineM_AVX2_1D(R.data(), r_size, GF, alpha);
			std::make_heap(R.begin(), R.end(), ComparePtr1D);
		}
		else {
			if (currM > Mmax) {
				if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
				else { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D); }
			}
			else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
		}
	}
	else {
		if (currM > Mmax) {
			if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
			else { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D); }
		}
		else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
	}
	R.back() = curr; std::push_heap(R.begin(), R.end(), ComparePtr1D);
	R.emplace_back(curr1); std::push_heap(R.begin(), R.end(), ComparePtr1D);
	const Interval1D* top = R.front();
	x_Rmax_1 = top->x1; x_Rmax_2 = top->x2; y_Rmax_1 = top->y1; y_Rmax_2 = top->y2;
}

}

extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_2D(const float N, const float global_iterations, const float a, const float b, const float c,
const float d, const float r, const bool mode, const float epsilon, const float seed,
float** const __restrict out_data, size_t* const __restrict out_len) noexcept
{
Slab* const __restrict slab = tls.local(); slab->current = slab->base;
int schetchick = 0, no_improve = 0; const int rank = g_world->rank(); const int world_size = g_world->size();
while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1));
const float x_addition = (b - a) * inv_divider, y_addition = (d - c) * inv_divider;
const float true_start = a + x_addition, true_end = b - x_addition;
float x_Rmax_1 = true_start, x_Rmax_2 = true_end;
const float initial_length = x_Rmax_2 - x_Rmax_1; float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length, inv_threshold_03 = 1.0f / threshold_03;
const float start_val = rank % 3 ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition);
float best_f = rank % 2 ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition);
float y_Rmax_1 = start_val, y_Rmax_2 = best_f;
std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval1D* __restrict, boost::alignment::aligned_allocator<Interval1D* __restrict, 64u>> R;
Extr.clear(); Extr.reserve(static_cast<size_t>(global_iterations) << 2u);
R.clear(); R.reserve(static_cast<size_t>(global_iterations) << 1u);
R.emplace_back(new Interval1D(true_start, true_end, start_val, best_f, 2.0f));
const Interval1D* __restrict top_ptr;
float Mmax = R.front()->M, m = r * Mmax;
while (true) {
const float interval_len = x_Rmax_2 - x_Rmax_1;
const bool stagnation = no_improve > 100 && schetchick > 270;
const float p = fmaf(-1.0f / initial_length, dmax, 1.0f);
while (g_world->iprobe(boost::mpi::any_source, 0)) {
CtrlMsg in; g_world->recv(boost::mpi::any_source, 0, in);
if (in.kind) {
if (!rank) {
Extr.emplace_back((float)schetchick); Extr.emplace_back(interval_len); *out_len = Extr.size();
out_data = reinterpret_cast<float __restrict>(CoTaskMemAlloc(sizeof(float) * (*out_len)));
memcpy(*out_data, Extr.data(), sizeof(float) * (out_len));
}
return;
}
const float sx = FindX2D_analytic(in.xchg.s_x1, in.xchg.s_x2);
const float ex = FindX2D_analytic(in.xchg.e_x1, in.xchg.e_x2);
Interval1D const __restrict injected = new Interval1D(sx, ex,
RastriginFunc(in.xchg.s_x1, in.xchg.s_x2), RastriginFunc(in.xchg.e_x1, in.xchg.e_x2), 2.0f);
injected->ChangeCharacteristic(m);
if (injected->R > 1.15f * top_ptr->R) {
const float k = stagnation ? fmaf(0.5819767068693265f, expm1f(p), 0.3f) : fmaf(0.3491860241215959f, expm1f(p), 0.6f);
injected->R = in.xchg.Rtop * k;
R.emplace_back(injected); std::push_heap(R.begin(), R.end(), ComparePtr1D);
}
}
const int T = (int)fmaf(-expm1f(p), 264.0f, 277.0f);
const bool want_term = interval_len < epsilon || schetchick == static_cast<int>(global_iterations);
if (!(++schetchick % T) || stagnation || want_term) {
CtrlMsg out; out.kind = want_term;
if (!out.kind) {
float s_x1, s_x2, e_x1, e_x2; HitTest2D_analytic(top_ptr->x1, s_x1, s_x2); HitTest2D_analytic(top_ptr->x2, e_x1, e_x2);
out.xchg = CrossMsg{ s_x1, s_x2, e_x1, e_x2, top_ptr->R };
}
for (int i = 0; i < world_size; ++i) if (i != rank) g_world->isend(i, 0, out);
if (out.kind) {
if (!rank) {
Extr.emplace_back((float)schetchick); Extr.emplace_back(interval_len); *out_len = Extr.size();
out_data = reinterpret_cast<float __restrict>(CoTaskMemAlloc(sizeof(float) * (out_len)));
memcpy(out_data, Extr.data(), sizeof(float) * (out_len));
}
return;
}
}
const float new_point = Shag(m, x_Rmax_1, x_Rmax_2, y_Rmax_1, y_Rmax_2, 2.0f, r);
float new_x1_val, new_x2_val; HitTest2D_analytic(new_point, new_x1_val, new_x2_val);
const float new_value = RastriginFunc(new_x1_val, new_x2_val);
if (new_value < best_f) { best_f = new_value; Extr.emplace_back(best_f); Extr.emplace_back(new_x1_val); Extr.emplace_back(new_x2_val); no_improve = 0; }
else { ++no_improve; }
std::pop_heap(R.begin(), R.end(), ComparePtr1D);
Interval1D const __restrict promej = R.back();
const float segment_x1 = promej->x1, segment_x2 = promej->x2;
const float len2 = segment_x2 - new_point, len1 = new_point - segment_x1;
Interval1D const __restrict curr = new Interval1D(segment_x1, new_point, promej->y1, new_value, 2.0f);
Interval1D const __restrict curr1 = new Interval1D(new_point, segment_x2, new_value, promej->y2, 2.0f);
const float currM = (std::max)(curr->M, curr1->M); const size_t r_size = R.size();
if (mode) {
if (len2 + len1 == dmax) { dmax = (std::max)(len1, len2); for (auto pI : R) { float L = pI->x2 - pI->x1; if (L > dmax) dmax = L; } }
if (threshold_03 > dmax && !(schetchick % 3) || 10.0f * dmax < initial_length) {
if (currM > Mmax) { Mmax = currM; m = r * Mmax; }
const float progress = fmaf(-inv_threshold_03, dmax, 1.0f);
const float alpha = progress * progress; const float betta = 2.0f - alpha;
const float MULTIPLIER = (1.0f / dmax) * Mmax;
const float global_coeff = fmaf(MULTIPLIER, r, -MULTIPLIER);
const float GLOBAL_FACTOR = betta * global_coeff;
curr->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len1, curr->M * alpha));
curr1->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, len2, curr1->M * alpha));
RecomputeR_AffineM_AVX2_1D(R.data(), r_size, GLOBAL_FACTOR, alpha);
std::make_heap(R.begin(), R.end(), ComparePtr1D);
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
else {
Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D);
}
}
else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
}
}
else {
if (currM > Mmax) {
if (currM < 1.15f * Mmax) { Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
else {
Mmax = currM; m = r * Mmax; curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m);
RecomputeR_ConstM_AVX2_1D(R.data(), r_size, m); std::make_heap(R.begin(), R.end(), ComparePtr1D);
}
}
else { curr->ChangeCharacteristic(m); curr1->ChangeCharacteristic(m); }
}
R.back() = curr; std::push_heap(R.begin(), R.end(), ComparePtr1D);
R.emplace_back(curr1); std::push_heap(R.begin(), R.end(), ComparePtr1D);
top_ptr = R.front(); x_Rmax_1 = top_ptr->x1; x_Rmax_2 = top_ptr->x2; y_Rmax_1 = top_ptr->y1; y_Rmax_2 = top_ptr->y2;
}
}

struct RunParams {
int32_t nSegments;
uint8_t varLen;
float minTheta;
float tx, ty;
int32_t levels;
int32_t maxIter;
float r;
uint8_t adaptive;
float eps;
uint32_t seed;

text
template<typename Archive>
void serialize(Archive& ar, const unsigned int) {
	ar& nSegments& varLen& minTheta& tx& ty
		& levels& maxIter& r& adaptive& eps& seed;
}

};

extern "C" __declspec(dllexport) __declspec(noalias)
void AgpStartManipND(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon, unsigned int seed) {
RunParams p; p.nSegments = nSegments; p.varLen = (uint8_t)variableLengths; p.minTheta = minTheta; p.tx = targetX; p.ty = targetY;
p.levels = peanoLevels; p.maxIter = maxIterPerBranch; p.r = r; p.adaptive = (uint8_t)adaptiveMode; p.eps = epsilon; p.seed = seed;
int i = 1; const int world = g_world->size();
while (i < world) { g_world->isend(i, 1, p); ++i; }
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun() {
RunParams p;
float* __restrict q; size_t qlen; float bx, by, bf;
while (true) {
if (g_world->iprobe(0, 1)) {
g_world->recv(0, 1, p);
AGP_Manip2D(p.nSegments, (bool)p.varLen, p.minTheta, p.tx, p.ty, p.levels, p.maxIter, p.r, (bool)p.adaptive, p.eps, p.seed, &q, &qlen, &bx, &by, &bf);
}
Sleep(0);
}
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpWaitStartAndRun2D() noexcept
{
int dummy;
float* __restrict buf;
size_t len;
while (true) {
if (g_world->iprobe(0, 1)) {
g_world->recv(0, 1, dummy);
AGP_2D(2.0f, 10000.0f, -2.2f, 1.8f, -2.2f, 1.8f, 2.5f, false, 0.00001f, GetTickCount(), &buf, &len);
}
Sleep(0);
}
}

extern "C" __declspec(dllexport) __declspec(noalias) void AgpStartWorkers() {
int i = 1; const int world = g_world->size();
while (i < world) { g_world->isend(i, 1, 0); ++i; }
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* p) { CoTaskMemFree(p); } кусок моего кода сейчас из dll, используется в управляемом коде: // MyForm.cpp — точка входа GUI, без функциональных изменений, просто подключение DLL и запуск формы

#include "MyForm.h"
#include <float.h>

using namespace System;
using namespace System::Windows::Forms;

typedef int(__cdecl* PInit)(int, float, float, float, float);
typedef void(__cdecl* PStartWorkers)();

[STAThread]
int main() {
HMODULE h = LoadLibraryW(L"TEST_FUNC.dll");
auto AgpInit = (PInit)GetProcAddress(h, "AgpInit");
auto AgpWaitStartAndRun = (PStartWorkers)GetProcAddress(h, "AgpWaitStartAndRun");

text
// Инициализация: глубина 12 как и была
const int rank = AgpInit(12, -2.2f, 1.8f, -2.2f, 1.8f);

if (!rank) {
	Application::EnableVisualStyles();
	Application::SetCompatibleTextRenderingDefault(false);
	Application::Run(gcnew TESTAGP::MyForm(h));
}
else {
	AgpWaitStartAndRun();
}
return 0;

}, #pragma once
#define WIN32_LEAN_AND_MEAN
#include <Windows.h>
#include <stdint.h>

using namespace System;
using namespace System::Drawing;
using namespace System::Windows::Forms;
using namespace System::Collections::Generic;
using namespace System::Drawing::Drawing2D;

typedef void(__cdecl* P_MANIP)(int, bool, float, float, float, int, int, float, bool, float, unsigned int, float**, size_t*, float*, float*, float*);
typedef void(__cdecl* P_FREE)(float*);
typedef void(__cdecl* P_START)(int, bool, float, float, float, int, int, float, bool, float, unsigned int);

namespace TESTAGP {
public ref class MyForm : public Form {
public:
MyForm(HMODULE hLib) : hLib(hLib) {
this->SetStyle(ControlStyles::AllPaintingInWmPaint | ControlStyles::UserPaint | ControlStyles::OptimizedDoubleBuffer, true);
this->Text = L"AGP Manipulator 2D"; this->ClientSize = System::Drawing::Size(1000, 700);
this->Resize += gcnew EventHandler(this, &MyForm::OnResize);
fManip = (P_MANIP)GetProcAddress(hLib, "AGP_Manip2D");
pFree = (P_FREE)GetProcAddress(hLib, "AGP_Free");
pStart = (P_START)GetProcAddress(hLib, "AgpStartManipND");
angles = gcnew List<float>(0);
lengths = gcnew List<float>(0);
InitUI();
ResetRandomConfig();
}
protected:
~MyForm() {}
private:
HMODULE hLib; P_MANIP fManip; P_FREE pFree; P_START pStart;
int nSegments; bool variableLengths; List<float>^ angles; List<float>^ lengths;
CheckBox^ cbVarLen; NumericUpDown^ nudMinTheta, ^ nudBaseLength, ^ nudStretchFactor, ^ nudTargetX, ^ nudTargetY, ^ nudLevels, ^ nudMaxIter;
CheckBox^ cbAdaptive; NumericUpDown^ nudR, ^ nudEps; Button^ btnAdd, ^ btnRem, ^ btnOptimize; Label^ lblInfo;
System::UInt32 rngState = 0xA5C39E0Du;

text
	void WireInvalidate(Control^ c) {
		// чтобы после любых изменений интерфейса можно было сразу жать «Оптимизировать»
		if (dynamic_cast<NumericUpDown^>(c)) ((NumericUpDown^)c)->ValueChanged += gcnew EventHandler(this, &MyForm::OnAnyChanged);
		else if (dynamic_cast<CheckBox^>(c)) ((CheckBox^)c)->CheckedChanged += gcnew EventHandler(this, &MyForm::OnAnyChanged);
	}

	void InitUI() {
		int y = 10, w = 180, h = 24, pad = 8, currentX = 10;
		Label^ L;

		// Мин угол (рад) — по умолчанию 1.5 рад (было 0.5)
		L = gcnew Label(); L->Text = L"Мин. угол (рад)"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudMinTheta = gcnew NumericUpDown(); nudMinTheta->Location = Point(currentX, y + h + 2); nudMinTheta->Width = w;
		nudMinTheta->DecimalPlaces = 3; nudMinTheta->Minimum = (Decimal)0.01; nudMinTheta->Maximum = (Decimal)3.14159; nudMinTheta->Value = (Decimal)1.5; this->Controls->Add(nudMinTheta);
		WireInvalidate(nudMinTheta);

		// Базовая длина — диапазон [0.5; 2.0]
		currentX += w + 20; L = gcnew Label(); L->Text = L"Базовая длина"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudBaseLength = gcnew NumericUpDown(); nudBaseLength->Location = Point(currentX, y + h + 2); nudBaseLength->Width = w;
		nudBaseLength->DecimalPlaces = 2; nudBaseLength->Minimum = (Decimal)0.50; nudBaseLength->Maximum = (Decimal)2.00; nudBaseLength->Value = (Decimal)1.00; this->Controls->Add(nudBaseLength);
		WireInvalidate(nudBaseLength);

		// Коэф. растяжения — [1.0; 1.5] с мелким шагом
		currentX += w + 20; L = gcnew Label(); L->Text = L"Коэф. растяжения"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudStretchFactor = gcnew NumericUpDown(); nudStretchFactor->Location = Point(currentX, y + h + 2); nudStretchFactor->Width = w;
		nudStretchFactor->DecimalPlaces = 2; nudStretchFactor->Minimum = (Decimal)1.00; nudStretchFactor->Maximum = (Decimal)1.50;
		nudStretchFactor->Increment = (Decimal)0.01; nudStretchFactor->Value = (Decimal)1.50; this->Controls->Add(nudStretchFactor);
		WireInvalidate(nudStretchFactor);

		// Переменные длины (флаг)
		currentX += w + 20; cbVarLen = gcnew CheckBox();
		cbVarLen->Text = L"Переменные длины"; cbVarLen->Location = Point(currentX, y + h + 2);
		cbVarLen->Width = w; cbVarLen->Checked = false; this->Controls->Add(cbVarLen);
		WireInvalidate(cbVarLen);

		// Следующая строка
		currentX = 10; y += h * 2 + pad + 10;

		// Цель X/Y
		L = gcnew Label(); L->Text = L"Цель X"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudTargetX = gcnew NumericUpDown(); nudTargetX->Location = Point(currentX, y + h + 2); nudTargetX->Width = w;
		nudTargetX->DecimalPlaces = 2; nudTargetX->Minimum = (Decimal)-10.0; nudTargetX->Maximum = (Decimal)10.0; nudTargetX->Value = (Decimal)3.5; this->Controls->Add(nudTargetX);
		WireInvalidate(nudTargetX);

		currentX += w + 20; L = gcnew Label(); L->Text = L"Цель Y"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudTargetY = gcnew NumericUpDown(); nudTargetY->Location = Point(currentX, y + h + 2); nudTargetY->Width = w;
		nudTargetY->DecimalPlaces = 2; nudTargetY->Minimum = (Decimal)-10.0; nudTargetY->Maximum = (Decimal)10.0; nudTargetY->Value = (Decimal)1.0; this->Controls->Add(nudTargetY);
		WireInvalidate(nudTargetY);

		// Глубина — минимум 7, максимум 20 (как просили)
		currentX += w + 20; L = gcnew Label(); L->Text = L"Глубина"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudLevels = gcnew NumericUpDown(); nudLevels->Location = Point(currentX, y + h + 2); nudLevels->Width = w;
		nudLevels->Minimum = 7; nudLevels->Maximum = 20; nudLevels->Value = 12; this->Controls->Add(nudLevels);
		WireInvalidate(nudLevels);

		// Надёжность r — [1.0; 20.0], по умолчанию 2.5
		currentX += w + 20; L = gcnew Label(); L->Text = L"Надежность (r)"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudR = gcnew NumericUpDown(); nudR->Location = Point(currentX, y + h + 2); nudR->Width = w;
		nudR->DecimalPlaces = 2; nudR->Minimum = (Decimal)1.00; nudR->Maximum = (Decimal)20.00; nudR->Value = (Decimal)2.50; this->Controls->Add(nudR);
		WireInvalidate(nudR);

		// Адаптивная схема (флаг)
		currentX += w + 20; cbAdaptive = gcnew CheckBox();
		cbAdaptive->Text = L"Адаптивная"; cbAdaptive->Location = Point(currentX, y + h + 2);
		cbAdaptive->Width = w; cbAdaptive->Checked = true; this->Controls->Add(cbAdaptive);
		WireInvalidate(cbAdaptive);

		// Следующая строка
		y += h * 2 + pad + 10; currentX = 10;

		// Точность — по умолчанию 0,00001; минимум 1e-9; максимум 0.1
		L = gcnew Label(); L->Text = L"Точность"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudEps = gcnew NumericUpDown(); nudEps->Location = Point(currentX, y + h + 2); nudEps->Width = w;
		// ВАЖНО: ставим 9 знаков, чтобы можно было выставить 1e-9; значение по умолчанию 0.00001
		nudEps->DecimalPlaces = 9; nudEps->Minimum = (Decimal)0.000000001; nudEps->Maximum = (Decimal)0.1; nudEps->Value = (Decimal)0.00001;
		// Примечание: WinForms будет показывать с 9 знаками; это нормально функционально.
		this->Controls->Add(nudEps);
		WireInvalidate(nudEps);

		// Макс. итераций — по умолчанию оставить 1000, но максимум поднять до 500000
		currentX += w + 20; L = gcnew Label(); L->Text = L"Макс. итераций"; L->Location = Point(currentX, y); L->Width = w; this->Controls->Add(L);
		nudMaxIter = gcnew NumericUpDown(); nudMaxIter->Location = Point(currentX, y + h + 2); nudMaxIter->Width = w;
		nudMaxIter->Minimum = 10; nudMaxIter->Maximum = 500000; nudMaxIter->Value = 1000; this->Controls->Add(nudMaxIter);
		WireInvalidate(nudMaxIter);

		// Кнопки и инфо
		currentX += 200; btnAdd = gcnew Button(); btnAdd->Text = L"+ Звено"; btnAdd->Location = Point(currentX, y + h + 2); btnAdd->Width = 80;
		btnAdd->Click += gcnew EventHandler(this, &MyForm::OnAddClick); this->Controls->Add(btnAdd);

		currentX += 85; btnRem = gcnew Button(); btnRem->Text = L"- Звено"; btnRem->Location = Point(currentX, y + h + 2); btnRem->Width = 80;
		btnRem->Click += gcnew EventHandler(this, &MyForm::OnRemClick); this->Controls->Add(btnRem);

		currentX += 125; btnOptimize = gcnew Button(); btnOptimize->Text = L"Оптимизировать"; btnOptimize->Location = Point(currentX, y + h + 2); btnOptimize->Width = 120;
		btnOptimize->Click += gcnew EventHandler(this, &MyForm::OnOptimizeClick); this->Controls->Add(btnOptimize);

		currentX += 125; lblInfo = gcnew Label(); lblInfo->Location = Point(currentX, y); lblInfo->Size = System::Drawing::Size(250, 60);
		lblInfo->BorderStyle = BorderStyle::FixedSingle; lblInfo->Text = L"Готов"; this->Controls->Add(lblInfo);
	}

	void ResetRandomConfig() {
		nSegments = 1;
		angles->Clear(); lengths->Clear();
		// По умолчанию первый сегмент смотрит строго вверх (π/2)
		angles->Add((float)(3.14159265358979323846 / 2.0));
		lengths->Add((float)nudBaseLength->Value);
		variableLengths = false;
		this->Invalidate();
	}

	float Rand01() { rngState ^= rngState << 13; rngState ^= rngState >> 17; rngState ^= rngState << 5; return (float)((double)(unsigned)rngState / 4294967296.0); }
	float RandAngle() { return (Rand01() * 6.28318530717958647692f) - 3.14159265358979323846f; }

	System::Void OnResize(System::Object^, System::EventArgs^) { this->Invalidate(); }
	System::Void OnAnyChanged(System::Object^, System::EventArgs^) { this->Invalidate(); }

	System::Void OnAddClick(System::Object^, System::EventArgs^) {
		++nSegments;
		angles->Add(RandAngle());
		lengths->Add((float)nudBaseLength->Value);
		this->Invalidate();
	}
	System::Void OnRemClick(System::Object^, System::EventArgs^) {
		if (nSegments > 1) { --nSegments; angles->RemoveAt(angles->Count - 1); lengths->RemoveAt(lengths->Count - 1); this->Invalidate(); }
	}

	System::Void OnOptimizeClick(System::Object^, System::EventArgs^) {
		variableLengths = cbVarLen->Checked;
		float minTheta = (float)nudMinTheta->Value;
		float tx = (float)nudTargetX->Value, ty = (float)nudTargetY->Value;
		int levels = (int)nudLevels->Value, maxIter = (int)nudMaxIter->Value;
		bool adaptive = cbAdaptive->Checked;
		float r_param = (float)nudR->Value, eps = (float)nudEps->Value;
		unsigned int seed = (unsigned)GetTickCount();

		pStart(nSegments, variableLengths, minTheta, tx, ty, levels, maxIter, r_param, adaptive, eps, seed);

		LARGE_INTEGER t0, t1, fq; QueryPerformanceCounter(&t0);
		float* bestQ = nullptr; size_t bestQLen = 0; float bestX = 0, bestY = 0, bestF = 0;
		fManip(nSegments, variableLengths, minTheta, tx, ty, levels, maxIter, r_param, adaptive, eps, seed,
			&bestQ, &bestQLen, &bestX, &bestY, &bestF);
		QueryPerformanceCounter(&t1); QueryPerformanceFrequency(&fq);
		double micros = 1e6 * (double)(t1.QuadPart - t0.QuadPart) / (double)fq.QuadPart;

		angles->Clear(); for (int i = 0; i < nSegments; ++i) angles->Add(bestQ[i]);
		if (variableLengths) { lengths->Clear(); for (int i = 0; i < nSegments; ++i) lengths->Add(bestQ[nSegments + i]); }
		else { lengths->Clear(); for (int i = 0; i < nSegments; ++i) lengths->Add((float)nudBaseLength->Value); }
		if (pFree) pFree(bestQ);
		lblInfo->Text = String::Format(L"Результат:\nЦель: {0:F5}\nТочка: ({1:F3}, {2:F3})\nВремя: {3:F0} мкс", bestF, bestX, bestY, micros);
		this->Invalidate();
	}

protected:
	virtual void OnPaint(PaintEventArgs^ e) override {
		Form::OnPaint(e); Graphics^ g = e->Graphics; g->SmoothingMode = System::Drawing::Drawing2D::SmoothingMode::HighQuality; g->Clear(this->BackColor);
		int topOffset = 150; System::Drawing::Rectangle drawArea = System::Drawing::Rectangle(0, topOffset, this->ClientSize.Width, this->ClientSize.Height - topOffset);
		g->FillRectangle(Brushes::White, drawArea); int leftWallX = drawArea.Left + this->ClientSize.Width * 25 / 100;
		Pen^ wallPen = gcnew Pen(Color::Black, 2); g->DrawLine(wallPen, leftWallX, drawArea.Top, leftWallX, drawArea.Bottom);
		HatchBrush^ hatchBrush = gcnew HatchBrush(HatchStyle::BackwardDiagonal, Color::LightGray, Color::White); int leftHatchWidth = 100;
		g->FillRectangle(hatchBrush, leftWallX - leftHatchWidth, drawArea.Top, leftHatchWidth, drawArea.Height);
		float targetX = (float)nudTargetX->Value; float targetY = (float)nudTargetY->Value; float scale = 160.0f; int baseX = leftWallX; int baseY = drawArea.Top + drawArea.Height / 2;
		float pixelTargetX = baseX + targetX * scale; float pixelTargetY = baseY - targetY * scale; int rightWallX = (int)(pixelTargetX + 8); rightWallX = Math::Min(rightWallX, drawArea.Right - 10);
		Pen^ dashedPen = gcnew Pen(Color::Black, 2); dashedPen->DashStyle = DashStyle::Dash; g->DrawLine(dashedPen, rightWallX, drawArea.Top, rightWallX, drawArea.Bottom);
		int rightHatchWidth = leftHatchWidth; g->FillRectangle(hatchBrush, rightWallX, drawArea.Top, rightHatchWidth, drawArea.Height);
		Pen^ targetPen = gcnew Pen(Color::Green, 1.5f); targetPen->DashStyle = DashStyle::Dot; g->DrawEllipse(targetPen, pixelTargetX - 8.0f, pixelTargetY - 8.0f, 16.0f, 16.0f);
		cli::array<PointF>^ pts = gcnew cli::array<PointF>(nSegments + 1); pts[0] = PointF((float)baseX, (float)baseY); float x = 0.0f, y = 0.0f, phi = 0.0f;
		array<float>^ anglesArray = angles->ToArray(); array<float>^ lengthsArray = lengths->ToArray();
		for (int i = 0; i < nSegments; ++i) {
			float theta = anglesArray[i]; float L = lengthsArray[i]; phi += theta; x += L * (float)Math::Cos((double)phi); y += L * (float)Math::Sin((double)phi);
			pts[i + 1] = PointF(baseX + x * scale, baseY - y * scale);
		}
		Pen^ penRod = gcnew Pen(Color::Red, 6.0f); for (int i = 0; i < nSegments; ++i) g->DrawLine(penRod, pts[i], pts[i + 1]);
		SolidBrush^ brJoint = gcnew SolidBrush(Color::Blue); const float R = 6.0f; for (int i = 0; i <= nSegments; ++i) g->FillEllipse(brJoint, pts[i].X - R, pts[i].Y - R, 2 * R, 2 * R);
		delete wallPen; delete dashedPen; delete targetPen; delete penRod; delete brJoint; delete hatchBrush;
	}
};

} - проблема в том, что получаю дэдлок, зависание или вылет: C:\Users\maxim>mpiexec -n 6 D:\TEST_AGP\x64\Release\TEST_AGP.exe

Необработанное исключение:

Необработанное исключение: System.Runtime.InteropServices.SEHException: Внешний компонент создал исключение.
в main() в D:\TEST_AGP\MyForm.cpp:строка 27
в mainCRTStartup()
System.Runtime.InteropServices.SEHException: Внешний компонент создал исключение.
в main() в D:\TEST_AGP\MyForm.cpp:строка 27
в mainCRTStartup()

job aborted:
[ranks] message

[0-1] terminated

[2] process exited without calling finalize

[3-5] terminated

---- error analysis -----

[2] on DESKTOP-C1VT6JF
D:\TEST_AGP\x64\Release\TEST_AGP.exe ended prematurely and may have crashed. exit code 0xc0000409

---- error analysis ----- - точно знаю что дело в неправильно логике MPI коммуникаций, или в слишком высокой нагрузке на MPI сеть, когда попросил тебя сделать надёжно работающее решение:
/// код до этого момента без изменений
enum : int {
TAG_STAGE_INT_BASE = 100,
TAG_STAGE_INT_BACK = 140,
TAG_STAGE_BEST_BASE = 180,
TAG_STAGE_BEST_BACK = 220
};

static inline int levels_for_world(int world) noexcept {
if (world <= 1) return 0;
unsigned v = (unsigned)(world - 1);
int L = 0; while (v) { ++L; v >>= 1; }
return L;
}

static inline bool stage_partner_xor(int level, int rank, int world,
int& partner, bool& am_sender, bool& is_last, int L) noexcept {
const int step = 1 << level;
partner = rank ^ step;
if (partner >= world) return false;
const int block = rank & ~((step << 1) - 1);
if (partner < block || partner >= block + (step << 1)) return false;
am_sender = ((rank & step) != 0);
is_last = (level == (L - 1));
return true;
}

static inline void fill_top_intervals_heap(std::vector<IntervalND*>& H, float out5x3[15], uint8_t& outCount) {
if (H.empty()) { outCount = 0; return; }
outCount = (uint8_t)((H.size() >= 3) ? 3 : H.size());
IntervalND* tops[3] = { H[0], (H.size() > 1 ? H[1] : H[0]), (H.size() > 2 ? H[2] : H.back()) };
float* dst = out5x3;
for (uint8_t i = 0; i < outCount; ++i) {
IntervalND* T = tops[i];
dst[0] = T->x1; dst[1] = 0.0f; dst[2] = T->x2; dst[3] = 0.0f; dst[4] = T->R;
dst += 5;
}
}

static inline void inject_received_intervals_p(
const MultiCrossMsg& mX,
const MortonND& map, const ManipCost& cost,
float r, float p, bool stagnation,
float& Mmax, float(&M_by_span)[12], std::vector<IntervalND*>& H,
const std::function<void(IntervalND*)>& update_pockets_and_Mmax,
const std::function<uint64_t(float)>& t_to_idx)
{
const uint8_t cnt = mX.count;
for (uint8_t ii = 0; ii < cnt; ++ii) {
const float* d = &mX.intervals[ii * 5];
float sx = d[0], ex = d[2];
if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f;
if (ex <= sx) continue;

text
    alignas(64) float tmp[32]; float tx, ty;
    map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty);
    map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty);
    IntervalND* inj = new IntervalND(sx, ex, y1i, y2i);
    inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex);
    inj->diam = map.block_diameter(inj->i1, inj->i2);
    inj->compute_span_level(map);

    inj->set_metric(inj->diam);

    if (update_pockets_and_Mmax) update_pockets_and_Mmax(inj);
    inj->ChangeCharacteristic(r * Mmax);

    if (!H.empty()) {
        IntervalND* topH = H.front();
        if (inj->R > 1.15f * topH->R) {
            const float k = stagnation
                ? fmaf(0.5819767068693265f, expm1f(p), 0.3f)
                : fmaf(0.3491860241215959f, expm1f(p), 0.6f);
            inj->R = d[4] * k;
        }
    }

    H.emplace_back(inj);
    std::push_heap(H.begin(), H.end(), ComparePtrND);
}

}

static inline void drain_stage_mailbox(
const MortonND& map, const ManipCost& cost,
float r, float p, bool stagnation,
float& Mmax, float(&M_by_span)[12], std::vector<IntervalND*>& H,
const std::function<void(IntervalND*)>& update_pockets_and_Mmax,
const std::function<uint64_t(float)>& t_to_idx)
{
using namespace boost;
while (g_world->iprobe(mpi::any_source, mpi::any_tag)) {
mpi::status st = g_world->probe(mpi::any_source, mpi::any_tag);
const int tag = st.tag();

text
    if (tag == 0) return;

    if (tag >= TAG_STAGE_INT_BASE && tag < TAG_STAGE_INT_BASE + 64) {
        CtrlMsgND in; g_world->recv(st.source(), tag, in);
        inject_received_intervals_p(in.multiXchg, map, cost, r, p, stagnation,
            Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx);
        continue;
    }

    if (tag >= TAG_STAGE_BEST_BASE && tag < TAG_STAGE_BEST_BASE + 64) {
        BestSolutionMsg sink; g_world->recv(st.source(), tag, sink);
        continue;
    }

    return;
}

}

static __forceinline void agp_run_branch_mpi(
const MortonND& map, const ManipCost& cost, int maxIter, float r, bool adaptive, float eps, unsigned seed,
std::vector<IntervalND*>& H, std::vector<float>& bestQ, float& bestF, float& bestX, float& bestY, float M_prior = 1e-3f)
{
const int n = cost.n;
const int dim = n + (cost.variableLen ? n : 0);

text
alignas(32) float M_by_span[12]; for (int i = 0; i < 12; ++i) M_by_span[i] = M_prior;
float Mmax = M_prior;

alignas(64) float q_local[32], phi[32], s_arr[32], c_arr[32], sum_s[32], sum_c[32], q_try[32];
bestQ.reserve(dim);
float x = 0.0f, y = 0.0f;
int no_improve = 0;

auto t_to_idx = [&](float t) -> uint64_t {
    float tt = (t <= 0.0f) ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t);
    uint64_t idx = (uint64_t)((double)tt * (double)map.scale);
    if (idx >= map.scale) idx = map.scale - 1ull;
    return idx;
    };

auto update_pockets_and_Mmax = [&](IntervalND* I) {
    const int k = I->span_level; if (I->M > M_by_span[k]) M_by_span[k] = I->M; if (M_by_span[k] > Mmax) Mmax = M_by_span[k];
    };

float a = 0.0f, b = 1.0f;

auto evalAt = [&](float t) -> float {
    map.map01ToPoint(t, q_local);
    float f = cost(q_local, x, y);

    if (f < bestF * 1.25f) {
        float acc = 0.0f; for (int i = 0; i < n; ++i) { acc += q_local[i]; phi[i] = acc; }
        FABE13_SINCOS(phi, s_arr, c_arr, n);
        float as = 0.0f, ac = 0.0f;
        for (int k = n - 1; k >= 0; --k) {
            const float Lk = cost.variableLen ? q_local[n + k] : 1.0f;
            as += Lk * s_arr[k]; ac += Lk * c_arr[k]; sum_s[k] = as; sum_c[k] = ac;
        }
        const float dx = x - cost.targetX, dy = y - cost.targetY;
        float dist = sqrtf(fmaf(dx, dx, dy * dy)) + 1e-8f;

        float eta = 0.125f;
        for (int stepI = 0; stepI < 3; ++stepI) {
            for (int i = 0; i < n; ++i) {
                float gpen = 0.0f;
                {
                    float ai = fabsf(q_local[i]);
                    float v = cost.minTheta - ai;
                    if (v > 0.0f) {
                        float scale = 2.0f / (cost.minTheta + 1e-6f);
                        float e = exp2f(scale * v);
                        float dpen_dtheta = cost.sharpW * (e * 0.69314718055994530941723212145818f * scale) * (-copysignf(1.0f, q_local[i]));
                        gpen += dpen_dtheta;
                    }
                }
                {
                    float tsg = -q_local[i] * cost.archBiasK;
                    float sig = 1.0f / (1.0f + expf(-tsg));
                    gpen += -(cost.archBiasW * cost.archBiasK) * sig;
                }

                float g = (dx * (-sum_s[i]) + dy * (sum_c[i])) / dist + gpen;
                q_try[i] = q_local[i] - eta * g;

                const float deg2rad = 3.14159265358979323846f / 180.0f;
                const float lo0 = -60.0f * deg2rad, hi0 = 150.0f * deg2rad;
                const float lo = -150.0f * deg2rad, hi = 150.0f * deg2rad;
                const float Lb = (i == 0) ? lo0 : lo;
                const float Hb = (i == 0) ? hi0 : hi;
                if (q_try[i] < Lb) q_try[i] = Lb; else if (q_try[i] > Hb) q_try[i] = Hb;
            }
            if (cost.variableLen) {
                for (int i = 0; i < n; ++i) {
                    float g = (dx * c_arr[i] + dy * s_arr[i]) / dist;
                    float v = q_local[n + i] - eta * g;
                    if (v < 0.5f) v = 0.5f; else if (v > 2.0f) v = 2.0f;
                    q_try[n + i] = v;
                }
            }
            float x2, y2; float f2 = cost(q_try, x2, y2);
            if (f2 < f) { memcpy(q_local, q_try, dim * sizeof(float)); f = f2; x = x2; y = y2; break; }
            eta *= 0.5f;
        }

        const int last = n - 1;
        const float deg2rad = 3.14159265358979323846f / 180.0f;
        const float lo = (last == 0) ? (-60.0f * deg2rad) : (-150.0f * deg2rad);
        const float hi = 150.0f * deg2rad;
        float bestLocF = f;
        float saved = q_local[last];
        for (float delta = 0.05f; delta >= 0.00625f; delta *= 0.5f) {
            for (int sgn = -1; sgn <= 1; sgn += 2) {
                float cand = saved + sgn * delta;
                if (cand < lo) cand = lo; else if (cand > hi) cand = hi;
                float backup = q_local[last]; q_local[last] = cand;
                float x2, y2; float f2 = cost(q_local, x2, y2);
                if (f2 < bestLocF) { bestLocF = f2; x = x2; y = y2; saved = cand; }
                q_local[last] = backup;
            }
        }
        if (bestLocF < f) { q_local[last] = saved; f = bestLocF; }
    }

    if (f < bestF) { bestF = f; bestQ.assign(q_local, q_local + dim); bestX = x; bestY = y; no_improve = 0; }
    else { ++no_improve; }
    return f;
    };

float f_a = evalAt(a), f_b = evalAt(b);
const int K = (std::min)((std::max)(2 * dim, 8), 128);

H.reserve((size_t)maxIter + K + 16);
const int rank = g_world->rank();
const int world = g_world->size();

alignas(64) float seeds[256 * 32];
const int seedCnt = generate_heuristic_seeds(cost, map, dim, seeds, 32, seed + rank * 7919u);

for (int i = 0; i < seedCnt; ++i) {
    const float* s = seeds + i * 32;
    float t_seed = map.pointToT(s);
    float interval_size = (i < 3) ? (0.0004f * (float)dim)
        : (0.00031f * (float)dim) * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(0.00025f / 0.00031f) * (float)(i - 3));
    float t1 = fmaxf(a, t_seed - interval_size), t2 = fminf(b, t_seed + interval_size);
    if (t2 <= t1) continue;
    alignas(64) float q1[32], q2[32]; float x1, y1, x2, y2;
    map.map01ToPoint(t1, q1); float f1 = cost(q1, x1, y1);
    map.map01ToPoint(t2, q2); float f2 = cost(q2, x2, y2);
    IntervalND* I = new IntervalND(t1, t2, f1, f2);
    I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2);
    I->diam = map.block_diameter(I->i1, I->i2);
    I->compute_span_level(map);
    I->set_metric(I->diam);
    update_pockets_and_Mmax(I);
    I->ChangeCharacteristic(r * Mmax);
    if (i < 3) I->R *= fmaf(0.01f, (float)dim, 0.85f);
    else {
        float start_mult = 0.214f * (float)dim;
        float end_mult = 0.174f * (float)dim;
        float mult = start_mult * exp2f((1.0f / (float)(seedCnt - 4)) * log2f(end_mult / start_mult) * (float)(i - 3));
        I->R *= mult;
    }
    H.emplace_back(I);
    std::push_heap(H.begin(), H.end(), ComparePtrND);
    if (f1 < bestF) { bestF = f1; bestQ.assign(q1, q1 + dim); bestX = x1; bestY = y1; }
    if (f2 < bestF) { bestF = f2; bestQ.assign(q2, q2 + dim); bestX = x2; bestY = y2; }
}

float prev_t = a, prev_f = f_a;
for (int k = 1; k <= K; ++k) {
    float t = a + (b - a) * ((float)k / (K + 1)) + (float)rank / (float)(world * (K + 1));
    float f = evalAt(t);
    IntervalND* I = new IntervalND(prev_t, t, prev_f, f);
    I->i1 = t_to_idx(prev_t); I->i2 = t_to_idx(t);
    I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map); I->set_metric(I->diam);
    update_pockets_and_Mmax(I);
    I->ChangeCharacteristic(r * Mmax);
    H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND);
    prev_t = t; prev_f = f;
}
IntervalND* tail = new IntervalND(prev_t, b, prev_f, f_b);
tail->i1 = t_to_idx(prev_t); tail->i2 = t_to_idx(b);
tail->diam = map.block_diameter(tail->i1, tail->i2); tail->compute_span_level(map); tail->set_metric(tail->diam);
update_pockets_and_Mmax(tail);
tail->ChangeCharacteristic(r * Mmax);
H.emplace_back(tail); std::push_heap(H.begin(), H.end(), ComparePtrND);

float dmax = b - a, initial_len = dmax, thr03 = 0.3f * initial_len, inv_thr03 = 1.0f / thr03;
int it = 0;

auto kickEveryByDim = [&](int dim) -> int {
    float z = 120.0f * exp2f(-0.05f * (float)dim);
    if (z < 60.0f) z = 60.0f;
    return (int)z;
    };

auto noImproveThrByDim = [&](int dim) -> int {
    float z = 80.0f * exp2f(-0.08f * (float)dim);
    if (z < 30.0f) z = 30.0f;
    return (int)z;
    };

while (it < maxIter) {
    if ((it % kickEveryByDim(dim)) == 0 && no_improve > noImproveThrByDim(dim)) {
        float t_best = map.pointToT(bestQ.data());
        for (int i = 0; i < 2; ++i) {
            float off = (i == 0) ? 0.01f : -0.01f;
            float t_seed = fminf(b, fmaxf(a, t_best + off));
            float f_seed = evalAt(t_seed);
            IntervalND* J = new IntervalND(t_seed - 0.005f, t_seed + 0.005f, f_seed, f_seed);
            J->i1 = t_to_idx(t_seed - 0.005f); J->i2 = t_to_idx(t_seed + 0.005f);
            J->diam = map.block_diameter(J->i1, J->i2); J->compute_span_level(map); J->set_metric(J->diam);
            update_pockets_and_Mmax(J);
            J->ChangeCharacteristic(r * Mmax); J->R *= 0.9f;
            H.emplace_back(J); std::push_heap(H.begin(), H.end(), ComparePtrND);
        }
        no_improve = 0;
    }

    const float p = fmaf(-1.0f / initial_len, dmax, 1.0f);
    const bool   stagnation = (no_improve > 100) && (it > 270);

    float A = 200.0f + 64.0f * exp2f(-0.06f * (float)dim);
    float B = 210.0f + 67.0f * exp2f(-0.06f * (float)dim);
    const int T = (int)fmaf(-expm1f(p), A, B);

    float  r_eff = fmaxf(1.0f, r * (0.7f + 0.3f * (1.0f - p)));

    std::pop_heap(H.begin(), H.end(), ComparePtrND);
    IntervalND* cur = H.back(); H.pop_back();

    const float x1 = cur->x1, x2 = cur->x2, y1 = cur->y1, y2 = cur->y2;
    float m = r_eff * Mmax;
    float tNew = Shag(m, x1, x2, y1, y2, (float)dim, r);
    tNew = fminf(fmaxf(tNew, a), b);
    float fNew = evalAt(tNew);

    IntervalND* L = new IntervalND(x1, tNew, y1, fNew);
    IntervalND* Rv = new IntervalND(tNew, x2, fNew, y2);

    L->i1 = t_to_idx(x1); L->i2 = t_to_idx(tNew);
    Rv->i1 = t_to_idx(tNew); Rv->i2 = t_to_idx(x2);
    L->diam = map.block_diameter(L->i1, L->i2); Rv->diam = map.block_diameter(Rv->i1, Rv->i2);
    L->compute_span_level(map); Rv->compute_span_level(map);
    L->set_metric(L->diam); Rv->set_metric(Rv->diam);

    float Mloc = (std::max)(L->M, Rv->M);
    update_pockets_and_Mmax(L); update_pockets_and_Mmax(Rv);

    const float prevMmax = Mmax; if (Mloc > Mmax) Mmax = Mloc;
    m = r_eff * Mmax;

    if (adaptive) {
        float len1 = tNew - x1, len2 = x2 - tNew;
        if (len1 + len2 == dmax) {
            dmax = (std::max)(len1, len2);
            for (auto pI : H) { float Ls = pI->x2 - pI->x1; if (Ls > dmax) dmax = Ls; }
        }
        if ((thr03 > dmax && !(it % 3)) || (10.0f * dmax < initial_len)) {
            const float progress = fmaf(-dmax, inv_thr03, 1.0f);
            const float alpha = progress * progress;
            const float beta = fmaf(-alpha, 1.0f, 2.0f);
            const float MULT = (1.0f / dmax) * Mmax;
            const float global_coeff = fmaf(MULT, r_eff, -MULT);
            const float GF = fmaf(beta, global_coeff, 0.0f);
            L->ChangeCharacteristic(fmaf(GF, len1, L->M * alpha));
            Rv->ChangeCharacteristic(fmaf(GF, len2, Rv->M * alpha));
            size_t sz = H.size();
            RecomputeR_AffineM_AVX2_ND(H.data(), sz, GF, alpha);
            std::make_heap(H.begin(), H.end(), ComparePtrND);
        }
        else {
            if (Mloc > prevMmax) {
                L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m);
                if (Mloc > 1.15f * prevMmax) {
                    size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND);
                }
            }
            else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); }
        }
    }
    else {
        if (Mloc > prevMmax) {
            L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m);
            if (Mloc > 1.15f * prevMmax) { size_t sz = H.size(); RecomputeR_ConstM_AVX2_ND(H.data(), sz, m); std::make_heap(H.begin(), H.end(), ComparePtrND); }
        }
        else { L->ChangeCharacteristic(m); Rv->ChangeCharacteristic(m); }
    }

    H.push_back(L);  std::push_heap(H.begin(), H.end(), ComparePtrND);
    H.push_back(Rv); std::push_heap(H.begin(), H.end(), ComparePtrND);
    if (!H.empty()) {
        _mm_prefetch((const char*)H[0], _MM_HINT_T0);
        if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0);
    }

    IntervalND* top = H.front();
    float interval_len = top->x2 - top->x1;

    bool want_term = (exp2f(log2f(interval_len) / (float)dim) < eps) || (it == maxIter - 1);
    if (!(it % T) || want_term) {
        if (want_term) {
            return;
        }

        drain_stage_mailbox(map, cost, r, p, stagnation, Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx);

        const int L = levels_for_world(world);
        int curLevel = (int)floorf(p * (float)L);
        if (curLevel < 0) curLevel = 0; else if (curLevel >= L) curLevel = L - 1;

        int partner; bool am_sender, is_last;
        if (stage_partner_xor(curLevel, rank, world, partner, am_sender, is_last, L)) {
            CtrlMsgND out; out.kind = 2;
            fill_top_intervals_heap(H, out.multiXchg.intervals, out.multiXchg.count);
            g_world->isend(partner, TAG_STAGE_INT_BASE + curLevel, out);
        }
    }

    if (!(it % 500)) {
        drain_stage_mailbox(map, cost, r, p, stagnation, Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx);

        const int Lb = levels_for_world(world);
        int curLevel = (int)floorf(p * (float)Lb);
        if (curLevel < 0) curLevel = 0; else if (curLevel >= Lb) curLevel = Lb - 1;

        BestSolutionMsg mine;
        mine.bestF = bestF; mine.bestX = bestX; mine.bestY = bestY;
        mine.dim = (uint8_t)bestQ.size();
        if (mine.dim) memcpy(mine.bestQ, bestQ.data(), mine.dim * sizeof(float));

        int partner; bool am_sender, is_last;
        if (stage_partner_xor(curLevel, rank, world, partner, am_sender, is_last, Lb)) {
            g_world->isend(partner, TAG_STAGE_BEST_BASE + curLevel, mine);
        }
    }

    while (g_world->iprobe(boost::mpi::any_source, boost::mpi::any_tag)) {
        auto st = g_world->probe(boost::mpi::any_source, boost::mpi::any_tag);
        const int tag = st.tag();
        if (tag == 0) {
            CtrlMsgND in; g_world->recv(st.source(), 0, in);
            if (in.kind == 0) { return; }
            else if (in.kind == 1) {
                float sx = in.xchg.s_x1, ex = in.xchg.e_x1;
                if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f;
                if (ex > sx) {
                    alignas(64) float tmp[32]; float tx, ty;
                    map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty);
                    map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty);
                    IntervalND* inj = new IntervalND(sx, ex, y1i, y2i);
                    inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex);
                    inj->diam = map.block_diameter(inj->i1, inj->i2);
                    inj->compute_span_level(map); inj->set_metric(inj->diam);
                    update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax);
                    if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); }
                    IntervalND* topH = H.front();
                    if (inj->R > 1.15f * topH->R) {
                        float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f);
                        float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f)
                            : fmaf(0.3491860241215959f, expm1f(p2), 0.6f);
                        inj->R = in.xchg.Rtop * k;
                        H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND);
                    }
                }
            }
            else if (in.kind == 2) {
                const MultiCrossMsg& mX = in.multiXchg;
                for (uint8_t ii = 0; ii < mX.count; ++ii) {
                    const float* d = &mX.intervals[ii * 5];
                    float sx = d[0], ex = d[2];
                    if (sx < 0.0f) sx = 0.0f; if (ex > 1.0f) ex = 1.0f;
                    if (ex > sx) {
                        alignas(64) float tmp[32]; float tx, ty;
                        map.map01ToPoint(sx, tmp); float y1i = cost(tmp, tx, ty);
                        map.map01ToPoint(ex, tmp); float y2i = cost(tmp, tx, ty);
                        IntervalND* inj = new IntervalND(sx, ex, y1i, y2i);
                        inj->i1 = t_to_idx(sx); inj->i2 = t_to_idx(ex);
                        inj->diam = map.block_diameter(inj->i1, inj->i2);
                        inj->compute_span_level(map); inj->set_metric(inj->diam);
                        update_pockets_and_Mmax(inj); inj->ChangeCharacteristic(r * Mmax);
                        if (!H.empty()) { _mm_prefetch((const char*)H[0], _MM_HINT_T0); if (H.size() > 1) _mm_prefetch((const char*)H[1], _MM_HINT_T0); }
                        IntervalND* topH = H.front();
                        if (inj->R > 1.15f * topH->R) {
                            float p2 = fmaf(-1.0f / initial_len, dmax, 1.0f);
                            float k = (no_improve > 100 && it > 270) ? fmaf(0.5819767068693265f, expm1f(p2), 0.3f)
                                : fmaf(0.3491860241215959f, expm1f(p2), 0.6f);
                            inj->R = d[4] * k;
                            H.emplace_back(inj); std::push_heap(H.begin(), H.end(), ComparePtrND);
                        }
                    }
                }
            }
            else if (in.kind == 3) {
                const BestSolutionMsg& bm = in.bestSol;
                if (bm.bestF < bestF * 1.15f) {
                    alignas(64) float tmp_q[32]; memcpy(tmp_q, bm.bestQ, bm.dim * sizeof(float));
                    float t_best = map.pointToT(tmp_q);
                    float t1 = fmaxf(a, t_best - 0.001f), t2 = fminf(b, t_best + 0.001f);
                    if (t2 > t1) {
                        alignas(64) float tq1[32], tq2[32]; float xx1, yy1, xx2, yy2;
                        map.map01ToPoint(t1, tq1); float f1 = cost(tq1, xx1, yy1);
                        map.map01ToPoint(t2, tq2); float f2 = cost(tq2, xx2, yy2);
                        IntervalND* I = new IntervalND(t1, t2, f1, f2);
                        I->i1 = t_to_idx(t1); I->i2 = t_to_idx(t2);
                        I->diam = map.block_diameter(I->i1, I->i2); I->compute_span_level(map);
                        I->set_metric(I->diam); update_pockets_and_Mmax(I);
                        I->ChangeCharacteristic(r * Mmax); I->R *= 0.90f;
                        H.emplace_back(I); std::push_heap(H.begin(), H.end(), ComparePtrND);
                    }
                    if (bm.bestF < bestF) { bestF = bm.bestF; bestX = bm.bestX; bestY = bm.bestY; bestQ.assign(bm.bestQ, bm.bestQ + bm.dim); }
                }
            }
            continue;
        }
        if (tag >= TAG_STAGE_INT_BASE && tag < TAG_STAGE_INT_BASE + 64) {
            CtrlMsgND in; g_world->recv(st.source(), tag, in);
            inject_received_intervals_p(in.multiXchg, map, cost, r, p, stagnation,
                Mmax, M_by_span, H, update_pockets_and_Mmax, t_to_idx);
            continue;
        }

        if (tag >= TAG_STAGE_BEST_BASE && tag < TAG_STAGE_BEST_BASE + 64) {
            BestSolutionMsg in; g_world->recv(st.source(), tag, in);
            if (in.bestF < bestF) {
                bestF = in.bestF; bestX = in.bestX; bestY = in.bestY;
                bestQ.assign(in.bestQ, in.bestQ + in.dim);
            }
            continue;
        }

        break;

    }
    ++it;
}

}

extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_Manip2D(int nSegments, bool variableLengths, float minTheta, float targetX, float targetY,
int peanoLevels, int maxIterPerBranch, float r, bool adaptiveMode, float epsilon,
unsigned int seed, float** out_bestQ, size_t* out_bestQLen, float* out_bestX,
float* out_bestY, float* out_bestF)
{
Slab* const __restrict slab = tls.local(); slab->current = slab->base;
while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
const int dim = nSegments + (variableLengths ? nSegments : 0);

text
g_mc.permCache.resize(dim); for (int i = 0; i < dim; ++i) g_mc.permCache[i] = i;
uint32_t s = g_mc.baseSeed;
for (int i = dim - 1; i > 0; --i) { s ^= s << 13; s ^= s >> 17; s ^= s << 5; uint32_t j = s % (uint32_t)(i + 1); std::swap(g_mc.permCache[i], g_mc.permCache[j]); }
g_mc.invMaskCache.resize(dim); for (int k = 0; k < dim; ++k) { s ^= s << 13; s ^= s >> 17; s ^= s << 5; g_mc.invMaskCache[k] = (uint64_t)s; }

const float deg2rad = 3.14159265358979323846f / 180.0f;
const float theta0Min = -60.0f * deg2rad, theta0Max = 150.0f * deg2rad;
const float thetaMin = -150.0f * deg2rad, thetaMax = 150.0f * deg2rad;
const float lenMin = 0.5f, lenMax = 2.0f;

std::vector<float> low;  low.reserve(dim);
std::vector<float> high; high.reserve(dim);
for (int i = 0; i < nSegments; ++i) { low.push_back(i == 0 ? theta0Min : thetaMin); high.push_back(i == 0 ? theta0Max : thetaMax); }
if (variableLengths) { for (int i = 0; i < nSegments; ++i) { low.push_back(lenMin); high.push_back(lenMax); } }

ManipCost cost(nSegments, variableLengths, targetX, targetY, minTheta);

const int rank = g_world->rank(), world = g_world->size();
std::vector<float> bestQ; float bestF = FLT_MAX, bx = 0.0f, by = 0.0f;

const int levels0 = (std::min)(peanoLevels, 8);
const int maxIter0 = (int)(maxIterPerBranch * 0.2f);
MortonND map0(dim, levels0, low.data(), high.data(), g_mc);

std::vector<IntervalND*> H_coarse; std::vector<float> bestQ_coarse;
float bestF_coarse = FLT_MAX, bx_coarse = 0.0f, by_coarse = 0.0f;

float M_prior = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << levels0));
if (variableLengths) M_prior += 1.41421356237f * (1.0f / (float)(1u << levels0));

agp_run_branch_mpi(map0, cost, maxIter0, r, adaptiveMode, epsilon, seed,
    H_coarse, bestQ_coarse, bestF_coarse, bx_coarse, by_coarse, M_prior);

if (bestF_coarse < bestF) { bestF = bestF_coarse; bestQ = bestQ_coarse; bx = bx_coarse; by = by_coarse; }

if (levels0 < peanoLevels) {
    while (g_world->iprobe(boost::mpi::any_source, 0)) { CtrlMsg dummy; g_world->recv(boost::mpi::any_source, 0, dummy); }
    MortonND map1(dim, peanoLevels, low.data(), high.data(), g_mc);
    std::vector<IntervalND*> H_fine;
    std::vector<float> bestQ_fine = bestQ;
    float bestF_fine = bestF, bx_fine = bx, by_fine = by;

    float M_prior_fine = (variableLengths ? 2.0f * nSegments : 2.0f * nSegments) * (1.0f / (float)(1u << peanoLevels));
    if (variableLengths) M_prior_fine += 1.41421356237f * (1.0f / (float)(1u << peanoLevels));

    std::sort(H_coarse.begin(), H_coarse.end(), [](const IntervalND* a, const IntervalND* b) { return a->R < b->R; });
    const size_t topCount = (size_t)(H_coarse.size() * 0.3f);
    auto t_to_idx_fine = [&](float t) -> uint64_t {
        float tt = (t <= 0.0f) ? 0.0f : (t >= 1.0f ? std::nextafterf(1.0f, 0.0f) : t);
        uint64_t idx = (uint64_t)((double)tt * (double)map1.scale);
        if (idx >= map1.scale) idx = map1.scale - 1ull;
        return idx;
        };
    for (size_t i = 0; i < topCount && i < H_coarse.size(); ++i) {
        const IntervalND* C = H_coarse[i];
        alignas(64) float q1[32], q2[32]; float x1, y1, x2, y2;
        map1.map01ToPoint(C->x1, q1); float f1 = cost(q1, x1, y1);
        map1.map01ToPoint(C->x2, q2); float f2 = cost(q2, x2, y2);
        IntervalND* I = new IntervalND(C->x1, C->x2, f1, f2);
        I->i1 = t_to_idx_fine(C->x1); I->i2 = t_to_idx_fine(C->x2);
        I->diam = map1.block_diameter(I->i1, I->i2);
        I->set_metric(I->diam);
        H_fine.push_back(I);
        if (f1 < bestF_fine) { bestF_fine = f1; bestQ_fine.assign(q1, q1 + dim); bx_fine = x1; by_fine = y1; }
        if (f2 < bestF_fine) { bestF_fine = f2; bestQ_fine.assign(q2, q2 + dim); bx_fine = x2; by_fine = y2; }
    }
    std::make_heap(H_fine.begin(), H_fine.end(), ComparePtrND);
    agp_run_branch_mpi(map1, cost, maxIterPerBranch - maxIter0, r, adaptiveMode, epsilon, seed,
        H_fine, bestQ_fine, bestF_fine, bx_fine, by_fine, M_prior_fine);

    if (bestF_fine < bestF) { bestF = bestF_fine; bestQ = bestQ_fine; bx = bx_fine; by = by_fine; }
}

BestSolutionMsg agg;
agg.bestF = bestF;
agg.bestX = bx;
agg.bestY = by;
agg.dim = (uint8_t)bestQ.size();
if (agg.dim) memcpy(agg.bestQ, bestQ.data(), agg.dim * sizeof(float));

const int L = levels_for_world(world);
for (int lvl = 0; lvl < L; ++lvl) {
    int partner; bool am_sender, is_last;
    if (!stage_partner_xor(lvl, rank, world, partner, am_sender, is_last, L)) continue;

    if (!is_last) {
        if (am_sender) {
            g_world->send(partner, 2 + lvl, agg);
            return;
        }
        else {
            BestSolutionMsg in;
            g_world->recv(partner, 2 + lvl, in);
            if (in.bestF < agg.bestF) agg = in;
        }
    }
    else {
        BestSolutionMsg in;
        g_world->sendrecv(partner, 2 + lvl + 128,
            agg, partner, 2 + lvl + 128, in);
        if (in.bestF < agg.bestF) agg = in;
    }
}

if (rank == 0) {
    *out_bestQLen = (size_t)agg.dim;
    *out_bestQ = (float*)CoTaskMemAlloc(sizeof(float) * (*out_bestQLen));
    memcpy(*out_bestQ, agg.bestQ, sizeof(float) * (*out_bestQLen));
    *out_bestX = agg.bestX;
    *out_bestY = agg.bestY;
    *out_bestF = agg.bestF;
}

}

///остальной код без изменений

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* p) { CoTaskMemFree(p); } - оно работает - но хуже по скорости моей реализации - твоя задача понять почему в моей реализации происходит дэдлок и исправить не меняя моей сути с древовидным обменом - оставиви в основном мою функциональность и логику обмена - правки должны быть минимальны чтобы устранить ошибку но не меняя сам смысл логики, известно что твоя реализация работала стабильно, думай долго и мне верни полный исправленный код основного файла dll, удовлетворяющий моим требованиям по устранению дэдлока, думай долго чтобы ответ удовлетворял моим требованиям