/*用c語言中的函式遞迴呼叫演算法實現n階矩陣的n次冪*/

#include

#include

#include

#include

//建立矩陣，矩陣用一維陣列儲存

double *matCreate(unsigned int m, unsigned int n)

{

double *p = (double *)malloc(sizeof(double) * m * n);

if (p == NULL) printf("建立矩陣失敗！\n");

return p;

}

//輸入矩陣元素

void matInput(double *a, unsigned int m, unsigned int n)

{

for (int i = 0; i

{

for (int j = 0; j

{

scanf("%f ", &a[i * n + j]);

}

}

return;

}

//隨機產生矩陣元素，均勻分佈於[from to]

void matInitRand(double *a, unsigned int m, unsigned int n, double from, double to)

{

if (a == NULL || m

double x;

srand(time(NULL));

for (int i = 0; i

{

for (int j = 0; j

{

x = (1.0 * rand() / RAND_MAX) * (to - from) + from;

a[i * n + j] = x;

}

}

return;

}

//轉置

void matTranspose(double *a, double *b, unsigned int m, unsigned int n)

{

for (int i = 0; i

{

for (int j = 0; j

{

b[j*n +i]=a[i * n + j] ;

}

}

}

//輸出矩陣

void matPrint(double *a, unsigned int m, unsigned int n)

{

for (int i = 0; i

{

for (int j = 0; j

{

printf("%8.4f ", a[i * n + j]);

}

putchar("\n");

}

return;

}

//矩陣乘法c=a*b

void matMul(double *a, double *b, double *c, unsigned int m, unsigned int n, unsigned int k)

{

if (a == NULL || b == NULL || c == NULL || m

double x = 0.0f;

for (int i = 0; i

{

for (int u = 0; u

{

x = 0.0f;

for (int j = 0; j

{

x += a[i * n + j] * b[j * k + u];

}

c[i * k + u] = x;

}

}

return;

}

//b=a^n, a:m*m階矩陣

void matFac(double *a, double *b, unsigned int n, unsigned int m)

{

double *c = (double *)malloc(sizeof(double) * m * m); //儲存臨時結果

if (n > 1)

{

matFac(a, c, n - 1, m);

matMul(a, c, b, m, m, m);

}

else

memcpy(b, a, sizeof(double)*m * m);

// printf("%d:\n",n);

// matPrint(b, m,m);

free(c); //回收記憶體

return ;

}

#define M 3

#define N 4

#define K N

int main(int argc, char const *argv[])

{

double *A, *B, *B1,*BT, *C;

A = matCreate(M, N);

B = matCreate(N, K);

B1 = matCreate(N, K);

BT = matCreate(K,N);

C = matCreate(M, K);

if (!A || !B || !B1 || !BT || !C) return -1;

matInitRand(A, M, N, 0.0f, 1.0f);

printf("A=\n");

matPrint(A, M, N);

matInitRand(B, N, K, 0.0f, 1.0f);

printf("B=\n");

matPrint(B, N, K);

matTranspose(B,BT,N,K);

printf("B"=\n");

matPrint(BT, K,N);

matMul(A, B, C, M, N, K);

printf("C=A*B\n");

matPrint(C, M, N);

matFac(B, B1, 4, N);

printf("B^4\n");

matPrint(B1, N, K);

return 0;

}