#include <iostream>
#include <algorithm>
#include <climits>
using namespace std;

// Индекс самого старшего выставленного бита, считая младший бит за №0:
// floor_log2(0b101000) == 5, floor_log2(1 << N) == N.
//              543210
int floor_log2(unsigned int p)
{
	return CHAR_BIT * sizeof(unsigned int) - 1 - __builtin_clz(p);
}

template<typename Item>
class matrix
{
public:
	matrix(size_t sx, size_t sy)
		: _sx(sx), _sy(sy), items(nullptr)
	{
		items = new Item[sx * sy];
		set_identity();
	}

	~matrix() { delete[] items; }
	size_t sx() const { return _sx; }
	size_t sy() const { return _sy; }

	const Item& operator()(size_t x, size_t y) const { return items[y * _sx + x]; }
	Item& operator()(size_t x, size_t y) { return items[y * _sx + x]; }

	void set_identity()
	{
		fill(items, items + _sx * _sy, 0);
		for (size_t xy = 0; xy < min(_sx, _sy); xy++) (*this)(xy, xy) = 1;
	}

	void set_copy(const matrix& b)
	{
		if (_sx != b._sx || _sy != b._sy)
			throw runtime_error("Неверные размеры для set_copy.");
		copy(b.items, b.items + b._sx * b._sy, items);
	}

	static void mul(matrix& r, const matrix& a, const matrix& b)
	{
		if (r._sx != b._sx || r._sy != a._sy || a._sx != b._sy)
			throw runtime_error("Неверные размеры для mul.");

		for (size_t ry = 0; ry < r._sy; ry++)
			for (size_t rx = 0; rx < r._sx; rx++)
			{
				// скалярное произведение (dot) ry-й строки a на rx-й столбец b. rc — сокращённо row-column.
				Item& rcdot = r(rx, ry) = a(0, ry) * b(rx, 0);
				for (size_t irc = 1; irc < a._sx; irc++) rcdot += a(irc, ry) * b(rx, irc);
			}
	}

	matrix operator *(const matrix& b)
	{
		matrix r(b._sx, _sy);
		mul(r, *this, b);
		return r;
	}

	static void pow(matrix& r, const matrix& m, unsigned int p)
	{
		if (r._sx != m._sx || r._sy != m._sy || m._sx != m._sy)
			throw runtime_error("Неверные размеры для pow.");

		if (p == 0) { r.set_identity(); return; }
		if (p == 1) { r.set_copy(m); return; }

		matrix temp(m._sx, m._sy);
		matrix *ping = &temp, *pong = &r;
		bool first = true; // на первой итерации умножаем m на m
		int ibit = floor_log2(p);
		if ((__builtin_popcount(p) + ibit) % 2 != 0) swap(ping, pong);

		while (--ibit >= 0)
		{
			mul(*pong, first ? m : *ping, first ? m : *ping);
			swap(ping, pong);
			first = false;

			if ((p & (1 << ibit)) != 0)
			{
				mul(*pong, *ping, m);
				swap(ping, pong);
			}
		}

		if (ping != &r) throw runtime_error("Этого не должно происходить.");
	}

	matrix pow(unsigned int p)
	{
		matrix r(_sx, _sy);
		pow(r, *this, p);
		return r;
	}

private:
	size_t _sx, _sy;
	Item* items;
};

template<typename Item>
ostream& operator <<(ostream& strm, const matrix<Item>& m)
{
	for (size_t y = 0; y < m.sy(); y++)
		for (size_t x = 0; x < m.sx(); x++)
			(x > 0 ? strm << " " : (y > 0 ? strm << endl : strm)) << m(x, y);
	return strm;
}

int main()
{
	matrix<double> m(4, 4);
	for (int i = 0; i < 16; i++) m(i%4, i/4) = i%3;
	cout << "Исходная матрица m:" << endl;
	cout << m << endl << endl;

	cout << "m * m * ...6 раз... * m:" << endl;
	cout << m*m*m*m*m*m << endl << endl;

	cout << "m.pow(6):" << endl;
	cout << m.pow(6) << endl << endl;

	cout << "m * m * ...7 раз... * m:" << endl;
	cout << m*m*m*m*m*m*m << endl << endl;

	cout << "m.pow(7):" << endl;
	cout << m.pow(7) << endl;
	return 0;
}