C++26 и SIMD: Data-Parallel Types

// Типичный код с интринсиками до C++26
__m128 a = _mm_set_ps(4.0f, 3.0f, 2.0f, 1.0f);
__m128 b = _mm_set_ps(8.0f, 7.0f, 6.0f, 5.0f);
__m128 c = _mm_add_ps(a, b);  // Векторное сложение

// Современный подход с data-parallel types в C++26
std::simd<float> a = {1.0f, 2.0f, 3.0f, 4.0f};
std::simd<float> b = {5.0f, 6.0f, 7.0f, 8.0f};
std::simd<float> c = a + b;  // То же самое векторное сложение, но портабельно!

std::vector<float> a = {1.0f, 2.0f, 3.0f, 4.0f};
std::vector<float> b = {5.0f, 6.0f, 7.0f, 8.0f};
std::vector<float> result(a.size());
 
for (size_t i = 0; i < a.size(); ++i) {
    result[i] = a[i] + b[i];  // Последовательная обработка
}

std::simd<float> a = {1.0f, 2.0f, 3.0f, 4.0f};
std::simd<float> b = {5.0f, 6.0f, 7.0f, 8.0f};
std::simd<float> result = a + b;  // Параллельная обработка всех элементов

#ifdef __AVX2__
    // Реализация для AVX2
#elif defined(__ARM_NEON)
    // Реализация для NEON
#else
    // Обычная реализация
#endif

template <typename T>
void process_array(T* data, size_t size) {
    using simd_t = std::simd<T>;
    constexpr size_t simd_size = simd_t::size();
    
    size_t i = 0;
    for (; i + simd_size <= size; i += simd_size) {
        simd_t values(data + i);
        values = some_operation(values);
        values.copy_to(data + i);
    }
    
    // Обработка остатка
    for (; i < size; ++i) {
        data[i] = some_operation_scalar(data[i]);
    }
}

for (int i = 0; i < size; ++i) {
    histogram[data[i]]++;
}

void add_arrays(float* a, float* b, float* result, int size) {
    for (int i = 0; i < size; ++i) {
        result[i] = a[i] + b[i];
    }
}

template <typename T>
void process_array(T* data, size_t size) {
    using simd_t = std::simd<T>;
    constexpr size_t simd_size = simd_t::size();
    
    // Обработка полных SIMD-векторов
    size_t i = 0;
    for (; i + simd_size <= size; i += simd_size) {
        simd_t values(data + i);
        // Обработка полного вектора
        values.copy_to(data + i);
    }
    
    // Обработка "хвоста" с помощью маски
    if (i < size) {
        size_t remaining = size - i;
        std::simd_mask<T> mask(false);
        for (size_t j = 0; j < remaining; ++j) {
            mask[j] = true;
        }
        
        simd_t values(data + i, mask);  // Загрузка с маской
        // Обработка с учетом маски
        values.copy_to(data + i, mask);  // Сохранение с маской
    }
}

// Выделение выровненной памяти (до C++17)
float* aligned_data = static_cast<float*>(
    _mm_malloc(size * sizeof(float), 32)  // Выравнивание на 32 байта (AVX)
);
 
// Современный подход с C++17
std::vector<float, aligned_allocator<float, 32>> data(size);

template <typename T>
void process_with_prefetch(T* data, size_t size) {
    using simd_t = std::simd<T>;
    constexpr size_t simd_size = simd_t::size();
    
    for (size_t i = 0; i < size; i += simd_size) {
        // Предварительная загрузка данных, которые понадобятся через несколько итераций
        if (i + simd_size * 16 < size) {
            __builtin_prefetch(data + i + simd_size * 16);
        }
        
        simd_t values(data + i);
        // Обработка...
        values.copy_to(data + i);
    }
}

std::vector<float> input = /* ... */;
std::vector<float> output(input.size());
 
std::transform(
    std::execution::par_unseq,  // Параллельное и векторизуемое выполнение
    input.begin(), input.end(),
    output.begin(),
    [](float x) { return std::sin(x); }
);

// Предполагаемый синтаксис в C++26
void f(std::vector<float>& data) {
    std::for_each(std::execution::simd, data.begin(), data.end(),
        [](auto& v) {
            v = std::sin(v);
        }
    );
}

// Создание вектора и заполнение одним значением
std::simd<float> a(3.14f);  // Все элементы равны 3.14
 
// Инициализация из массива
float data[] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f};
std::simd<float> b(data);  // Загрузка из массива
 
// Инициализация списком (если количество элементов совпадает с шириной вектора)
std::simd<int> c = {1, 2, 3, 4};  // Работает только если std::simd<int>::size() == 4
 
// Инициализация с использованием генератора
std::simd<double> d = std::experimental::simd_cast<double>(
    std::iota(std::experimental::simd_index<std::simd<double>>)
);  // Заполнение последовательными числами [0, 1, 2, ...]

std::simd<float> a = {1.0f, 2.0f, 3.0f, 4.0f};
std::simd<float> b = {5.0f, 6.0f, 7.0f, 8.0f};
 
// Арифметические операции
std::simd<float> sum = a + b;      // {6.0f, 8.0f, 10.0f, 12.0f}
std::simd<float> diff = a - b;     // {-4.0f, -4.0f, -4.0f, -4.0f}
std::simd<float> prod = a * b;     // {5.0f, 12.0f, 21.0f, 32.0f}
std::simd<float> div = a / b;      // {0.2f, 0.333f, 0.428f, 0.5f}
 
// Смешанные операции со скалярами
std::simd<float> scaled = a * 2.0f;  // {2.0f, 4.0f, 6.0f, 8.0f}
 
// Встроенные функции
std::simd<float> sqr = std::sqrt(a); // {1.0f, 1.414f, 1.732f, 2.0f}
std::simd<float> e_pow = std::exp(a); // {2.718f, 7.389f, 20.086f, 54.598f}

std::simd<float> values = {1.0f, 2.0f, 3.0f, 4.0f};
std::simd<float> threshold(2.5f);
 
// Создание маски: true для элементов, где условие выполняется
std::simd_mask<float> mask = values > threshold;  // {false, false, true, true}
 
// Условное присваивание: заменить только элементы, где маска true
std::simd<float> result = values;
where(mask, result) = result * 2.0f;  // {1.0f, 2.0f, 6.0f, 8.0f}
 
// Условный выбор между двумя векторами
std::simd<float> a = {10.0f, 20.0f, 30.0f, 40.0f};
std::simd<float> b = {1.0f, 2.0f, 3.0f, 4.0f};
std::simd<float> selected = std::experimental::choose(mask, a, b);
// {1.0f, 2.0f, 30.0f, 40.0f}

std::simd<float> angles = {0.0f, 0.5f, 1.0f, 1.5f};
 
// Тригонометрические функции
std::simd<float> sines = std::sin(angles);
std::simd<float> cosines = std::cos(angles);
std::simd<float> tangents = std::tan(angles);
 
// Экспоненциальные и логарифмические функции
std::simd<float> logs = std::log(angles + 1.0f);  // +1.0f чтобы избежать log(0)
std::simd<float> exp_values = std::exp(angles);
std::simd<float> pow_values = std::pow(angles + 1.0f, 2.0f);
 
// Округление и абсолютное значение
std::simd<float> rounded = std::round(angles);
std::simd<float> abs_values = std::abs(angles - 1.0f);

void process_array(float* input, float* output, size_t size) {
    using simd_t = std::simd<float>;
    constexpr size_t simd_size = simd_t::size();
    
    // Обработка основной части массива блоками по simd_size элементов
    size_t i = 0;
    for (; i + simd_size <= size; i += simd_size) {
        simd_t block(input + i);
        
        // Выполнение каких-то вычислений над блоком
        block = std::sqrt(block) + simd_t(1.0f);
        
        // Сохранение результатов обратно в массив
        block.copy_to(output + i);
    }
    
    // Обработка оставшихся элементов
    for (; i < size; ++i) {
        output[i] = std::sqrt(input[i]) + 1.0f;
    }
}

void process_array(std::vector<float>& data) {
    std::transform(
        std::execution::simd,
        data.begin(), data.end(),
        data.begin(),
        [](float x) { return std::sqrt(x) + 1.0f; }
    );
}

void blur_image(uint8_t* input, uint8_t* output, int width, int height) {
    using simd_t = std::simd<float>;
    constexpr size_t simd_size = simd_t::size();
    
    for (int y = 1; y < height - 1; ++y) {
        for (int x = 1; x < width - 1; x += simd_size) {
            // Ограничиваем ширину блока, чтобы не выйти за границы
            int block_size = std::min<int>(simd_size, width - 1 - x);
            std::simd_mask<float> mask = std::simd_mask<float>(true);
            for (int i = block_size; i < simd_size; ++i) {
                mask[i] = false;
            }
            
            // Загружаем 9 блоков пикселей (матрица 3x3 вокруг текущего пикселя)
            simd_t top_left(input + (y-1)*width + (x-1), mask);
            simd_t top(input + (y-1)*width + x, mask);
            simd_t top_right(input + (y-1)*width + (x+1), mask);
            
            simd_t left(input + y*width + (x-1), mask);
            simd_t center(input + y*width + x, mask);
            simd_t right(input + y*width + (x+1), mask);
            
            simd_t bottom_left(input + (y+1)*width + (x-1), mask);
            simd_t bottom(input + (y+1)*width + x, mask);
            simd_t bottom_right(input + (y+1)*width + (x+1), mask);
            
            // Вычисляем среднее значение (простейший фильтр размытия)
            simd_t result = (top_left + top + top_right +
                            left + center + right +
                            bottom_left + bottom + bottom_right) / 9.0f;
            
            // Сохраняем результат
            result.copy_to(output + y*width + x, mask);
        }
    }
}

// Вычисление полинома a*x^3 + b*x^2 + c*x + d
std::vector<float> evaluate_polynomial(
    const std::vector<float>& x_values,
    float a, float b, float c, float d
) {
    std::vector<float> results(x_values.size());
    
    using simd_t = std::simd<float>;
    constexpr size_t simd_size = simd_t::size();
    
    simd_t a_vec(a), b_vec(b), c_vec(c), d_vec(d);
    
    for (size_t i = 0; i + simd_size <= x_values.size(); i += simd_size) {
        simd_t x(x_values.data() + i);
        
        // Вычисление полинома методом Горнера
        simd_t result = a_vec;
        result = result * x + b_vec;
        result = result * x + c_vec;
        result = result * x + d_vec;
        
        result.copy_to(results.data() + i);
    }
    
    // Обработка оставшихся элементов
    for (size_t i = (x_values.size() / simd_size) * simd_size; i < x_values.size(); ++i) {
        float x = x_values[i];
        results[i] = ((a * x + b) * x + c) * x + d;
    }
    
    return results;
}

void compress_audio(float* samples, size_t count, float threshold, float ratio) {
    using simd_t = std::simd<float>;
    constexpr size_t simd_size = simd_t::size();
    
    simd_t threshold_vec(threshold);
    simd_t ratio_vec(ratio);
    simd_t one(1.0f);
    
    for (size_t i = 0; i + simd_size <= count; i += simd_size) {
        simd_t input(samples + i);
        simd_t abs_input = std::abs(input);
        
        // Создаем маску для значений выше порога
        simd_mask<float> mask = abs_input > threshold_vec;
        
        // Применяем компрессию только к значениям выше порога
        simd_t compressed = input;
        where(mask, compressed) = threshold_vec + 
            (abs_input - threshold_vec) / ratio_vec * 
            (input / abs_input);
        
        compressed.copy_to(samples + i);
    }
    
    // Обработка оставшихся элементов скалярным кодом
    for (size_t i = (count / simd_size) * simd_size; i < count; ++i) {
        float input = samples[i];
        float abs_input = std::abs(input);
        if (abs_input > threshold) {
            float sign = (input > 0.0f) ? 1.0f : -1.0f;
            samples[i] = sign * (threshold + (abs_input - threshold) / ratio);
        }
    }
}

template <typename T, typename Func>
bool verify_simd_computation(const std::vector<T>& input, Func scalar_func) {
    // Вычисление с использованием SIMD
    std::vector<T> simd_result(input.size());
    {
        using simd_t = std::simd<T>;
        constexpr size_t simd_size = simd_t::size();
        
        for (size_t i = 0; i + simd_size <= input.size(); i += simd_size) {
            simd_t x(input.data() + i);
            simd_t result = scalar_func(x);  // Векторизованная функция
            result.copy_to(simd_result.data() + i);
        }
        
        // Обработка оставшихся элементов
        for (size_t i = (input.size() / simd_size) * simd_size; i < input.size(); ++i) {
            simd_result[i] = scalar_func(input[i]);  // Скалярная функция
        }
    }
    
    // Вычисление с использованием скалярных операций
    std::vector<T> scalar_result(input.size());
    for (size_t i = 0; i < input.size(); ++i) {
        scalar_result[i] = scalar_func(input[i]);
    }
    
    // Сравнение результатов с допустимой погрешностью
    constexpr T epsilon = std::numeric_limits<T>::epsilon() * 10;  // Допустимая погрешность
    for (size_t i = 0; i < input.size(); ++i) {
        if (std::abs(simd_result[i] - scalar_result[i]) > epsilon * std::abs(scalar_result[i])) {
            std::cerr << "Несоответствие в позиции " << i 
                      << ": SIMD=" << simd_result[i] 
                      << ", скаляр=" << scalar_result[i] << std::endl;
            return false;
        }
    }
    
    return true;
}

// Векторизованная сортировка битонной последовательности
void bitonic_sort_simd(float* data, size_t size) {
    using simd_t = std::simd<float>;
    constexpr size_t simd_size = simd_t::size();
    
    // Предполагаем, что size кратно simd_size и является степенью 2
    for (size_t k = 2; k <= size; k *= 2) {
        for (size_t j = k / 2; j > 0; j /= 2) {
            for (size_t i = 0; i < size; i += j * 2) {
                for (size_t l = 0; l < j; l += simd_size) {
                    simd_t a(data + i + l);
                    simd_t b(data + i + j + l);
                    
                    // Сравнение и обмен векторов
                    simd_t min_vals = std::min(a, b);
                    simd_t max_vals = std::max(a, b);
                    
                    min_vals.copy_to(data + i + l);
                    max_vals.copy_to(data + i + j + l);
                }
            }
        }
    }
}

// Векторизованный бинарный поиск
template <typename T>
std::vector<size_t> simd_binary_search(
    const std::vector<T>& sorted_array,
    const std::vector<T>& keys
) {
    using simd_t = std::simd<T>;
    using index_t = std::simd<int>;
    constexpr size_t simd_size = simd_t::size();
    
    std::vector<size_t> results(keys.size(), -1);
    
    for (size_t k = 0; k + simd_size <= keys.size(); k += simd_size) {
        simd_t key_vec(keys.data() + k);
        
        // Начальные границы поиска
        index_t left(0);
        index_t right(sorted_array.size() - 1);
        
        // Векторизованный бинарный поиск
        simd_mask<T> not_found = left <= right;
        while (any_of(not_found)) {
            index_t mid = (left + right) / 2;
            
            // Загружаем средние элементы для всех ключей
            simd_t mid_vals;
            for (size_t i = 0; i < simd_size; ++i) {
                if (not_found[i]) {
                    mid_vals[i] = sorted_array[mid[i]];
                }
            }
            
            // Обновляем границы поиска
            simd_mask<T> is_less = key_vec < mid_vals;
            where(is_less && not_found, right) = mid - 1;
            where(!is_less && not_found, left) = mid + 1;
            
            // Нашли точное совпадение?
            simd_mask<T> found = mid_vals == key_vec;
            
            // Сохраняем результаты
            for (size_t i = 0; i < simd_size; ++i) {
                if (found[i]) {
                    results[k + i] = mid[i];
                    not_found[i] = false;
                }
            }
            
            // Выход из цикла, если все не найдены и границы поиска пересеклись
            not_found &= left <= right;
        }
    }
    
    // Обработка оставшихся ключей скалярным методом
    // ...
    
    return results;
}

// Векторизованное умножение разреженной матрицы на вектор (формат CSR)
void sparse_matrix_vector_mul_simd(
    const std::vector<float>& values,     // Ненулевые элементы
    const std::vector<int>& col_indices,  // Индексы столбцов
    const std::vector<int>& row_ptr,      // Указатели на начало строк
    const std::vector<float>& x,          // Вектор, на который умножаем
    std::vector<float>& result            // Результат
) {
    using simd_t = std::simd<float>;
    constexpr size_t simd_size = simd_t::size();
    
    for (size_t i = 0; i < row_ptr.size() - 1; ++i) {
        int row_start = row_ptr[i];
        int row_end = row_ptr[i+1];
        
        simd_t sum(0.0f);
        
        for (int j = row_start; j + simd_size <= row_end; j += simd_size) {
            simd_t val_vec(values.data() + j);
            
            // Собираем значения из вектора x по индексам столбцов
            simd_t x_vec;
            for (size_t k = 0; k < simd_size; ++k) {
                x_vec[k] = x[col_indices[j + k]];
            }
            
            sum += val_vec * x_vec;
        }
        
        // Суммируем элементы вектора
        float row_result = reduce(sum, std::plus<float>());
        
        // Обработка оставшихся элементов скалярным способом
        for (int j = row_end - ((row_end - row_start) % simd_size); j < row_end; ++j) {
            row_result += values[j] * x[col_indices[j]];
        }
        
        result[i] = row_result;
    }
}

// Выделение выровненной памяти с C++17
std::vector<float, std::aligned_allocator<float, 32>> data(size);
 
// В C++26 с Data-Parallel Types можно использовать:
constexpr size_t alignment = alignof(std::simd<float>);
std::vector<float, std::aligned_allocator<float, alignment>> data(size);

// Этот код будет работать везде, но с разной эффективностью
std::simd<float> values = /* ... */;
float sum = reduce(values, std::plus<float>());

// Скалярный код - прост для понимания
float result = 0.0f;
for (size_t i = 0; i < data.size(); ++i) {
    if (data[i] > threshold) {
        result += std::sin(data[i]);
    }
}
 
// Векторизованный код - менее очевиден
std::simd<float> sum(0.0f);
for (size_t i = 0; i + std::simd<float>::size() <= data.size(); i += std::simd<float>::size()) {
    std::simd<float> values(data.data() + i);
    std::simd_mask<float> mask = values > std::simd<float>(threshold);
    sum += std::sin(values) & mask;
}
float result = reduce(sum, std::plus<float>());
// Плюс еще код обработки "хвоста"...

// Такой код сложно эффективно векторизовать
for (size_t i = 0; i < data.size(); ++i) {
    if (data[i] < 0)
        result[i] = -data[i];
    else if (data[i] < threshold)
        result[i] = data[i] * factor1;
    else if (data[i] < upper_bound)
        result[i] = std::pow(data[i], exponent);
    else
        result[i] = max_value;
}

// Функция для верификации SIMD-вычислений
template <typename T, typename Func>
bool verify_simd_vs_scalar(const std::vector<T>& input, Func func) {
    std::vector<T> simd_results = compute_with_simd(input, func);
    std::vector<T> scalar_results = compute_scalar(input, func);
    
    // Сравнение с допустимой погрешностью
    for (size_t i = 0; i < input.size(); ++i) {
        if (!almost_equal(simd_results[i], scalar_results[i])) {
            std::cerr << "Mismatch at " << i << ": " 
                      << simd_results[i] << " vs " << scalar_results[i] << std::endl;
            return false;
        }
    }
    return true;
}

// Векторизованная реализация фильтра Гаусса с использованием std::simd
void apply_gaussian_blur_simd(Image& image, float sigma) {
    using simd_float = std::simd<float>;
    constexpr size_t simd_size = simd_float::size();
    
    // Создаём ядро фильтра Гаусса
    const int kernel_size = static_cast<int>(ceil(sigma * 3) * 2 + 1);
    const int kernel_radius = kernel_size / 2;
    std::vector<float> kernel(kernel_size);
    
    // Заполняем ядро значениями гауссианы
    float sum = 0.0f;
    for (int i = 0; i < kernel_size; ++i) {
        int x = i - kernel_radius;
        kernel[i] = std::exp(-(x * x) / (2 * sigma * sigma));
        sum += kernel[i];
    }
    
    // Нормализуем ядро
    for (int i = 0; i < kernel_size; ++i) {
        kernel[i] /= sum;
    }
    
    // Создаём временное изображение для промежуточных результатов
    Image temp(image.width, image.height);
    
    // Горизонтальная свёртка (по оси X)
    for (int y = 0; y < image.height; ++y) {
        for (int x = 0; x + simd_size <= image.width; x += simd_size) {
            // Для каждого канала (R, G, B)
            for (int c = 0; c < 3; ++c) {
                simd_float sum(0.0f);
                
                // Применяем ядро свёртки
                for (int k = -kernel_radius; k <= kernel_radius; ++k) {
                    // Вычисляем позицию с учётом граничных условий
                    simd_float values;
                    for (size_t i = 0; i < simd_size; ++i) {
                        int pos_x = std::clamp(x + i + k, 0, image.width - 1);
                        values[i] = image.pixel(pos_x, y)[c];
                    }
                    sum += values * kernel[k + kernel_radius];
                }
                
                // Сохраняем результаты
                for (size_t i = 0; i < simd_size; ++i) {
                    temp.pixel(x + i, y)[c] = std::clamp(sum[i], 0.0f, 255.0f);
                }
            }
        }
        
        // Обработка оставшихся пикселей в строке (хвост)
        // ...
    }
    
    // Вертикальная свёртка (по оси Y) - аналогично
    // ...
    
    // Копируем результат обратно в исходное изображение
    image = std::move(temp);
}