// E04: two products using the complete old vector at each step.
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
using Vec = std::array<double, 2>;
using Mat = std::array<double, 4>;

Vec matvec(const Mat& A, const Vec& x) {
    Vec y{};  // Separate result: do not overwrite x while reading it.
    for (int i = 0; i < 2; ++i) {
        double sum = 0.0;
        for (int j = 0; j < 2; ++j) {
            sum += A[i * 2 + j] * x[j];
        }
        y[i] = sum;
    }
    return y;
}
int main() {
    const Mat A = {1, 2, 3, 4};
    const Vec x0 = {1, 2};
    const Vec x1 = matvec(A, x0);
    const Vec x2 = matvec(A, x1);
    const double n1 = std::sqrt(x1[0]*x1[0] + x1[1]*x1[1]);
    const double n2 = std::sqrt(x2[0]*x2[0] + x2[1]*x2[1]);
    std::cout << std::setprecision(12);
    std::cout << "x1 = " << x1[0] << ' ' << x1[1] << '\n';
    std::cout << "x2 = " << x2[0] << ' ' << x2[1] << '\n';
    std::cout << "norms = " << n1 << ' ' << n2 << '\n';
    if (x1 != Vec{5,11} || x2 != Vec{27,59} || x0 != Vec{1,2}) return 1;
    if (!std::isfinite(n1) || !std::isfinite(n2)) return 1;
    if (std::abs(n1 - std::sqrt(146.0)) > 1e-12 ||
        std::abs(n2 - std::sqrt(4210.0)) > 1e-12) return 1;
    return 0;
}