lab3_integral/lab3_integral.cpp - master - Konst_And/Labs_cpp

lab3_integral/lab3_integral.cpp
248 строк9 KB
And Konst
add lab 3
07 май 2026, 15:17
07 май 2026, 15:174d12eb3
#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>

// 1. Наивное суммирование

float mean_naive(const float psi[], const float pdf[], float dv, unsigned size) {
    float sum = 0.0f;
    for (unsigned i = 0; i < size; ++i)
        sum += psi[i] * pdf[i];
    return dv * sum;
}

// 2. Рекурсивное суммирование (деление пополам)

static float recursive_sum(const float arr[], unsigned l, unsigned r) {
    if (l + 1 == r) return arr[l];
    unsigned m = (l + r) / 2;
    return recursive_sum(arr, l, m) + recursive_sum(arr, m, r);
}

float mean_recursive(const float psi[], const float pdf[], float dv, unsigned size) {
    float* prod = new float[size];
    for (unsigned i = 0; i < size; ++i)
        prod[i] = psi[i] * pdf[i];
    float sum = recursive_sum(prod, 0, size);
    delete[] prod;
    return dv * sum;
}

// 3. Итеративное попарное суммирование (без рекурсии)

float mean_pairwise(const float psi[], const float pdf[], float dv, unsigned size) {
    float* arr = new float[size];
    for (unsigned i = 0; i < size; ++i)
        arr[i] = psi[i] * pdf[i];
    unsigned n = size;
    while (n > 1) {
        unsigned half = n / 2;
        for (unsigned i = 0; i < half; ++i)
            arr[i] = arr[2*i] + arr[2*i+1];
        if (n % 2 == 1) {
            arr[half] = arr[n-1];
            n = half + 1;
        } else {
            n = half;
        }
    }
    float sum = arr[0];
    delete[] arr;
    return dv * sum;
}

// 4. Суммирование Кэхена

float mean_kahan(const float psi[], const float pdf[], float dv, unsigned size) {
    float sum = 0.0f;
    float comp = 0.0f;
    for (unsigned i = 0; i < size; ++i) {
        float x = psi[i] * pdf[i];
        float y = x - comp;
        float t = sum + y;
        comp = (t - sum) - y;
        sum = t;
    }
    return dv * sum;
}

// 5. Суммирование с FMA

void Split(float x, float &x_high, float &x_low) {
    const float C = 4097.0f;   // 2^12 + 1
    float tmp = C * x;
    x_high = tmp - (tmp - x);
    x_low = x - x_high;
}

void TwoMult(float a, float b, float &s, float &t) {
    s = a * b;
    float a_high, a_low, b_high, b_low;
    Split(a, a_high, a_low);
    Split(b, b_high, b_low);
    t = -s + a_high * b_high;
    t += a_high * b_low;
    t += a_low * b_high;
    t += a_low * b_low;
}

void TwoSum(float a, float b, float &s, float &t) {
    s = a + b;
    float z = s - a;
    t = (a - (s - z)) + (b - z);
}

float my_fma(float a, float b, float c) {
    float prod, err_prod;
    TwoMult(a, b, prod, err_prod);
    float sum, err_sum;
    TwoSum(prod, c, sum, err_sum);
    return sum + (err_prod + err_sum);
}

float mean_fma(const float psi[], const float pdf[], float dv, unsigned size) {
    float sum = 0.0f;
    for (unsigned i = 0; i < size; ++i)
        sum = my_fma(psi[i], pdf[i], sum);
    return dv * sum;
}

// 5b. Отдельное сложение и умножение

float mean_fma_separate(const float psi[], const float pdf[], float dv, unsigned size) {
    float sum = 0.0f;
    for (unsigned i = 0; i < size; ++i) {
        float prod = psi[i] * pdf[i];
        sum = sum + prod;
    }
    return dv * sum;
}

// 5c. Аппаратный FMA
float mean_std_fma(const float psi[], const float pdf[], float dv, unsigned size) {
    float sum = 0.0f;
    for (unsigned i = 0; i < size; ++i) {
        sum = std::fma(psi[i], pdf[i], sum);
    }
    return dv * sum;
}

// 6. Суммирование в double

double mean_double(const float psi[], const float pdf[], float dv, unsigned size) {
    double sum = 0.0;
    for (unsigned i = 0; i < size; ++i)
        sum += static_cast<double>(psi[i]) * static_cast<double>(pdf[i]);
    return dv * sum;
}

// main

int main() {
    // Параметры расчёта
    const float T_vals[] = {0.1f, 1.0f, 100.0f, 10000.0f};
    const unsigned nT = sizeof(T_vals) / sizeof(T_vals[0]);

    // Различные размеры сетки 
    const unsigned sizes[] = {10, 50, 100, 1000, 10000, 2000000};
    const unsigned nSizes = sizeof(sizes) / sizeof(sizes[0]);

    std::ofstream out("results.csv");
    if (!out.is_open()) {
        std::cerr << "Ошибка открытия results.csv\n";
        return 1;
    }

    // 7 знаков после запятой
    out << std::scientific << std::setprecision(7);
    out << "# T size method mean_abs rel_err_abs mean_sq rel_err_sq\n";

    for (unsigned t_idx = 0; t_idx < nT; ++t_idx) {
        float T = T_vals[t_idx];
        double abs_analit = std::sqrt(T / M_PI);
        double sq_analit  = static_cast<double>(T) / 2.0;

        for (unsigned sz_idx = 0; sz_idx < nSizes; ++sz_idx) {
            unsigned size = sizes[sz_idx];
            // Хвосты (10σ)
            float v_max = 10.0f * std::sqrt(T);
            float dv = 2.0f * v_max / (size - 1);

  
            float* psi_abs = new float[size];
            float* psi_sq  = new float[size];
            float* pdf     = new float[size];

            // Заполнение массивов
            float inv_sqrt_pi_T = 1.0f / std::sqrt(static_cast<float>(M_PI) * T);
            for (unsigned i = 0; i < size; ++i) {
                float x = -v_max + i * dv;
                psi_abs[i] = std::fabs(x);
                psi_sq[i]  = x * x;
                pdf[i]     = std::exp(-x * x / T) * inv_sqrt_pi_T;
            }

            // Вычисление среднего |v|
            float abs_naive   = mean_naive(psi_abs, pdf, dv, size);
            float abs_rec     = mean_recursive(psi_abs, pdf, dv, size);
            float abs_pair    = mean_pairwise(psi_abs, pdf, dv, size);
            float abs_kahan   = mean_kahan(psi_abs, pdf, dv, size);
            float abs_fma     = mean_fma(psi_abs, pdf, dv, size);
            float abs_fma_sep = mean_fma_separate(psi_abs, pdf, dv, size);
            double abs_double = mean_double(psi_abs, pdf, dv, size);
            float abs_std_fma = mean_std_fma(psi_abs, pdf, dv, size);

            // Вычисление среднего v^2
            float sq_naive   = mean_naive(psi_sq, pdf, dv, size);
            float sq_rec     = mean_recursive(psi_sq, pdf, dv, size);
            float sq_pair    = mean_pairwise(psi_sq, pdf, dv, size);
            float sq_kahan   = mean_kahan(psi_sq, pdf, dv, size);
            float sq_fma     = mean_fma(psi_sq, pdf, dv, size);
            float sq_fma_sep = mean_fma_separate(psi_sq, pdf, dv, size);
            double sq_double = mean_double(psi_sq, pdf, dv, size);
            float sq_std_fma = mean_std_fma(psi_sq, pdf, dv, size);

            // Расчёт относительной ошибки
            auto rel_err = [](double val, double analit) -> double {
                return std::fabs(val - analit) / analit;
            };

            out << T << " " << size << " naive "
                << abs_naive << " " << rel_err(abs_naive, abs_analit) << " "
                << sq_naive  << " " << rel_err(sq_naive, sq_analit) << "\n";
            out << T << " " << size << " recursive "
                << abs_rec << " " << rel_err(abs_rec, abs_analit) << " "
                << sq_rec  << " " << rel_err(sq_rec, sq_analit) << "\n";
            out << T << " " << size << " pairwise "
                << abs_pair << " " << rel_err(abs_pair, abs_analit) << " "
                << sq_pair  << " " << rel_err(sq_pair, sq_analit) << "\n";
            out << T << " " << size << " kahan "
                << abs_kahan << " " << rel_err(abs_kahan, abs_analit) << " "
                << sq_kahan  << " " << rel_err(sq_kahan, sq_analit) << "\n";
            out << T << " " << size << " fma "
                << abs_fma << " " << rel_err(abs_fma, abs_analit) << " "
                << sq_fma  << " " << rel_err(sq_fma, sq_analit) << "\n";
            out << T << " " << size << " fma_separate "
                << abs_fma_sep << " " << rel_err(abs_fma_sep, abs_analit) << " "
                << sq_fma_sep  << " " << rel_err(sq_fma_sep, sq_analit) << "\n";
            out << T << " " << size << " std_fma "
                << abs_std_fma << " " << rel_err(abs_std_fma, abs_analit) << " "
                << sq_std_fma  << " " << rel_err(sq_std_fma, sq_analit) << "\n";
            out << T << " " << size << " double "
                << abs_double << " " << rel_err(abs_double, abs_analit) << " "
                << sq_double  << " " << rel_err(sq_double, sq_analit) << "\n";

            // Освобождение памяти
            delete[] psi_abs;
            delete[] psi_sq;
            delete[] pdf;
        }
    }

    out.close();

    // компилировать с -O0 !!!

    return 0;
}