#include<iostream.h>
#include<stdlib.h>
#define BREAK 2
#define a11 a->p[0]
#define a12 a->p[1]
#define a21 a->p[2]
#define a22 a->p[3]
#define b11 b->p[0]
#define b12 b->p[1]
#define b21 b->p[2]
#define b22 b->p[3]
#define c11 c->p[0]
#define c12 c->p[1]
#define c21 c->p[2]
#define c22 c->p[3]
#define d11 d->p[0]
#define d12 d->p[1]
#define d21 d->p[2]
#define d22 d->p[3]
typedef double **matrix;
typedef union _strassen_matrix
{
matrix d;
union _strassen_matrix **p;
} *strassen_matrix;
matrix new_matrix(int);
strassen_matrix new_strassen(int);
void normal_to_strassen(matrix, strassen_matrix, int);
void strassen_to_normal(strassen_matrix, matrix, int);
matrix strassen_submatrix(strassen_matrix, int, int, int);
void copy_matrix(matrix, matrix, int);
void add_matrix(matrix, matrix, matrix, int);
void sub_matrix(matrix, matrix, matrix, int);
void copy_strassen(strassen_matrix, strassen_matrix, int);
void add_strassen(strassen_matrix, strassen_matrix, strassen_matrix, int);
void sub_strassen(strassen_matrix, strassen_matrix, strassen_matrix, int);
void mul_matrix(matrix, matrix, matrix, int);
void mul_strassen(strassen_matrix, strassen_matrix, strassen_matrix, strassen_matrix, int);
void print_matrix(matrix, int);
int least_power_of_two(int);
matrix new_matrix(int n)
{
matrix a = (matrix) malloc(sizeof(double *) * n);
for (int j = 0; j < n; j++)
a[j] = (double *) malloc(sizeof(double) * n);
return(a);
}
strassen_matrix new_strassen(int n)
{
strassen_matrix a;
a = (strassen_matrix)malloc(sizeof(*a));
if (n <= BREAK)
a->d = (matrix ) new_matrix(n);
else
{
register int m = n/2;
a->p = (strassen_matrix *)malloc(4*sizeof(strassen_matrix));
a11 = new_strassen(m);
a12 = new_strassen(m);
a21 = new_strassen(m);
a22 = new_strassen(m);
}
return a;
}
matrix strassen_submatrix(strassen_matrix a,
int i,
int j,
int n
)
{
if (n <= BREAK)
return(a->d);
else
{
int cur_bit, bit_num;
strassen_matrix cur_ptr = a;
bit_num = least_power_of_two(n)-1;
cur_bit = n/2;
while (cur_bit >= BREAK)
{
cur_ptr = cur_ptr->p[(((j & cur_bit) | ((i & cur_bit)*2)) >> bit_num)];
cur_bit >>= 1;
bit_num--;
}
return (cur_ptr->d);
}
}
void normal_to_strassen(matrix a,strassen_matrix b,int n)
{
if (n <= BREAK)
copy_matrix(a,b->d,n);
else
{
int i,j,ii,jj;
matrix sub;
for (i=0; i<n; i += BREAK)
{
for (j=0; j<n; j += BREAK)
{
sub = strassen_submatrix(b,i,j,n);
for (ii=0; ii<BREAK; ii++)
for (jj=0; jj<BREAK; jj++)
sub[ii][jj] = a[i+ii][j+jj];
}
}
}
}
void strassen_to_normal(strassen_matrix a,matrix b,int n)
{
if (n <= BREAK)
copy_matrix(a->d,b, n);
else
{
matrix sub;
for (int i=0; i<n; i += BREAK)
{
for (int j=0; j<n; j += BREAK)
{
sub = strassen_submatrix(a,i,j,n);
for (int ii=0; ii<BREAK; ii++)
for (int jj=0; jj<BREAK; jj++)
b[i+ii][j+jj] = sub[ii][jj];
}
}
}
}
void copy_matrix(
matrix a,
matrix b,
int n
)
{
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
b[i][j] = a[i][j];
}
void add_matrix(
matrix a,
matrix b,
matrix c,
int n
)
{
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
c[i][j] = b[i][j] + a[i][j];
}
void sub_matrix(
matrix a,
matrix b,
matrix c,
int n
)
{
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
c[i][j] = a[i][j] - b[i][j];
}
void add_strassen(
strassen_matrix a,
strassen_matrix b,
strassen_matrix c,
int n
)
{
if (n <= BREAK)
add_matrix(a->d, b->d, c->d, n);
else
{
int m=n/2;
add_strassen(a11, b11, c11, m);
add_strassen(a12, b12, c12, m);
add_strassen(a21, b21, c21, m);
add_strassen(a22, b22, c22, m);
}
}
void sub_strassen(
strassen_matrix a,
strassen_matrix b,
strassen_matrix c,
int n
)
{
if (n <= BREAK)
{
sub_matrix(a->d, b->d, c->d, n);
}
else
{
int m = n/2;
sub_strassen(a11, b11, c11, m);
sub_strassen(a12, b12, c12, m);
sub_strassen(a21, b21, c21, m);
sub_strassen(a22, b22, c22, m);
}
}
void mul_matrix(
matrix a,
matrix b,
matrix c,
int n
)
{
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
{
c[i][j] = 0.0;
for(int k=0; k<n; k++)
c[i][j] += a[i][k] * b[k][j];
}
}
void mul_strassen(
strassen_matrix a,
strassen_matrix b,
strassen_matrix c,
strassen_matrix d,
int n
)
{
if (n <= BREAK)
mul_matrix(a->d,b->d,c->d,n);
else
{
int m = n/2;
sub_strassen(a12, a22, d11, m);
add_strassen(b21, b22, d12, m);
mul_strassen(d11, d12, c11, d21, m);
sub_strassen(a21, a11, d11, m);
add_strassen(b11, b12, d12, m);
mul_strassen(d11, d12, c22, d21, m);
add_strassen(a11, a12, d11, m);
mul_strassen(d11, b22, c12, d12, m);
sub_strassen(c11, c12, c11, m);
sub_strassen(b21, b11, d11, m);
mul_strassen(a22, d11, c21, d12, m);
add_strassen(c21, c11, c11, m);
sub_strassen(b12, b22, d11, m);
mul_strassen(a11, d11, d12, d21, m);
add_strassen(d12, c12, c12, m);
add_strassen(d12, c22, c22, m);
add_strassen(a21, a22, d11, m);
mul_strassen(d11, b11, d12, d21, m);
add_strassen(d12, c21, c21, m);
sub_strassen(c22, d12, c22, m);
add_strassen(a11, a22, d11, m);
add_strassen(b11, b22, d12, m);
mul_strassen(d11, d12, d21, d22, m);
add_strassen(d21, c11, c11, m);
add_strassen(d21, c22, c22, m);
}
}
void print_matrix(matrix a,int n)
{
for(int i=0;i<n;i++){
for(int j=0;j<n;j++)
cout<<a[i][j]<<"\t";
cout<<endl;
}
}
int least_power_of_two(int n )
{
int i = 1, k = 1;
if (n==1)
return (0);
while ((k <<= 1) < n)
i++;
return(i);
}
void readMatrix(matrix a,int n)
{
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
cin>>a[i][j];
}
void main()
{
int n;