qualia_core_db/solvers/calculus/
analysis.rs

1//! Finite-dimensional analysis foundations used by native calculus.
2
3#[derive(Debug, Clone, Copy, PartialEq)]
4pub enum AnalysisError {
5    NonFinite,
6    Degenerate,
7    Singular,
8    InvalidDomain,
9    NotCertified,
10    IterationLimit { residual: f64 },
11    OutputBufferFull,
12}
13
14#[repr(C)]
15#[derive(Debug, Clone, Copy, PartialEq)]
16pub struct Vector<const N: usize> {
17    pub data: [f64; N],
18}
19
20impl<const N: usize> Vector<N> {
21    pub const fn new(data: [f64; N]) -> Self {
22        Self { data }
23    }
24
25    pub const fn zero() -> Self {
26        Self { data: [0.0; N] }
27    }
28
29    pub fn validate(&self) -> Result<(), AnalysisError> {
30        if self.data.iter().all(|value| value.is_finite()) {
31            Ok(())
32        } else {
33            Err(AnalysisError::NonFinite)
34        }
35    }
36
37    pub fn dot(self, other: Self) -> Result<f64, AnalysisError> {
38        self.validate()?;
39        other.validate()?;
40        Ok(self
41            .data
42            .iter()
43            .zip(other.data)
44            .map(|(left, right)| left * right)
45            .sum())
46    }
47
48    pub fn norm(self) -> Result<f64, AnalysisError> {
49        Ok(self.dot(self)?.sqrt())
50    }
51
52    pub fn distance(self, other: Self) -> Result<f64, AnalysisError> {
53        self.sub(other).norm()
54    }
55
56    pub fn add(self, other: Self) -> Self {
57        let mut result = [0.0; N];
58        for (index, value) in result.iter_mut().enumerate() {
59            *value = self.data[index] + other.data[index];
60        }
61        Self::new(result)
62    }
63
64    pub fn sub(self, other: Self) -> Self {
65        let mut result = [0.0; N];
66        for (index, value) in result.iter_mut().enumerate() {
67            *value = self.data[index] - other.data[index];
68        }
69        Self::new(result)
70    }
71
72    pub fn scale(self, scalar: f64) -> Result<Self, AnalysisError> {
73        if !scalar.is_finite() {
74            return Err(AnalysisError::NonFinite);
75        }
76        let mut result = self;
77        for value in &mut result.data {
78            *value *= scalar;
79        }
80        result.validate()?;
81        Ok(result)
82    }
83
84    pub fn project_onto(self, direction: Self) -> Result<Self, AnalysisError> {
85        let denominator = direction.dot(direction)?;
86        if denominator <= f64::MIN_POSITIVE {
87            return Err(AnalysisError::Degenerate);
88        }
89        direction.scale(self.dot(direction)? / denominator)
90    }
91}
92
93#[repr(C)]
94#[derive(Debug, Clone, Copy, PartialEq)]
95pub struct LinearMap<const R: usize, const C: usize> {
96    pub coefficients: [[f64; C]; R],
97}
98
99impl<const R: usize, const C: usize> LinearMap<R, C> {
100    pub const fn new(coefficients: [[f64; C]; R]) -> Self {
101        Self { coefficients }
102    }
103
104    pub fn validate(&self) -> Result<(), AnalysisError> {
105        if self
106            .coefficients
107            .iter()
108            .flatten()
109            .all(|value| value.is_finite())
110        {
111            Ok(())
112        } else {
113            Err(AnalysisError::NonFinite)
114        }
115    }
116
117    pub fn apply(&self, input: Vector<C>) -> Result<Vector<R>, AnalysisError> {
118        self.validate()?;
119        input.validate()?;
120        let mut output = [0.0; R];
121        for (row, value) in output.iter_mut().enumerate() {
122            *value = self.coefficients[row]
123                .iter()
124                .zip(input.data)
125                .map(|(coefficient, component)| coefficient * component)
126                .sum();
127        }
128        Ok(Vector::new(output))
129    }
130
131    pub fn transpose(&self) -> LinearMap<C, R> {
132        let mut result = [[0.0; R]; C];
133        for (row, coefficients) in self.coefficients.iter().enumerate() {
134            for (column, coefficient) in coefficients.iter().enumerate() {
135                result[column][row] = *coefficient;
136            }
137        }
138        LinearMap::new(result)
139    }
140}
141
142impl<const N: usize> LinearMap<N, N> {
143    pub const fn identity() -> Self {
144        let mut coefficients = [[0.0; N]; N];
145        let mut index = 0;
146        while index < N {
147            coefficients[index][index] = 1.0;
148            index += 1;
149        }
150        Self::new(coefficients)
151    }
152
153    pub fn trace(&self) -> Result<f64, AnalysisError> {
154        self.validate()?;
155        Ok((0..N).map(|index| self.coefficients[index][index]).sum())
156    }
157
158    pub fn determinant(&self) -> Result<f64, AnalysisError> {
159        self.validate()?;
160        let mut matrix = self.coefficients;
161        let mut determinant = 1.0;
162        let mut sign = 1.0;
163        for column in 0..N {
164            let mut pivot = column;
165            for row in column + 1..N {
166                if matrix[row][column].abs() > matrix[pivot][column].abs() {
167                    pivot = row;
168                }
169            }
170            if matrix[pivot][column].abs() <= f64::EPSILON {
171                return Ok(0.0);
172            }
173            if pivot != column {
174                matrix.swap(pivot, column);
175                sign = -sign;
176            }
177            let diagonal = matrix[column][column];
178            determinant *= diagonal;
179            for row in column + 1..N {
180                let factor = matrix[row][column] / diagonal;
181                for trailing in column + 1..N {
182                    matrix[row][trailing] -= factor * matrix[column][trailing];
183                }
184            }
185        }
186        Ok(sign * determinant)
187    }
188
189    pub fn solve(&self, rhs: Vector<N>) -> Result<Vector<N>, AnalysisError> {
190        self.validate()?;
191        rhs.validate()?;
192        let mut matrix = self.coefficients;
193        let mut values = rhs.data;
194        for column in 0..N {
195            let mut pivot = column;
196            for row in column + 1..N {
197                if matrix[row][column].abs() > matrix[pivot][column].abs() {
198                    pivot = row;
199                }
200            }
201            if matrix[pivot][column].abs() <= 64.0 * f64::EPSILON {
202                return Err(AnalysisError::Singular);
203            }
204            matrix.swap(pivot, column);
205            values.swap(pivot, column);
206            for row in column + 1..N {
207                let factor = matrix[row][column] / matrix[column][column];
208                matrix[row][column] = 0.0;
209                for trailing in column + 1..N {
210                    matrix[row][trailing] -= factor * matrix[column][trailing];
211                }
212                values[row] -= factor * values[column];
213            }
214        }
215        let mut solution = [0.0; N];
216        for row in (0..N).rev() {
217            let mut value = values[row];
218            for column in row + 1..N {
219                value -= matrix[row][column] * solution[column];
220            }
221            solution[row] = value / matrix[row][row];
222        }
223        let result = Vector::new(solution);
224        result.validate()?;
225        Ok(result)
226    }
227}
228
229#[repr(C)]
230#[derive(Debug, Clone, Copy, PartialEq)]
231pub struct Basis<const N: usize> {
232    /// Basis vectors are stored as columns.
233    pub matrix: LinearMap<N, N>,
234}
235
236impl<const N: usize> Basis<N> {
237    pub fn new(matrix: LinearMap<N, N>) -> Result<Self, AnalysisError> {
238        if matrix.determinant()?.abs() <= 64.0 * f64::EPSILON {
239            return Err(AnalysisError::Singular);
240        }
241        Ok(Self { matrix })
242    }
243
244    pub fn from_coordinates(&self, coordinates: Vector<N>) -> Result<Vector<N>, AnalysisError> {
245        self.matrix.apply(coordinates)
246    }
247
248    pub fn to_coordinates(&self, vector: Vector<N>) -> Result<Vector<N>, AnalysisError> {
249        self.matrix.solve(vector)
250    }
251}
252
253#[repr(C)]
254#[derive(Debug, Clone, Copy, PartialEq)]
255pub struct Interval {
256    pub lower: f64,
257    pub upper: f64,
258}
259
260impl Interval {
261    pub fn new(lower: f64, upper: f64) -> Result<Self, AnalysisError> {
262        if !lower.is_finite() || !upper.is_finite() || lower > upper {
263            return Err(AnalysisError::InvalidDomain);
264        }
265        Ok(Self { lower, upper })
266    }
267
268    pub fn contains(self, value: f64) -> bool {
269        value.is_finite() && value >= self.lower && value <= self.upper
270    }
271
272    pub fn uniform_cover(
273        self,
274        radius: f64,
275        out_centers: &mut [f64],
276    ) -> Result<usize, AnalysisError> {
277        if !radius.is_finite() || radius <= 0.0 {
278            return Err(AnalysisError::InvalidDomain);
279        }
280        let required = (((self.upper - self.lower) / (2.0 * radius)).ceil() as usize).max(1);
281        if out_centers.len() < required {
282            return Err(AnalysisError::OutputBufferFull);
283        }
284        let width = (self.upper - self.lower) / required as f64;
285        for (index, center) in out_centers[..required].iter_mut().enumerate() {
286            *center = self.lower + (index as f64 + 0.5) * width;
287        }
288        Ok(required)
289    }
290}
291
292#[repr(C)]
293#[derive(Debug, Clone, Copy, PartialEq)]
294pub struct FixedPointCertificate<const N: usize> {
295    pub point: Vector<N>,
296    pub contraction_factor: f64,
297    pub iterations: u32,
298    pub residual: f64,
299    pub a_posteriori_error_bound: f64,
300}
301
302pub fn contraction_fixed_point<const N: usize, F>(
303    map: F,
304    initial: Vector<N>,
305    contraction_factor: f64,
306    tolerance: f64,
307    max_iterations: u32,
308) -> Result<FixedPointCertificate<N>, AnalysisError>
309where
310    F: Fn(Vector<N>) -> Vector<N>,
311{
312    if !contraction_factor.is_finite()
313        || !(0.0..1.0).contains(&contraction_factor)
314        || !tolerance.is_finite()
315        || tolerance <= 0.0
316    {
317        return Err(AnalysisError::NotCertified);
318    }
319    initial.validate()?;
320    let mut current = initial;
321    for iteration in 1..=max_iterations {
322        let next = map(current);
323        next.validate()?;
324        let residual = next.distance(current)?;
325        let error_bound = contraction_factor * residual / (1.0 - contraction_factor);
326        if error_bound <= tolerance {
327            return Ok(FixedPointCertificate {
328                point: next,
329                contraction_factor,
330                iterations: iteration,
331                residual,
332                a_posteriori_error_bound: error_bound,
333            });
334        }
335        current = next;
336    }
337    let residual = map(current).distance(current)?;
338    Err(AnalysisError::IterationLimit { residual })
339}
340
341#[repr(C)]
342#[derive(Debug, Clone, Copy, PartialEq)]
343pub struct Complex64 {
344    pub re: f64,
345    pub im: f64,
346}
347
348impl Complex64 {
349    pub const I: Self = Self { re: 0.0, im: 1.0 };
350
351    pub const fn new(re: f64, im: f64) -> Self {
352        Self { re, im }
353    }
354
355    pub fn norm(self) -> f64 {
356        self.re.hypot(self.im)
357    }
358
359    pub fn add(self, other: Self) -> Self {
360        Self::new(self.re + other.re, self.im + other.im)
361    }
362
363    pub fn mul(self, other: Self) -> Self {
364        Self::new(
365            self.re * other.re - self.im * other.im,
366            self.re * other.im + self.im * other.re,
367        )
368    }
369
370    pub fn exp(self) -> Self {
371        let magnitude = self.re.exp();
372        Self::new(magnitude * self.im.cos(), magnitude * self.im.sin())
373    }
374
375    /// Principal logarithm with argument in `(-pi, pi]`.
376    pub fn principal_log(self) -> Result<Self, AnalysisError> {
377        let norm = self.norm();
378        if !norm.is_finite() || norm == 0.0 {
379            return Err(AnalysisError::InvalidDomain);
380        }
381        Ok(Self::new(norm.ln(), self.im.atan2(self.re)))
382    }
383
384    /// Principal square root with non-negative real part.
385    pub fn principal_sqrt(self) -> Result<Self, AnalysisError> {
386        if !self.re.is_finite() || !self.im.is_finite() {
387            return Err(AnalysisError::NonFinite);
388        }
389        let magnitude = self.norm();
390        let re = ((magnitude + self.re) * 0.5).max(0.0).sqrt();
391        let im_magnitude = ((magnitude - self.re) * 0.5).max(0.0).sqrt();
392        let im = if self.im < 0.0 {
393            -im_magnitude
394        } else {
395            im_magnitude
396        };
397        Ok(Self::new(re, im))
398    }
399}
400
401#[cfg(test)]
402mod tests {
403    use super::*;
404
405    #[test]
406    fn basis_coordinate_round_trip_and_singular_rejection() {
407        let basis = Basis::new(LinearMap::new([[1.0, 1.0], [0.0, 2.0]])).unwrap();
408        let coordinates = Vector::new([3.0, -1.0]);
409        let vector = basis.from_coordinates(coordinates).unwrap();
410        let recovered = basis.to_coordinates(vector).unwrap();
411        assert_eq!(recovered, coordinates);
412        assert_eq!(
413            Basis::new(LinearMap::new([[1.0, 2.0], [2.0, 4.0]])),
414            Err(AnalysisError::Singular)
415        );
416    }
417
418    #[test]
419    fn projection_is_idempotent_and_residual_is_orthogonal() {
420        let vector = Vector::new([2.0, 3.0, 4.0]);
421        let direction = Vector::new([1.0, -1.0, 0.0]);
422        let projection = vector.project_onto(direction).unwrap();
423        let second = projection.project_onto(direction).unwrap();
424        assert!(projection.distance(second).unwrap() < 1e-14);
425        assert!(vector.sub(projection).dot(direction).unwrap().abs() < 1e-14);
426    }
427
428    #[test]
429    fn contraction_certificate_contains_analytic_fixed_point() {
430        let certificate = contraction_fixed_point(
431            |x: Vector<1>| Vector::new([0.5 * x.data[0] + 1.0]),
432            Vector::new([0.0]),
433            0.5,
434            1e-10,
435            128,
436        )
437        .unwrap();
438        assert!((certificate.point.data[0] - 2.0).abs() <= 1e-10);
439        assert!(certificate.a_posteriori_error_bound <= 1e-10);
440    }
441
442    #[test]
443    fn complex_principal_branches_are_explicit_and_consistent() {
444        let z = Complex64::new(-3.0, 4.0);
445        let root = z.principal_sqrt().unwrap();
446        let squared = root.mul(root);
447        assert!((squared.re - z.re).abs() < 1e-14);
448        assert!((squared.im - z.im).abs() < 1e-14);
449
450        let value = Complex64::new(0.3, -0.7);
451        let round_trip = value.exp().principal_log().unwrap();
452        assert!((round_trip.re - value.re).abs() < 1e-14);
453        assert!((round_trip.im - value.im).abs() < 1e-14);
454    }
455
456    #[test]
457    fn interval_cover_reports_capacity() {
458        let interval = Interval::new(0.0, 1.0).unwrap();
459        let mut too_small = [0.0; 2];
460        assert_eq!(
461            interval.uniform_cover(0.1, &mut too_small),
462            Err(AnalysisError::OutputBufferFull)
463        );
464        let mut centers = [0.0; 5];
465        assert_eq!(interval.uniform_cover(0.1, &mut centers), Ok(5));
466    }
467}
qualia_core_db/solvers/calculus/analysis.rs

qualia_core_db/solvers/calculus/
analysis.rs
