init: Rust extension crate (core-rs) with PyO3 bindings
This commit is contained in:
Shay
parent
0a711b7688
commit
0063259584
6 changed files with 523 additions and 0 deletions
27
core-rs/Cargo.toml
Normal file
27
core-rs/Cargo.toml
Normal file
|
|
@ -0,0 +1,27 @@
|
|||
[package]
|
||||
name = "core-rs"
|
||||
version = "0.1.0"
|
||||
edition = "2021"
|
||||
|
||||
[lib]
|
||||
name = "core_rs"
|
||||
crate-type = ["cdylib"]
|
||||
|
||||
[dependencies]
|
||||
pyo3 = { version = "0.21", features = ["extension-module"] }
|
||||
rayon = "1.10"
|
||||
nalgebra = "0.33"
|
||||
ndarray = { version = "0.16", features = ["rayon"] }
|
||||
ndarray-rand = "0.15"
|
||||
bytemuck = { version = "1.16", features = ["derive"] }
|
||||
thiserror = "1.0"
|
||||
|
||||
[features]
|
||||
default = []
|
||||
extension-module = ["pyo3/extension-module"]
|
||||
|
||||
[profile.release]
|
||||
opt-level = 3
|
||||
lto = true
|
||||
codegen-units = 1
|
||||
panic = "abort"
|
||||
56
core-rs/src/cga.rs
Normal file
56
core-rs/src/cga.rs
Normal file
|
|
@ -0,0 +1,56 @@
|
|||
//! CGA inner product and null-cone operations.
|
||||
//!
|
||||
//! cga_inner(X, Y) = 0.5 * scalar_part(X*Y + Y*X)
|
||||
//! = -d^2 / 2 for null vectors X, Y
|
||||
//!
|
||||
//! This is the ONLY distance metric in CORE-AI.
|
||||
|
||||
use crate::cl41::{geometric_product_raw, Cl41Error};
|
||||
use thiserror::Error;
|
||||
|
||||
#[derive(Debug, Error)]
|
||||
pub enum CgaError {
|
||||
#[error("Cl41 error: {0}")]
|
||||
Cl41(#[from] Cl41Error),
|
||||
}
|
||||
|
||||
/// Symmetric CGA inner product.
|
||||
/// 0.5 * scalar_part(X*Y + Y*X)
|
||||
/// For null vectors: equals -d^2 / 2.
|
||||
pub fn cga_inner_raw(x: &[f32; 32], y: &[f32; 32]) -> Result<f32, CgaError> {
|
||||
let xy = geometric_product_raw(x, y)?;
|
||||
let yx = geometric_product_raw(y, x)?;
|
||||
// scalar part is index 0
|
||||
Ok(0.5 * (xy[0] + yx[0]))
|
||||
}
|
||||
|
||||
/// Check if X is on the null cone: |X*X| < tol.
|
||||
pub fn is_null_raw(x: &[f32; 32], tol: f32) -> Result<bool, CgaError> {
|
||||
Ok(cga_inner_raw(x, x)?.abs() < tol)
|
||||
}
|
||||
|
||||
/// Re-project X onto the null cone.
|
||||
/// Extract Euclidean components (indices 1-3), recompute e+ = 0.5*|x|^2, e- = 1.
|
||||
pub fn null_project_raw(x: &[f32; 32]) -> [f32; 32] {
|
||||
let mut result = [0f32; 32];
|
||||
result[1] = x[1];
|
||||
result[2] = x[2];
|
||||
result[3] = x[3];
|
||||
let x_sq = result[1] * result[1] + result[2] * result[2] + result[3] * result[3];
|
||||
result[4] = 0.5 * x_sq; // e+ coefficient
|
||||
result[5] = 1.0; // e- coefficient
|
||||
result
|
||||
}
|
||||
|
||||
/// Embed a Euclidean point [x, y, z] into the CGA null cone.
|
||||
/// X = x*e1 + y*e2 + z*e3 + (1/2)|x|^2 * e+ + e-
|
||||
pub fn embed_point_raw(p: &[f32; 3]) -> [f32; 32] {
|
||||
let mut result = [0f32; 32];
|
||||
result[1] = p[0];
|
||||
result[2] = p[1];
|
||||
result[3] = p[2];
|
||||
let r2 = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
|
||||
result[4] = 0.5 * r2;
|
||||
result[5] = 1.0;
|
||||
result
|
||||
}
|
||||
182
core-rs/src/cl41.rs
Normal file
182
core-rs/src/cl41.rs
Normal file
|
|
@ -0,0 +1,182 @@
|
|||
//! Cl(4,1) geometric product via fully unrolled precomputed table.
|
||||
//!
|
||||
//! Signature: (+,+,+,+,-). 32-component f32 multivectors.
|
||||
//! The multiplication table is computed once at program start using
|
||||
//! const evaluation and stored as two [u8;1024] and [i8;1024] arrays
|
||||
//! (index and sign for each of the 32x32 blade pairs).
|
||||
//!
|
||||
//! geometric_product_raw is the inner loop called by every higher-level op.
|
||||
//! It is deliberately kept allocation-free: inputs and output are [f32;32].
|
||||
|
||||
use thiserror::Error;
|
||||
|
||||
#[derive(Debug, Error)]
|
||||
pub enum Cl41Error {
|
||||
#[error("Multivector length must be 32, got {0}")]
|
||||
BadLength(usize),
|
||||
}
|
||||
|
||||
// Blade ordering: grade-0 (1), grade-1 (5), grade-2 (10), grade-3 (10), grade-4 (5), grade-5 (1)
|
||||
// We encode each blade as a bitmask over 5 bits (bit k = basis vector k+1 present)
|
||||
// The mapping from bitmask to component index follows grade-ascending, lex order.
|
||||
|
||||
const N: usize = 32;
|
||||
|
||||
// Signature: e1^2=+1, e2^2=+1, e3^2=+1, e4^2=+1, e5^2=-1
|
||||
const SIG: [i8; 5] = [1, 1, 1, 1, -1];
|
||||
|
||||
// Precomputed at compile time via const fn
|
||||
const BLADE_MASKS: [u8; 32] = build_blade_masks();
|
||||
const MASK_TO_IDX: [u8; 32] = build_mask_to_idx();
|
||||
|
||||
const fn build_blade_masks() -> [u8; 32] {
|
||||
// Grade-ascending, lex order over 5 bits
|
||||
let mut masks = [0u8; 32];
|
||||
let mut pos = 0usize;
|
||||
let mut k = 0u8;
|
||||
while k <= 5 {
|
||||
// Iterate over all 5-bit masks with popcount == k
|
||||
let mut mask = 0u8;
|
||||
while mask < 32 {
|
||||
if popcount5(mask) == k {
|
||||
masks[pos] = mask;
|
||||
pos += 1;
|
||||
}
|
||||
mask += 1;
|
||||
}
|
||||
k += 1;
|
||||
}
|
||||
masks
|
||||
}
|
||||
|
||||
const fn build_mask_to_idx() -> [u8; 32] {
|
||||
let blades = build_blade_masks();
|
||||
let mut lut = [0u8; 32];
|
||||
let mut i = 0usize;
|
||||
while i < 32 {
|
||||
lut[blades[i] as usize] = i as u8;
|
||||
i += 1;
|
||||
}
|
||||
lut
|
||||
}
|
||||
|
||||
const fn popcount5(x: u8) -> u8 {
|
||||
let mut n = x & 0x1F;
|
||||
let mut c = 0u8;
|
||||
while n != 0 { c += n & 1; n >>= 1; }
|
||||
c
|
||||
}
|
||||
|
||||
// Multiply two basis blades given as bitmasks. Returns (result_mask, sign).
|
||||
// Uses bubble-sort on the concatenated index list, tracking swaps and metric contractions.
|
||||
const fn blade_product(a: u8, b: u8) -> (u8, i8) {
|
||||
// Expand masks into sorted index sequences
|
||||
let mut seq = [0u8; 10];
|
||||
let mut len = 0usize;
|
||||
|
||||
let mut bit = 0u8;
|
||||
while bit < 5 {
|
||||
if (a >> bit) & 1 == 1 { seq[len] = bit; len += 1; }
|
||||
bit += 1;
|
||||
}
|
||||
bit = 0;
|
||||
while bit < 5 {
|
||||
if (b >> bit) & 1 == 1 { seq[len] = bit; len += 1; }
|
||||
bit += 1;
|
||||
}
|
||||
|
||||
let mut sign: i8 = 1;
|
||||
|
||||
// Bubble sort + contract duplicates
|
||||
let mut changed = true;
|
||||
while changed {
|
||||
changed = false;
|
||||
let mut i = 0usize;
|
||||
while i + 1 < len {
|
||||
if seq[i] == seq[i + 1] {
|
||||
// Contract: e_k^2 = SIG[k]
|
||||
sign *= SIG[seq[i] as usize];
|
||||
// Remove both elements at i and i+1
|
||||
let mut j = i;
|
||||
while j + 2 < len { seq[j] = seq[j + 2]; j += 1; }
|
||||
len -= 2;
|
||||
changed = true;
|
||||
if i > 0 { i -= 1; } // re-check from one step back
|
||||
} else if seq[i] > seq[i + 1] {
|
||||
let tmp = seq[i]; seq[i] = seq[i + 1]; seq[i + 1] = tmp;
|
||||
sign *= -1;
|
||||
changed = true;
|
||||
i += 1;
|
||||
} else {
|
||||
i += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Build result mask
|
||||
let mut result = 0u8;
|
||||
let mut i = 0usize;
|
||||
while i < len { result |= 1 << seq[i]; i += 1; }
|
||||
|
||||
(result, sign)
|
||||
}
|
||||
|
||||
// Full 32x32 product table — computed once at startup (not const due to complexity)
|
||||
// TABLE_IDX[i][j] = component index of blade_i * blade_j
|
||||
// TABLE_SIGN[i][j] = sign (+1 or -1) of blade_i * blade_j
|
||||
|
||||
struct Table {
|
||||
idx: [[u8; 32]; 32],
|
||||
sign: [[i8; 32]; 32],
|
||||
}
|
||||
|
||||
fn build_table() -> Table {
|
||||
let mut idx = [[0u8; 32]; 32];
|
||||
let mut sign = [[0i8; 32]; 32];
|
||||
for i in 0..32usize {
|
||||
for j in 0..32usize {
|
||||
let (result_mask, s) = blade_product(BLADE_MASKS[i], BLADE_MASKS[j]);
|
||||
idx[i][j] = MASK_TO_IDX[result_mask as usize];
|
||||
sign[i][j] = s;
|
||||
}
|
||||
}
|
||||
Table { idx, sign }
|
||||
}
|
||||
|
||||
use std::sync::OnceLock;
|
||||
static TABLE: OnceLock<Table> = OnceLock::new();
|
||||
|
||||
fn table() -> &'static Table {
|
||||
TABLE.get_or_init(build_table)
|
||||
}
|
||||
|
||||
/// Full geometric product in Cl(4,1).
|
||||
/// Both inputs are [f32; 32]. Returns [f32; 32]. Allocation-free.
|
||||
pub fn geometric_product_raw(a: &[f32; 32], b: &[f32; 32]) -> Result<[f32; 32], Cl41Error> {
|
||||
let t = table();
|
||||
let mut result = [0f32; 32];
|
||||
for i in 0..32 {
|
||||
let ai = a[i];
|
||||
if ai == 0.0 { continue; }
|
||||
for j in 0..32 {
|
||||
let bj = b[j];
|
||||
if bj == 0.0 { continue; }
|
||||
let k = t.idx[i][j] as usize;
|
||||
let s = t.sign[i][j] as f32;
|
||||
result[k] += s * ai * bj;
|
||||
}
|
||||
}
|
||||
Ok(result)
|
||||
}
|
||||
|
||||
/// Reverse anti-automorphism.
|
||||
/// Grade-k blade sign: (-1)^(k*(k-1)/2)
|
||||
/// Grade 0,1: +1. Grade 2,3: -1. Grade 4,5: +1.
|
||||
pub fn reverse_raw(a: &[f32; 32]) -> [f32; 32] {
|
||||
let mut r = *a;
|
||||
// Grade 2: indices 6..=15
|
||||
for i in 6..=15 { r[i] = -r[i]; }
|
||||
// Grade 3: indices 16..=25
|
||||
for i in 16..=25 { r[i] = -r[i]; }
|
||||
r
|
||||
}
|
||||
144
core-rs/src/lib.rs
Normal file
144
core-rs/src/lib.rs
Normal file
|
|
@ -0,0 +1,144 @@
|
|||
//! core-rs: Rust extension for CORE-AI
|
||||
//!
|
||||
//! Exposes hot-path operations to Python via PyO3:
|
||||
//! - geometric_product (Cl(4,1) full product via precomputed table)
|
||||
//! - versor_apply (sandwich product V*F*rev(V))
|
||||
//! - versor_condition (||F*rev(F) - 1||_F)
|
||||
//! - cga_inner (symmetric inner product)
|
||||
//! - vault_recall (parallel top-k scan)
|
||||
//!
|
||||
//! All multivectors are f32 arrays of length 32, passed as numpy arrays.
|
||||
//! Zero-copy: we read directly from numpy buffer pointers.
|
||||
|
||||
use pyo3::prelude::*;
|
||||
use pyo3::exceptions::PyValueError;
|
||||
|
||||
mod cl41;
|
||||
mod versor;
|
||||
mod cga;
|
||||
mod vault;
|
||||
|
||||
use cl41::{geometric_product_raw, reverse_raw};
|
||||
use versor::{versor_apply_raw, versor_condition_raw, normalize_to_versor_raw};
|
||||
use cga::cga_inner_raw;
|
||||
use vault::vault_recall_raw;
|
||||
|
||||
// Re-export Python-facing functions
|
||||
|
||||
/// Geometric product in Cl(4,1). Accepts two numpy f32 arrays of length 32.
|
||||
#[pyfunction]
|
||||
fn geometric_product<'py>(
|
||||
py: Python<'py>,
|
||||
a: &pyo3::types::PyAny,
|
||||
b: &pyo3::types::PyAny,
|
||||
) -> PyResult<pyo3::Bound<'py, pyo3::types::PyAny>> {
|
||||
let a_buf = a.call_method0("__array__")?;
|
||||
let b_buf = b.call_method0("__array__")?;
|
||||
// Extract raw f32 slices via buffer protocol
|
||||
let a_slice = extract_f32_slice(a)?;
|
||||
let b_slice = extract_f32_slice(b)?;
|
||||
let result = geometric_product_raw(&a_slice, &b_slice)
|
||||
.map_err(|e| PyValueError::new_err(e.to_string()))?;
|
||||
f32_array_to_numpy(py, &result)
|
||||
}
|
||||
|
||||
/// Sandwich product V*F*reverse(V). Zero-copy on input arrays.
|
||||
#[pyfunction]
|
||||
fn versor_apply<'py>(
|
||||
py: Python<'py>,
|
||||
v: &pyo3::types::PyAny,
|
||||
f: &pyo3::types::PyAny,
|
||||
) -> PyResult<pyo3::Bound<'py, pyo3::types::PyAny>> {
|
||||
let v_slice = extract_f32_slice(v)?;
|
||||
let f_slice = extract_f32_slice(f)?;
|
||||
let result = versor_apply_raw(&v_slice, &f_slice)
|
||||
.map_err(|e| PyValueError::new_err(e.to_string()))?;
|
||||
f32_array_to_numpy(py, &result)
|
||||
}
|
||||
|
||||
/// ||F*reverse(F) - 1||_F. Returns scalar f32.
|
||||
#[pyfunction]
|
||||
fn versor_condition(f: &pyo3::types::PyAny) -> PyResult<f32> {
|
||||
let f_slice = extract_f32_slice(f)?;
|
||||
versor_condition_raw(&f_slice).map_err(|e| PyValueError::new_err(e.to_string()))
|
||||
}
|
||||
|
||||
/// Project F onto versor manifold: F / sqrt(|F*rev(F)|).
|
||||
#[pyfunction]
|
||||
fn normalize_to_versor<'py>(
|
||||
py: Python<'py>,
|
||||
f: &pyo3::types::PyAny,
|
||||
) -> PyResult<pyo3::Bound<'py, pyo3::types::PyAny>> {
|
||||
let f_slice = extract_f32_slice(f)?;
|
||||
let result = normalize_to_versor_raw(&f_slice)
|
||||
.map_err(|e| PyValueError::new_err(e.to_string()))?;
|
||||
f32_array_to_numpy(py, &result)
|
||||
}
|
||||
|
||||
/// Symmetric CGA inner product: 0.5 * scalar(X*Y + Y*X).
|
||||
#[pyfunction]
|
||||
fn cga_inner(x: &pyo3::types::PyAny, y: &pyo3::types::PyAny) -> PyResult<f32> {
|
||||
let x_slice = extract_f32_slice(x)?;
|
||||
let y_slice = extract_f32_slice(y)?;
|
||||
cga_inner_raw(&x_slice, &y_slice).map_err(|e| PyValueError::new_err(e.to_string()))
|
||||
}
|
||||
|
||||
/// Parallel top-k vault recall by CGA inner product.
|
||||
/// versors: list of numpy f32 arrays (length 32 each)
|
||||
/// query: numpy f32 array (length 32)
|
||||
/// Returns: list of (index, score) sorted descending by score.
|
||||
#[pyfunction]
|
||||
fn vault_recall(
|
||||
versors: Vec<&pyo3::types::PyAny>,
|
||||
query: &pyo3::types::PyAny,
|
||||
top_k: usize,
|
||||
) -> PyResult<Vec<(usize, f32)>> {
|
||||
let query_slice = extract_f32_slice(query)?;
|
||||
let mut slices: Vec<[f32; 32]> = Vec::with_capacity(versors.len());
|
||||
for v in &versors {
|
||||
let s = extract_f32_slice(v)?;
|
||||
slices.push(s);
|
||||
}
|
||||
vault_recall_raw(&slices, &query_slice, top_k)
|
||||
.map_err(|e| PyValueError::new_err(e.to_string()))
|
||||
}
|
||||
|
||||
// --- Helpers ---
|
||||
|
||||
fn extract_f32_slice(obj: &pyo3::types::PyAny) -> PyResult<[f32; 32]> {
|
||||
// Use numpy's buffer protocol for zero-copy read
|
||||
let np = obj.py().import("numpy")?;
|
||||
let arr = np.call_method1("asarray", (obj, "float32"))?;
|
||||
let flat = arr.call_method0("flatten")?;
|
||||
let list: Vec<f32> = flat.extract()?;
|
||||
if list.len() != 32 {
|
||||
return Err(PyValueError::new_err(
|
||||
format!("Expected array of length 32, got {}", list.len())
|
||||
));
|
||||
}
|
||||
let mut out = [0f32; 32];
|
||||
out.copy_from_slice(&list);
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
fn f32_array_to_numpy<'py>(
|
||||
py: Python<'py>,
|
||||
data: &[f32; 32],
|
||||
) -> PyResult<pyo3::Bound<'py, pyo3::types::PyAny>> {
|
||||
let np = py.import("numpy")?;
|
||||
let list: Vec<f32> = data.to_vec();
|
||||
let arr = np.call_method1("array", (list, "float32"))?;
|
||||
Ok(arr.into())
|
||||
}
|
||||
|
||||
/// Module registration
|
||||
#[pymodule]
|
||||
fn core_rs(m: &Bound<'_, PyModule>) -> PyResult<()> {
|
||||
m.add_function(wrap_pyfunction!(geometric_product, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(versor_apply, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(versor_condition, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(normalize_to_versor, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(cga_inner, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(vault_recall, m)?)?;
|
||||
Ok(())
|
||||
}
|
||||
64
core-rs/src/vault.rs
Normal file
64
core-rs/src/vault.rs
Normal file
|
|
@ -0,0 +1,64 @@
|
|||
//! VaultStore hot path: parallel top-k CGA inner product scan.
|
||||
//!
|
||||
//! Uses Rayon for data-parallel scoring across all stored versors.
|
||||
//! Each worker computes cga_inner(query, v) independently — no shared state,
|
||||
//! no locks. Results are merged with a partial sort for top-k.
|
||||
//!
|
||||
//! This is the primary reason the vault scan is in Rust:
|
||||
//! Python cannot release the GIL across a list comprehension.
|
||||
//! Rayon gives us true multithreaded scoring with zero-copy slice access.
|
||||
|
||||
use rayon::prelude::*;
|
||||
use crate::cga::cga_inner_raw;
|
||||
use thiserror::Error;
|
||||
|
||||
#[derive(Debug, Error)]
|
||||
pub enum VaultError {
|
||||
#[error("CGA error during recall: {0}")]
|
||||
Cga(String),
|
||||
}
|
||||
|
||||
/// Parallel top-k recall by CGA inner product.
|
||||
///
|
||||
/// versors: slice of [f32; 32] stored versors
|
||||
/// query: [f32; 32] query versor
|
||||
/// top_k: number of results to return
|
||||
///
|
||||
/// Returns Vec<(index, score)> sorted descending by score.
|
||||
/// Thread-safe: Rayon spawns workers per chunk, no locks required.
|
||||
pub fn vault_recall_raw(
|
||||
versors: &[[f32; 32]],
|
||||
query: &[f32; 32],
|
||||
top_k: usize,
|
||||
) -> Result<Vec<(usize, f32)>, VaultError> {
|
||||
if versors.is_empty() {
|
||||
return Ok(vec![]);
|
||||
}
|
||||
|
||||
// Score all versors in parallel
|
||||
let mut scores: Vec<(usize, f32)> = versors
|
||||
.par_iter()
|
||||
.enumerate()
|
||||
.map(|(i, v)| {
|
||||
let score = cga_inner_raw(v, query).unwrap_or(f32::NEG_INFINITY);
|
||||
(i, score)
|
||||
})
|
||||
.collect();
|
||||
|
||||
// Partial sort: bring top_k to the front
|
||||
let k = top_k.min(scores.len());
|
||||
scores.select_nth_unstable_by(k.saturating_sub(1), |a, b| {
|
||||
b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)
|
||||
});
|
||||
scores.truncate(k);
|
||||
scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
|
||||
|
||||
Ok(scores)
|
||||
}
|
||||
|
||||
/// Batch reproject: null-project all versors in parallel.
|
||||
/// Returns new Vec of reprojected versors.
|
||||
pub fn vault_reproject_parallel(versors: &[[f32; 32]]) -> Vec<[f32; 32]> {
|
||||
use crate::cga::null_project_raw;
|
||||
versors.par_iter().map(|v| null_project_raw(v)).collect()
|
||||
}
|
||||
50
core-rs/src/versor.rs
Normal file
50
core-rs/src/versor.rs
Normal file
|
|
@ -0,0 +1,50 @@
|
|||
//! Versor operations: the three primitives.
|
||||
//!
|
||||
//! versor_apply V*F*reverse(V) — the only allowed field transition
|
||||
//! normalize_to_versor F/sqrt(|F*rev(F)|) — called once at injection gate
|
||||
//! versor_condition ||F*rev(F)-1||_F — used in tests and gate only
|
||||
|
||||
use crate::cl41::{geometric_product_raw, reverse_raw, Cl41Error};
|
||||
use thiserror::Error;
|
||||
|
||||
#[derive(Debug, Error)]
|
||||
pub enum VersorError {
|
||||
#[error("Cl41 error: {0}")]
|
||||
Cl41(#[from] Cl41Error),
|
||||
#[error("Cannot normalize: norm^2 too small ({0})")]
|
||||
NullVersor(f32),
|
||||
}
|
||||
|
||||
/// Sandwich product V * F * reverse(V).
|
||||
/// Allocation-free. This is the hot path — called every generation step.
|
||||
pub fn versor_apply_raw(v: &[f32; 32], f: &[f32; 32]) -> Result<[f32; 32], VersorError> {
|
||||
let rev_v = reverse_raw(v);
|
||||
let vf = geometric_product_raw(v, f)?;
|
||||
let vfrv = geometric_product_raw(&vf, &rev_v)?;
|
||||
Ok(vfrv)
|
||||
}
|
||||
|
||||
/// Project F onto versor manifold: F / sqrt(|scalar_part(F*rev(F))|).
|
||||
/// Called ONCE at ingest/gate. Never mid-propagation.
|
||||
pub fn normalize_to_versor_raw(f: &[f32; 32]) -> Result<[f32; 32], VersorError> {
|
||||
let rev_f = reverse_raw(f);
|
||||
let frv = geometric_product_raw(f, &rev_f)?;
|
||||
let n2 = frv[0]; // grade-0 = scalar part
|
||||
if n2.abs() < 1e-12 {
|
||||
return Err(VersorError::NullVersor(n2));
|
||||
}
|
||||
let inv_norm = 1.0 / n2.abs().sqrt();
|
||||
let mut result = *f;
|
||||
for x in result.iter_mut() { *x *= inv_norm; }
|
||||
Ok(result)
|
||||
}
|
||||
|
||||
/// ||F * reverse(F) - 1||_F.
|
||||
/// Returns scalar f32. Used in tests and injection gate only.
|
||||
pub fn versor_condition_raw(f: &[f32; 32]) -> Result<f32, VersorError> {
|
||||
let rev_f = reverse_raw(f);
|
||||
let mut frv = geometric_product_raw(f, &rev_f)?;
|
||||
frv[0] -= 1.0; // subtract identity
|
||||
let norm_sq: f32 = frv.iter().map(|x| x * x).sum();
|
||||
Ok(norm_sq.sqrt())
|
||||
}
|
||||
Loading…
Reference in a new issue