1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
|
#include <iostream>
#include <iomanip>
#include <vector>
#include <random>
#include <cassert>
#include <cmath>
#include <ctime>
using namespace std;
using vec = vector<double>;
using matrix = vector<vec>;
//==========================================================
mt19937 gen( time(0) );
uniform_real_distribution<double> uni( 0.0, 1.0 );
double RNG( double a, double b ){ return a + ( b - a ) * uni( gen ); }
//==========================================================
ostream & operator << ( ostream &out, const matrix &M )
{
for ( auto &row : M )
{
for ( auto e : row ) out << setw( 15 ) << e << " ";
out << '\n';
}
return out;
}
//==========================================================
matrix operator *( const matrix& A, const matrix &B )
{
assert( A[0].size() == B.size() ); // check for compatibility
matrix M( A.size(), vec( B[0].size(), 0.0 ) );
for ( int i = 0; i < A.size(); i++ )
{
for ( int j = 0; j < B[0].size(); j++ )
{
for ( int k = 0; k < A[0].size(); k++ ) M[i][j] += A[i][k] * B[k][j];
}
}
return M;
}
//==========================================================
matrix identity( unsigned n )
{
matrix M( n, vec( n, 0.0 ) );
for ( int i = 0; i < n; i++ ) M[i][i] = 1;
return M;
}
//==========================================================
matrix randMatrix( unsigned n, double a, double b )
{
matrix M( n, vec( n ) );
for ( auto &row : M )
for ( auto &e : row ) e = RNG( a, b );
return M;
}
//==========================================================
matrix pow( const matrix& A, unsigned n ) // cheap and cheerful by repeated multiply (quicker by powers of 2)
{
assert( A.size() == A[0].size() ); // check square matrix
if ( n == 0 ) return identity( A.size() );
if ( n == 1 ) return A;
return A * pow( A, n - 1 );
}
//==========================================================
int main()
{
matrix M = randMatrix( 3, -10.0, 10.0 );
cout << "M=\n" << M << "\n\n";
for ( unsigned i = 0; i < 5; i++ ) cout << "pow(M," << i << ") = \n" << pow( M, i ) << "\n\n";
}
|