fork download
copy 
  1. #include <iostream>
  2. #include <algorithm>
  3. #include <climits>
  4. using namespace std;
  5.  
  6. // Индекс самого старшего выставленного бита, считая младший бит за №0:
  7. // floor_log2(0b101000) == 5, floor_log2(1 << N) == N.
  8. //              543210
  9. int floor_log2(unsigned int p)
  10. {
  11. 	return CHAR_BIT * sizeof(unsigned int) - 1 - __builtin_clz(p);
  12. }
  13.  
  14. template<typename Item>
  15. class matrix
  16. {
  17. public:
  18. 	matrix(size_t sx, size_t sy)
  19. 		: _sx(sx), _sy(sy), items(nullptr)
  20. 	{
  21. 		items = new Item[sx * sy];
  22. 		set_identity();
  23. 	}
  24.  
  25. 	~matrix() { delete[] items; }
  26. 	size_t sx() const { return _sx; }
  27. 	size_t sy() const { return _sy; }
  28.  
  29. 	const Item& operator()(size_t x, size_t y) const { return items[y * _sx + x]; }
  30. 	Item& operator()(size_t x, size_t y) { return items[y * _sx + x]; }
  31.  
  32. 	void set_identity()
  33. 	{
  34. 		fill(items, items + _sx * _sy, 0);
  35. 		for (size_t xy = 0; xy < min(_sx, _sy); xy++) (*this)(xy, xy) = 1;
  36. 	}
  37.  
  38. 	void set_copy(const matrix& b)
  39. 	{
  40. 		if (_sx != b._sx || _sy != b._sy)
  41. 			throw runtime_error("Неверные размеры для set_copy.");
  42. 		copy(b.items, b.items + b._sx * b._sy, items);
  43. 	}
  44.  
  45. 	static void mul(matrix& r, const matrix& a, const matrix& b)
  46. 	{
  47. 		if (r._sx != b._sx || r._sy != a._sy || a._sx != b._sy)
  48. 			throw runtime_error("Неверные размеры для mul.");
  49.  
  50. 		for (size_t ry = 0; ry < r._sy; ry++)
  51. 			for (size_t rx = 0; rx < r._sx; rx++)
  52. 			{
  53. 				// скалярное произведение (dot) ry-й строки a на rx-й столбец b. rc — сокращённо row-column.
  54. 				Item& rcdot = r(rx, ry) = a(0, ry) * b(rx, 0);
  55. 				for (size_t irc = 1; irc < a._sx; irc++) rcdot += a(irc, ry) * b(rx, irc);
  56. 			}
  57. 	}
  58.  
  59. 	matrix operator *(const matrix& b)
  60. 	{
  61. 		matrix r(b._sx, _sy);
  62. 		mul(r, *this, b);
  63. 		return r;
  64. 	}
  65.  
  66. 	static void pow(matrix& r, const matrix& m, unsigned int p)
  67. 	{
  68. 		if (r._sx != m._sx || r._sy != m._sy || m._sx != m._sy)
  69. 			throw runtime_error("Неверные размеры для pow.");
  70.  
  71. 		if (p == 0) { r.set_identity(); return; }
  72. 		if (p == 1) { r.set_copy(m); return; }
  73.  
  74. 		matrix temp(m._sx, m._sy);
  75. 		matrix *ping = &temp, *pong = &r;
  76. 		bool first = true; // на первой итерации умножаем m на m
  77. 		int ibit = floor_log2(p);
  78. 		if ((__builtin_popcount(p) + ibit) % 2 != 0) swap(ping, pong);
  79.  
  80. 		while (--ibit >= 0)
  81. 		{
  82. 			mul(*pong, first ? m : *ping, first ? m : *ping);
  83. 			swap(ping, pong);
  84. 			first = false;
  85.  
  86. 			if ((p & (1 << ibit)) != 0)
  87. 			{
  88. 				mul(*pong, *ping, m);
  89. 				swap(ping, pong);
  90. 			}
  91. 		}
  92.  
  93. 		if (ping != &r) throw runtime_error("Этого не должно происходить.");
  94. 	}
  95.  
  96. 	matrix pow(unsigned int p)
  97. 	{
  98. 		matrix r(_sx, _sy);
  99. 		pow(r, *this, p);
  100. 		return r;
  101. 	}
  102.  
  103. private:
  104. 	size_t _sx, _sy;
  105. 	Item* items;
  106. };
  107.  
  108. template<typename Item>
  109. ostream& operator <<(ostream& strm, const matrix<Item>& m)
  110. {
  111. 	for (size_t y = 0; y < m.sy(); y++)
  112. 		for (size_t x = 0; x < m.sx(); x++)
  113. 			(x > 0 ? strm << " " : (y > 0 ? strm << endl : strm)) << m(x, y);
  114. 	return strm;
  115. }
  116.  
  117. int main()
  118. {
  119. 	matrix<double> m(4, 4);
  120. 	for (int i = 0; i < 16; i++) m(i%4, i/4) = i%3;
  121. 	cout << "Исходная матрица m:" << endl;
  122. 	cout << m << endl << endl;
  123.  
  124. 	cout << "m * m * ...6 раз... * m:" << endl;
  125. 	cout << m*m*m*m*m*m << endl << endl;
  126.  
  127. 	cout << "m.pow(6):" << endl;
  128. 	cout << m.pow(6) << endl << endl;
  129.  
  130. 	cout << "m * m * ...7 раз... * m:" << endl;
  131. 	cout << m*m*m*m*m*m*m << endl << endl;
  132.  
  133. 	cout << "m.pow(7):" << endl;
  134. 	cout << m.pow(7) << endl;
  135. 	return 0;
  136. }
Success #stdin #stdout 0s 4812KB
comments ()
stdin
copy 
Standard input is empty
stdout
copy 
Исходная матрица m:
0 1 2 0
1 2 0 1
2 0 1 2
0 1 2 0

m * m * ...6 раз... * m:
752 656 806 752
656 792 1012 656
806 1012 1401 806
752 656 806 752

m.pow(6):
752 656 806 752
656 792 1012 656
806 1012 1401 806
752 656 806 752

m * m * ...7 раз... * m:
2268 2816 3814 2268
2816 2896 3636 2816
3814 3636 4625 3814
2268 2816 3814 2268

m.pow(7):
2268 2816 3814 2268
2816 2896 3636 2816
3814 3636 4625 3814
2268 2816 3814 2268
https://ideone.com/wKp7YY
language:
C++14 (gcc 8.3)
created:
4 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!