Raising a square array to the power of N in C

Question

I'm writing a C program to compute a square array raised to the power of N.

I'd like to write these in separate functions to keep it clean. I have a function

arrMult(int M, int R[M][M], int X[M][M], int out[M][M])

which essentially does out = R*X. To raise it to some Nth power, I wrote this function:

int arrNthPow(int M, int N, int R[M][M], int out[M][M])
{
    // Initialize
    arrPrint(M, M, out);
    arrMult(M, R, R, out);
    for (int i = 1; i < N - 1; i++){
        arrPrint(M, M, out);
        arrMult(M, R, out, out);
    }
    return 0;
}

However, I realise that this doesn't give me the result I expect. Instead, if I do

arrMult(M, R, R, out1);
arrMult(M, out1, R, out2);
arrMult(M, out2, R, out3);

out3 gets me the answer I'm looking for. I'm assuming that I have to somehow use pointers to copy the value of each array output after arrMult before using arrMult again for the next iteration but I'm not sure how.

I'm not very familiar with pointers however and have no clue on how to write this without copying the values to a temp variable and using additional for loops.

Answer 1

There are multiple problems in your code:

• there is an extra argument in arrMult(int M, int M, int R[M][M], int X[M][M], int out[M][M])
• the out matrix should be set to the unity matrix if N == 0.
• the out matrix should receive a copy of R if N == 1.
• otherwise out starts at R², but it is unlikely you can pass the same matrix as input and output to arrMult. You should use a temporary matrix to store the result and copy it to out either inside arrMult or after each matrix multiplication.
• there is a more efficient algorithm to compute powers in log N steps rather than N steps as posted: https://en.wikipedia.org/wiki/Exponentiation_by_squaring

Here is a simple implementation:

void arrIdentity(int M, int out[M][M]) {
    /* initialize array to identity matrix. */
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < M; j++) {
            out[i][j] = (i == j);
        }
    }
}    

void arrCopy(int M, const int mat[M][M], int out[M][M]) {
    /* copy a matrix. */
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < M; j++) {
            out[i][j] = mat[i][j];
        }
    }
}    

void arrMult(int M, const int R[M][M], const int X[M][M], int out[M][M]) {
    /* compute matrix multiplication: out = R * X. */
    int temp[M][M];

    for (int i = 0; i < M; i++) {
        for (int j = 0; j < M; j++) {
            int sum = 0;
            for (int k = 0; k < M; k++) {
                sum += R[i][k] * X[k][j];
            }
            temp[i][j] = sum;
        }
        arrCopy(M, temp, out);
    }
}

int arrNthPow(int M, int N, const int R[M][M], int out[M][M]) {
    /* Naive matrix exponentiation */
    if (N < 0)
        return -1;
    if (N == 0) {
        arrIdentity(M, out);
        return 0;
    } else {
        arrCopy(M, R, out);
        for (int i = 1; i < N; i++) {
            arrMult(M, R, out, out);
        }
    }
    return 0;
}

Here is a more efficient alternative using the method explained in the Wikipedia article:

int arrNthPow2(int M, int N, const int R[M][M], int out[M][M]) {
    /* Matrix exponentiation by squaring */
    int X[M][M];

    if (N < 0)
        return -1;
    if (N == 0) {
        arrIdentity(M, out);
        return 0;
    }
    if (N == 1) {
        arrCopy(M, R, out);
        return 0;
    }
    if (N == 2) {
        arrMult(M, R, R, out);
        return 0;
    }
    arrIdentity(M, out);
    arrCopy(M, R, X);
    while (N > 1) {
        if (N & 1)
            arrMult(M, X, out, out);
        arrMult(M, X, X, X);
        N >>= 1;
    }
    arrMult(M, X, out, out);
    return 0;
}