#include "pch.h" #define XOR_RAND(state, result_v...

Question

#include "pch.h"

#define XOR_RAND(state, result_var)
do {
int s = state;
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
state = s;
result_var = state * 0x1.0p-32f;
} while(0)

#define XOR_RAND_GRSH(state, result_var)
do {
int s = state;
s ^= s << 13;
s ^= s >> 17;
s ^= s << 5;
state = s;
result_var = fmaf(state, 0x1.0p-31f, -1.0f);
} while(0)

#define FABE13_COS(x, result_var)
do {
const float abs_val = fabsf(x);
float reduced = fmodf(abs_val, 6.28318530718f);
if(reduced > 3.14159265359f) {
reduced = 6.28318530718f - reduced;
}
if(reduced < 1.57079632679f) {
const float val2 = reduced * reduced;
const float val4 = val2 * val2;
result_var = fmaf(val4, fmaf(val2, -0.0013888889f, 0.0416666667f), fmaf(val2, -0.5f, 1.0f));
} else {
reduced = 3.14159265359f - reduced;
const float val2 = reduced * reduced;
const float val4 = val2 * val2;
result_var = -fmaf(val4, fmaf(val2, -0.0013888889f, 0.0416666667f), fmaf(val2, -0.5f, 1.0f));
}
} while(0)

#define FABE13_SINCOS(in, sin_out, cos_out, n)
do {
int i = 0;
const int limit = n & ~7;
if(n >= 8) {
static __declspec(align(32)) const __m256 VEC_TWOPI = _mm256_set1_ps(6.28318530718f);
static __declspec(align(32)) const __m256 VEC_PI = _mm256_set1_ps(3.14159265359f);
static __declspec(align(32)) const __m256 VEC_PI_2 = _mm256_set1_ps(1.57079632679f);
static __declspec(align(32)) const __m256 INV_TWOPI = _mm256_set1_ps(0.15915494309189535f);
static __declspec(align(32)) const __m256 BIAS = _mm256_set1_ps(12582912.0f);
static __declspec(align(32)) const __m256 VEC_COS_P5 = _mm256_set1_ps(-0.0013888889f);
static __declspec(align(32)) const __m256 VEC_COS_P3 = _mm256_set1_ps(0.0416666667f);
static __declspec(align(32)) const __m256 VEC_COS_P1 = _mm256_set1_ps(-0.5f);
static __declspec(align(32)) const __m256 VEC_COS_P0 = _mm256_set1_ps(1.0f);
static __declspec(align(32)) const __m256 VEC_SIN_P5 = _mm256_set1_ps(-0.0001984127f);
static __declspec(align(32)) const __m256 VEC_SIN_P3 = _mm256_set1_ps(0.0083333333f);
static __declspec(align(32)) const __m256 VEC_SIN_P1 = _mm256_set1_ps(-0.16666666f);
static __declspec(align(32)) const __m256 VEC_SIN_P0 = _mm256_set1_ps(1.0f);
static __declspec(align(32)) const __m256 VEC_ZERO = _mm256_setzero_ps();
while(i < limit) {
const __m256 vx = _mm256_load_ps(&in[i]);
const __m256 vax = _mm256_andnot_ps(_mm256_set1_ps(-0.0f), vx);
__m256 q = _mm256_fmadd_ps(vax, INV_TWOPI, BIAS);
q = _mm256_sub_ps(q, BIAS);
const __m256 r = _mm256_fnmadd_ps(VEC_TWOPI, q, vax);
const __m256 r1 = _mm256_min_ps(r, _mm256_sub_ps(VEC_TWOPI, r));
const __m256 r2 = _mm256_min_ps(r1, _mm256_sub_ps(VEC_PI, r1));
const __m256 t2 = _mm256_mul_ps(r2, r2);
const __m256 cosv = _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_COS_P5, VEC_COS_P3), VEC_COS_P1), VEC_COS_P0);
const __m256 sinv = _mm256_mul_ps(_mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, _mm256_fmadd_ps(t2, VEC_SIN_P5, VEC_SIN_P3), VEC_SIN_P1), VEC_SIN_P0), r2);
const __m256 cflip = _mm256_cmp_ps(r1, VEC_PI_2, _CMP_GT_OQ);
const __m256 sflip = _mm256_xor_ps(_mm256_cmp_ps(vx, VEC_ZERO, _CMP_LT_OQ), _mm256_cmp_ps(r, VEC_PI, _CMP_GT_OQ));
_mm256_store_ps(&cos_out[i], _mm256_blendv_ps(cosv, _mm256_sub_ps(VEC_ZERO, cosv), cflip));
_mm256_store_ps(&sin_out[i], _mm256_blendv_ps(sinv, _mm256_sub_ps(VEC_ZERO, sinv), sflip));
i += 8;
}
}
while(i < n) {
const float x = in[i];
const float ax = fabsf(x);
float q = fmaf(ax, 0.15915494309189535f, 12582912.0f);
q -= 12582912.0f;
float r = fmaf(-6.28318530718f, q, ax);
const bool sflip = r > 3.14159265359f;
if(sflip) r = 6.28318530718f - r;
const bool cflip = r > 1.57079632679f;
if(cflip) r = 3.14159265359f - r;
const float t2 = r * r;
const float c = fmaf(t2, fmaf(t2, fmaf(t2, -0.0013888889f, 0.0416666667f), -0.5f), 1.0f);
const float s = fmaf(t2, fmaf(t2, fmaf(t2, -0.0001984127f, 0.0083333333f), -0.16666666f), 1.0f) * r;
cos_out[i] = cflip ? -c : c;
sin_out[i] = ((x < 0.0f) ^ sflip) ? -s : s;
++i;
}
} while(0)

__declspec(align(16)) struct Slab final {
char* const __restrict base;
char* __restrict current;
char* const __restrict end;

text
__declspec(noalias) __forceinline Slab(void* const __restrict memory, const size_t usable_size) noexcept
	: base(static_cast<char* __restrict>(memory))
	, current(base)
	, end(base + (usable_size & ~static_cast<size_t>(63u)))
{
}

};

static tbb::enumerable_thread_specific<Slab*> tls( noexcept {
void* const __restrict memory = _aligned_malloc(16777216u, 16u);
Slab* const __restrict slab = static_cast<Slab*>(_aligned_malloc(32u, 16u));
new (slab) Slab(memory, 16777216u);

char* __restrict p = slab->base;

#pragma loop(ivdep)
while (p < slab->end) {
*p = 0u;
p += 4096u;
}
return slab;
}());

__declspec(align(64)) struct Interval final {
const float x1;
const float x2;
const float y1;
const float y2;
const float delta_y;
const float ordinate_factor;
const float N_factor;
const float quadratic_term;
const float M;
float R;

text
__declspec(noalias) __forceinline void* operator new(const size_t) noexcept {
	Slab* const __restrict s = tls.local();
	char* const __restrict result = s->current;
	s->current += 64u;
	return result;
}

__declspec(noalias) __forceinline Interval(const float _x1, const float _x2,
	const float _y1, const float _y2, const float _N) noexcept
	: x1(_x1), x2(_x2), y1(_y1), y2(_y2)
	, delta_y(_y2 - _y1)
	, ordinate_factor(-(y1 + y2) * 2.0f)
	, N_factor(_N == 1.0f ? _x2 - _x1 : _N == 2.0f ? sqrtf(_x2 - _x1) : powf(_x2 - _x1, 1.0f / _N))
	, quadratic_term((1.0f / N_factor)* delta_y* delta_y)
	, M((1.0f / N_factor)* fabsf(delta_y))
{
}

__declspec(noalias) __forceinline void ChangeCharacteristic(const float _m) noexcept
{
	R = fmaf(1.0f / _m, quadratic_term, fmaf(_m, N_factor, ordinate_factor));
}

};

enum List : uint8_t { Top = 0b00u, Down = 0b01u, Left = 0b10u, Right = 0b11u };

__declspec(align(4)) struct Step final {
const uint8_t next;
const uint8_t dx;
const uint8_t dy;
};

__declspec(align(4)) struct InvStep final {
const uint8_t q;
const uint8_t next;
};

__declspec(align(64)) static const Step g_step_tbl[4][4] = {
{ {Right,0u,0u}, {Top,0u,1u}, {Top,1u,1u}, {Left,1u,0u} },
{ {Left,1u,1u}, {Down,1u,0u}, {Down,0u,0u}, {Right,0u,1u} },
{ {Down,1u,1u}, {Left,0u,1u}, {Left,0u,0u}, {Top,1u,0u} },
{ {Top,0u,0u}, {Right,1u,0u}, {Right,1u,1u}, {Down,0u,1u} }
};

__declspec(align(64)) static const InvStep g_inv_tbl[4][4] = {
{ {0u,Right}, {1u,Top}, {3u,Left}, {2u,Top} },
{ {2u,Down}, {3u,Right}, {1u,Down}, {0u,Left} },
{ {2u,Left}, {1u,Left}, {3u,Top}, {0u,Down} },
{ {0u,Top}, {3u,Down}, {1u,Right}, {2u,Right} }
};

__declspec(align(16)) struct Peano2DMap final {
const int levels;
const float a, b, c, d;
const float lenx, leny;
const float inv_lenx;
const uint32_t scale;
const uint8_t start;

text
__declspec(noalias) __forceinline Peano2DMap(
	const int L,
	const float _a,
	const float _b,
	const float _c,
	const float _d,
	const uint8_t startType
) noexcept
	: levels(L)
	, a(_a), b(_b), c(_c), d(_d)
	, lenx(_b - _a)
	, leny(_d - _c)
	, inv_lenx(1.0f / (_b - _a))
	, scale(static_cast<uint32_t>(1u) << (L << 1))
	, start(startType)
{
}

};

static Peano2DMap gActiveMap(0, 0.0f, 0.0f, 0.0f, 0.0f, 0b00u);

static const boost::mpi::environment* __restrict g_env;
static const boost::mpi::communicator* __restrict g_world;

__declspec(align(16)) struct CrossMsg final {
float s_x1, s_x2;
float e_x1, e_x2;
float Rtop;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept { ar& s_x1& s_x2& e_x1& e_x2& Rtop; }
};

__declspec(align(16)) struct CtrlMsg final {
bool kind;
CrossMsg xchg;
template<typename Archive>
__declspec(noalias) __forceinline void serialize(Archive& ar, const unsigned int) noexcept {
ar& kind& xchg;
}
};

static __declspec(noalias) __forceinline bool ComparePtr(const Interval* const __restrict a, const Interval* const __restrict b) noexcept
{
return a->R < b->R;
}

static __declspec(noalias) __forceinline float ShekelFunc(const float x, const float seed) noexcept
{
int i = 0;
float current_state = seed, current_res, current_res2, res = 0.0f;
while (i < 10) {
XOR_RAND(current_state, current_res);
const float x_part = fmaf(-current_res, 10.0f, x);
XOR_RAND(current_state, current_res);
XOR_RAND(current_state, current_res2);
float delimiter = fmaf(fmaf(current_res, 20.0f, 5.0f), x_part * x_part, fmaf(current_res2, 0.2f, 1.0f));
delimiter = copysignf(fmaxf(fabsf(delimiter), FLT_MIN), delimiter);
res -= 1.0f / delimiter;
++i;
}
return res;
}

static __declspec(noalias) __forceinline float RastriginFunc(const float x1, const float x2) noexcept
{
const float term1 = fmaf(x1, x1, x2 * x2);
float cos1, cos2;
FABE13_COS(6.28318530717958647692f * x1, cos1);
FABE13_COS(6.28318530717958647692f * x2, cos2);
return (term1 - fmaf(cos1 + cos2, 10.0f, -14.6f)) * fmaf(-term1, 0.25f, 18.42f);
}

static __declspec(noalias) __forceinline float HillFunc(const float x, const float seed) noexcept
{
int j = 0;
__declspec(align(32)) float angles[14u];
const float start_angle = 6.28318530717958647692f * x;
#pragma loop(ivdep)
while (j < 14) {
angles[j] = start_angle * static_cast<float>(j + 1);
++j;
}
__declspec(align(32)) float sin_vals[14u];
__declspec(align(32)) float cos_vals[14u];
FABE13_SINCOS(angles, sin_vals, cos_vals, 14u);
float current_state = seed, current_res, current_res2;
XOR_RAND(current_state, current_res);
float res = fmaf(current_res, 2.0f, -1.1f);
--j;
while (j >= 0) {
XOR_RAND(current_state, current_res);
XOR_RAND(current_state, current_res2);
res += fmaf(fmaf(current_res, 2.0f, -1.1f), sin_vals[j], fmaf(current_res2, 2.0f, -1.1f) * cos_vals[j]);
--j;
}
return res;
}

static __declspec(noalias) __forceinline float GrishaginFunc(const float x1, const float x2, const float seed) noexcept
{
int j = 0;
__declspec(align(32)) float angles_j[8u];
__declspec(align(32)) float angles_k[8u];
#pragma loop(ivdep)
while (j < 8) {
const float pj_mult = 3.14159265358979323846f * static_cast<float>(j + 1);
angles_j[j] = pj_mult * x1;
angles_k[j] = pj_mult * x2;
++j;
}
__declspec(align(32)) float sin_j[8u], cos_j[8u];
__declspec(align(32)) float sin_k[8u], cos_k[8u];
FABE13_SINCOS(angles_j, sin_j, cos_j, 8u);
FABE13_SINCOS(angles_k, sin_k, cos_k, 8u);
--j;
float part1 = 0.0f;
float part2 = 0.0f;
float current_state = seed, current_res, current_res2;
while (j >= 0) {
size_t k = 0u;
while (k < 8u) {
const float sin_term = sin_j[j] * sin_j[j];
const float cos_term = cos_k[k] * cos_k[k];
XOR_RAND_GRSH(current_state, current_res);
XOR_RAND_GRSH(current_state, current_res2);
part1 = fmaf(current_res, sin_term, fmaf(current_res2, cos_term, part1));
XOR_RAND_GRSH(current_state, current_res);
XOR_RAND_GRSH(current_state, current_res2);
part2 = fmaf(-current_res, cos_term, fmaf(current_res2, sin_term, part2));
++k;
}
--j;
}
return -sqrtf(fmaf(part1, part1, part2 * part2));
}

static __declspec(noalias) __forceinline float Shag(const float _m, const float x1, const float x2, const float y1,
const float y2, const float _N, const float _r) noexcept
{
const float diff = y2 - y1;
const float sign_mult = _mm_cvtss_f32(_mm_castsi128_ps(_mm_set1_epi32(
0x3F800000u | ((reinterpret_cast<const uint32_t>(&diff) & 0x80000000u) ^ 0x80000000u))));
return _N == 1.0f
? fmaf(-(1.0f / _m), diff, x1 + x2) * 0.5f
: _N == 2.0f
? fmaf(sign_mult / (_m * _m), diff * diff * _r, x1 + x2) * 0.5f
: fmaf(sign_mult / powf(_m, _N), powf(diff, _N) * _r, x1 + x2) * 0.5f;
}

static __declspec(noalias) __forceinline void RecomputeR_ConstM_AVX2(Interval* const* const __restrict arr, const size_t n, const float m) noexcept {
const __m256 vm = _mm256_set1_ps(m);

text
__m256 vinvm = _mm256_rcp_ps(vm);
vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));

size_t i = 0;
const int limit = static_cast<int>(n & ~7u);
if (n >= 8u) {
	while (i < static_cast<size_t>(limit)) {
		__declspec(align(32)) float q[8], nf[8], ord[8];
		int k = 0;

#pragma loop(ivdep)
while (k < 8) {
const Interval* const __restrict p = arr[i + k];
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
ord[k] = p->ordinate_factor;
++k;
}
const __m256 vq = _mm256_load_ps(q);
const __m256 vnf = _mm256_load_ps(nf);
const __m256 vod = _mm256_load_ps(ord);

text
		const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
		const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);

		__declspec(align(32)) float out[8];
		_mm256_store_ps(out, res);
		k = 0;

#pragma loop(ivdep)
while (k < 8) {
arr[i + k]->R = out[k];
++k;
}
i += 8u;
}
}
while (i < n) {
arr[i]->ChangeCharacteristic(m);
++i;
}
}

static __declspec(noalias) __forceinline void RecomputeR_AffineM_AVX2(Interval* const* const __restrict arr, const size_t n, const float GLOBAL_FACTOR, const float alpha) noexcept {
const __m256 vGF = _mm256_set1_ps(GLOBAL_FACTOR);
const __m256 va = _mm256_set1_ps(alpha);

text
size_t i = 0;
const int limit = static_cast<int>(n & ~7u);
if (n >= 8u) {
	while (i < static_cast<size_t>(limit)) {
		__declspec(align(32)) float len[8], Mv[8], q[8], nf[8], ord[8];
		int k = 0;

#pragma loop(ivdep)
while (k < 8) {
const Interval* const p = arr[i + k];
len[k] = p->x2 - p->x1;
Mv[k] = p->M;
q[k] = p->quadratic_term;
nf[k] = p->N_factor;
ord[k] = p->ordinate_factor;
++k;
}
const __m256 vlen = _mm256_load_ps(len);
const __m256 vM = _mm256_load_ps(Mv);
const __m256 vq = _mm256_load_ps(q);
const __m256 vnf = _mm256_load_ps(nf);
const __m256 vod = _mm256_load_ps(ord);

text
		const __m256 vm = _mm256_fmadd_ps(vGF, vlen, _mm256_mul_ps(va, vM));

		__m256 vinvm = _mm256_rcp_ps(vm);
		vinvm = _mm256_mul_ps(vinvm, _mm256_fnmadd_ps(vm, vinvm, _mm256_set1_ps(2.0f)));

		const __m256 t = _mm256_fmadd_ps(vm, vnf, vod);
		const __m256 res = _mm256_fmadd_ps(vq, vinvm, t);

		__declspec(align(32)) float out[8];
		_mm256_store_ps(out, res);
		k = 0;

#pragma loop(ivdep)
while (k < 8) {
arr[i + k]->R = out[k];
++k;
}
i += 8u;
}
}
while (i < n) {
const Interval* const __restrict p = arr[i];
const float mi = fmaf(GLOBAL_FACTOR, p->x2 - p->x1, p->M * alpha);
arr[i]->R = fmaf(1.0f / mi, p->quadratic_term, fmaf(mi, p->N_factor, p->ordinate_factor));
++i;
}
}

__declspec(noalias) __forceinline void HitTest2D_analytic(const float x_param, float& out_x1, float& out_x2) noexcept
{
const float a = gActiveMap.a;
const float inv_lenx = gActiveMap.inv_lenx;
const uint32_t scale = gActiveMap.scale;
const uint32_t scale_minus_1 = scale - 1u;
const float lenx = gActiveMap.lenx;
const float leny = gActiveMap.leny;
const float c = gActiveMap.c;
const uint8_t start = gActiveMap.start;
const int levels = gActiveMap.levels;

text
float norm = (x_param - a) * inv_lenx;
norm = fminf(fmaxf(norm, 0.0f), 0x1.fffffep-1f);

uint32_t idx = static_cast<uint32_t>(norm * static_cast<float>(scale));
idx = idx > scale_minus_1 ? scale_minus_1 : idx;

float sx = lenx, sy = leny;
float x1 = a, x2 = c;
uint8_t type = start;

int l = levels - 1;

#pragma loop(ivdep)
while (l >= 0) {
const uint32_t q = (idx >> (l * 2)) & 3u;
const Step s = g_step_tbl[type][q];
type = s.next;
sx *= 0.5f; sy *= 0.5f;
x1 += s.dx ? sx : 0.0f;
x2 += s.dy ? sy : 0.0f;
--l;
}
out_x1 = x1 + sx * 0.5f;
out_x2 = x2 + sy * 0.5f;
}

__declspec(noalias) __forceinline float FindX2D_analytic(const float px, const float py) noexcept
{
const float a = gActiveMap.a;
const float b = gActiveMap.b;
const float c = gActiveMap.c;
const float d = gActiveMap.d;
const float lenx = gActiveMap.lenx;
const float leny = gActiveMap.leny;
const uint32_t scale = gActiveMap.scale;
const uint8_t start = gActiveMap.start;
const int levels = gActiveMap.levels;

text
const float clamped_px = fminf(fmaxf(px, a), b);
const float clamped_py = fminf(fmaxf(py, c), d);

float sx = lenx, sy = leny;
float x0 = a, y0 = c;
uint32_t idx = 0u;
uint8_t type = start;

int l = 0;

#pragma loop(ivdep)
while (l < levels) {
sx *= 0.5f; sy *= 0.5f;
const float mx = x0 + sx;
const float my = y0 + sy;

text
	const uint32_t tr = static_cast<uint32_t>((clamped_px > mx) & (clamped_py > my));
	const uint32_t tl = static_cast<uint32_t>((clamped_px < mx) & (clamped_py > my));
	const uint32_t dl = static_cast<uint32_t>((clamped_px < mx) & (clamped_py < my));
	const uint32_t none = static_cast<uint32_t>(1u ^ (tr | tl | dl));

	const uint32_t dd = (tr << 1) | tr | tl | (none << 1);

	const InvStep inv = g_inv_tbl[type][dd];
	type = inv.next;
	idx = (idx << 2) | inv.q;

	const uint32_t dx = dd >> 1;
	const uint32_t dy = dd & 1u;
	x0 += dx ? sx : 0.0f;
	y0 += dy ? sy : 0.0f;
	++l;
}

const float scale_reciprocal = 1.0f / static_cast<float>(scale);
return fmaf(static_cast<float>(idx) * scale_reciprocal, lenx, a);

}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline int AgpInit(const int peanoLevel, const float a, const float b, const float c, const float d) noexcept
{
g_env = new boost::mpi::environment();
g_world = new boost::mpi::communicator();
const int rank = g_world->rank();

text
new(&gActiveMap) Peano2DMap(peanoLevel, a, b, c, d, rank ? static_cast<uint8_t>(Down) : static_cast<uint8_t>(Top));
return rank;

}

extern "C" __declspec(dllexport) __declspec(noalias)
void AGP_1D(const float global_iterations,
const float a, const float b, 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;

text
int schetchick = 0;
const float initial_length = b - a;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float 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;
float x_Rmax_2 = b;
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;

std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __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 Interval(a, b, start_val, best_f, 1.0f));

float Mmax = R.front()->M;
float m = r * Mmax;

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(), ComparePtr);
	const Interval* const __restrict promejutochny_otrezok = R.back();

	const float new_x1 = promejutochny_otrezok->x1;
	const float new_x2 = promejutochny_otrezok->x2;
	const float len2 = new_x2 - new_point;
	const float len1 = new_point - new_x1;
	const float interval_len = len1 < len2 ? len1 : len2;

	if (interval_len < epsilon || ++schetchick == static_cast<int>(global_iterations)) {
		Extr.emplace_back(static_cast<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;
	}

	Interval* const __restrict curr = new Interval(new_x1, new_point, promejutochny_otrezok->y1, new_value, 1.0f);
	Interval* const __restrict curr1 = new Interval(new_point, new_x2, new_value, promejutochny_otrezok->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;
			size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}

text
		if (threshold_03 > dmax && schetchick % 3 == 0 || 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 = fmaf(progress, progress, 1.0f);
			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));

			if (r_size < 64u) {
				size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}

text
	R.back() = curr;
	std::push_heap(R.begin(), R.end(), ComparePtr);
	R.emplace_back(curr1);
	std::push_heap(R.begin(), R.end(), ComparePtr);

	const Interval* const __restrict 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;
}

}

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;

text
int schetchick = 0;
const int rank = g_world->rank();
const int partner = rank ^ 1;
int dummy;

const float inv_divider = ldexpf(1.0f, -((gActiveMap.levels << 1) + 1));
const float x_addition = (b - a) * inv_divider;
const float y_addition = (d - c) * inv_divider;
const float true_start = a + x_addition;
const float true_end = b - x_addition;

float x_Rmax_1 = true_start;
float x_Rmax_2 = true_end;
const float initial_length = true_end - true_start;
float dmax = initial_length;
const float threshold_03 = 0.3f * initial_length;
const float inv_threshold_03 = 1.0f / threshold_03;
const float start_val = rank ? RastriginFunc(true_end, d - y_addition) : RastriginFunc(true_start, c + y_addition);
float best_f = rank ? RastriginFunc(true_start, d - y_addition) : RastriginFunc(true_end, c + y_addition);
float y_Rmax_1 = start_val;
float y_Rmax_2 = best_f;

std::vector<float, boost::alignment::aligned_allocator<float, 16u>> Extr;
std::vector<Interval* __restrict, boost::alignment::aligned_allocator<Interval* __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 Interval(true_start, true_end, start_val, best_f, 2.0f));
const Interval* __restrict top_ptr;

float Mmax = R.front()->M;
float m = r * Mmax;

while (true) {
	const float cooling = ldexpf(1.0f, -(1.0f / 138.63f) * ++schetchick);
	const int T = static_cast<int>(fmaf(20.0f, cooling, 10.0f));
	const float k = fmaf(0.2f, cooling, 0.7f);
	const float interval_len = x_Rmax_2 - x_Rmax_1;

	const bool want_xchg = schetchick % T == 0;
	const bool want_term = interval_len < epsilon || schetchick == static_cast<int>(global_iterations);

	if (want_xchg || want_term || g_world->iprobe(partner, 0)) {
		CtrlMsg out{};
		CtrlMsg in{};
		out.kind = want_term ? false : true;
		if (out.kind == true) {
			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 };
		}
		g_world->sendrecv(partner, 0, out, partner, 0, in);

		if (in.kind == false || out.kind == false) {
			if (partner) {
				Extr.emplace_back(static_cast<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);

		Interval* const __restrict injected = new Interval(sx, ex, RastriginFunc(in.xchg.s_x1, in.xchg.s_x2), RastriginFunc(in.xchg.e_x1, in.xchg.e_x2), 2.0f);
		injected->R = in.xchg.Rtop * k;
		R.emplace_back(injected);
		std::push_heap(R.begin(), R.end(), ComparePtr);
	}

	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);
	}

	std::pop_heap(R.begin(), R.end(), ComparePtr);
	Interval* const __restrict promejutochny_otrezok = R.back();

	const float segment_x1 = promejutochny_otrezok->x1;
	const float segment_x2 = promejutochny_otrezok->x2;
	const float len2 = segment_x2 - new_point;
	const float len1 = new_point - segment_x1;

	Interval* const __restrict curr = new Interval(segment_x1, new_point, promejutochny_otrezok->y1, new_value, 2.0f);
	Interval* const __restrict curr1 = new Interval(new_point, segment_x2, new_value, promejutochny_otrezok->y2, 2.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;
			size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
const float len_item = R[i]->x2 - R[i]->x1;
if (len_item > dmax) dmax = len_item;
++i;
}
}

text
		if (threshold_03 > dmax && schetchick % 3 == 0 || 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 = fmaf(progress, progress, 1.0f);
			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));

			if (r_size < 64u) {
				size_t i = 0u;

#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(fmaf(GLOBAL_FACTOR, R[i]->x2 - R[i]->x1, R[i]->M * alpha));
++i;
}
}
else {
RecomputeR_AffineM_AVX2(R.data(), r_size, GLOBAL_FACTOR, alpha);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}
}
else {
if (currM > Mmax) {
if (currM - Mmax < Mmax * 0.15f) {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
else {
Mmax = currM;
m = r * Mmax;
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
if (r_size < 64u) {
size_t i = 0u;
#pragma loop(ivdep)
while (i < r_size) {
R[i]->ChangeCharacteristic(m);
++i;
}
}
else {
RecomputeR_ConstM_AVX2(R.data(), r_size, m);
}
std::make_heap(R.begin(), R.end(), ComparePtr);
}
}
else {
curr->ChangeCharacteristic(m);
curr1->ChangeCharacteristic(m);
}
}

text
	R.back() = curr;
	std::push_heap(R.begin(), R.end(), ComparePtr);
	R.emplace_back(curr1);
	std::push_heap(R.begin(), R.end(), ComparePtr);

	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;
}

}

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

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AgpStartWorkers() noexcept
{
g_world->isend(1, 0, 0);
}

extern "C" __declspec(dllexport) __declspec(noalias) __forceinline void AGP_Free(float* const __restrict p) noexcept
{
CoTaskMemFree(p);
} рассмотри этот код - в целом он работает быстро и корректно - но нужно подправить механизм кросс-инъекций - сейчас они происходят раз в T итераций, и со временем этот период плавно уменьшается так как в начале поисковая информация может быть не так надёжна и в начале можно делать более редкие обмены - тогда как на заключительной стадии поиска нужны более частые обмены, но сейчас сделано так, что на заключительной части поиска вклад характеристик при обмене снижается и мне кажется это не правильно, сейчас вклад характеристик обмениваемых интервалов определяется через const float k = fmaf(0.2f, cooling, 0.7f); и этот коэффициент плавно снижается со временем - но мне кажется это неправильно, наоборот со временем поисковая информация становится более надёжной и вклад должен увеличиваться, так как это алгоритм глобального поиска - то поищи информацию о том какие стратегии обмена поисковой информации между процессами существуют - какие из них считаются самыми продвинутыми - и что из них я могу применить в своём коде - возможно стоит делать обмены не по периоду итераций а по длине максимального интервала dmax - если dmax достаточно мал это значит что уже идёт заключительная часть поиска и можно обмениваться чаще - такой подход может быть более надёжным - в любом случае предложи подход который обеспечит максимально быструю сходимость для моего случая - для тестовых функций с тригонометрическими коэффициентами

Answer 1

классный кусок — читается как «параллельный Strongin/Sergeyev» с Пеано-развёрткой и характеристиками. По сути у тебя две развилки:

когда обмениваться (частота),
как «доверять» пришедшей характеристике (вес k).

Сейчас и то и другое завязано на псевдотемпературу cooling, из-за чего в финале поиска обмены ещё более редкие, а вклад k падает — это действительно против логики глобального поиска: ближе к финалу информация становится надёжнее и её стоит внедрять агрессивнее.

Ниже — практичный план апгрейда + минимальные патчи прямо в твой код.

Что сделать концептуально

1) Триггер обмена: не по итерациям, а по «прогрессу»

Самая устойчивая метрика — нормированный прогресс по укрупнённой сетке:

p \;=\; 1 - \frac{d_{\max}}{d_{\text{init}}}\;\in[0,1]

Где dmax — текущая максимальная длина интервала в очереди, d_init — стартовая. Тогда:

Ранний этап (p≈0): редкие обмены (экономим трафик, инфа «сырая»).
Поздний этап (p→1): частые обмены (инфа стабильная, нужна интенсивная кооперация).

Практика «островных» (island) моделей показывает, что адаптивные интервалы миграции и асинхронные обмены стабильно выигрывают у фиксированного периода, особенно на неоднородных/многомодальных функциях. См. обзоры и эксперименты по адаптивному интервалу миграций и асинхронным топологиям. (CMap)

Надстройки к триггеру

Добавь ещё два «событийных» условия:

Застой: если нет улучшения best_f ≥ Icrit итераций — форсируем обмен.
Финишная стадия: если текущий «топовый» интервал стал короче порога (например, < 0.1 * d_init) — также форсируем.

Такие правила в island-моделях обычно дают лучшую динамику, чем чисто периодические обмены. (CMap)

2) Вес пришедшей характеристики: пересчитывать локально и нарастать к финалу

Сейчас ты ставишь injected->R = in.xchg.Rtop * k с уменьшающимся k. Лучше:

Сначала пересчитать R локально через твою ChangeCharacteristic(m) (иначе несогласованность из-за разных m у процессов; это классическая рекомендация для «методов по характеристикам» Strongin/Sergeyev). (unina2.on-line.it)
Потом умножить на возрастающий по прогрессу коэффициент доверия: $k(p)=k_{\min}+(k_{\max}-k_{\min})\cdot p^{\alpha}$ где, например, k_min=0.8f, k_max=1.15f, α∈[1,1.5].

Это отражает идею: ближе к концу чужая «верхняя характеристика» действительно более надёжна. В island-литературе аналог — fitness/performance-based migration weight, который усиливают при стабильной динамике «донора». (CMap)

Почему важно пересчитывать локально?
В твоём методе характеристика $R=\frac{Q}{m}+mN + \text{ord}$ зависит от локального Lipschitz-оценщика m=r*Mmax. Передача «сырых» $R$ без пересчёта ломает сопоставимость при разном m. Strongin/Sergeyev прямо подчёркивают пересчёт характеристик под локальные константы. (unina2.on-line.it)

3) (Опционально) Что именно обменивать

Сейчас передаёшь один «топовый» сегмент. Для триггера «финиш/застой» я бы высылал K=2–3 неперекрывающихся лучших промежутка (сжатое описание: [x1,x2],R), это часто ускоряет схождение на «гребёнчатых» триг-функциях (Rastrigin/Grishagin). Но это потребует расширить CrossMsg. Если не хочется — оставь один, но включи триггеры из пункта 1.

Патч (минимальные правки)

Внутри AGP_2D замени участок с cooling, T и k, а также место, где ты присваиваешь injected->R.

1) Вводим константы (до основного цикла):

cpp
int last_best_iter = 0;              // трекаем улучшения
const int   T_min   = 6,  T_max = 48;
const float k_min   = 0.80f, k_max = 1.15f;
const float alphaT  = 1.0f;          // кривизна частоты обменов
const float alphaK  = 1.2f;          // кривизна роста доверия
const int   Icrit   = 150;           // застойный порог по итерациям

2) В начале каждой итерации цикла (while (true)) считаем прогресс и частоту:

cpp
const float progress = fminf(fmaxf((initial_length - dmax) / initial_length, 0.0f), 1.0f);
const int   T        = T_min + (int)((T_max - T_min) * (1.0f - powf(progress, alphaT)));
const float k        = k_min + (k_max - k_min) * powf(progress, alphaK);

const float interval_len = x_Rmax_2 - x_Rmax_1;
const bool want_xchg = (schetchick % T == 0) || ((schetchick - last_best_iter) >= Icrit) || (interval_len < 0.10f * initial_length);

// прежний блок с iprobe оставляем как есть, только используем want_xchg
if (want_xchg || g_world->iprobe(partner, 0)) {
    CtrlMsg out{}, in{};
    out.kind = true;
    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 };
    g_world->sendrecv(partner, 0, out, partner, 0, in);

    if (in.kind == false) { /* как у тебя */ }
    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);
    Interval* __restrict injected =
        new Interval(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);
    injected->R *= k;   // усиливаем вклад по мере роста progress

    R.emplace_back(injected);
    std::push_heap(R.begin(), R.end(), ComparePtr);
}

3) Сброс счётчика застоя при улучшении (в том же цикле, сразу после блока с new_value):

cpp
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);
    last_best_iter = schetchick;     // сбрасываем застой
}

Всё остальное (адаптация m, перерасчёт R, AVX-ветки) — не трогаем.

Почему это сработает именно на «тригонометрических» тестах

Rastrigin/Grishagin/«холмистые» функции имеют много регулярных, но быстро затухающего масштаба ложных минимумов. На ранней фазе перекрёстная информация «шумновата» — редкие обмены логичны. Ближе к финалу локальные Lipschitz-оценки стабилизируются, и «чужие» узкие/мощные интервалы помогают быстро добить глобальник. Это ровно то, что делают адаптивные и асинхронные миграции в островных моделях. (CMap)
Пересчёт $R$ локально устраняет систематическую ошибку сравнения интервалов при разных m, о чём пишут Strongin/Sergeyev для методов по характеристикам и Пеано-развёрток. (unina2.on-line.it)

Что ещё можно (по желанию)

Отправлять K=2–3 интервала при want_xchg (требует расширить CrossMsg). В island-моделях обмен несколькими «неперекрывающимися» лучшими особями ускоряет сходимость на гребёнках. (ScienceDirect)
Для >2 процессов — кольцо/гиперкуб и асинхронные обмены без барьеров (см. ESA/ACT и недавние обзоры). (European Space Agency)

На что опираюсь

Адаптивные интервалы миграции и события «застоя/интенсивизации» в островных моделях ЭА/МОЭА. (CMap)
Асинхронные, производительные топологии обменов. (European Space Agency)
Параллельные детерминированные глобальные методы по характеристикам с Пеано-развёртками и рекомендацией локального пересчёта характеристик. (unina2.on-line.it)
Динамическая маршрутизация мигрантов к «перспективным» островам (аналог «взвешивания доверия»). (ertello.github.io)

Если хочешь, могу сразу прислать вариант CrossMsg на несколько интервалов и аккуратный sendrecv без дедлоков для K>1 (с фиксированным максимумом, чтобы размер сообщения был статический).

#include "pch.h" #define XOR_RAND(state, result_v...

Frage

Gedanke

Proposing optimization rules

Improving strategy with reliability and weight adjustments

Simplifying strategy for exchange and adaptive periods

Simplifying the approach for handling injections

Code and formula adjustments

Antwort