Example: Strassen's Matrix Multiplication in C++

// Strassen's Matrix Multiplication
#include <iostream>
#include <vector>
#include <cmath> // For std::pow, std::ceil, std::log2

using namespace std;

// Function to print a matrix
void printMatrix(const vector<vector<int>>& matrix) {
    for (const auto& row : matrix) {
        for (int val : row) {
            cout << val << "\t";
        }
        cout << endl;
    }
}

// Function to add two matrices
vector<vector<int>> addMatrices(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> C(n, vector<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];
        }
    }
    return C;
}

// Function to subtract two matrices
vector<vector<int>> subtractMatrices(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> C(n, vector<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];
        }
    }
    return C;
}

// Strassen's matrix multiplication function
vector<vector<int>> strassenMultiply(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();

    // Base case: If matrix is 1x1, perform scalar multiplication
    if (n == 1) {
        return {{A[0][0] * B[0][0]}};
    }

    // Divide matrices into 4 sub-matrices
    int k = n / 2;
    vector<vector<int>> A11(k, vector<int>(k));
    vector<vector<int>> A12(k, vector<int>(k));
    vector<vector<int>> A21(k, vector<int>(k));
    vector<vector<int>> A22(k, vector<int>(k));

    vector<vector<int>> B11(k, vector<int>(k));
    vector<vector<int>> B12(k, vector<int>(k));
    vector<vector<int>> B21(k, vector<int>(k));
    vector<vector<int>> B22(k, vector<int>(k));

    for (int i = 0; i < k; ++i) {
        for (int j = 0; j < k; ++j) {
            A11[i][j] = A[i][j];
            A12[i][j] = A[i][j + k];
            A21[i][j] = A[i + k][j];
            A22[i][j] = A[i + k][j + k];

            B11[i][j] = B[i][j];
            B12[i][j] = B[i][j + k];
            B21[i][j] = B[i + k][j];
            B22[i][j] = B[i + k][j + k];
        }
    }

    // Calculate 7 products (P1 to P7) recursively
    vector<vector<int>> P1 = strassenMultiply(addMatrices(A11, A22), addMatrices(B11, B22));
    vector<vector<int>> P2 = strassenMultiply(addMatrices(A21, A22), B11);
    vector<vector<int>> P3 = strassenMultiply(A11, subtractMatrices(B12, B22));
    vector<vector<int>> P4 = strassenMultiply(A22, subtractMatrices(B21, B11));
    vector<vector<int>> P5 = strassenMultiply(addMatrices(A11, A12), B22);
    vector<vector<int>> P6 = strassenMultiply(subtractMatrices(A21, A11), addMatrices(B11, B12));
    vector<vector<int>> P7 = strassenMultiply(subtractMatrices(A12, A22), addMatrices(B21, B22));

    // Combine products to form the result sub-matrices
    vector<vector<int>> C11 = addMatrices(subtractMatrices(addMatrices(P1, P4), P5), P7);
    vector<vector<int>> C12 = addMatrices(P3, P5);
    vector<vector<int>> C21 = addMatrices(P2, P4);
    vector<vector<int>> C22 = addMatrices(subtractMatrices(addMatrices(P1, P3), P2), P6);

    // Assemble the final result matrix C
    vector<vector<int>> C(n, vector<int>(n));
    for (int i = 0; i < k; ++i) {
        for (int j = 0; j < k; ++j) {
            C[i][j] = C11[i][j];
            C[i][j + k] = C12[i][j];
            C[i + k][j] = C21[i][j];
            C[i + k][j + k] = C22[i][j];
        }
    }

    return C;
}

// Function to get the next power of 2
int getNextPowerOf2(int n) {
    if (n == 0) return 1;
    return pow(2, ceil(log2(n)));
}

int main() {
    // Step 1: Define input matrices
    vector<vector<int>> A_orig = {{1, 2}, {3, 4}};
    vector<vector<int>> B_orig = {{5, 6}, {7, 8}};

    // vector<vector<int>> A_orig = {{1, 1, 1, 1},
    //                               {2, 2, 2, 2},
    //                               {3, 3, 3, 3},
    //                               {4, 4, 4, 4}};
    // vector<vector<int>> B_orig = {{1, 1, 1, 1},
    //                               {2, 2, 2, 2},
    //                               {3, 3, 3, 3},
    //                               {4, 4, 4, 4}};

    cout << "Matrix A:" << endl;
    printMatrix(A_orig);
    cout << "\nMatrix B:" << endl;
    printMatrix(B_orig);

    // Step 2: Determine actual size and padded size
    int n_orig = A_orig.size();
    int padded_size = getNextPowerOf2(n_orig);

    // Step 3: Pad matrices with zeros if necessary
    vector<vector<int>> A_padded(padded_size, vector<int>(padded_size, 0));
    vector<vector<int>> B_padded(padded_size, vector<int>(padded_size, 0));

    for (int i = 0; i < n_orig; ++i) {
        for (int j = 0; j < n_orig; ++j) {
            A_padded[i][j] = A_orig[i][j];
            B_padded[i][j] = B_orig[i][j];
        }
    }

    // Step 4: Perform Strassen's multiplication on padded matrices
    vector<vector<int>> C_padded = strassenMultiply(A_padded, B_padded);

    // Step 5: Extract the original sized result from the padded result
    vector<vector<int>> C_result(n_orig, vector<int>(n_orig));
    for (int i = 0; i < n_orig; ++i) {
        for (int j = 0; j < n_orig; ++j) {
            C_result[i][j] = C_padded[i][j];
        }
    }

    cout << "\nResult C (using Strassen's algorithm):" << endl;
    printMatrix(C_result);

    return 0;
}