Example: Strassen's Matrix Multiplication in C++

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

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

// Function to add two matrices
std::vector<std::vector<int>> addMatrices(
    const std::vector<std::vector<int>>& A,
    const std::vector<std::vector<int>>& B) {
    int n = A.size();
    std::vector<std::vector<int>> C(n, std::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
std::vector<std::vector<int>> subtractMatrices(
    const std::vector<std::vector<int>>& A,
    const std::vector<std::vector<int>>& B) {
    int n = A.size();
    std::vector<std::vector<int>> C(n, std::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 for standard matrix multiplication (base case for Strassen)
std::vector<std::vector<int>> baseCaseMultiply(
    const std::vector<std::vector<int>>& A,
    const std::vector<std::vector<int>>& B) {
    int n = A.size();
    std::vector<std::vector<int>> C(n, std::vector<int>(n));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            for (int k = 0; k < n; ++k) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}


// Strassen's Matrix Multiplication
std::vector<std::vector<int>> strassenMultiply(
    const std::vector<std::vector<int>>& A,
    const std::vector<std::vector<int>>& B) {

    int n = A.size();

    // Base case: If matrix size is small, use standard multiplication
    if (n <= 16) { // Threshold can be tuned for performance
        return baseCaseMultiply(A, B);
    }

    // Step 1: Divide matrices into 4 sub-matrices
    int newSize = n / 2;
    std::vector<std::vector<int>> A11(newSize, std::vector<int>(newSize));
    std::vector<std::vector<int>> A12(newSize, std::vector<int>(newSize));
    std::vector<std::vector<int>> A21(newSize, std::vector<int>(newSize));
    std::vector<std::vector<int>> A22(newSize, std::vector<int>(newSize));

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

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

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

    // Step 2: Compute the 7 products recursively
    std::vector<std::vector<int>> P1 = strassenMultiply(A11, subtractMatrices(B12, B22));
    std::vector<std::vector<int>> P2 = strassenMultiply(addMatrices(A11, A12), B22);
    std::vector<std::vector<int>> P3 = strassenMultiply(addMatrices(A21, A22), B11);
    std::vector<std::vector<int>> P4 = strassenMultiply(A22, subtractMatrices(B21, B11));
    std::vector<std::vector<int>> P5 = strassenMultiply(addMatrices(A11, A22), addMatrices(B11, B22));
    std::vector<std::vector<int>> P6 = strassenMultiply(subtractMatrices(A12, A22), addMatrices(B21, B22));
    std::vector<std::vector<int>> P7 = strassenMultiply(subtractMatrices(A11, A21), addMatrices(B11, B12));

    // Step 3: Compute the 4 resulting sub-matrices
    std::vector<std::vector<int>> C11 = addMatrices(subtractMatrices(addMatrices(P5, P4), P2), P6);
    std::vector<std::vector<int>> C12 = addMatrices(P1, P2);
    std::vector<std::vector<int>> C21 = addMatrices(P3, P4);
    std::vector<std::vector<int>> C22 = subtractMatrices(subtractMatrices(addMatrices(P5, P1), P3), P7);

    // Step 4: Combine sub-matrices into the final result
    std::vector<std::vector<int>> C(n, std::vector<int>(n));
    for (int i = 0; i < newSize; ++i) {
        for (int j = 0; j < newSize; ++j) {
            C[i][j] = C11[i][j];
            C[i][j + newSize] = C12[i][j];
            C[i + newSize][j] = C21[i][j];
            C[i + newSize][j + newSize] = C22[i][j];
        }
    }
    return C;
}


int main() {
    // Example matrices must have dimensions that are powers of 2 for this simplified implementation
    // For general matrices, padding with zeros to the next power of 2 is required.
    std::vector<std::vector<int>> A = {{1, 2, 1, 2},
                                        {3, 4, 3, 4},
                                        {1, 2, 1, 2},
                                        {3, 4, 3, 4}};

    std::vector<std::vector<int>> B = {{5, 6, 5, 6},
                                        {7, 8, 7, 8},
                                        {5, 6, 5, 6},
                                        {7, 8, 7, 8}};

    std::cout << "Matrix A (4x4):" << std::endl;
    printMatrix(A);
    std::cout << "\nMatrix B (4x4):" << std::endl;
    printMatrix(B);

    std::vector<std::vector<int>> C = strassenMultiply(A, B);

    std::cout << "\nResult of Strassen's Multiplication (C = A * B):" << std::endl;
    printMatrix(C);

    // Another example (2x2)
    std::vector<std::vector<int>> A_small = {{1, 2}, {3, 4}};
    std::vector<std::vector<int>> B_small = {{5, 6}, {7, 8}};

    std::cout << "\nMatrix A (2x2):" << std::endl;
    printMatrix(A_small);
    std::cout << "\nMatrix B (2x2):" << std::endl;
    printMatrix(B_small);

    std::vector<std::vector<int>> C_small = strassenMultiply(A_small, B_small);

    std::cout << "\nResult of Strassen's Multiplication (C = A * B) for 2x2:" << std::endl;
    printMatrix(C_small);

    return 0;
}