qfall_math/integer/mat_poly_over_z/
tensor.rs

1// Copyright © 2023 Marvin Beckmann
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//! This module contains the implementation of the `tensor` product.
10
11use super::MatPolyOverZ;
12use crate::{
13    integer::PolyOverZ,
14    traits::{MatrixDimensions, MatrixGetEntry, Tensor},
15};
16use flint_sys::{fmpz_poly::fmpz_poly_mul, fmpz_poly_mat::fmpz_poly_mat_entry};
17
18impl Tensor for MatPolyOverZ {
19    /// Computes the tensor product of `self` with `other`.
20    ///
21    /// Parameters:
22    /// - `other`: the value with which the tensor product is computed.
23    ///
24    /// Returns the tensor product of `self` with `other`.
25    ///
26    /// # Examples
27    /// ```
28    /// use qfall_math::integer::MatPolyOverZ;
29    /// use qfall_math::traits::Tensor;
30    /// use std::str::FromStr;
31    ///
32    /// let mat_1 = MatPolyOverZ::from_str("[[1  1, 2  1 1]]").unwrap();
33    /// let mat_2 = MatPolyOverZ::from_str("[[1  1, 1  2]]").unwrap();
34    ///
35    /// let mat_ab = mat_1.tensor_product(&mat_2);
36    /// let mat_ba = mat_2.tensor_product(&mat_1);
37    ///
38    /// let res_ab = "[[1  1, 1  2, 2  1 1, 2  2 2]]";
39    /// let res_ba = "[[1  1, 2  1 1, 1  2, 2  2 2]]";
40    /// assert_eq!(mat_ab, MatPolyOverZ::from_str(res_ab).unwrap());
41    /// assert_eq!(mat_ba, MatPolyOverZ::from_str(res_ba).unwrap());
42    /// ```
43    fn tensor_product(&self, other: &Self) -> Self {
44        let mut out = MatPolyOverZ::new(
45            self.get_num_rows() * other.get_num_rows(),
46            self.get_num_columns() * other.get_num_columns(),
47        );
48
49        for i in 0..self.get_num_rows() {
50            for j in 0..self.get_num_columns() {
51                let entry = unsafe { self.get_entry_unchecked(i, j) };
52
53                if !entry.is_zero() {
54                    unsafe { set_matrix_window_mul(&mut out, i, j, entry, other) }
55                }
56            }
57        }
58
59        out
60    }
61}
62
63/// This function sets a specific window of the provided matrix `out`
64/// according to the `scalar` multiple of `matrix`.
65///
66/// Sets the entries
67/// `[i*rows_other, j*columns_other]` up till `[row_left*(row_other +1), column_upper*(columns_other + 1)]`
68///
69/// Parameters:
70/// - `out`: the matrix in which the result is saved
71/// - `row_left`: defines the leftmost row of the set window
72/// - `column_upper`: defines the highest column of the set window
73/// - `scalar`: defines the value with which the part of the tensor product
74///   is calculated
75/// - `matrix`: the matrix with which the scalar is multiplied
76///   before setting the entries in `out`
77///
78/// Implicitly sets the entries of the matrix according to the definition
79/// of the tensor product.
80///
81/// # Security
82/// This function accesses memory directly without checking whether the memory is
83/// actually obtained by the matrix out.
84/// This means that this function should only be called wisely.
85/// If `row_left` or `row_upper` together with the length of the matrix exceeds the
86/// range of the matrix other memory could be overwritten.
87/// We included asserts to check whether this occurs, but we advise careful usage.
88unsafe fn set_matrix_window_mul(
89    out: &mut MatPolyOverZ,
90    row_left: i64,
91    column_upper: i64,
92    scalar: PolyOverZ,
93    matrix: &MatPolyOverZ,
94) {
95    let columns_other = matrix.get_num_columns();
96    let rows_other = matrix.get_num_rows();
97
98    assert!(row_left >= 0 && row_left + rows_other <= out.get_num_rows());
99    assert!(column_upper >= 0 && column_upper + columns_other <= out.get_num_columns());
100
101    for i_other in 0..rows_other {
102        for j_other in 0..columns_other {
103            unsafe {
104                fmpz_poly_mul(
105                    fmpz_poly_mat_entry(
106                        &out.matrix,
107                        row_left * rows_other + i_other,
108                        column_upper * columns_other + j_other,
109                    ),
110                    &scalar.poly,
111                    fmpz_poly_mat_entry(&matrix.matrix, i_other, j_other),
112                )
113            }
114        }
115    }
116}
117
118#[cfg(test)]
119mod test_tensor {
120    use crate::{
121        integer::MatPolyOverZ,
122        traits::{MatrixDimensions, Tensor},
123    };
124    use std::str::FromStr;
125
126    /// Ensure that the dimensions of the tensor product are taken over correctly.
127    #[test]
128    fn dimensions_fit() {
129        let mat_1 = MatPolyOverZ::new(17, 13);
130        let mat_2 = MatPolyOverZ::new(3, 4);
131
132        let mat_3 = mat_1.tensor_product(&mat_2);
133
134        assert_eq!(51, mat_3.get_num_rows());
135        assert_eq!(52, mat_3.get_num_columns());
136    }
137
138    /// Ensure that the tensor works correctly with identity.
139    #[test]
140    fn identity() {
141        let identity = MatPolyOverZ::identity(2, 2);
142        let mat_1 = MatPolyOverZ::from_str(&format!(
143            "[[1  1, 1  {}, 1  1],[0, 1  {}, 1  -1]]",
144            u64::MAX,
145            i64::MIN
146        ))
147        .unwrap();
148
149        let mat_2 = identity.tensor_product(&mat_1);
150        let mat_3 = mat_1.tensor_product(&identity);
151
152        let cmp_mat_2 = MatPolyOverZ::from_str(&format!(
153            "[[1  1, 1  {}, 1  1, 0, 0, 0],[0, 1  {}, 1  -1, 0, 0, 0],[0, 0, 0, 1  1, 1  {}, 1  1],[0, 0, 0, 0, 1  {}, 1  -1]]",
154            u64::MAX,
155            i64::MIN,
156            u64::MAX,
157            i64::MIN
158        ))
159        .unwrap();
160        let cmp_mat_3 = MatPolyOverZ::from_str(&format!(
161            "[[1  1, 0, 1  {}, 0, 1  1, 0],[0, 1  1, 0, 1  {}, 0, 1  1],[0, 0, 1  {}, 0, 1  -1, 0],[0, 0, 0, 1  {}, 0, 1  -1]]",
162            u64::MAX,
163            u64::MAX,
164            i64::MIN,
165            i64::MIN
166        ))
167        .unwrap();
168
169        assert_eq!(cmp_mat_2, mat_2);
170        assert_eq!(cmp_mat_3, mat_3);
171    }
172
173    /// Ensure the tensor product works where one is a vector and the other is a matrix.
174    #[test]
175    fn vector_matrix() {
176        let vector = MatPolyOverZ::from_str("[[1  1],[1  -1]]").unwrap();
177        let mat_1 = MatPolyOverZ::from_str(&format!(
178            "[[1  1, 1  {}, 1  1],[0, 1  {}, 1  -1]]",
179            u64::MAX,
180            i64::MAX
181        ))
182        .unwrap();
183
184        let mat_2 = vector.tensor_product(&mat_1);
185        let mat_3 = mat_1.tensor_product(&vector);
186
187        let cmp_mat_2 = MatPolyOverZ::from_str(&format!(
188            "[[1  1, 1  {}, 1  1],[0, 1  {}, 1  -1],[1  -1, 1  -{}, 1  -1],[0, 1  -{}, 1  1]]",
189            u64::MAX,
190            i64::MAX,
191            u64::MAX,
192            i64::MAX
193        ))
194        .unwrap();
195        let cmp_mat_3 = MatPolyOverZ::from_str(&format!(
196            "[[1  1, 1  {}, 1  1],[1  -1, 1  -{}, 1  -1],[0, 1  {}, 1  -1],[0, 1  -{}, 1  1]]",
197            u64::MAX,
198            u64::MAX,
199            i64::MAX,
200            i64::MAX
201        ))
202        .unwrap();
203
204        assert_eq!(cmp_mat_2, mat_2);
205        assert_eq!(cmp_mat_3, mat_3);
206    }
207
208    /// Ensure that the tensor product works correctly with two vectors.
209    #[test]
210    fn vector_vector() {
211        let vec_1 = MatPolyOverZ::from_str("[[1  2],[1  1]]").unwrap();
212        let vec_2 = MatPolyOverZ::from_str(&format!(
213            "[[1  {}],[1  {}]]",
214            (u64::MAX - 1) / 2,
215            i64::MIN / 2
216        ))
217        .unwrap();
218
219        let vec_3 = vec_1.tensor_product(&vec_2);
220        let vec_4 = vec_2.tensor_product(&vec_1);
221
222        let cmp_vec_3 = MatPolyOverZ::from_str(&format!(
223            "[[1  {}],[1  {}],[1  {}],[1  {}]]",
224            u64::MAX - 1,
225            i64::MIN,
226            (u64::MAX - 1) / 2,
227            i64::MIN / 2
228        ))
229        .unwrap();
230        let cmp_vec_4 = MatPolyOverZ::from_str(&format!(
231            "[[1  {}],[1  {}],[1  {}],[1  {}]]",
232            u64::MAX - 1,
233            (u64::MAX - 1) / 2,
234            i64::MIN,
235            i64::MIN / 2
236        ))
237        .unwrap();
238
239        assert_eq!(cmp_vec_3, vec_3);
240        assert_eq!(cmp_vec_4, vec_4);
241    }
242
243    /// Ensures that the tensor product works for higher degree polynomials.
244    #[test]
245    fn higher_degree() {
246        let higher_degree = MatPolyOverZ::from_str("[[1  1, 2  0 1, 2  1 1]]").unwrap();
247        let mat_1 =
248            MatPolyOverZ::from_str(&format!("[[1  1, 1  {}, 2  1 {}]]", u64::MAX, i64::MIN))
249                .unwrap();
250
251        let mat_2 = higher_degree.tensor_product(&mat_1);
252
253        let cmp_mat_2 = MatPolyOverZ::from_str(&format!(
254            "[[1  1, 1  {}, 2  1 {}, 2  0 1, 2  0 {}, 3  0 1 {}, 2  1 1, 2  {} {}, 3  1 {} {}]]",
255            u64::MAX,
256            i64::MIN,
257            u64::MAX,
258            i64::MIN,
259            u64::MAX,
260            u64::MAX,
261            i64::MIN + 1,
262            i64::MIN,
263        ))
264        .unwrap();
265
266        assert_eq!(cmp_mat_2, mat_2);
267    }
268}
qfall_math/integer/mat_poly_over_z/tensor.rs

qfall_math/integer/mat_poly_over_z/
tensor.rs
