Skeptical about that parallel speedup by 23 times?!?!

• L3 Cache (30 MiB total in both CPU's): This is huge. It's likely why your matrix inversion is running so well; the Xeon can keep large chunks of your BlockMatrix close to the silicon instead of going out to the slower RAM.
• AVX2 Support: Your Rust code is likely using "SIMD" (Single Instruction, Multiple Data) instructions automatically. This means that inside each of those 24 threads, the CPU is processing multiple numbers in a single clock cycle.

The Rayon Parallelism Engine

--- ENCRYPTION KEY (E) ---
[[[[296694, 390531, 243675], [320449, 970431, 98244], [627164, 953395, 916306]], 
  [[430087, 647583, 421775], [475138, 345706, 880216], [733525, 368995, 692281]]], 
 [[[637979, 335126, 854517], [71159, 285713, 495290], [214226, 716943, 923249]], 
  [[137641, 656211, 147803], [735157, 506297, 348918], [455121, 649579, 248797]]]]

--- DECRYPTION KEY (D) ---
[[[[813163, 768189, 61393], [826065, 60489, 275868], [261426, 297233, 455745]], 
  [[249617, 528552, 440182], [520685, 959134, 550012], [804146, 877024, 121624]]], 
 [[[361085, 589785, 121899], [113114, 125523, 957137], [150194, 96425, 201282]], 
  [[46970, 404510, 885054], [678588, 155444, 608800], [25970, 885063, 979838]]]]

Technical Whitepaper: The Matrix-with-Matrix Cryptosystem

1. Executive Summary

The Matrix-within-Matrix Cryptosystem is an advanced evolution of classical polyalphabetic ciphers. By treating every element of the encryption process as a sub-matrix rather than a single integer, the system creates a massive computational barrier. It operates within a non-commutative ring of matrices over a finite field, ensuring that every byte of data is subjected to quadratic diffusion.

2. Algebraic Glossary

• Finite Field / Galois Field / GF(P): A closed mathematical system where all results stay within a fixed range (the modulus).
  See Wikipedia for more information about Finite Fields.


• Matrix Ring: In abstract algebra, a matrix ring is defined as a collection of matrices with entries from a ring R that satisfies the properties of a ring under matrix addition and multiplication. The set of all n×n matrices over a ring R is denoted as M_n(R). A structure where elements are square matrices. It is a "Ring" because singular matrices exist which do not have an inverse.
  See Wikipedia for more information about Matrix Rings.


• Non-Commutative: A property where A × B ≠ B × A, preventing simple algebraic division attacks.
  See Wikipedia for more information about Matrix Multiplication and non-commutativity.


• Block Identity Matrix: The target matrix for verification, featuring identity sub-matrices on the diagonal and zero sub-matrices elsewhere. A block identity matrix is a special type of block-diagonal matrix where one or more diagonal blocks are identity matrices and all off-diagonal blocks are zero matrices.
  See Wikipedia for more information about Block Matrices, especially the discussion of "block diagonal matrix".

3. Mathematical Work Multiplier

The "Strength" of the system is measured by the increase in operations compared to a standard Hill cipher. While a baseline cipher uses O(N) operations per element, the Matrix-within-Matrix system scales at O(SUB_DIM²) per element.

In a configuration where SUB_DIM is 20, each byte of plaintext is processed by 400 times more mathematical operations than a standard matrix cipher. At SUB_DIM 48, this multiplier increases to 2,304.

4. Performance Benchmarking

Empirical testing demonstrates that the system scales predictably without hitting a "performance wall," even when key sizes exceed typical L3 cache limits.

Data derived from hardware testing on 2026-02-22 and 2026-02-23.

5. Security & Post-Quantum (PQ) Assessment

The system provides a "Lattice-Like" defense. Unlike RSA or ECC, it does not rely on factoring or discrete logarithms, making it resistant to Shor’s Algorithm. To remain secure against Grover’s Algorithm (quantum brute force), the total keyspace is expanded through the high DIM and SUB_DIM parameters.

6. Resilience to Known-Plaintext Attacks (KPA)

The system mitigates KPA through Dimensionality Explosion. At a 48x48 configuration, an attacker must solve for over 5.3 million unknown variables. The non-commutative nature of the ring ensures that standard linear solver shortcuts are inapplicable.

Deep Dive: The Quadratic Work Multiplier & Matrix Ring Security

1. The Geometry of Diffusion

In a standard Hill cipher, a single byte of plaintext is multiplied by a single scalar in the key matrix. In the Matrix-within-Matrix Matrix Cryptosystem, a single byte of plaintext is conceptually expanded into a coordinate within a sub-matrix.

Because every "element" in our DIM × DIM grid is actually a SUB_DIM × SUB_DIM matrix, the system achieves Horizontal and Vertical Diffusion simultaneously:

• Horizontal Diffusion: The byte interacts with an entire row of the inner sub-matrix.
• Vertical Diffusion: The resulting product is then distributed across the entire outer block structure.

2. Understanding the "Math Multiplier"

The primary benefit of this architecture is the Quadratic Work Multiplier. We calculate the cost of processing a single byte of data as follows:

In our system, for every block-level interaction, we perform a matrix-matrix multiplication. For an inner matrix of size N (where N = SUB_DIM), this requires N³ multiplications and additions. However, when we average this work across the total elements in the block, the cost per byte scales at O(SUB_DIM²).

Note: This multiplier represents the "computational tax" an attacker must pay for every single byte they attempt to analyze or brute-force.

3. The Non-Commutative Advantage

One of the most significant benefits is operating in a Non-Commutative Ring. In standard modular arithmetic, A × B = B × A. This allows attackers to use simple division to isolate keys.

In our matrix ring, the order of operations is immutable (A × B ≠ B × A). This renders standard algebraic "shortcut" attacks useless, as the attacker cannot simply "divide out" the plaintext from the ciphertext without knowing the exact orientation and sequence of the sub-matrix multiplications.

4. Performance Efficiency at Scale

Despite the massive increase in mathematical density, the system remains efficient due to Block Gauss-Jordan Elimination. Our benchmarking shows that even when the complexity increases 1,630x, the hardware handles the operations predictably:

• Predictable Memory Access: The algorithm processes sub-matrices in contiguous memory blocks, which is highly compatible with CPU pre-fetching.
• Cache Resilience: By working on SUB_DIM blocks, the "hot" data often fits within the L1/L2 cache even when the total Key File is as large as 84MB.

5. Security Conclusion

The "Quadratic Nature" of the system means that for every small increase in user effort (parameter size), the attacker's burden grows exponentially. This makes the Matrix-within-Matrix Matrix Cryptosystem an ideal candidate for high-security environments where mathematical "weight" is a priority.

The Engineering Behind Hill-GF-4D: Why Rust?

1. Memory Management: Stack vs. Heap

In the Hill-GF-4D system, performance is driven by how Rust manages data storage. Small, fixed-size variables like the modulus (MOD) or loop counters are stored on the Stack, which is incredibly fast but limited in size. However, the massive 4D matrix structures—which can reach 84MB in size—are allocated on the Heap.

By using the Vec<Vec<u64>> type, the program requests large contiguous blocks of memory from the system's RAM. This allows the Hill-GF-4D engine to handle millions of elements without crashing, as the Heap can expand to use all available system memory while the Stack stays light and responsive.

2. Zero-Cost Abstractions and Ownership

One of the primary reasons for the system's speed (such as the 0.101s inversion for 10x20 keys) is Rust's Ownership model. In many languages, a "Garbage Collector" must periodically stop the program to clean up memory, causing "stutters" during heavy math operations. Rust eliminates this.

When the block_invert function finishes, Rust automatically deallocates the temporary "Augmented" matrices. This "deterministic destruction" ensures that the program's memory footprint doesn't grow uncontrollably.

3. Explicit Cloning and Data Integrity

In the source code, you will notice the .clone() command frequently used before matrix multiplications. Because matrices live on the Heap, Rust's safety rules prevent multiple parts of the program from changing the same data at the same time.

By explicitly cloning the factor matrix (line 97), the program creates a dedicated "workspace" for that specific row-reduction step. This ensures that the original Encryption Key remains untouched and corruption-free while the Decryption Key is being mathematically derived.

4. Optimized Mathematical Throughput

By leveraging unsigned 64-bit integers (u64), the program bypasses the complexity of signed-bit logic. Modern CPUs can process u64 operations in parallel more efficiently than signed counterparts. This architectural choice is why the relationship between key size and timing remains remarkably consistent, allowing the system to scale to 5.3 million elements without hitting a performance wall.

5. Cache-Optimized Design (The "Conveyor Belt" Effect)

Modern CPUs have a tiny, lightning-fast memory called the L3 Cache. The Hill-GF-4D system is highly efficient because it stores data in "contiguous blocks"—meaning all the numbers for a sub-matrix are lined up right next to each other in RAM.

This allows the CPU to treat the data like a conveyor belt: it predicts what the program needs next and pre-loads it into the L3 cache before the math even starts. Even though your largest key files (84MB) are bigger than most caches, the program stays fast (0.101s to 164s) because it only ever "looks" at one small, perfectly-sized block at a time, never forcing the CPU to wait for slow system RAM.

6. Parallel Potential (The "Factory Floor")

The math in Hill-GF-4D is "embarrassingly parallel". In the main.rs code, when we multiply the outer DIM × DIM blocks, each of those smaller SUB_DIM multiplications is completely independent of the others.

• Keyspace Expansion: Moving to GF(p^n) means the number of possible "elements" in the matrix grows from p to p^n. Even for small degrees like 8, the search space for an attacker becomes astronomical.
• Non-Linearity: While the matrix operations are linear over the extension field, the relationship between the individual integer coefficients becomes highly complex due to the polynomial reduction step (the modulo of the irreducible polynomial).
• Efficiency: As the timings showed, even though the math becomes much more sophisticated, the actual encryption of data remains incredibly fast because it still relies on the fundamental efficiency of matrix-vector multiplication.

// 1. Initialize random Encryption Matrix (E)
    let mut encr_data = vec![vec![vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM]; DIM]; DIM];
    for i in 0..DIM {
        for j in 0..DIM {
            for r in 0..SUB_DIM {
                for c in 0..SUB_DIM {
                    let mut coeffs = vec![0; DEGREE];
                    for k in 0..DEGREE { coeffs[k] = rng.gen_range(0..MOD); }
                    encr_data[i][j][r][c] = Poly { coeffs };
                }
            }
        }
    }

use rand::Rng;
use std::time::Instant; 
use rayon::prelude::*;
use std::fs::File;
use std::io::{Write, BufWriter, Read}; 

const MOD: u64 = 997727;
const DEGREE: usize = 8;
const DIM: usize = 8;        
const SUB_DIM: usize = 8;   

lazy_static::lazy_static! {
    static ref IRREDUCIBLE: Vec<u64> = {
        let mut v = vec![0; DEGREE + 1];
        v[DEGREE] = 1; 
        v[1] = 1;      
        v[0] = 1;      
        v
    };
}

#[derive(Clone, Debug, PartialEq, Eq)]
struct Poly {
    coeffs: Vec<u64>,
}

impl Poly {
    fn zero() -> Self { Poly { coeffs: vec![0; DEGREE] } }
    fn identity() -> Self {
        let mut p = Self::zero();
        p.coeffs[0] = 1;
        p
    }
    fn add(&self, other: &Self) -> Self {
        let mut res = vec![0; DEGREE];
        for i in 0..DEGREE { res[i] = (self.coeffs[i] + other.coeffs[i]) % MOD; }
        Poly { coeffs: res }
    }
    fn sub(&self, other: &Self) -> Self {
        let mut res = vec![0; DEGREE];
        for i in 0..DEGREE { res[i] = (self.coeffs[i] + MOD - other.coeffs[i]) % MOD; }
        Poly { coeffs: res }
    }
    fn mul(&self, other: &Self) -> Self {
        let mut prod = vec![0u64; 2 * DEGREE - 1];
        for i in 0..DEGREE {
            for j in 0..DEGREE {
                prod[i + j] = (prod[i + j] + (self.coeffs[i] * other.coeffs[j])) % MOD;
            }
        }
        for i in (DEGREE..prod.len()).rev() {
            let factor = prod[i];
            for j in 0..DEGREE {
                let term = (factor * IRREDUCIBLE[j]) % MOD;
                prod[i - DEGREE + j] = (prod[i - DEGREE + j] + MOD - term) % MOD;
            }
        }
        Poly { coeffs: prod[0..DEGREE].to_vec() }
    }
}

// --- Math Helpers ---
fn mod_inv(a: u64, m: u64) -> u64 {
    let mut mn = (m as i64, a as i64);
    let mut xy = (0, 1);
    while mn.1 != 0 {
        xy = (xy.1, xy.0 - (mn.0 / mn.1) * xy.1);
        mn = (mn.1, mn.0 % mn.1);
    }
    let mut res = xy.0;
    if res < 0 { res += m as i64; }
    res as u64
}

fn poly_inv(poly: &Poly) -> Option<Poly> {
    let mut r0 = IRREDUCIBLE.clone();
    let mut r1 = poly.coeffs.clone();
    while r1.len() > 0 && r1[r1.len()-1] == 0 { r1.pop(); }
    if r1.is_empty() { return None; }
    let mut t0 = vec![0];
    let mut t1 = vec![1];
    while !r1.is_empty() {
        let mut q = vec![0; r0.len() - r1.len() + 1];
        let inv_lead = mod_inv(r1[r1.len() - 1], MOD);
        let mut temp_r = r0.clone();
        for i in (0..q.len()).rev() {
            q[i] = (temp_r[r1.len() - 1 + i] * inv_lead) % MOD;
            for j in 0..r1.len() {
                let sub = (q[i] * r1[j]) % MOD;
                temp_r[i + j] = (temp_r[i + j] + MOD - sub) % MOD;
            }
        }
        r0 = r1;
        while temp_r.len() > 0 && temp_r[temp_r.len()-1] == 0 { temp_r.pop(); }
        r1 = temp_r;
        let mut q_t1 = vec![0; q.len() + t1.len() - 1];
        for i in 0..q.len() {
            for j in 0..t1.len() { q_t1[i+j] = (q_t1[i+j] + q[i] * t1[j]) % MOD; }
        }
        let mut next_t = t0.clone();
        if next_t.len() < q_t1.len() { next_t.resize(q_t1.len(), 0); }
        for i in 0..q_t1.len() { next_t[i] = (next_t[i] + MOD - q_t1[i]) % MOD; }
        t0 = t1; t1 = next_t;
    }
    let scale = mod_inv(r0[0], MOD);
    let mut final_coeffs = vec![0; DEGREE];
    for i in 0..t0.len().min(DEGREE) { final_coeffs[i] = (t0[i] * scale) % MOD; }
    Some(Poly { coeffs: final_coeffs })
}

type Matrix = Vec<Vec<Poly>>;
struct BlockMatrix { data: Vec<Vec<Matrix>> }

fn get_identity_sub() -> Matrix {
    let mut m = vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM];
    for i in 0..SUB_DIM { m[i][i] = Poly::identity(); }
    m
}
fn get_zero_sub() -> Matrix { vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM] }

fn sub_mat_mul(a: &Matrix, b: &Matrix) -> Matrix {
    let mut res = vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM];
    for i in 0..SUB_DIM {
        for j in 0..SUB_DIM {
            let mut sum = Poly::zero();
            for k in 0..SUB_DIM { sum = sum.add(&a[i][k].mul(&b[k][j])); }
            res[i][j] = sum;
        }
    }
    res
}

fn block_mat_mul(a: &BlockMatrix, b: &BlockMatrix) -> BlockMatrix {
    let mut res_data = vec![vec![get_zero_sub(); DIM]; DIM];
    for i in 0..DIM {
        for j in 0..DIM {
            let mut sum_mat = get_zero_sub();
            for k in 0..DIM {
                let prod = sub_mat_mul(&a.data[i][k], &b.data[k][j]);
                for r in 0..SUB_DIM {
                    for c in 0..SUB_DIM {
                        sum_mat[r][c] = sum_mat[r][c].add(&prod[r][c]);
                    }
                }
            }
            res_data[i][j] = sum_mat;
        }
    }
    BlockMatrix { data: res_data }
}

fn sub_mat_inv(a: &Matrix) -> Option<Matrix> {
    let mut aug = vec![vec![Poly::zero(); 2 * SUB_DIM]; SUB_DIM];
    for i in 0..SUB_DIM {
        for j in 0..SUB_DIM { aug[i][j] = a[i][j].clone(); }
        aug[i][SUB_DIM + i] = Poly::identity();
    }
    for i in 0..SUB_DIM {
        let mut pivot = i;
        while pivot < SUB_DIM && aug[pivot][i] == Poly::zero() { pivot += 1; }
        if pivot == SUB_DIM { return None; }
        aug.swap(i, pivot);
        let inv = poly_inv(&aug[i][i])?;
        for j in 0..2 * SUB_DIM { aug[i][j] = aug[i][j].mul(&inv); }
        for k in 0..SUB_DIM {
            if k != i {
                let factor = aug[k][i].clone();
                for j in 0..2 * SUB_DIM {
                    let sub = factor.mul(&aug[i][j]);
                    aug[k][j] = aug[k][j].sub(&sub);
                }
            }
        }
    }
    let mut res = vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM];
    for i in 0..SUB_DIM {
        for j in 0..SUB_DIM { res[i][j] = aug[i][SUB_DIM + j].clone(); }
    }
    Some(res)
}

fn block_invert(bm: &BlockMatrix) -> BlockMatrix {
    let mut aug = vec![vec![vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM]; 2 * DIM]; DIM];
    for i in 0..DIM {
        for j in 0..DIM { aug[i][j] = bm.data[i][j].clone(); }
        for j in DIM..2 * DIM { aug[i][j] = if j - DIM == i { get_identity_sub() } else { get_zero_sub() }; }
    }
    for i in 0..DIM {
        let p_inv = sub_mat_inv(&aug[i][i]).expect("Singular Block");
        aug[i].par_iter_mut().for_each(|block| *block = sub_mat_mul(&p_inv, block));
        let pivot_row = aug[i].clone();
        aug.par_iter_mut().enumerate().for_each(|(k, row)| {
            if k != i {
                let factor = row[i].clone();
                for j in 0..2 * DIM {
                    let prod = sub_mat_mul(&factor, &pivot_row[j]);
                    for r in 0..SUB_DIM {
                        for c in 0..SUB_DIM { row[j][r][c] = row[j][r][c].sub(&prod[r][c]); }
                    }
                }
            }
        });
    }
    let mut inv_data = vec![vec![get_zero_sub(); DIM]; DIM];
    for i in 0..DIM { for j in 0..DIM { inv_data[i][j] = aug[i][DIM + j].clone(); } }
    BlockMatrix { data: inv_data }
}

fn process_block(p_vec: &Vec<Vec<Poly>>, key: &BlockMatrix) -> Vec<Vec<Poly>> {
    let mut res_vec = vec![vec![Poly::zero(); SUB_DIM]; DIM];
    for j in 0..DIM {
        for sc in 0..SUB_DIM {
            let mut sum = Poly::zero();
            for i in 0..DIM {
                for sr in 0..SUB_DIM {
                    sum = sum.add(&p_vec[i][sr].mul(&key.data[i][j][sr][sc]));
                }
            }
            res_vec[j][sc] = sum;
        }
    }
    res_vec
}

// --- IO Logic ---

fn save_nested_data(data: &Vec<Vec<Poly>>, filename: &str, header: &str) -> std::io::Result<()> {
    let file = File::create(filename)?;
    let mut writer = BufWriter::new(file);
    writeln!(writer, "{}", header)?;
    writeln!(writer, "MOD: {}, DIM: {}, SUB_DIM: {}, DEGREE: {}", MOD, DIM, SUB_DIM, DEGREE)?;
    for i in 0..DIM {
        write!(writer, "Block-Row [{}]: [", i)?;
        for j in 0..SUB_DIM {
            let poly = &data[i][j];
            write!(writer, "(")?;
            for k in 0..DEGREE {
                write!(writer, "{}", poly.coeffs[k])?;
                if k < DEGREE - 1 { write!(writer, ",")?; }
            }
            write!(writer, ")")?;
            if j < SUB_DIM - 1 { write!(writer, ", ")?; }
        }
        writeln!(writer, "]")?;
    }
    Ok(())
}

fn save_key_pair_file(encr: &BlockMatrix, decr: &BlockMatrix, filename: &str) -> std::io::Result<()> {
    let file = File::create(filename)?;
    let mut writer = BufWriter::new(file);
    writeln!(writer, "--- HILL-GF-POLY FULL KEY PAIR ---")?;
    writeln!(writer, "MOD: {}, DIM: {}, SUB_DIM: {}, DEGREE: {}", MOD, DIM, SUB_DIM, DEGREE)?;
    writeln!(writer, "\n=== ENCRYPTION MATRIX (E) ==="")?;
    for i in 0..DIM {
        for j in 0..DIM {
            write!(writer, "E[{},{}]: [", i, j)?;
            for r in 0..SUB_DIM {
                for c in 0..SUB_DIM {
                    let poly = &encr.data[i][j][r][c];
                    write!(writer, "(")?;
                    for k in 0..DEGREE {
                        write!(writer, "{}", poly.coeffs[k])?;
                        if k < DEGREE - 1 { write!(writer, ",")?; }
                    }
                    write!(writer, ")")?;
                    if r == SUB_DIM - 1 && c == SUB_DIM - 1 { } else { write!(writer, ", ")?; }
                }
            }
            writeln!(writer, "]")?;
        }
    }
    writeln!(writer, "\n=== DECRYPTION MATRIX (D) ==="")?;
    for i in 0..DIM {
        for j in 0..DIM {
            write!(writer, "D[{},{}]: [", i, j)?;
            for r in 0..SUB_DIM {
                for c in 0..SUB_DIM {
                    let poly = &decr.data[i][j][r][c];
                    write!(writer, "(")?;
                    for k in 0..DEGREE {
                        write!(writer, "{}", poly.coeffs[k])?;
                        if k < DEGREE - 1 { write!(writer, ",")?; }
                    }
                    write!(writer, ")")?;
                    if r == SUB_DIM - 1 && c == SUB_DIM - 1 { } else { write!(writer, ", ")?; }
                }
            }
            writeln!(writer, "]")?;
        }
    }
    Ok(())
}

fn main() {
    let mut rng = rand::thread_rng();
    println!("--- Hill-GF-Poly 5D Cipher Engine ---");
    println!("CONFIG: MOD={}, DIM={}, SUB_DIM={}, DEGREE={}", MOD, DIM, SUB_DIM, DEGREE);

    // 1. Initialize random Encryption Matrix (E)
    let mut encr_data = vec![vec![vec![vec![Poly::zero(); SUB_DIM]; SUB_DIM]; DIM]; DIM];
    for i in 0..DIM {
        for j in 0..DIM {
            for r in 0..SUB_DIM {
                for c in 0..SUB_DIM {
                    let mut coeffs = vec![0; DEGREE];
                    for k in 0..DEGREE { coeffs[k] = rng.gen_range(0..MOD); }
                    encr_data[i][j][r][c] = Poly { coeffs };
                }
            }
        }
    }
    let encr = BlockMatrix { data: encr_data };
    println!("Step 1: Encryption Matrix E initialized.");

    // 2. Calculate functional Inverse (D)
    let start_keys = Instant::now();
    let decr = block_invert(&encr);
    println!("Step 2: Decryption Matrix D calculated in {:?}", start_keys.elapsed());

    // 3. Write Key Pair to disk
    let key_filename = format!("POLY-M{}-D{}-S{}-DEG{}.key", MOD, DIM, SUB_DIM, DEGREE);
    save_key_pair_file(&encr, &decr, &key_filename).expect("Failed to save keys");
    println!("Step 3: Key file [{}] written.", key_filename);

    // 4. EXPLICIT CHECK: E * D = I
    println!("Step 4: Performing explicit Matrix Check (E * D)...");
    let start_check = Instant::now();
    let result_mat = block_mat_mul(&encr, &decr);
    println!("Step 4: Matrix multiplication completed in {:?}", start_check.elapsed());

    // 5. Verify Identity
    let mut is_identity = true;
    for i in 0..DIM {
        for j in 0..DIM {
            let expected = if i == j { get_identity_sub() } else { get_zero_sub() };
            if result_mat.data[i][j] != expected { is_identity = false; }
        }
    }
    if is_identity {
        println!("✅ Step 5: SUCCESS! E * D = Identity Matrix.");
    } else {
        println!("❌ Step 5: FAILURE! Matrices are not inverses.");
    }

    // 6. Load and display Plaintext
    let mut f = File::open("plaintext.txt").expect("Could not find plaintext.txt");
    let mut raw_bytes = Vec::new();
    f.read_to_end(&mut raw_bytes).unwrap();
    let block_bytes = DIM * SUB_DIM * DEGREE;
    let input_content = if raw_bytes.len() >= block_bytes {
        raw_bytes[0..block_bytes].to_vec()
    } else {
        let mut padded = raw_bytes.clone();
        while padded.len() < block_bytes { padded.push(b' '); }
        padded
    };
    println!("\nStep 6: Original Plaintext Block:");
    println!("\"{}\"\n", String::from_utf8_lossy(&input_content));

    // 7. Convert text to Polynomial Vector
    let mut p_vec = vec![vec![Poly::zero(); SUB_DIM]; DIM];
    let mut b_idx = 0;
    for i in 0..DIM {
        for j in 0..SUB_DIM {
            let mut c = vec![0; DEGREE];
            for k in 0..DEGREE { c[k] = input_content[b_idx] as u64; b_idx += 1; }
            p_vec[i][j] = Poly { coeffs: c };
        }
    }
    println!("Step 7: Plaintext converted to polynomial vector.");

    // 8. Encrypt
    let enc_start = Instant::now();
    let cipher_vec = process_block(&p_vec, &encr);
    save_nested_data(&cipher_vec, "ciphertext.txt", "POLY-HILL-5D-CIPHERTEXT").unwrap();
    println!("Step 8: Encryption complete in {:?}. Saved to ciphertext.txt", enc_start.elapsed());

    // 9. Decrypt (Functional Verification)
    let dec_start = Instant::now();
    let recovered_vec = process_block(&cipher_vec, &decr);
    println!("Step 9: Decryption complete in {:?}.", dec_start.elapsed());

    // 10. Recover text and display
    let mut recovered_bytes = Vec::new();
    for i in 0..DIM {
        for j in 0..SUB_DIM {
            for k in 0..DEGREE { recovered_bytes.push(recovered_vec[i][j].coeffs[k] as u8); }
        }
    }
    println!("\nStep 10: Recovered Plaintext Block:");
    println!("\"{}\"\n", String::from_utf8_lossy(&recovered_bytes));

    // 11. Final Integrity check
    if String::from_utf8_lossy(&input_content).trim_end() == String::from_utf8_lossy(&recovered_bytes).trim_end() {
        println!("✅ Step 11: SUCCESS! Full data cycle verified.");
    } else {
        println!("❌ Step 11: FAILURE! Data mismatch detected.");
    }
}

Plaintext: The Preamble to the Constitution of the United States

The Encryption using 16x16x16x16 4D matrix with polynomials in GF(997727^8)

# cargo run --bin parallel-poly --release
   Compiling matrix-within-matrix v0.1.0 (/var/www/html/crypto/GF997727/matrix-within-matrix)
    Finished `release` profile [optimized] target(s) in 2.81s
     Running `target/release/parallel-poly`
--- Hill-GF-Poly 5D Cipher Engine ---
CONFIG: MOD=997727, DIM=16, SUB_DIM=16, DEGREE=8
Step 1: Encryption Matrix E initialized.
Step 2: Decryption Matrix D calculated in 2.7023054s
Step 3: Key file [POLY-M997727-D16-S16-DEG8.key] written.
Step 4: Performing explicit Matrix Check (E * D)...
Step 4: Matrix multiplication completed in 10.421212572s
✅ Step 5: SUCCESS! E * D = Identity Matrix.
Step 6: Original Plaintext Block:
"We the People of the United States, in Order to form a more perfect Union, establish Justice, insure domestic Tranquility, provide for the common defence, 
promote the general Welfare, and secure the Blessings of Liberty to ourselves and our Posterity, do ordain and establish this Constitution for the United States of America."
Step 7: Plaintext converted to polynomial vector.
Step 8: Encryption complete in 45.473363ms. Saved to ciphertext.txt
Step 9: Decryption complete in 45.029343ms.
Step 10: Recovered Plaintext Block:
"We the People of the United States, in Order to form a more perfect Union, establish Justice, insure domestic Tranquility, provide for the common defence,
 promote the general Welfare, and secure the Blessings of Liberty to ourselves and our Posterity, do ordain and establish this Constitution for the United States of America."
✅ Step 11: SUCCESS! Full data cycle verified.

Ciphertext using 16x16 matrices inside every cell of a 16x16 matrix, filled with 8-coefficient polynomials in GF(997727^8)