#include<iostream>
#include<cstdlib>
#include<ctime>
class matrix;
class matrix
{
private:
static const int threshold;
static bool random;
bool erase;
int* mat;
int r_size, c_size;
int r, c;
void assign( const matrix& a, const int& rr, const int& cc, const int& k = 1 );
void sum( const matrix& a, const int& rr, const int& cc, const int& k = 1 );
void sum( const matrix& a, const matrix& b, const int& ak = 1, const int& bk = 1 );
void mul_classic( const matrix& a, const matrix& b );
void mul_strassen( const matrix& a, const matrix& b );
public:
matrix( const int& r_size, const int& c_size );
matrix( const int& r, const int& c, const int& top_num );
matrix( const int* m, const int& r, const int& c );
matrix( const matrix& m, const int& r_init, const int& c_init, const int& r, const int& c );
~matrix();
bool operator==( const matrix& a );
matrix operator*( const matrix& a );
matrix classic_mul( const matrix &a );
void show() const;
};
#define M 1000
#define N 1000
#define P 1000
int main()
{
matrix a( N, N, 2 );
// a.show();
// std::cout << std::endl;
matrix b( N, N, 2 );
// b.show();
// std::cout << std::endl;
matrix c = a * b;
std::cout << "Hecho! multiplicacion por Strassen." << std::endl;
// c.show();
// std::cout << std::endl;
matrix d = a.classic_mul( b );
std::cout << "Hecho! multiplicacion por la forma clasica." << std::endl;
// d.show();
// std::cout << std::endl;
return 0;
}
/* matrix */
bool matrix::random = false;
const int matrix::threshold = 48;
// asignar memoria para una matriz. Mucho cuidado, no se inicializan sus dimensiones.
matrix::matrix( const int& r_size, const int& c_size )
: r_size(r_size), c_size(c_size), erase(true)
{
mat = new int[ r_size*c_size ];
}
// crear matriz r x c con valores aleatorios entre -top_num y top_num
matrix::matrix( const int& r, const int& c, const int& top_num )
: r_size(r), c_size(c), r(r), c(c), erase(true)
{
mat = new int[ r_size*c_size ];
if( not random )
{
std::srand( time(NULL) );
random = true;
}
for( int i = 0; i < r*c; ++i )
mat[i] = -top_num + std::rand() % (2*top_num+1);
}
// crear matriz r x c a partir de un array
matrix::matrix( const int* m, const int& r, const int& c )
: r_size(r), c_size(c), r(r), c(c), erase(true)
{
mat = new int[ r * c ];
for( int i = 0; i < r*c; ++i )
mat[i] = m[i];
}
// crear submatriz en funcion de otra matriz. Coste O(1)
matrix::matrix( const matrix& m, const int& r_init, const int& c_init, const int& r, const int& c )
: mat(m.mat + r_init * m.c_size + c_init), r_size(m.r_size), c_size(m.c_size), r(r), c(c), erase(false)
{ }
matrix::~matrix()
{
if( erase ) delete [] mat;
}
matrix matrix::operator*( const matrix& a )
{
matrix b( r, a.c );
b.r = r;
b.c = a.c;
if( c != a.r )
std::cout << "error: las matrices no se pueden multiplicar" << std::endl;
else
b.mul_strassen( *this, a );
return b;
}
matrix matrix::classic_mul( const matrix& a )
{
matrix b( r, a.c );
b.r = r;
b.c = a.c;
if( c != a.r )
std::cout << "error: las matrices no se pueden multiplicar" << std::endl;
else
b.mul_classic( *this, a );
return b;
}
bool matrix::operator==( const matrix& a )
{
if( c != a.c || r != a.r )
return false;
int i;
for( i = 0; i < r*c && mat[i] == a.mat[i] ; ++i );
return i >= r*c;
}
// this = k * a( rr x cc )
void matrix::assign( const matrix& a, const int& rr, const int& cc, const int& k )
{
int i, j;
for( i = 0; i < rr; ++i )
for( j = 0; j < cc; ++j )
mat[ i * c_size + j ] = k * a.mat[ i * a.c_size + j ];
}
// this += k * a( rr x cc)
void matrix::sum( const matrix& a, const int& rr, const int& cc, const int& k )
{
int i, j;
for( i = 0; i < rr; ++i )
for( j = 0; j < cc; ++j )
mat[ i * c_size + j ] += k * a.mat[ i * a.c_size + j ];
}
// this = ak*a + bk*b
void matrix::sum( const matrix& a, const matrix& b, const int& ak, const int& bk )
{
if( erase )
{
r = a.r;
c = a.c;
assign( a, a.r, a.c, ak );
sum ( b, b.r, b.c, bk );
}
else
{
assign( a, r, c, ak );
sum ( b, r, c, bk );
}
}
// this = a x b
void matrix::mul_classic( const matrix& a, const matrix& b )
{
int i, j, k;
r = a.r;
c = b.c;
for( i = 0; i < r; ++i )
for( j = 0; j < c; ++j )
mat[ i * c_size + j ] = 0;
for( i = 0; i < r; ++i )
for( j = 0; j < c; ++j )
for( k = 0; k < b.r; ++k ) // cambio hace poco
mat[ i * c_size + j ] += a.mat[ i * a.c_size + k ] * b.mat[ k * b.c_size + j ];
}
void matrix::show() const
{
for( int i = 0; i < r; ++i )
{
for( int j = 0; j < c; ++j )
{
std::cout.width(3);
std::cout << mat[ i*c_size + j ];
}
std::cout << std::endl;
}
}
void matrix::mul_strassen( const matrix& a, const matrix& b )
{
if( a.r <= threshold ||
a.c <= threshold ||
b.r <= threshold ||
b.c <= threshold )
{
mul_classic( a, b );
return;
}
r = a.r;
c = b.c;
int arm = r / 2;
int acm = a.c / 2;
int brm = b.r / 2;
int bcm = c / 2;
matrix xa(r - arm, a.c - acm );
matrix ya(r - arm, a.c - acm );
matrix A00( a, 0, 0, arm , acm );
matrix A01( a, 0,acm, arm , a.c-acm );
matrix A10( a,arm, 0, r-arm, acm );
matrix A11( a,arm,acm, r-arm, a.c-acm );
matrix xb(b.r - brm, c - bcm );
matrix yb(b.r - brm, c - bcm );
matrix B00( b, 0, 0, brm , bcm );
matrix B01( b, 0,bcm, brm , c-bcm );
matrix B10( b,brm, 0, b.r-brm, bcm );
matrix B11( b,brm,bcm, b.r-brm, c-bcm );
matrix m( xa.r_size, xb.c_size );
matrix n( xa.r_size, xb.c_size );
matrix o( xa.r_size, xb.c_size );
matrix C00( *this, 0, 0, arm , bcm );
matrix C01( *this, 0,bcm, arm , c-bcm );
matrix C10( *this,arm, 0, r-arm, bcm );
matrix C11( *this,arm,bcm, r-arm, c-bcm );
// m = A00 * B00
m.mul_strassen( A00, B00 );
// n = A01 * B10
n.mul_strassen( A01, B10 );
// C00 = m + n
C00.sum( m, n );
// xa = -A10 + A00
xa.sum( A10, A00, -1, 1 );
// xb = B11 - B01
xb.sum( B11, B01, 1, -1 );
// ya = A11 - xa
ya.sum( A11, xa, 1, -1 );
// yb = xb + B00
yb.sum( xb, B00 );
// o = ya * yb
o.mul_strassen( ya, yb );
// o += m
o.sum( m, m.r, m.c );
// redimensionar xa en caso necesario
if( xa.c < xb.r ){
for( int i = 0; i < xa.r; ++i )
xa.mat[ i*xa.c_size + xa.c ] = 0;
++xa.c;
}
// n = xa * xb
n.mul_strassen( xa, xb );
// n += o
n.sum( o, o.r, o.c );
// xa = A11 + A10
xa.sum( A11, A10 );
// xb = B10 - B00
xb.sum( B01, B00, 1, -1 );
// redimensionar xb en caso necesario
if( xb.r < xa.c ){
for( int j = 0; j < xb.c; ++j )
xb.mat[ xb.r*xb.c_size + j ] = 0;
++xb.r;
}
// m = xa * xb
m.mul_strassen( xa, xb ); // m' x m
// C11 = n + m
C11.sum( n, m );
// m += o
m.sum( o, o.r, o.c ); // m x m'''
// xa = -ya + A01;
xa.sum( ya, A01, -1, 1 );
// o = xa * B11
o.mul_strassen( xa, B11 );
// C01 = m + o
C01.sum( m, o );
// xb = yb - B10;
xb.sum( yb, B10, 1, -1 );
// o = A11 * xb
o.mul_strassen( A11, xb );
// C10 = n - o
C10.sum( n, o, 1, -1 );
} |
F.A.Q.
Conéctese para revisar sus mensajes
.