I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Strassen's algorithm

From Rosetta Code
Strassen's algorithm is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Description

In linear algebra, the Strassen algorithm, named after Volker Strassen, is an algorithm for matrix multiplication. It is faster than the standard matrix multiplication algorithm and is useful in practice for large matrices, but would be slower than the fastest known algorithms for extremely large matrices.

Task

Write a routine, function, procedure etc. in your language to implement the Strassen algorithm for matrix multiplication.

While practical implementations of Strassen's algorithm usually switch to standard methods of matrix multiplication for small enough submatrices (currently anything < 512x512 according to wp), for the purposes of this task you should not switch until reaching a size of 1 or 2.

Related task
See also


Go[edit]

Translation of: Wren

Rather than use a library such as gonum, we create a simple Matrix type which is adequate for this task.

package main
 
import (
"fmt"
"log"
"math"
)
 
type Matrix [][]float64
 
func (m Matrix) rows() int { return len(m) }
func (m Matrix) cols() int { return len(m[0]) }
 
func (m Matrix) add(m2 Matrix) Matrix {
if m.rows() != m2.rows() || m.cols() != m2.cols() {
log.Fatal("Matrices must have the same dimensions.")
}
c := make(Matrix, m.rows())
for i := 0; i < m.rows(); i++ {
c[i] = make([]float64, m.cols())
for j := 0; j < m.cols(); j++ {
c[i][j] = m[i][j] + m2[i][j]
}
}
return c
}
 
func (m Matrix) sub(m2 Matrix) Matrix {
if m.rows() != m2.rows() || m.cols() != m2.cols() {
log.Fatal("Matrices must have the same dimensions.")
}
c := make(Matrix, m.rows())
for i := 0; i < m.rows(); i++ {
c[i] = make([]float64, m.cols())
for j := 0; j < m.cols(); j++ {
c[i][j] = m[i][j] - m2[i][j]
}
}
return c
}
 
func (m Matrix) mul(m2 Matrix) Matrix {
if m.cols() != m2.rows() {
log.Fatal("Cannot multiply these matrices.")
}
c := make(Matrix, m.rows())
for i := 0; i < m.rows(); i++ {
c[i] = make([]float64, m2.cols())
for j := 0; j < m2.cols(); j++ {
for k := 0; k < m2.rows(); k++ {
c[i][j] += m[i][k] * m2[k][j]
}
}
}
return c
}
 
func (m Matrix) toString(p int) string {
s := make([]string, m.rows())
pow := math.Pow(10, float64(p))
for i := 0; i < m.rows(); i++ {
t := make([]string, m.cols())
for j := 0; j < m.cols(); j++ {
r := math.Round(m[i][j]*pow) / pow
t[j] = fmt.Sprintf("%g", r)
if t[j] == "-0" {
t[j] = "0"
}
}
s[i] = fmt.Sprintf("%v", t)
}
return fmt.Sprintf("%v", s)
}
 
func params(r, c int) [4][6]int {
return [4][6]int{
{0, r, 0, c, 0, 0},
{0, r, c, 2 * c, 0, c},
{r, 2 * r, 0, c, r, 0},
{r, 2 * r, c, 2 * c, r, c},
}
}
 
func toQuarters(m Matrix) [4]Matrix {
r := m.rows() / 2
c := m.cols() / 2
p := params(r, c)
var quarters [4]Matrix
for k := 0; k < 4; k++ {
q := make(Matrix, r)
for i := p[k][0]; i < p[k][1]; i++ {
q[i-p[k][4]] = make([]float64, c)
for j := p[k][2]; j < p[k][3]; j++ {
q[i-p[k][4]][j-p[k][5]] = m[i][j]
}
}
quarters[k] = q
}
return quarters
}
 
func fromQuarters(q [4]Matrix) Matrix {
r := q[0].rows()
c := q[0].cols()
p := params(r, c)
r *= 2
c *= 2
m := make(Matrix, r)
for i := 0; i < c; i++ {
m[i] = make([]float64, c)
}
for k := 0; k < 4; k++ {
for i := p[k][0]; i < p[k][1]; i++ {
for j := p[k][2]; j < p[k][3]; j++ {
m[i][j] = q[k][i-p[k][4]][j-p[k][5]]
}
}
}
return m
}
 
func strassen(a, b Matrix) Matrix {
if a.rows() != a.cols() || b.rows() != b.cols() || a.rows() != b.rows() {
log.Fatal("Matrices must be square and of equal size.")
}
if a.rows() == 0 || (a.rows()&(a.rows()-1)) != 0 {
log.Fatal("Size of matrices must be a power of two.")
}
if a.rows() == 1 {
return a.mul(b)
}
qa := toQuarters(a)
qb := toQuarters(b)
p1 := strassen(qa[1].sub(qa[3]), qb[2].add(qb[3]))
p2 := strassen(qa[0].add(qa[3]), qb[0].add(qb[3]))
p3 := strassen(qa[0].sub(qa[2]), qb[0].add(qb[1]))
p4 := strassen(qa[0].add(qa[1]), qb[3])
p5 := strassen(qa[0], qb[1].sub(qb[3]))
p6 := strassen(qa[3], qb[2].sub(qb[0]))
p7 := strassen(qa[2].add(qa[3]), qb[0])
var q [4]Matrix
q[0] = p1.add(p2).sub(p4).add(p6)
q[1] = p4.add(p5)
q[2] = p6.add(p7)
q[3] = p2.sub(p3).add(p5).sub(p7)
return fromQuarters(q)
}
 
func main() {
a := Matrix{{1, 2}, {3, 4}}
b := Matrix{{5, 6}, {7, 8}}
c := Matrix{{1, 1, 1, 1}, {2, 4, 8, 16}, {3, 9, 27, 81}, {4, 16, 64, 256}}
d := Matrix{{4, -3, 4.0 / 3, -1.0 / 4}, {-13.0 / 3, 19.0 / 4, -7.0 / 3, 11.0 / 24},
{3.0 / 2, -2, 7.0 / 6, -1.0 / 4}, {-1.0 / 6, 1.0 / 4, -1.0 / 6, 1.0 / 24}}
e := Matrix{{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 16}}
f := Matrix{{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}
fmt.Println("Using 'normal' matrix multiplication:")
fmt.Printf(" a * b = %v\n", a.mul(b))
fmt.Printf(" c * d = %v\n", c.mul(d).toString(6))
fmt.Printf(" e * f = %v\n", e.mul(f))
fmt.Println("\nUsing 'Strassen' matrix multiplication:")
fmt.Printf(" a * b = %v\n", strassen(a, b))
fmt.Printf(" c * d = %v\n", strassen(c, d).toString(6))
fmt.Printf(" e * f = %v\n", strassen(e, f))
}
Output:
Using 'normal' matrix multiplication:
  a * b = [[19 22] [43 50]]
  c * d = [[1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1]]
  e * f = [[1 2 3 4] [5 6 7 8] [9 10 11 12] [13 14 15 16]]

Using 'Strassen' matrix multiplication:
  a * b = [[19 22] [43 50]]
  c * d = [[1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1]]
  e * f = [[1 2 3 4] [5 6 7 8] [9 10 11 12] [13 14 15 16]]

Julia[edit]

With dynamic padding[edit]

Because Julia uses column major in matrices, sometimes the code uses the adjoint of a matrix in order to match examples as written.

"""
Strassen's matrix multiplication algorithm.
Use dynamic padding in order to reduce required auxiliary memory.
"""
function strassen(x::Matrix, y::Matrix)
# Check that the sizes of these matrices match.
(r1, c1) = size(x)
(r2, c2) = size(y)
if c1 != r2
error("Multiplying $r1 x $c1 and $r2 x $c2 matrix: dimensions do not match.")
end
 
# Put a matrix into the top left of a matrix of zeros.
# `rows` and `cols` are the dimensions of the output matrix.
function embed(mat, rows, cols)
# If the matrix is already of the right dimensions, don't allocate new memory.
(r, c) = size(mat)
if (r, c) == (rows, cols)
return mat
end
 
# Pad the matrix with zeros to be the right size.
out = zeros(Int, rows, cols)
out[1:r, 1:c] = mat
out
end
 
# Make sure both matrices are the same size.
# This is exclusively for simplicity:
# this algorithm can be implemented with matrices of different sizes.
r = max(r1, r2); c = max(c1, c2)
x = embed(x, r, c)
y = embed(y, r, c)
 
# Our recursive multiplication function.
function block_mult(a, b, rows, cols)
# For small matrices, resort to naive multiplication.
# if rows <= 128 || cols <= 128
if rows == 1 && cols == 1
# if rows == 2 && cols == 2
return a * b
end
 
# Apply dynamic padding.
if rows % 2 == 1 && cols % 2 == 1
a = embed(a, rows + 1, cols + 1)
b = embed(b, rows + 1, cols + 1)
elseif rows % 2 == 1
a = embed(a, rows + 1, cols)
b = embed(b, rows + 1, cols)
elseif cols % 2 == 1
a = embed(a, rows, cols + 1)
b = embed(b, rows, cols + 1)
end
 
half_rows = Int(size(a, 1) / 2)
half_cols = Int(size(a, 2) / 2)
 
# Subdivide input matrices.
a11 = a[1:half_rows, 1:half_cols]
b11 = b[1:half_rows, 1:half_cols]
 
a12 = a[1:half_rows, half_cols+1:size(a, 2)]
b12 = b[1:half_rows, half_cols+1:size(a, 2)]
 
a21 = a[half_rows+1:size(a, 1), 1:half_cols]
b21 = b[half_rows+1:size(a, 1), 1:half_cols]
 
a22 = a[half_rows+1:size(a, 1), half_cols+1:size(a, 2)]
b22 = b[half_rows+1:size(a, 1), half_cols+1:size(a, 2)]
 
# Compute intermediate values.
multip(x, y) = block_mult(x, y, half_rows, half_cols)
m1 = multip(a11 + a22, b11 + b22)
m2 = multip(a21 + a22, b11)
m3 = multip(a11, b12 - b22)
m4 = multip(a22, b21 - b11)
m5 = multip(a11 + a12, b22)
m6 = multip(a21 - a11, b11 + b12)
m7 = multip(a12 - a22, b21 + b22)
 
# Combine intermediate values into the output.
c11 = m1 + m4 - m5 + m7
c12 = m3 + m5
c21 = m2 + m4
c22 = m1 - m2 + m3 + m6
 
# Crop output to the desired size (undo dynamic padding).
out = [c11 c12; c21 c22]
out[1:rows, 1:cols]
end
 
block_mult(x, y, r, c)
end
 
const A = [[1, 2] [3, 4]]
const B = [[5, 6] [7, 8]]
const C = [[1, 1, 1, 1] [2, 4, 8, 16] [3, 9, 27, 81] [4, 16, 64, 256]]
const D = [[4, -3, 4/3, -1/4] [-13/3, 19/4, -7/3, 11/24] [3/2, -2, 7/6, -1/4] [-1/6, 1/4, -1/6, 1/24]]
const E = [[1, 2, 3, 4] [5, 6, 7, 8] [9, 10, 11, 12] [13, 14, 15, 16]]
const F = [[1, 0, 0, 0] [0, 1, 0, 0] [0, 0, 1, 0] [0, 0, 0, 1]]
 
""" For pretty printing: change matrix to integer if it is within 0.00000001 of an integer """
intprint(s, mat) = println(s, map(x -> Int(round(x, digits=8)), mat)')
 
intprint("Regular multiply: ", A' * B')
intprint("Strassen multiply: ", strassen(Matrix(A'), Matrix(B')))
intprint("Regular multiply: ", C * D)
intprint("Strassen multiply: ", strassen(C, D))
intprint("Regular multiply: ", E * F)
intprint("Strassen multiply: ", strassen(E, F))
 
const r = sqrt(2)/2
const R = [[r, r] [-r, r]]
 
intprint("Regular multiply: ", R * R)
intprint("Strassen multiply: ", strassen(R,R))
 
Output:
Regular multiply: [19 43; 22 50]
Strassen multiply: [19 43; 22 50]
Regular multiply: [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
Strassen multiply: [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
Regular multiply: [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
Strassen multiply: [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
Regular multiply: [0 1; -1 0]
Strassen multiply: [0 1; -1 0]

Recursive[edit]

Output is the same as the dynamically padded version.

function Strassen(A, B)
n = size(A, 1)
if n == 1
return A * B
end
@views A11 = A[1:n÷2, 1:n÷2]
@views A12 = A[1:n÷2, n÷2+1:n]
@views A21 = A[n÷2+1:n, 1:n÷2]
@views A22 = A[n÷2+1:n, n÷2+1:n]
@views B11 = B[1:n÷2, 1:n÷2]
@views B12 = B[1:n÷2, n÷2+1:n]
@views B21 = B[n÷2+1:n, 1:n÷2]
@views B22 = B[n÷2+1:n, n÷2+1:n]
 
P1 = Strassen(A12 - A22, B21 + B22)
P2 = Strassen(A11 + A22, B11 + B22)
P3 = Strassen(A11 - A21, B11 + B12)
P4 = Strassen(A11 + A12, B22)
P5 = Strassen(A11, B12 - B22)
P6 = Strassen(A22, B21 - B11)
P7 = Strassen(A21 + A22, B11)
 
C11 = P1 + P2 - P4 + P6
C12 = P4 + P5
C21 = P6 + P7
C22 = P2 - P3 + P5 - P7
 
return [C11 C12; C21 C22]
end
 
const A = [[1, 2] [3, 4]]
const B = [[5, 6] [7, 8]]
const C = [[1, 1, 1, 1] [2, 4, 8, 16] [3, 9, 27, 81] [4, 16, 64, 256]]
const D = [[4, -3, 4/3, -1/4] [-13/3, 19/4, -7/3, 11/24] [3/2, -2, 7/6, -1/4] [-1/6, 1/4, -1/6, 1/24]]
const E = [[1, 2, 3, 4] [5, 6, 7, 8] [9, 10, 11, 12] [13, 14, 15, 16]]
const F = [[1, 0, 0, 0] [0, 1, 0, 0] [0, 0, 1, 0] [0, 0, 0, 1]]
 
intprint(s, mat) = println(s, map(x -> Int(round(x, digits=8)), mat)')
intprint("Regular multiply: ", A' * B')
intprint("Strassen multiply: ", Strassen(Matrix(A'), Matrix(B')))
intprint("Regular multiply: ", C * D)
intprint("Strassen multiply: ", Strassen(C, D))
intprint("Regular multiply: ", E * F)
intprint("Strassen multiply: ", Strassen(E, F))
 
const r = sqrt(2)/2
const R = [[r, r] [-r, r]]
 
intprint("Regular multiply: ", R * R)
intprint("Strassen multiply: ", Strassen(R,R))
 

Phix[edit]

As noted on wp, you could pad with zeroes, and strip them on exit, instead of crashing for non-square 2n matrices.

function strassen(sequence a, b)
integer l = length(a)
if length(a[1])!=l
or length(b)!=l
or length(b[1])!=l then
crash("two equal square matrices only")
end if
if l=1 then return sq_mul(a,b) end if
if remainder(l,1) then
crash("2^n matrices only")
end if
integer h = l/2
sequence {a11,a12,a21,a22,b11,b12,b21,b22} @= repeat(repeat(0,h),h)
for i=1 to h do
for j=1 to h do
a11[i][j] = a[i][j]
a12[i][j] = a[i][j+h]
a21[i][j] = a[i+h][j]
a22[i][j] = a[i+h][j+h]
b11[i][j] = b[i][j]
b12[i][j] = b[i][j+h]
b21[i][j] = b[i+h][j]
b22[i][j] = b[i+h][j+h]
end for
end for
sequence p1 = strassen(sq_sub(a12,a22), sq_add(b21,b22)),
p2 = strassen(sq_add(a11,a22), sq_add(b11,b22)),
p3 = strassen(sq_sub(a11,a21), sq_add(b11,b12)),
p4 = strassen(sq_add(a11,a12), b22),
p5 = strassen(a11, sq_sub(b12,b22)),
p6 = strassen(a22, sq_sub(b21,b11)),
p7 = strassen(sq_add(a21,a22), b11),
 
c11 = sq_add(sq_sub(sq_add(p1,p2),p4),p6),
c12 = sq_add(p4,p5),
c21 = sq_add(p6,p7),
c22 = sq_sub(sq_add(sq_sub(p2,p3),p5),p7),
c = repeat(repeat(0,l),l)
for i=1 to h do
for j=1 to h do
c[i][j] = c11[i][j]
c[i][j+h] = c12[i][j]
c[i+h][j] = c21[i][j]
c[i+h][j+h] = c22[i][j]
end for
end for
return c
end function
 
ppOpt({pp_Nest,1,pp_IntFmt,"%3d",pp_FltFmt,"%3.0f",pp_IntCh,false})
 
constant A = {{1,2},
{3,4}},
B = {{5,6},
{7,8}}
pp(strassen(A,B))
 
constant C = { { 1, 1, 1, 1 },
{ 2, 4, 8, 16 },
{ 3, 9, 27, 81 },
{ 4, 16, 64, 256 }},
D = { { 4, -3, 4/3, -1/ 4 },
{-13/3, 19/4, -7/3, 11/24 },
{ 3/2, -2, 7/6, -1/ 4 },
{ -1/6, 1/4, -1/6, 1/24 }}
pp(strassen(C,D))
 
constant E = {{ 1, 2, 3, 4},
{ 5, 6, 7, 8},
{ 9,10,11,12},
{13,14,15,16}},
F = {{1, 0, 0, 0},
{0, 1, 0, 0},
{0, 0, 1, 0},
{0, 0, 0, 1}}
pp(strassen(E,F))
 
constant r = sqrt(2)/2,
R = {{ r,r},
{-r,r}}
pp(strassen(R,R))
Output:

Matches that of Matrix_multiplication#Phix, when given the same inputs

{{ 19, 22},
 { 43, 50}}
{{  1,  0,  0,  0},
 {  0,  1,  0,  0},
 {  0,  0,  1,  0},
 {  0,  0,  0,  1}}
{{  1,  2,  3,  4},
 {  5,  6,  7,  8},
 {  9, 10, 11, 12},
 { 13, 14, 15, 16}}
{{  0,  1},
 { -1,  0}}

Wren[edit]

Library: Wren-fmt

Wren doesn't currently have a matrix module so I've written a rudimentary Matrix class with sufficient functionality to complete this task.

I've used the Phix entry's examples to test the Strassen algorithm implementation.

import "/fmt" for Fmt
 
class Matrix {
construct new(a) {
if (a.type != List || a.count == 0 || a[0].type != List || a[0].count == 0 || a[0][0].type != Num) {
Fiber.abort("Argument must be a non-empty two dimensional list of numbers.")
}
_a = a
}
 
rows { _a.count }
cols { _a[0].count }
 
+(b) {
if (b.type != Matrix) Fiber.abort("Argument must be a matrix.")
if ((this.rows != b.rows) || (this.cols != b.cols)) {
Fiber.abort("Matrices must have the same dimensions.")
}
var c = List.filled(rows, null)
for (i in 0...rows) {
c[i] = List.filled(cols, 0)
for (j in 0...cols) c[i][j] = _a[i][j] + b[i, j]
}
return Matrix.new(c)
}
 
- { this * -1 }
 
-(b) { this + (-b) }
 
*(b) {
var c = List.filled(rows, null)
if (b is Num) {
for (i in 0...rows) {
c[i] = List.filled(cols, 0)
for (j in 0...cols) c[i][j] = _a[i][j] * b
}
} else if (b is Matrix) {
if (this.cols != b.rows) Fiber.abort("Cannot multiply these matrices.")
for (i in 0...rows) {
c[i] = List.filled(b.cols, 0)
for (j in 0...b.cols) {
for (k in 0...b.rows) c[i][j] = c[i][j] + _a[i][k] * b[k, j]
}
}
} else {
Fiber.abort("Argument must be a matrix or a number.")
}
return Matrix.new(c)
}
 
[i] { _a[i].toList }
 
[i, j] { _a[i][j] }
 
toString { _a.toString }
 
// rounds all elements to 'p' places
toString(p) {
var s = List.filled(rows, "")
var pow = 10.pow(p)
for (i in 0...rows) {
var t = List.filled(cols, "")
for (j in 0...cols) {
var r = (_a[i][j]*pow).round / pow
t[j] = r.toString
if (t[j] == "-0") t[j] = "0"
}
s[i] = t.toString
}
return s
}
}
 
var params = Fn.new { |r, c|
return [
[0...r, 0...c, 0, 0],
[0...r, c...2*c, 0, c],
[r...2*r, 0...c, r, 0],
[r...2*r, c...2*c, r, c]
]
}
 
var toQuarters = Fn.new { |m|
var r = (m.rows/2).floor
var c = (m.cols/2).floor
var p = params.call(r, c)
var quarters = []
for (k in 0..3) {
var q = List.filled(r, null)
for (i in p[k][0]) {
q[i - p[k][2]] = List.filled(c, 0)
for (j in p[k][1]) q[i - p[k][2]][j - p[k][3]] = m[i, j]
}
quarters.add(Matrix.new(q))
}
return quarters
}
 
var fromQuarters = Fn.new { |q|
var r = q[0].rows
var c = q[0].cols
var p = params.call(r, c)
r = r * 2
c = c * 2
var m = List.filled(r, null)
for (i in 0...c) m[i] = List.filled(c, 0)
for (k in 0..3) {
for (i in p[k][0]) {
for (j in p[k][1]) m[i][j] = q[k][i - p[k][2], j - p[k][3]]
}
}
return Matrix.new(m)
}
 
var strassen // recursive
strassen = Fn.new { |a, b|
if (a.rows != a.cols || b.rows != b.cols || a.rows != b.rows) {
Fiber.abort("Matrices must be square and of equal size.")
}
if (a.rows == 0 || (a.rows & (a.rows - 1)) != 0) {
Fiber.abort("Size of matrices must be a power of two.")
}
if (a.rows == 1) return a * b
var qa = toQuarters.call(a)
var qb = toQuarters.call(b)
var p1 = strassen.call(qa[1] - qa[3], qb[2] + qb[3])
var p2 = strassen.call(qa[0] + qa[3], qb[0] + qb[3])
var p3 = strassen.call(qa[0] - qa[2], qb[0] + qb[1])
var p4 = strassen.call(qa[0] + qa[1], qb[3])
var p5 = strassen.call(qa[0], qb[1] - qb[3])
var p6 = strassen.call(qa[3], qb[2] - qb[0])
var p7 = strassen.call(qa[2] + qa[3], qb[0])
var q = List.filled(4, null)
q[0] = p1 + p2 - p4 + p6
q[1] = p4 + p5
q[2] = p6 + p7
q[3] = p2 - p3 + p5 - p7
return fromQuarters.call(q)
}
 
var a = Matrix.new([ [1,2], [3, 4] ])
var b = Matrix.new([ [5,6], [7, 8] ])
var c = Matrix.new([ [1, 1, 1, 1], [2, 4, 8, 16], [3, 9, 27, 81], [4, 16, 64, 256] ])
var d = Matrix.new([ [4, -3, 4/3, -1/4], [-13/3, 19/4, -7/3, 11/24],
[3/2, -2, 7/6, -1/4], [-1/6, 1/4, -1/6, 1/24] ])
var e = Matrix.new([ [1, 2, 3, 4], [5, 6, 7, 8], [9,10,11,12], [13,14,15,16] ])
var f = Matrix.new([ [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1] ])
System.print("Using 'normal' matrix multiplication:")
System.print(" a * b = %(a * b)")
System.print(" c * d = %((c * d).toString(6))")
System.print(" e * f = %(e * f)")
System.print("\nUsing 'Strassen' matrix multiplication:")
System.print(" a * b = %(strassen.call(a, b))")
System.print(" c * d = %(strassen.call(c, d).toString(6))")
System.print(" e * f = %(strassen.call(e, f))")
Output:
Using 'normal' matrix multiplication:
  a * b = [[19, 22], [43, 50]]
  c * d = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
  e * f = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]

Using 'Strassen' matrix multiplication:
  a * b = [[19, 22], [43, 50]]
  c * d = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
  e * f = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]