// E08: a bounded timing demonstration with a separate analytic oracle.
#include <algorithm>
#include <chrono>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

void multiply(const std::vector<double>& A, const std::vector<double>& B,
              std::vector<double>& C, int n) {
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            double sum = 0.0;
            for (int k = 0; k < n; ++k) sum += A[i*n+k] * B[k*n+j];
            C[i*n+j] = sum;
        }
    }
}
int main() {
    const int n = 96;
    std::vector<double> A(n*n), B(n*n, 1.0), C(n*n);
    std::vector<double> times;
    std::cout << std::setprecision(12);
    for (int t = 0; t <= 5; ++t) {  // t=0 is the unreported warm-up.
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) A[i*n+j] = i+j+1+t;
        const auto start = std::chrono::steady_clock::now();
        multiply(A, B, C, n);
        const auto stop = std::chrono::steady_clock::now();
        const double seconds = std::chrono::duration<double>(stop-start).count();
        if (!std::isfinite(seconds) || seconds < 0) return 1;
        double checksum = 0.0;
        for (int i = 0; i < n; ++i) {
            const double expected = n*(i+1.0+t) + n*(n-1)/2.0;
            for (int j = 0; j < n; ++j) {
                const double value = C[i*n+j];
                if (!std::isfinite(value) || std::abs(value-expected) > 1e-9)
                    return 1;
                checksum += value;
            }
        }
        if (t > 0) {
            times.push_back(seconds);
            std::cout << "run=" << t << " seconds=" << seconds
                      << " checksum=" << checksum << '\n';
        }
    }
    std::sort(times.begin(), times.end());
    std::cout << "median_seconds=" << times[2]
              << " min=" << times.front() << " max=" << times.back() << '\n';
    std::cout << "All entries checked in warm-up and five timed runs\n";
    return 0;
}