Lab3.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <chrono>
#include <random>
#include <numeric>
#include <map>
#include <algorithm>
#include <functional>
#include <iomanip>
#include <boost/math/distributions/chi_squared.hpp>

using namespace std;
using namespace chrono;

// Модифицированный Middle Square
vector<uint32_t> generate_middle_square_modified(uint64_t seed, size_t n) {
    vector<uint32_t> result(n);
    uint64_t x = seed;
    uint64_t w = 0;
    uint64_t s = 8527;
    for (size_t i = 0; i < n; ++i) {
        x *= x;
        x += (w += s);
        x = (x >> 13) | (x << 13);
        x %= 10000;
        result[i] = x;
    }
    return result;
}

// Модифицированный XorShift
vector<uint32_t> generate_xor_shift_modified(uint32_t seed, size_t n) {
    vector<uint32_t> result(n);
    uint32_t x = seed;
    for (int i = 0; i < n; ++i) {
        x ^= x << 13;
        x ^= x >> 17;
        x ^= x << 5;
        x %= 10000;
        result[i] = x;
    }
    return result;
}

// Модифицированный LCG
vector<uint32_t> generate_lcd_modifier(uint32_t seed, size_t n) {
    vector<uint32_t> result(n);
    uint32_t x = seed;
    for (size_t i = 0; i < n; ++i) {
        x = ((8963 * x + 7933) ^ (x << 5) ^ (x >> 3));
        x %= 10000;
        result[i] = x;
    }
    return result;
}

// Сохранения сгенерированной выборки в файл
void save_sample_to_file(const vector<uint32_t>& data, int sample_size, const string& filename) {
    ofstream out(filename);
    if (!out) return;

    for (int i = 0; i < sample_size; ++i) {
        out << data[i] << '\n';
    }
    out.close();
}

// Статистика: среднее, стандартное отклонение, коэффициент вариации
tuple<double, double, double> compute_stats(const vector<uint32_t>& data) {
    double mean = accumulate(data.begin(), data.end(), 0.0) / data.size();
    double variance = 0;
    for (auto val : data)
        variance += (val - mean) * (val - mean);
    variance /= data.size();
    double stddev = sqrt(variance);
    double cv = (mean != 0) ? stddev / mean : 0;
    return {mean, stddev, cv};
}

// Хи-квадрат тест равномерности
double chi_squared_test(const vector<uint32_t>& data) {
    if (data.empty()) return 0.0;
    int n = static_cast<int>(data.size());
    int k = static_cast<int>(1 + 3.322 * log10(n));
    if (k < 2) k = 2;
    vector<double> normalized(n);
    uint32_t min_val = *min_element(data.begin(), data.end());
    uint32_t max_val = *max_element(data.begin(), data.end());
    if (min_val == max_val) return 0.0;
    for (int i = 0; i < n; ++i) {
        normalized[i] = (static_cast<double>(data[i]) - min_val) / (max_val - min_val);
    }
    vector<int> counts(k, 0);
    for (int i = 0; i < n; ++i) {
        int index = static_cast<int>(normalized[i] * k);
        if (index == k) index = k - 1;
        counts[index]++;
    }
    double expected = static_cast<double>(n) / k;
    double chi2 = 0.0;
    for (int i = 0; i < k; ++i) {
        double diff = counts[i] - expected;
        chi2 += (diff * diff) / expected;
    }
    return chi2;
}

// Сохранение chi^2
void save_chi_squared_summary(const string& filename,
                              const string& sample_name,
                              double chi2_exp) {
    ofstream out(filename, ios::app);
    if (!out) return;
    const double alphas[] = {0.01, 0.1, 0.5};
    out << setw(25) << left << sample_name
        << setw(14) << fixed << setprecision(4) << chi2_exp;
    map<double, double> chi2_table = {{0.01, 21.6660}, {0.1, 14.6837}, {0.5, 8.3428}};
    for (int i = 0; i < 3; ++i) {
        double alpha = alphas[i];
        double chi2_crit = chi2_table.at(alpha);
        string decision = (chi2_exp <= chi2_crit) ? "H0" : "H1";
        out << setw(12) << chi2_crit
            << setw(6)  << decision;
    }
    out << '\n';
}

// Сохранение статистики
void save_stats_to_file(const string& stats_filename, const string& sample_name, double mean, double stddev, double variation) {
    ofstream out(stats_filename, ios::app);
    if (!out) return;
    out << setw(25) << left << sample_name
        << fixed    << setprecision(4)
        << setw(20) << mean
        << setw(20) << stddev
        << setw(20) << variation << '\n';
    out.close();
}

// Извлечение битов из выборки
vector<bool> extract_bits(const vector<uint32_t>& data) {
    vector<bool> bits;
    bits.reserve(data.size() * 13);

    for (uint32_t x : data) {
        for (int b = 12; b >= 0; --b) {
            bits.push_back((x >> b) & 1);
        }
    }
    return bits;
}

// Тест 1. Frequency (Monobit) Test
// Проверяет, насколько близко количество 1 и 0 в последовательности к равному числу.
// Основан на статистике суммы знаков (+1 за 1, −1 за 0), нормированной на sqrt(n).
// Предполагается, что в случайной последовательности единицы и нули равновероятны.
bool frequency_monobit_test(const vector<bool>& bits, double& p_value) {
    int n = static_cast<int>(bits.size());
    if (n == 0) {
        p_value = 0.0;
        return false;
    }
    int sum = 0;
    for (int i = 0; i < n; ++i) {
        sum += bits[i] ? 1 : -1;
    }
    double s_obs = abs(sum) / sqrt(n);
    p_value = erfc(s_obs / sqrt(2.0));
    // p_value это вероятность того, что отклонения возникли случайно, при условии что исходная последовательность случайна
    return p_value >= 0.01;
}

// Тест 2. Block Frequency Test
// Делит битовую последовательность на блоки фиксированного размера.
// Для каждого блока считается доля единиц и сравнивается с теоретической 0.5.
// Вычисляется xi2-статистика по всем блокам для оценки отклонения.
bool block_frequency_test(const vector<bool>& bits, int block_size, double& p_value) {
    int n = static_cast<int>(bits.size());
    if (block_size <= 0 || n < block_size) {
        p_value = 0.0;
        return false;
    }
    int num_blocks = n / block_size;
    if (num_blocks == 0) {
        p_value = 0.0;
        return false;
    }
    double chi2 = 0.0;
    for (int i = 0; i < num_blocks; ++i) {
        int count_ones = 0;
        for (int j = 0; j < block_size; ++j) {
            if (bits[i * block_size + j]) count_ones++;
        }
        double pi = static_cast<double>(count_ones) / block_size;
        chi2 += pow(pi - 0.5, 2);
    }
    chi2 *= 4.0 * block_size;
    boost::math::chi_squared_distribution<> dist(num_blocks);
    p_value = boost::math::cdf(boost::math::complement(dist, chi2));
    return p_value >= 0.01;
}

// Тест 3. Runs Test
// Проверяет число серий (подряд идущих одинаковых битов) в последовательности.
// При равномерной случайности число таких серий должно быть определённым при данной плотности единиц.
// Тест применим только если плотность близка к 0.5.
bool runs_test(const vector<bool>& bits, double& p_value) {
    int n = static_cast<int>(bits.size());
    if (n < 100) {
        p_value = 0.0;
        return false;
    }
    int ones = 0;
    for (int i = 0; i < n; ++i) {
        if (bits[i]) ones++;
    }
    double pi = static_cast<double>(ones) / n;
    if (abs(pi - 0.5) >= (2.0 / sqrt(n))) {
        p_value = 0.0;
        return false; 
    }
    int V = 1;
    for (int i = 1; i < n; ++i) {
        if (bits[i] != bits[i - 1]) {
            V++;
        }
    }
    double numerator = abs(V - 2.0 * n * pi * (1.0 - pi));
    double denominator = 2.0 * sqrt(2.0 * n) * pi * (1.0 - pi);
    double z = numerator / denominator;
    p_value = erfc(z / sqrt(2.0));
    return p_value >= 0.01;
}

// Тест 4. Longest Run of Ones in a Block Test
// Делит последовательность на блоки и измеряет максимальную длину серии единиц в каждом.
// Сравнивает частоты таких максимальных серий с теоретическими значениями.
// Использует категориальные интервалы и χ²-критерий для оценки отклонения.
bool longest_run_ones_in_block_test(const vector<bool>& bits, double& p_value) {
    int n = static_cast<int>(bits.size());
    const int M = 8;  // Размер блока
    const int K = 4;  // Количество категорий значений(из документации)
    int N = n / M;
    vector<int> counts(K, 0);
    // Ожидаемые вероятности (из документации)
    const double pi[K] = {0.2148, 0.3672, 0.2305, 0.1875};
    for (int i = 0; i < N; ++i) {
        int max_run = 0;
        int run = 0;
        for (int j = 0; j < M; ++j) {
            bool bit = bits[i * M + j];
            if (bit) {
                run++;
                if (run > max_run) max_run = run;
            } else {
                run = 0;
            }
        }
        if (max_run <= 1) counts[0]++;
        else if (max_run == 2) counts[1]++;
        else if (max_run == 3) counts[2]++;
        else counts[3]++;
    }
    double chi2 = 0.0;
    for (int i = 0; i < K; ++i) {
        double expected = pi[i] * N;
        chi2 += pow(counts[i] - expected, 2) / expected;
    }
    boost::math::chi_squared_distribution<> dist(K - 1);
    p_value = boost::math::cdf(boost::math::complement(dist, chi2));
    return p_value >= 0.01;
}

// Тест 5. Cumulative Sums Test (Forward и Backward)
// Преобразует биты в +1/-1 и считает накопленную сумму слева направо (или справа).
// Ищет максимальное отклонение от 0 — это мера "смещения" от равновесия.
// Ожидается, что сумма будет колебаться около нуля для случайной последовательности.
bool cumulative_sums_test(const vector<bool>& bits, bool forward, double& p_value) {
    int n = static_cast<int>(bits.size());
    if (n == 0) {
        p_value = 0.0;
        return false;
    }
    int sum = 0;
    int max_z = 0;
    for (int i = 0; i < n; ++i) {
        int bit = forward ? bits[i] : bits[n - 1 - i];
        int val = bit ? 1 : -1;
        sum += val;
        if (abs(sum) > max_z) {
            max_z = abs(sum);
        }
    }
    double z = max_z / sqrt(n);
    p_value = erfc(z / sqrt(2.0));
    return p_value >= 0.01;
}

// Запуск всех nist тестов
void run_all_nist_tests(const string& sample_name,
                        const vector<uint32_t>& data,
                        const string& output_file)
{
    vector<bool> bits = extract_bits(data);
    ofstream out(output_file, ios::app);
    if (!out) return;
    out << setw(25) << left << sample_name;
    double p;
    bool ok;
    ok = frequency_monobit_test(bits, p);
    out << setw(10) << fixed << setprecision(4) << p
        << setw(6)  << (ok ? "OK" : "FAIL");
    ok = block_frequency_test(bits, 128, p);
    out << setw(10) << p << setw(6) << (ok ? "OK" : "FAIL");
    ok = runs_test(bits, p);
    out << setw(10) << p << setw(6) << (ok ? "OK" : "FAIL");
    ok = longest_run_ones_in_block_test(bits, p);
    out << setw(10) << p << setw(6) << (ok ? "OK" : "FAIL");
    ok = cumulative_sums_test(bits, true, p);
    out << setw(10) << p << setw(6) << (ok ? "OK" : "FAIL");
    ok = cumulative_sums_test(bits, false, p);
    out << setw(10) << p << setw(6) << (ok ? "OK" : "FAIL");
    out.close();
}


int main() {
    const int num_samples = 20;
    const int sample_size = 1000;

    map<string, function<vector<uint32_t>(uint32_t, size_t)>> generators = {
        {"LCG_Modified", [](uint32_t seed, size_t n) { return generate_lcd_modifier(seed, n); }},
        {"MiddleSquare_Modified", [](uint32_t seed, size_t n) { return generate_middle_square_modified(seed, n); }},
        {"XorShift_Modified", [](uint32_t seed, size_t n) { return generate_xor_shift_modified(seed, n); }}
    };

    ofstream stats_file("results/stats.txt");
    stats_file << setw(25) << left       << "Sample"
               << setw(20) << "Mean"
               << setw(20) << "StdDev"
               << setw(20) << "VarCoeff" << "\n";
    stats_file.close();

    ofstream xi_out("results/xi.txt");
    xi_out << setw(25) << left        << "Sample"
           << setw(14) << "Chi2_exp"
           << setw(12) << "Chi2_0.01" << setw(6) << "H"
           << setw(12) << "Chi2_0.1"  << setw(6) << "H"
           << setw(12) << "Chi2_0.5"  << setw(6) << "H" << "\n";
    xi_out.close();

    ofstream nist_file("results/nist.txt");
    nist_file << setw(25) << left << "Sample"
              << setw(10) << "Freq"    << setw(6) << ""
              << setw(10) << "Block"   << setw(6) << ""
              << setw(10) << "Runs"    << setw(6) << ""
              << setw(10) << "Run1s"   << setw(6) << ""
              << setw(10) << "Cusum+"  << setw(6) << ""
              << setw(10) << "Cusum-"  << "\n";
    nist_file.close();

    for (const auto& [name, gen] : generators) {
        for (int i = 0; i < num_samples; ++i) {
            uint32_t seed = rand() % 10000 + 1;
            auto data = gen(seed, sample_size);
            auto [mean, stddev, cv] = compute_stats(data);
            auto chi2 = chi_squared_test(data);
            save_stats_to_file("results/stats.txt", name, mean, stddev, cv);
            save_sample_to_file(data, sample_size, "out_data/" + name + to_string(i) + ".txt");
            save_chi_squared_summary("results/xi.txt", name, chi2);
            run_all_nist_tests(name, data, "results/nist.txt");
        }
    }

    for (int i = 0; i < num_samples; ++i) {
        vector<uint32_t> data(sample_size);
        for (int j = 0; j < sample_size; ++j)
            data[j] = rand() % 10000;
        auto [mean, stddev, cv] = compute_stats(data);
        auto chi2 = chi_squared_test(data);
        save_stats_to_file("results/stats.txt", "C++ rand", mean, stddev, cv);
        save_sample_to_file(data, sample_size, "out_data/C++ rand" + to_string(i) + ".txt");
        save_chi_squared_summary("results/xi.txt", "C++ rand", chi2);
        run_all_nist_tests("C++ rand", data, "results/nist.txt");
    }

    vector<string> timing_lines = {"Method,Size,Time_ms"};
    vector<int> sizes = {1000, 2500, 5000, 7500, 10000, 25000, 50000, 75000, 100000, 250000, 500000, 750000, 1000000};

    for (const auto& [name, gen] : generators) {
        for (auto size : sizes) {
            auto start = high_resolution_clock::now();
            gen(12345, size);
            auto end = high_resolution_clock::now();
            auto ms = duration_cast<microseconds>(end - start).count();
            timing_lines.push_back(name + "," + to_string(size) + "," + to_string(ms));
        }
    }

    for (auto size : sizes) {
        auto start = high_resolution_clock::now();
        vector<int> v(size);
        for (int i = 0; i < size; ++i)
            v[i] = rand();
        auto end = high_resolution_clock::now();
        auto ms = duration_cast<microseconds>(end - start).count();
        timing_lines.push_back("C++ rand," + to_string(size) + "," + to_string(ms));
    }

    ofstream file("timing.csv");
    for (const auto& line : timing_lines)
        file << line << '\n';
    file.close();
    return 0;
}