Example: Strassen's Matrix Multiplication in C

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

// Helper function to print a matrix
void printMatrix(const std::vector<std::vector<int>>& matrix, int n) {
    // Step 1: Iterate through each row up to the actual size 'n'.
    for (int i = 0; i < n; ++i) {
        // Step 2: Iterate through each element in the current row up to 'n'.
        for (int j = 0; j < n; ++j) {
            // Step 3: Print the element followed by a space.
            std::cout << matrix[i][j] << " ";
        }
        // Step 4: Move to the next line after printing all elements in a row.
        std::cout << std::endl;
    }
}

// Helper function to add two matrices
std::vector<std::vector<int>> add(const std::vector<std::vector<int>>& A,
                                   const std::vector<std::vector<int>>& B, int size) {
    // Step 1: Initialize result matrix C with dimensions 'size' x 'size'.
    std::vector<std::vector<int>> C(size, std::vector<int>(size));
    // Step 2: Iterate through rows and columns.
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j) {
            // Step 3: Perform element-wise addition.
            C[i][j] = A[i][j] + B[i][j];
        }
    }
    // Step 4: Return the sum matrix.
    return C;
}

// Helper function to subtract two matrices
std::vector<std::vector<int>> subtract(const std::vector<std::vector<int>>& A,
                                       const std::vector<std::vector<int>>& B, int size) {
    // Step 1: Initialize result matrix C with dimensions 'size' x 'size'.
    std::vector<std::vector<int>> C(size, std::vector<int>(size));
    // Step 2: Iterate through rows and columns.
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j) {
            // Step 3: Perform element-wise subtraction.
            C[i][j] = A[i][j] - B[i][j];
        }
    }
    // Step 4: Return the difference matrix.
    return C;
}

// Helper function to split a matrix into four sub-matrices
void split(const std::vector<std::vector<int>>& P,
           std::vector<std::vector<int>>& P11, std::vector<std::vector<int>>& P12,
           std::vector<std::vector<int>>& P21, std::vector<std::vector<int>>& P22, int size) {
    int halfSize = size / 2;
    // Step 1: Iterate through the rows of the parent matrix P.
    for (int i = 0; i < halfSize; ++i) {
        // Step 2: Iterate through the columns of the parent matrix P.
        for (int j = 0; j < halfSize; ++j) {
            // Step 3: Assign elements to the top-left sub-matrix P11.
            P11[i][j] = P[i][j];
            // Step 4: Assign elements to the top-right sub-matrix P12.
            P12[i][j] = P[i][j + halfSize];
            // Step 5: Assign elements to the bottom-left sub-matrix P21.
            P21[i][j] = P[i + halfSize][j];
            // Step 6: Assign elements to the bottom-right sub-matrix P22.
            P22[i][j] = P[i + halfSize][j + halfSize];
        }
    }
}

// Helper function to join four sub-matrices into a single matrix
void join(const std::vector<std::vector<int>>& C11, const std::vector<std::vector<int>>& C12,
          const std::vector<std::vector<int>>& C21, const std::vector<std::vector<int>>& C22,
          std::vector<std::vector<int>>& C, int size) {
    int halfSize = size / 2;
    // Step 1: Iterate through the rows of the result matrix C.
    for (int i = 0; i < halfSize; ++i) {
        // Step 2: Iterate through the columns of the result matrix C.
        for (int j = 0; j < halfSize; ++j) {
            // Step 3: Assign C11 elements to the top-left quadrant of C.
            C[i][j] = C11[i][j];
            // Step 4: Assign C12 elements to the top-right quadrant of C.
            C[i][j + halfSize] = C12[i][j];
            // Step 5: Assign C21 elements to the bottom-left quadrant of C.
            C[i + halfSize][j] = C21[i][j];
            // Step 6: Assign C22 elements to the bottom-right quadrant of C.
            C[i + halfSize][j + halfSize] = C22[i][j];
        }
    }
}

// Recursive Strassen's Matrix Multiplication
std::vector<std::vector<int>> strassen_multiply(const std::vector<std::vector<int>>& A,
                                                 const std::vector<std::vector<int>>& B, int size) {
    // Step 1: Base case for recursion: if matrix size is 1x1, perform direct multiplication.
    if (size == 1) {
        std::vector<std::vector<int>> C(1, std::vector<int>(1));
        C[0][0] = A[0][0] * B[0][0];
        return C;
    }

    // Step 2: Calculate half the current matrix size.
    int halfSize = size / 2;

    // Step 3: Initialize sub-matrices for A and B.
    std::vector<std::vector<int>> A11(halfSize, std::vector<int>(halfSize));
    std::vector<std::vector<int>> A12(halfSize, std::vector<int>(halfSize));
    std::vector<std::vector<int>> A21(halfSize, std::vector<int>(halfSize));
    std::vector<std::vector<int>> A22(halfSize, std::vector<int>(halfSize));

    std::vector<std::vector<int>> B11(halfSize, std::vector<int>(halfSize));
    std::vector<std::vector<int>> B12(halfSize, std::vector<int>(halfSize));
    std::vector<std::vector<int>> B21(halfSize, std::vector<int>(halfSize));
    std::vector<std::vector<int>> B22(halfSize, std::vector<int>(halfSize));

    // Step 4: Split A and B into their respective sub-matrices.
    split(A, A11, A12, A21, A22, size);
    split(B, B11, B12, B21, B22, size);

    // Step 5: Compute the 7 products (M1 to M7) recursively.
    // M1 = (A11 + A22) * (B11 + B22)
    std::vector<std::vector<int>> M1 = strassen_multiply(add(A11, A22, halfSize), add(B11, B22, halfSize), halfSize);
    // M2 = (A21 + A22) * B11
    std::vector<std::vector<int>> M2 = strassen_multiply(add(A21, A22, halfSize), B11, halfSize);
    // M3 = A11 * (B12 - B22)
    std::vector<std::vector<int>> M3 = strassen_multiply(A11, subtract(B12, B22, halfSize), halfSize);
    // M4 = A22 * (B21 - B11)
    std::vector<std::vector<int>> M4 = strassen_multiply(A22, subtract(B21, B11, halfSize), halfSize);
    // M5 = (A11 + A12) * B22
    std::vector<std::vector<int>> M5 = strassen_multiply(add(A11, A12, halfSize), B22, halfSize);
    // M6 = (A21 - A11) * (B11 + B12)
    std::vector<std::vector<int>> M6 = strassen_multiply(subtract(A21, A11, halfSize), add(B11, B12, halfSize), halfSize);
    // M7 = (A12 - A22) * (B21 + B22)
    std::vector<std::vector<int>> M7 = strassen_multiply(subtract(A12, A22, halfSize), add(B21, B22, halfSize), halfSize);

    // Step 6: Compute the four sub-matrices of the result C (C11, C12, C21, C22).
    // C11 = M1 + M4 - M5 + M7
    std::vector<std::vector<int>> C11 = add(subtract(add(M1, M4, halfSize), M5, halfSize), M7, halfSize);
    // C12 = M3 + M5
    std::vector<std::vector<int>> C12 = add(M3, M5, halfSize);
    // C21 = M2 + M4
    std::vector<std::vector<int>> C21 = add(M2, M4, halfSize);
    // C22 = M1 - M2 + M3 + M6
    std::vector<std::vector<int>> C22 = add(subtract(add(M1, M3, halfSize), M2, halfSize), M6, halfSize);

    // Step 7: Initialize the final result matrix C.
    std::vector<std::vector<int>> C(size, std::vector<int>(size));
    // Step 8: Join the four computed sub-matrices into C.
    join(C11, C12, C21, C22, C, size);

    // Step 9: Return the final product matrix C.
    return C;
}

// Function to get the next power of 2 for padding
int getNextPowerOf2(int n) {
    if (n == 0) return 1;
    if ((n & (n - 1)) == 0) return n; // Already a power of 2
    return std::pow(2, std::ceil(std::log2(n)));
}

// Wrapper function to handle padding for Strassen's algorithm
std::vector<std::vector<int>> strassen_wrapper(const std::vector<std::vector<int>>& A,
                                                const std::vector<std::vector<int>>& B,
                                                int original_n) {
    int new_n = getNextPowerOf2(original_n);

    // Step 1: Create padded matrices Ap and Bp.
    std::vector<std::vector<int>> Ap(new_n, std::vector<int>(new_n, 0));
    std::vector<std::vector<int>> Bp(new_n, std::vector<int>(new_n, 0));

    // Step 2: Copy original matrices into padded matrices.
    for (int i = 0; i < original_n; ++i) {
        for (int j = 0; j < original_n; ++j) {
            Ap[i][j] = A[i][j];
            Bp[i][j] = B[i][j];
        }
    }

    // Step 3: Perform Strassen's multiplication on padded matrices.
    std::vector<std::vector<int>> C_padded = strassen_multiply(Ap, Bp, new_n);

    // Step 4: Extract the result of the original dimensions.
    std::vector<std::vector<int>> C_original(original_n, std::vector<int>(original_n));
    for (int i = 0; i < original_n; ++i) {
        for (int j = 0; j < original_n; ++j) {
            C_original[i][j] = C_padded[i][j];
        }
    }

    // Step 5: Return the result for the original matrix size.
    return C_original;
}

int main() {
    // Step 1: Define matrices A and B (2x2 example).
    std::vector<std::vector<int>> A = {{1, 2}, {3, 4}};
    std::vector<std::vector<int>> B = {{5, 6}, {7, 8}};
    int n = A.size();

    if (n == 0 || A[0].size() != n || B.size() != n || B[0].size() != n) {
        std::cout << "Invalid matrix dimensions for square matrices." << std::endl;
        return 1;
    }

    // Step 2: Print original matrices.
    std::cout << "Matrix A:" << std::endl;
    printMatrix(A, n);
    std::cout << "\nMatrix B:" << std::endl;
    printMatrix(B, n);

    // Step 3: Perform Strassen's multiplication using the wrapper to handle padding.
    std::vector<std::vector<int>> C = strassen_wrapper(A, B, n);

    // Step 4: Print the result matrix.
    std::cout << "\nResult of Strassen's Multiplication (C = A * B):" << std::endl;
    printMatrix(C, n);

    // Example with 3x3 matrix (will be padded to 4x4)
    std::cout << "\n--- 3x3 Matrix Example (padded to 4x4 for Strassen) ---" << std::endl;
    std::vector<std::vector<int>> A_3x3 = {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}
    };
    std::vector<std::vector<int>> B_3x3 = {
        {9, 8, 7},
        {6, 5, 4},
        {3, 2, 1}
    };
    int n_3x3 = A_3x3.size();

    std::cout << "Matrix A (3x3):" << std::endl;
    printMatrix(A_3x3, n_3x3);
    std::cout << "\nMatrix B (3x3):" << std::endl;
    printMatrix(B_3x3, n_3x3);

    std::vector<std::vector<int>> C_3x3_strassen = strassen_wrapper(A_3x3, B_3x3, n_3x3);
    std::cout << "\nResult of Strassen's Multiplication (3x3):" << std::endl;
    printMatrix(C_3x3_strassen, n_3x3);

    // Verify with naive for 3x3
    std::vector<std::vector<int>> C_3x3_naive = multiplyNaive(A_3x3, B_3x3); // Needs multiplyNaive to handle non-power of 2, which it does.
    std::cout << "\nResult of Naive Multiplication (3x3 for comparison):" << std::endl;
    printMatrix(C_3x3_naive, n_3x3);


    return 0;
}