qfall_math/rational/mat_q/
gso.rs

1// Copyright © 2023 Marcel Luca Schmidt
2//
3// This file is part of qFALL-math.
4//
5// qFALL-math is free software: you can redistribute it and/or modify it under
6// the terms of the Mozilla Public License Version 2.0 as published by the
7// Mozilla Foundation. See <https://mozilla.org/en-US/MPL/2.0/>.
8
9//! Implementation of the Gram-Schmidt Orthogonalization for [`MatQ`] matrices.
10
11use super::MatQ;
12use crate::traits::MatrixDimensions;
13use flint_sys::fmpq_mat::fmpq_mat_gso;
14
15impl MatQ {
16    /// Computes the Gram-Schmidt Orthogonalization of the matrix and returns a [`MatQ`] with the corresponding matrix.
17    ///
18    /// # Examples
19    /// ```
20    /// use qfall_math::rational::MatQ;
21    /// use std::str::FromStr;
22    ///
23    /// let mat = MatQ::from_str("[[1/2, 1],[2, 5]]").unwrap();
24    /// let mat_gso = mat.gso();
25    /// ```
26    pub fn gso(&self) -> Self {
27        let mut out = MatQ::new(self.get_num_rows(), self.get_num_columns());
28        unsafe {
29            fmpq_mat_gso(&mut out.matrix, &self.matrix);
30        }
31        out
32    }
33}
34
35#[cfg(test)]
36mod test_gso {
37    use crate::{
38        rational::{MatQ, Q},
39        traits::{MatrixGetEntry, MatrixGetSubmatrix},
40    };
41    use std::str::FromStr;
42
43    /// Ensure that the generated vectors from the gso are orthogonal
44    #[test]
45    fn gso_works() {
46        let mat = MatQ::from_str(
47        "[[-1, 2/7, 3/9, 4, 5/2],[-123/1000, 235/5, 123, 643/7172721, 123],[124/8981, 212, 452/2140, 12/5, 1],[0, 0, 0, 1, 1],[1, 2, 3, 4/3, 1]]",
48    )
49    .unwrap();
50
51        let mat_gso = mat.gso();
52
53        let cmp = Q::ZERO;
54        for i in 0..5 {
55            let vec_i = mat_gso.get_column(i).unwrap();
56            assert_ne!(cmp, vec_i.dot_product(&vec_i).unwrap());
57            for j in i + 1..5 {
58                let vec_j = mat_gso.get_column(j).unwrap();
59                assert_eq!(cmp, vec_i.dot_product(&vec_j).unwrap());
60            }
61        }
62    }
63
64    /// Ensure that the generated vectors have the expected values
65    #[test]
66    fn gso_correct_values() {
67        let mat = MatQ::from_str("[[1, -1, 1],[1, 0, 1],[1, 1, 2]]").unwrap();
68
69        let mat_gso = mat.gso();
70
71        assert_eq!(
72            MatQ::from_str("[[1, -1, 1/6],[1, 0, -1/3],[1, 1, 1/6]]").unwrap(),
73            mat_gso
74        );
75    }
76
77    /// Ensure that gso works with independent vectors (more columns than rows)
78    #[test]
79    fn gso_dependent_columns() {
80        let mat = MatQ::from_str("[[1, 2, 3, 4, 4],[1, 2, 3, 4, 5]]").unwrap();
81
82        let mat_gso = mat.gso();
83
84        let cmp = Q::ZERO;
85        for i in 0..5 {
86            let vec_i = mat_gso.get_column(i).unwrap();
87            for j in i + 1..5 {
88                let vec_j = mat_gso.get_column(j).unwrap();
89                assert_eq!(cmp, vec_i.dot_product(&vec_j).unwrap());
90            }
91        }
92        let vec_1 = mat_gso.get_column(0).unwrap();
93        assert_ne!(cmp, vec_1.dot_product(&vec_1).unwrap());
94    }
95
96    /// Ensure that gso works with more rows than columns
97    #[test]
98    fn gso_dependent_rows() {
99        let mat = MatQ::from_str("[[1, 2/7],[1, 2/7],[10, -2],[0, 4],[0, 0]]").unwrap();
100
101        let mat_gso = mat.gso();
102
103        let cmp = Q::ZERO;
104        for i in 0..2 {
105            let vec_i = mat_gso.get_column(i).unwrap();
106            assert_ne!(cmp, vec_i.dot_product(&vec_i).unwrap());
107            for j in i + 1..2 {
108                let vec_j = mat_gso.get_column(j).unwrap();
109                assert_eq!(cmp, vec_i.dot_product(&vec_j).unwrap());
110            }
111        }
112    }
113
114    /// Ensure that gso works with large values
115    #[test]
116    fn gso_large_values() {
117        let mat = MatQ::from_str(&format!(
118            "[[1, {}/7, 2],[1, 2/{}, 10],[10, -2, 8]]",
119            i64::MAX,
120            i64::MAX
121        ))
122        .unwrap();
123
124        let mat_gso = mat.gso();
125
126        let cmp = Q::ZERO;
127        for i in 0..3 {
128            let vec_i = mat_gso.get_column(i).unwrap();
129            assert_ne!(cmp, vec_i.dot_product(&vec_i).unwrap());
130            for j in i + 1..3 {
131                let vec_j = mat_gso.get_column(j).unwrap();
132                assert_eq!(cmp, vec_i.dot_product(&vec_j).unwrap());
133            }
134        }
135    }
136
137    /// Ensure that gso works on edge cases
138    #[test]
139    fn gso_edge_cases() {
140        let mat_1 = MatQ::from_str("[[1]]").unwrap();
141        let mat_2 = MatQ::from_str("[[1, 2/2, 3/7]]").unwrap();
142        let mat_3 = MatQ::from_str("[[1],[2/2],[3/7]]").unwrap();
143
144        let mat_1_gso = mat_1.gso();
145        let mat_2_gso = mat_2.gso();
146        let mat_3_gso = mat_3.gso();
147
148        assert_eq!(Q::ONE, mat_1_gso.get_entry(0, 0).unwrap());
149
150        let cmp = Q::ZERO;
151        for i in 0..3 {
152            for j in i + 1..3 {
153                let vec_i = mat_2_gso.get_column(i).unwrap();
154                let vec_j = mat_2_gso.get_column(j).unwrap();
155                assert_eq!(cmp, vec_i.dot_product(&vec_j).unwrap());
156            }
157        }
158        let vec_1 = mat_2_gso.get_column(0).unwrap();
159        assert_ne!(cmp, vec_1.dot_product(&vec_1).unwrap());
160
161        assert_eq!(mat_3, mat_3_gso);
162    }
163}
qfall_math/rational/mat_q/gso.rs

qfall_math/rational/mat_q/
gso.rs
