Rust, score≈7.9

Rust port of @Ander Kaseorg's C++ code.

Cargo.toml

[package]
name = "fast_matrix_multiplicator_evaluator"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
rand_pcg = "0.3.1"
rand = "0.8.0"


[profile.release]
lto = true
panic = "abort"

src/main.rs

use std::env;
use std::time::{Duration, Instant};

// Import necessary traits and structs from the rand crate
use rand::Rng;
use rand::SeedableRng;
use rand_pcg::Pcg64;

// Function to initialize a vector with pseudo-random values using a fixed seed
fn initialize_vector(size: usize) -> Vec<f64> {
    // Initialize the RNG with a fixed seed for reproducibility
    let mut rng = Pcg64::seed_from_u64(0);
    let mut vec = Vec::with_capacity(size);
    for _ in 0..size {
        // Generate a random f64 between 0.0 and 1.0
        vec.push(rng.gen::<f64>());
    }
    vec
}

// The `error` function translated from C
fn error_func(gamma: usize, sigma: usize, A: &[f64], B: &[f64], C: &[f64]) -> f64 {
    let mut E = 0.0;
    for r in 0..sigma {
        for i in 0..gamma {
            for j in 0..gamma {
                let mut x = 0.0;
                for l in 0..gamma {
                    let rkl = r * gamma * gamma + j * gamma + l;
                    let rnm = r * gamma * gamma + i * gamma + l;
                    x += B[rkl] * C[rnm];
                }
                let rij = r * gamma * gamma + i * gamma + j;
                E += A[rij] * x;
            }
        }
    }
    E = (gamma as f64).powi(3) - 2.0 * E;
    for r in 0..sigma {
        for s in r..sigma {
            let mut xa = 0.0;
            let mut xb = 0.0;
            let mut xc = 0.0;
            for ij in 0..(gamma * gamma) {
                let rij = r * gamma * gamma + ij;
                let sij = s * gamma * gamma + ij;
                xa += A[rij] * A[sij];
                xb += B[rij] * B[sij];
                xc += C[rij] * C[sij];
            }
            E += (1.0 + if r != s { 1.0 } else { 0.0 }) * xa * xb * xc;
        }
    }
    E
}

// The `reference` function translated from C
fn reference_func(gamma: usize, sigma: usize, A: &[f64], B: &[f64], C: &[f64]) -> f64 {
    let mut E = 0.0;
    for i in 0..gamma {
        for j in 0..gamma {
            for k in 0..gamma {
                for l in 0..gamma {
                    for m in 0..gamma {
                        for n in 0..gamma {
                            let mut s = 0.0;
                            for r in 0..sigma {
                                let rij = r * gamma * gamma + i * gamma + j;
                                let rkl = r * gamma * gamma + k * gamma + l;
                                let rnm = r * gamma * gamma + n * gamma + m;
                                s += A[rij] * B[rkl] * C[rnm];
                            }
                            let e = s - if n == i && j == k && l == m {
                                1.0
                            } else {
                                0.0
                            };
                            E += e * e;
                        }
                    }
                }
            }
        }
    }
    E
}

fn main() {
    // Parse command-line arguments
    let args: Vec<String> = env::args().collect();
    if args.len() != 3 {
        eprintln!("Usage: {} <gamma> <sigma>", args[0]);
        std::process::exit(1);
    }
    let gamma: usize = args[1].parse().expect("Invalid gamma");
    let sigma: usize = args[2].parse().expect("Invalid sigma");

    // Allocate and initialize A, B, C
    let size = sigma * gamma * gamma;
    let A = initialize_vector(size);
    let B = initialize_vector(size);
    let C = initialize_vector(size);

    // Measure execution time for `reference_func`
    let reference_start = Instant::now();
    let mut reference_count = 0;
    let mut reference_E = 0.0;
    while reference_start.elapsed() < Duration::from_secs(1) {
        reference_E = reference_func(gamma, sigma, &A, &B, &C);
        reference_count += 1;
    }
    let reference_time = reference_start.elapsed().as_secs_f64() / reference_count as f64;
    println!("reference: E={} time={} s", reference_E, reference_time);

    // Measure execution time for `error_func`
    let error_start = Instant::now();
    let mut error_count = 0;
    let mut error_E = 0.0;
    while error_start.elapsed() < Duration::from_secs(1) {
        error_E = error_func(gamma, sigma, &A, &B, &C);
        error_count += 1;
    }
    let error_time = error_start.elapsed().as_secs_f64() / error_count as f64;
    println!("me:        E={} time={} s", error_E, error_time);

    // Calculate and print the score
    println!("score={}", reference_time / error_time);
}

Run

$ cargo.exe run --package fast_matrix_multiplicator_evaluator --bin fast_matrix_multiplicator_evaluator --release 
$ .\target\release\fast_matrix_multiplicator_evaluator.exe 6 159
reference: E=18104555.908983983 time=0.004142786363636363 s
me:        E=18104555.90898394 time=0.0005224228720626632 s
score=7.929948295103448

C (gcc), score ≈ 6.7×, 8.5×, 13.6x, 19.5×

For the four categories (3, 22), (4, 48), (5, 97), (6, 159).

This algorithm runs in \$O(γ^3σ + γ^2σ^2)\$ rather than \$O(γ^6σ)\$, using the rearranged formula:

\begin{multline*} \textrm{error}(A, B, C) = γ^3 - 2\sum_{r=1}^σ \sum_{i,j,l=1}^γ A_{ij}^rB_{jl}^rC_{il}^r \\ + \sum_{r,s=1}^σ\left(\sum_{i,j=1}^γ A_{ij}^rA_{ij}^s\right)\left(\sum_{k,l=1}^γ B_{kl}^rB_{kl}^s\right)\left(\sum_{n,m=1}^γ C_{nm}^rC_{nm}^s\right). \end{multline*}

double error(size_t gamma, size_t sigma, double *A, double *B, double *C) {
  double E = 0.0;
  for (int r = 0; r < sigma; ++r) {
    for (int i = 0; i < gamma; ++i) {
      for (int j = 0; j < gamma; ++j) {
        double x = 0.0;
        for (int l = 0; l < gamma; ++l) {
          int rkl = r * gamma * gamma + j * gamma + l;
          int rnm = r * gamma * gamma + i * gamma + l;
          x += B[rkl] * C[rnm];
        }
        int rij = r * gamma * gamma + i * gamma + j;
        E += A[rij] * x;
      }
    }
  }
  E = gamma * gamma * gamma - 2 * E;
  for (int r = 0; r < sigma; ++r) {
    for (int s = r; s < sigma; ++s) {
      double xa = 0.0, xb = 0.0, xc = 0.0;
      for (int ij = 0; ij < gamma * gamma; ++ij) {
        int rij = r * gamma * gamma + ij;
        int sij = s * gamma * gamma + ij;
        xa += A[rij] * A[sij];
        xb += B[rij] * B[sij];
        xc += C[rij] * C[sij];
      }
      E += (1 + (r != s)) * xa * xb * xc;
    }
  }
  return E;
}

Try it online!