lattice_qcd_rs/
su3.rs

1//! Module for SU(3) matrices and su(3) (that is the generators of SU(3) )
2//!
3//! The module defines the SU(3) generator we use the same matrices as on
4//! [wikipedia](https://en.wikipedia.org/w/index.php?title=Gell-Mann_matrices&oldid=988659438#Matrices)
5//! **divided by two** such that `Tr(T^a T^b) = \delta^{ab} /2 `.
6
7use std::iter::FusedIterator;
8
9use na::{base::allocator::Allocator, ComplexField, DefaultAllocator, OMatrix};
10use rand::Rng;
11use rand_distr::{Distribution, Uniform};
12
13use super::{
14    field::Su3Adjoint, su2, utils, utils::FactorialNumber, CMatrix2, CMatrix3, Complex, Real, I,
15    ONE, ZERO,
16};
17
18/// SU(3) generator
19/// ```textrust
20/// 0    0.5  0
21/// 0.5  0    0
22/// 0    0    0
23/// ```
24pub const GENERATOR_1: CMatrix3 = CMatrix3::new(
25    ZERO,
26    Complex::new(0.5_f64, 0_f64),
27    ZERO,
28    // ---
29    Complex::new(0.5_f64, 0_f64),
30    ZERO,
31    ZERO,
32    // ---
33    ZERO,
34    ZERO,
35    ZERO,
36);
37
38/// SU(3) generator
39/// ```textrust
40/// 0   -i/2   0
41/// i/2  0     0
42/// 0    0     0
43/// ```
44pub const GENERATOR_2: CMatrix3 = CMatrix3::new(
45    ZERO,
46    Complex::new(0_f64, -0.5_f64),
47    ZERO,
48    // ---
49    Complex::new(0_f64, 0.5_f64),
50    ZERO,
51    ZERO,
52    // ---
53    ZERO,
54    ZERO,
55    ZERO,
56);
57
58/// SU(3) generator
59/// ```textrust
60/// 0.5  0    0
61/// 0   -0.5  0
62/// 0    0    0
63/// ```
64pub const GENERATOR_3: CMatrix3 = CMatrix3::new(
65    Complex::new(0.5_f64, 0_f64),
66    ZERO,
67    ZERO,
68    // ---
69    ZERO,
70    Complex::new(-0.5_f64, 0_f64),
71    ZERO,
72    // ---
73    ZERO,
74    ZERO,
75    ZERO,
76);
77
78/// SU(3) generator
79/// ```textrust
80/// 0    0    0.5
81/// 0    0    0
82/// 0.5  0    0
83/// ```
84pub const GENERATOR_4: CMatrix3 = CMatrix3::new(
85    ZERO,
86    ZERO,
87    Complex::new(0.5_f64, 0_f64),
88    // ---
89    ZERO,
90    ZERO,
91    ZERO,
92    // ---
93    Complex::new(0.5_f64, 0_f64),
94    ZERO,
95    ZERO,
96);
97
98/// SU(3) generator
99/// ```textrust
100/// 0    0   -i/2
101/// 0    0    0
102/// i/2  0    0
103/// ```
104pub const GENERATOR_5: CMatrix3 = CMatrix3::new(
105    ZERO,
106    ZERO,
107    Complex::new(0_f64, -0.5_f64),
108    // ---
109    ZERO,
110    ZERO,
111    ZERO,
112    // ---
113    Complex::new(0_f64, 0.5_f64),
114    ZERO,
115    ZERO,
116);
117
118/// SU(3) generator
119/// ```textrust
120/// 0    0    0
121/// 0    0    0.5
122/// 0    0.5  0
123/// ```
124pub const GENERATOR_6: CMatrix3 = CMatrix3::new(
125    ZERO,
126    ZERO,
127    ZERO,
128    // ---
129    ZERO,
130    ZERO,
131    Complex::new(0.5_f64, 0_f64),
132    // ---
133    ZERO,
134    Complex::new(0.5_f64, 0_f64),
135    ZERO,
136);
137
138/// SU(3) generator
139/// ```textrust
140/// 0    0    0
141/// 0    0   -i/2
142/// 0    i/2  0
143/// ```
144pub const GENERATOR_7: CMatrix3 = CMatrix3::new(
145    ZERO,
146    ZERO,
147    ZERO,
148    // ---
149    ZERO,
150    ZERO,
151    Complex::new(0_f64, -0.5_f64),
152    // ---
153    ZERO,
154    Complex::new(0_f64, 0.5_f64),
155    ZERO,
156);
157
158/// SU(3) generator
159/// ```textrust
160/// 0.5/sqrt(3)  0            0
161/// 0            0.5/sqrt(3)  0
162/// 0            0           -1/sqrt(3)
163/// ```
164pub const GENERATOR_8: CMatrix3 = CMatrix3::new(
165    Complex::new(ONE_OVER_2_SQRT_3, 0_f64),
166    ZERO,
167    ZERO,
168    // ---
169    ZERO,
170    Complex::new(ONE_OVER_2_SQRT_3, 0_f64),
171    ZERO,
172    // ---
173    ZERO,
174    ZERO,
175    Complex::new(MINUS_ONE_OVER_SQRT_3, 0_f64),
176);
177
178/// -1/sqrt(3), used for [`GENERATOR_8`].
179const MINUS_ONE_OVER_SQRT_3: f64 = -0.577_350_269_189_625_8_f64;
180/// 2/sqrt(3), used for [`GENERATOR_8`].
181const ONE_OVER_2_SQRT_3: f64 = 0.288_675_134_594_812_9_f64;
182
183/// list of SU(3) generators
184/// they are normalize such that `Tr(T^a T^b) = \frac{1}{2}\delta^{ab}`
185pub const GENERATORS: [&CMatrix3; 8] = [
186    &GENERATOR_1,
187    &GENERATOR_2,
188    &GENERATOR_3,
189    &GENERATOR_4,
190    &GENERATOR_5,
191    &GENERATOR_6,
192    &GENERATOR_7,
193    &GENERATOR_8,
194];
195
196/// Exponential of matrices.
197///
198/// Note prefer using [`su3_exp_r`] and [`su3_exp_i`] when possible.
199#[deprecated(since = "0.2.0", note = "Please use nalgebra exp instead")]
200pub trait MatrixExp<T = Self> {
201    /// Return the exponential of the matrix
202    fn exp(self) -> T;
203}
204
205/// Basic implementation of matrix exponential for complex matrices.
206/// It does it by first diagonalizing the matrix then exponentiate the diagonal
207/// and retransforms it back to the original basis.
208#[allow(deprecated)]
209impl<T, D> MatrixExp<OMatrix<T, D, D>> for OMatrix<T, D, D>
210where
211    T: ComplexField + Copy,
212    D: na::DimName + na::DimSub<na::U1>,
213    DefaultAllocator: Allocator<T, D, na::DimDiff<D, na::U1>>
214        + Allocator<T, na::DimDiff<D, na::U1>>
215        + Allocator<T, D, D>
216        + Allocator<T, D>,
217    OMatrix<T, D, D>: Clone,
218{
219    fn exp(self) -> OMatrix<T, D, D> {
220        let decomposition = self.schur();
221        // a complex matrix is always diagonalizable
222        let eigens = decomposition.eigenvalues().unwrap();
223        let new_matrix = Self::from_diagonal(&eigens.map(|el| el.exp()));
224        let (q, _) = decomposition.unpack();
225        // q is always invertible
226        q.clone() * new_matrix * q.try_inverse().unwrap()
227    }
228}
229
230/// Create a [`super::CMatrix3`] of the form `[v1, v2, v1* x v2*]` where v1* is the conjugate of v1 and
231/// `x` is the cross product.
232///
233/// # Example
234/// ```
235/// # use nalgebra::{Vector3, Complex, Matrix3};
236/// # use lattice_qcd_rs::{assert_eq_matrix};
237/// # use lattice_qcd_rs::su3::create_matrix_from_2_vector;
238/// let v1 = Vector3::new(
239///     Complex::new(1_f64, 0_f64),
240///     Complex::new(0_f64, 0_f64),
241///     Complex::new(0_f64, 0_f64),
242/// );
243/// let v2 = Vector3::new(
244///     Complex::new(0_f64, 0_f64),
245///     Complex::new(1_f64, 0_f64),
246///     Complex::new(0_f64, 0_f64),
247/// );
248/// assert_eq_matrix!(
249///     Matrix3::identity(),
250///     create_matrix_from_2_vector(v1, v2),
251///     0.000_000_1_f64
252/// );
253///
254/// let v1 = Vector3::new(
255///     Complex::new(0_f64, 1_f64),
256///     Complex::new(0_f64, 0_f64),
257///     Complex::new(0_f64, 0_f64),
258/// );
259/// let v2 = Vector3::new(
260///     Complex::new(0_f64, 0_f64),
261///     Complex::new(1_f64, 0_f64),
262///     Complex::new(0_f64, 0_f64),
263/// );
264/// let m = Matrix3::new(
265///     Complex::new(0_f64, 1_f64),
266///     Complex::new(0_f64, 0_f64),
267///     Complex::new(0_f64, 0_f64),
268///     // ---
269///     Complex::new(0_f64, 0_f64),
270///     Complex::new(1_f64, 0_f64),
271///     Complex::new(0_f64, 0_f64),
272///     // ---
273///     Complex::new(0_f64, 0_f64),
274///     Complex::new(0_f64, 0_f64),
275///     Complex::new(0_f64, -1_f64),
276/// );
277/// assert_eq_matrix!(m, create_matrix_from_2_vector(v1, v2), 0.000_000_1_f64);
278/// ```
279pub fn create_matrix_from_2_vector(
280    v1: na::Vector3<Complex>,
281    v2: na::Vector3<Complex>,
282) -> na::Matrix3<Complex> {
283    // TODO find a better way
284    let cross_vec: na::Vector3<Complex> = v1.conjugate().cross(&v2.conjugate());
285    let iter = v1.iter().chain(v2.iter()).chain(cross_vec.iter()).copied();
286    na::Matrix3::<Complex>::from_iterator(iter)
287}
288
289/// get an orthonormalize matrix from two vector.
290fn ortho_matrix_from_2_vector(v1: na::Vector3<Complex>, v2: na::Vector3<Complex>) -> CMatrix3 {
291    let v1_new = v1.try_normalize(f64::EPSILON).unwrap_or(v1);
292    let v2_temp = v2 - v1_new * v1_new.conjugate().dot(&v2);
293    let v2_new = v2_temp.try_normalize(f64::EPSILON).unwrap_or(v2_temp);
294    create_matrix_from_2_vector(v1_new, v2_new)
295}
296
297/// Try orthonormalize the given matrix.
298// TODO example
299pub fn orthonormalize_matrix(matrix: &CMatrix3) -> CMatrix3 {
300    let v1 = na::Vector3::from_iterator(matrix.column(0).iter().copied());
301    let v2 = na::Vector3::from_iterator(matrix.column(1).iter().copied());
302    ortho_matrix_from_2_vector(v1, v2)
303}
304
305/// Orthonormalize the given matrix by mutating its content.
306// TODO example
307pub fn orthonormalize_matrix_mut(matrix: &mut CMatrix3) {
308    *matrix = orthonormalize_matrix(matrix);
309}
310
311/// Generate Uniformly distributed SU(3) matrix.
312///
313/// # Example
314/// ```
315/// # use lattice_qcd_rs::{assert_matrix_is_su_3,su3::random_su3};
316/// # use rand::SeedableRng;
317/// # let mut rng = rand::rngs::StdRng::seed_from_u64(0);
318/// for _ in 0..10 {
319///     assert_matrix_is_su_3!(random_su3(&mut rng), 9_f64 * f64::EPSILON);
320/// }
321/// ```
322pub fn random_su3<Rng>(rng: &mut Rng) -> CMatrix3
323where
324    Rng: rand::Rng + ?Sized,
325{
326    rand_su3_with_dis(rng, &rand::distributions::Uniform::new(-1_f64, 1_f64))
327}
328
329/// Get a random SU3 with the given distribution.
330///
331/// The given distribution can be quite opaque on the distribution of the SU(3) matrix.
332/// For a matrix Uniformly distributed among SU(3) use [`get_random_su3`].
333/// For a matrix close to unity use [`get_random_su3_close_to_unity`]
334fn rand_su3_with_dis<Rng>(rng: &mut Rng, d: &impl rand_distr::Distribution<Real>) -> CMatrix3
335where
336    Rng: rand::Rng + ?Sized,
337{
338    let mut v1 = random_vec_3(rng, d);
339    while v1.norm() <= f64::EPSILON {
340        v1 = random_vec_3(rng, d);
341    }
342    let mut v2 = random_vec_3(rng, d);
343    while v1.dot(&v2).modulus() <= f64::EPSILON {
344        v2 = random_vec_3(rng, d);
345    }
346    ortho_matrix_from_2_vector(v1, v2)
347}
348
349/// get a random [`na::Vector3<Complex>`].
350fn random_vec_3<Rng>(rng: &mut Rng, d: &impl rand_distr::Distribution<Real>) -> na::Vector3<Complex>
351where
352    Rng: rand::Rng + ?Sized,
353{
354    na::Vector3::from_fn(|_, _| Complex::new(d.sample(rng), d.sample(rng)))
355}
356
357/// Get a radom SU(3) matrix close to `[get_r] (+/- 1) * [get_s] (+/- 1) * [get_t] (+/- 1)`.
358///
359/// Note that it diverges from SU(3) slightly.
360/// `spread_parameter` should be between between 0 and 1 both excluded to generate valid data.
361/// outside this bound it will not panic but can have unexpected results.
362///
363/// # Example
364/// ```
365/// # use lattice_qcd_rs::{assert_matrix_is_su_3,su3::random_su3_close_to_unity};
366/// # use rand::SeedableRng;
367/// # let mut rng = rand::rngs::StdRng::seed_from_u64(0);
368/// for _ in 0..10 {
369///     assert_matrix_is_su_3!(
370///         random_su3_close_to_unity(0.000_000_001_f64, &mut rng),
371///         0.000_000_1_f64
372///     );
373/// }
374/// ```
375/// but it will be not close to SU(3) up to [`f64::EPSILON`].
376/// ```should_panic
377/// # use lattice_qcd_rs::{assert_matrix_is_su_3,su3::random_su3_close_to_unity};
378/// # use rand::SeedableRng;
379/// # let mut rng = rand::rngs::StdRng::seed_from_u64(0);
380/// assert_matrix_is_su_3!(
381///     random_su3_close_to_unity(0.000_000_001_f64, &mut rng),
382///     f64::EPSILON * 180_f64
383/// );
384/// ```
385pub fn random_su3_close_to_unity<R>(spread_parameter: Real, rng: &mut R) -> CMatrix3
386where
387    R: rand::Rng + ?Sized,
388{
389    let r = get_r(su2::random_su2_close_to_unity(spread_parameter, rng));
390    let s = get_s(su2::random_su2_close_to_unity(spread_parameter, rng));
391    let t = get_t(su2::random_su2_close_to_unity(spread_parameter, rng));
392    let distribution = rand::distributions::Bernoulli::new(0.5_f64).unwrap();
393    let mut x = r * s * t;
394    if distribution.sample(rng) {
395        x = x.adjoint();
396    }
397    x
398}
399
400/// Embed a Matrix2 inside Matrix3 leaving the last row and column be the same as identity.
401///
402/// # Example
403/// ```
404/// # use lattice_qcd_rs::{CMatrix3, CMatrix2, su3::get_r, Complex, assert_eq_matrix};
405/// let m1 = CMatrix2::new(
406///     Complex::from(1_f64),
407///     Complex::from(2_f64),
408///     // ---
409///     Complex::from(3_f64),
410///     Complex::from(4_f64),
411/// );
412///
413/// let m2 = CMatrix3::new(
414///     Complex::from(1_f64),
415///     Complex::from(2_f64),
416///     Complex::from(0_f64),
417///     // ---
418///     Complex::from(3_f64),
419///     Complex::from(4_f64),
420///     Complex::from(0_f64),
421///     // ---
422///     Complex::from(0_f64),
423///     Complex::from(0_f64),
424///     Complex::from(1_f64),
425/// );
426/// assert_eq_matrix!(get_r(m1), m2, f64::EPSILON);
427/// ```
428pub fn get_r(m: CMatrix2) -> CMatrix3 {
429    CMatrix3::new(
430        m[(0, 0)],
431        m[(0, 1)],
432        ZERO,
433        // ---
434        m[(1, 0)],
435        m[(1, 1)],
436        ZERO,
437        // ---
438        ZERO,
439        ZERO,
440        ONE,
441    )
442}
443
444/// Embed a Matrix2 inside Matrix3 leaving the second row and column be the same as identity.
445///
446/// # Example
447/// ```
448/// # use lattice_qcd_rs::{CMatrix3, CMatrix2, su3::get_s, Complex, assert_eq_matrix};
449/// let m1 = CMatrix2::new(
450///     Complex::from(1_f64),
451///     Complex::from(2_f64),
452///     // ---
453///     Complex::from(3_f64),
454///     Complex::from(4_f64),
455/// );
456///
457/// let m2 = CMatrix3::new(
458///     Complex::from(1_f64),
459///     Complex::from(0_f64),
460///     Complex::from(2_f64),
461///     // ---
462///     Complex::from(0_f64),
463///     Complex::from(1_f64),
464///     Complex::from(0_f64),
465///     // ---
466///     Complex::from(3_f64),
467///     Complex::from(0_f64),
468///     Complex::from(4_f64),
469/// );
470/// assert_eq_matrix!(get_s(m1), m2, f64::EPSILON);
471/// ```
472pub fn get_s(m: CMatrix2) -> CMatrix3 {
473    CMatrix3::new(
474        m[(0, 0)],
475        ZERO,
476        m[(0, 1)],
477        // ---
478        ZERO,
479        ONE,
480        ZERO,
481        // ---
482        m[(1, 0)],
483        ZERO,
484        m[(1, 1)],
485    )
486}
487
488/// Embed a Matrix2 inside Matrix3 leaving the first row and column be the same as identity.
489///
490/// # Example
491/// ```
492/// # use lattice_qcd_rs::{CMatrix3, CMatrix2, su3::get_t, Complex, assert_eq_matrix};
493/// let m1 = CMatrix2::new(
494///     Complex::from(1_f64),
495///     Complex::from(2_f64),
496///     // ---
497///     Complex::from(3_f64),
498///     Complex::from(4_f64),
499/// );
500///
501/// let m2 = CMatrix3::new(
502///     Complex::from(1_f64),
503///     Complex::from(0_f64),
504///     Complex::from(0_f64),
505///     // ---
506///     Complex::from(0_f64),
507///     Complex::from(1_f64),
508///     Complex::from(2_f64),
509///     // ---
510///     Complex::from(0_f64),
511///     Complex::from(3_f64),
512///     Complex::from(4_f64),
513/// );
514/// assert_eq_matrix!(get_t(m1), m2, f64::EPSILON);
515/// ```
516pub fn get_t(m: CMatrix2) -> CMatrix3 {
517    CMatrix3::new(
518        ONE,
519        ZERO,
520        ZERO,
521        // ---
522        ZERO,
523        m[(0, 0)],
524        m[(0, 1)],
525        // ---
526        ZERO,
527        m[(1, 0)],
528        m[(1, 1)],
529    )
530}
531
532/// get the [`CMatrix2`] sub block corresponding to [`get_r`]
533///
534/// # Example
535/// ```
536/// # use lattice_qcd_rs::{CMatrix3, CMatrix2, su3::get_sub_block_r, Complex, assert_eq_matrix};
537///  let m1 = CMatrix3::new(
538///     Complex::from(1_f64),
539///     Complex::from(2_f64),
540///     Complex::from(3_f64),
541///     // ---
542///     Complex::from(4_f64),
543///     Complex::from(5_f64),
544///     Complex::from(6_f64),
545///     // ---
546///     Complex::from(7_f64),
547///     Complex::from(8_f64),
548///     Complex::from(9_f64),
549/// );
550/// let m2 = CMatrix2::new(
551///     Complex::from(1_f64),
552///     Complex::from(2_f64),
553///     // ---
554///     Complex::from(4_f64),
555///     Complex::from(5_f64),
556/// );
557/// assert_eq_matrix!(get_sub_block_r(m1), m2, f64::EPSILON);
558/// ```
559pub fn get_sub_block_r(m: CMatrix3) -> CMatrix2 {
560    CMatrix2::new(
561        m[(0, 0)],
562        m[(0, 1)],
563        // ---
564        m[(1, 0)],
565        m[(1, 1)],
566    )
567}
568
569/// get the [`CMatrix2`] sub block corresponding to [`get_s`]
570///
571/// # Example
572/// ```
573/// # use lattice_qcd_rs::{CMatrix3, CMatrix2, su3::get_sub_block_s, Complex, assert_eq_matrix};
574///  let m1 = CMatrix3::new(
575///     Complex::from(1_f64),
576///     Complex::from(2_f64),
577///     Complex::from(3_f64),
578///     // ---
579///     Complex::from(4_f64),
580///     Complex::from(5_f64),
581///     Complex::from(6_f64),
582///     // ---
583///     Complex::from(7_f64),
584///     Complex::from(8_f64),
585///     Complex::from(9_f64),
586/// );
587/// let m2 = CMatrix2::new(
588///     Complex::from(1_f64),
589///     Complex::from(3_f64),
590///     // ---
591///     Complex::from(7_f64),
592///     Complex::from(9_f64),
593/// );
594/// assert_eq_matrix!(get_sub_block_s(m1), m2, f64::EPSILON);
595/// ```
596pub fn get_sub_block_s(m: CMatrix3) -> CMatrix2 {
597    CMatrix2::new(
598        m[(0, 0)],
599        m[(0, 2)],
600        // ---
601        m[(2, 0)],
602        m[(2, 2)],
603    )
604}
605
606/// Get the [`CMatrix2`] sub block corresponding to [`get_t`]
607///
608/// # Example
609/// ```
610/// # use lattice_qcd_rs::{CMatrix3, CMatrix2, su3::get_sub_block_t, Complex, assert_eq_matrix};
611///  let m1 = CMatrix3::new(
612///     Complex::from(1_f64),
613///     Complex::from(2_f64),
614///     Complex::from(3_f64),
615///     // ---
616///     Complex::from(4_f64),
617///     Complex::from(5_f64),
618///     Complex::from(6_f64),
619///     // ---
620///     Complex::from(7_f64),
621///     Complex::from(8_f64),
622///     Complex::from(9_f64),
623/// );
624/// let m2 = CMatrix2::new(
625///     Complex::from(5_f64),
626///     Complex::from(6_f64),
627///     // ---
628///     Complex::from(8_f64),
629///     Complex::from(9_f64),
630/// );
631/// assert_eq_matrix!(get_sub_block_t(m1), m2, f64::EPSILON);
632/// ```
633pub fn get_sub_block_t(m: CMatrix3) -> CMatrix2 {
634    CMatrix2::new(
635        m[(1, 1)],
636        m[(1, 2)],
637        // ---
638        m[(2, 1)],
639        m[(2, 2)],
640    )
641}
642
643/// Get the unormalize SU(2) sub matrix of an SU(3) matrix corresponding to the "r" sub block see
644/// [`get_sub_block_r`] and [`su2::project_to_su2_unorm`].
645pub fn get_su2_r_unorm(input: CMatrix3) -> CMatrix2 {
646    su2::project_to_su2_unorm(get_sub_block_r(input))
647}
648
649/// Get the unormalize SU(2) sub matrix of an SU(3) matrix corresponding to the "s" sub block see
650/// [`get_sub_block_s`] and [`su2::project_to_su2_unorm`].
651pub fn get_su2_s_unorm(input: CMatrix3) -> CMatrix2 {
652    su2::project_to_su2_unorm(get_sub_block_s(input))
653}
654
655/// Get the unormalize SU(2) sub matrix of an SU(3) matrix corresponding to the "t" sub block see
656/// [`get_sub_block_t`] and [`su2::project_to_su2_unorm`].
657pub fn get_su2_t_unorm(input: CMatrix3) -> CMatrix2 {
658    su2::project_to_su2_unorm(get_sub_block_t(input))
659}
660
661/// Get the three unormalize sub SU(2) matrix of the given SU(3) matrix, ordered `r, s, t`
662/// see [`get_sub_block_r`], [`get_sub_block_s`] and [`get_sub_block_t`]
663pub fn extract_su2_unorm(m: CMatrix3) -> [CMatrix2; 3] {
664    [get_su2_r_unorm(m), get_su2_s_unorm(m), get_su2_t_unorm(m)]
665}
666
667/// Return the matrix
668/// ```textrust
669/// U'_{ij} =| -U_{ij} if i not equal to j
670///          | U_{ii}  if i = j
671/// ```
672/// # Example
673/// ```
674/// # use lattice_qcd_rs::{CMatrix3, su3::reverse, Complex};
675/// assert_eq!(CMatrix3::identity(), reverse(CMatrix3::identity()));
676/// let m1 = CMatrix3::new(
677///     Complex::from(1_f64),
678///     Complex::from(2_f64),
679///     Complex::from(3_f64),
680///     // ---
681///     Complex::from(4_f64),
682///     Complex::from(5_f64),
683///     Complex::from(6_f64),
684///     // ---
685///     Complex::from(7_f64),
686///     Complex::from(8_f64),
687///     Complex::from(9_f64),
688/// );
689/// let m2 = CMatrix3::new(
690///     Complex::from(1_f64),
691///     -Complex::from(2_f64),
692///     -Complex::from(3_f64),
693///     // ---
694///     -Complex::from(4_f64),
695///     Complex::from(5_f64),
696///     -Complex::from(6_f64),
697///     // ---
698///     -Complex::from(7_f64),
699///     -Complex::from(8_f64),
700///     Complex::from(9_f64),
701/// );
702/// assert_eq!(m2, reverse(m1));
703/// assert_eq!(m1, reverse(m2));
704/// ```
705pub fn reverse(input: CMatrix3) -> CMatrix3 {
706    input.map_with_location(|i, j, el| {
707        if i == j {
708            el
709        }
710        else {
711            -el
712        }
713    })
714}
715
716/// Return N such that `1/(N-7)!` < [`f64::EPSILON`].
717///
718/// This number is needed for the computation of exponential matrix
719#[allow(clippy::as_conversions)] // no try into for f64
720pub fn factorial_size_for_exp() -> usize {
721    let mut n: usize = 7;
722    let mut factorial_value = 1;
723    while 1_f64 / (factorial_value as f64) >= Real::EPSILON {
724        n += 1;
725        factorial_value *= n - 7;
726    }
727    n
728}
729
730/// Size of the array for FactorialStorageStatic
731const FACTORIAL_STORAGE_STAT_SIZE: usize = utils::MAX_NUMBER_FACTORIAL + 1;
732
733/// Static store for factorial number.
734#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq, PartialOrd, Ord)]
735struct FactorialStorageStatic {
736    data: [FactorialNumber; FACTORIAL_STORAGE_STAT_SIZE],
737}
738
739macro_rules! set_factorial_storage {
740    ($data:ident, 0) => {
741        $data[0] = 1;
742    };
743    ($data:ident, $e:expr) => {
744        $data[$e] = $data[$e - 1] * $e as FactorialNumber;
745    };
746}
747
748impl FactorialStorageStatic {
749    /// compile time evaluation of all 34 factorial numbers
750    #[allow(clippy::as_conversions)] // constant function cant use try into
751    pub const fn new() -> Self {
752        let mut data: [FactorialNumber; FACTORIAL_STORAGE_STAT_SIZE] =
753            [1; FACTORIAL_STORAGE_STAT_SIZE];
754        let mut i = 1;
755        while i < FACTORIAL_STORAGE_STAT_SIZE {
756            // still not for loop in const fn
757            set_factorial_storage!(data, i);
758            i += 1;
759        }
760        Self { data }
761    }
762
763    /// access in O(1). Return None if `value` is bigger than 34.
764    pub fn try_get_factorial(&self, value: usize) -> Option<&FactorialNumber> {
765        self.data.get(value)
766    }
767
768    /// Get an iterator over the factorial number form `0!` up to `34!`.
769    pub fn iter(
770        &self,
771    ) -> impl Iterator<Item = &FactorialNumber> + ExactSizeIterator + FusedIterator {
772        self.data.iter()
773    }
774
775    /// Get the slice of factorial number.
776    pub const fn as_slice(&self) -> &[FactorialNumber; FACTORIAL_STORAGE_STAT_SIZE] {
777        &self.data
778    }
779}
780
781impl Default for FactorialStorageStatic {
782    fn default() -> Self {
783        Self::new()
784    }
785}
786
787impl std::fmt::Display for FactorialStorageStatic {
788    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
789        write!(f, "factorial Storage : ")?;
790        for (i, n) in self.iter().enumerate() {
791            write!(f, "{}! = {}", i, n)?;
792            if i < self.data.len() - 1 {
793                write!(f, ", ")?;
794            }
795        }
796        Ok(())
797    }
798}
799
800impl<'a> IntoIterator for &'a FactorialStorageStatic {
801    type IntoIter = <&'a [u128; FACTORIAL_STORAGE_STAT_SIZE] as IntoIterator>::IntoIter;
802    type Item = &'a u128;
803
804    fn into_iter(self) -> Self::IntoIter {
805        self.data.iter()
806    }
807}
808
809impl AsRef<[FactorialNumber; FACTORIAL_STORAGE_STAT_SIZE]> for FactorialStorageStatic {
810    fn as_ref(&self) -> &[FactorialNumber; FACTORIAL_STORAGE_STAT_SIZE] {
811        self.as_slice()
812    }
813}
814
815/// factorial number storage in order to find the exponential in O(1) for a set storage
816/// the set if for all number `N` such that `\frac{1}{(N-7)!} >= \mathrm{f64::EPSILON}`
817const FACTORIAL_STORAGE_STAT: FactorialStorageStatic = FactorialStorageStatic::new();
818
819/// number of step for the computation of matrix exponential using the Cayley–Hamilton theorem.
820const N: usize = 26;
821
822/// give the SU3 matrix from the adjoint rep, i.e compute `exp(i v^a T^a )`
823///
824/// The algorithm use is much more efficient the diagonalization method.
825/// It use the Cayley–Hamilton theorem. If you wish to find more about it you can read the
826/// [OpenQCD](https://luscher.web.cern.ch/luscher/openQCD/) documentation that can be found
827/// [here](https://github.com/sa2c/OpenQCD-AVX512/blob/master/doc/su3_fcts.pdf) or by downloading a release.
828/// Note that the documentation above explain the algorithm for exp(X) here it is a modified version for
829/// exp(i X).
830#[inline]
831#[allow(clippy::as_conversions)] // no try into for f64
832pub fn su3_exp_i(su3_adj: Su3Adjoint) -> CMatrix3 {
833    // todo optimize even more using f64 to reduce the number of operation using complex that might be useless
834    const N_LOOP: usize = N - 1;
835    let mut q0: Complex =
836        Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(N_LOOP).unwrap() as f64);
837    let mut q1: Complex = Complex::from(0_f64);
838    let mut q2: Complex = Complex::from(0_f64);
839    let d: Complex = su3_adj.d();
840    let t: Complex = su3_adj.t().into();
841    for i in (0..N_LOOP).rev() {
842        let q0_n =
843            Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(i).unwrap() as f64)
844                + d * q2;
845        let q1_n = I * (q0 - t * q2);
846        let q2_n = I * q1;
847
848        q0 = q0_n;
849        q1 = q1_n;
850        q2 = q2_n;
851    }
852
853    let m = su3_adj.to_matrix();
854    CMatrix3::from_diagonal_element(q0) + m * q1 + m * m * q2
855}
856
857/// the input must be a su(3) matrix (generator of SU(3)), gives the SU3 matrix from the adjoint rep,
858/// i.e compute `exp(i v^a T^a )`
859///
860/// The algorithm use is much more efficient the diagonalization method.
861/// It use the Cayley–Hamilton theorem. If you wish to find more about it you can read the
862/// [OpenQCD](https://luscher.web.cern.ch/luscher/openQCD/) documentation that can be found
863/// [here](https://github.com/sa2c/OpenQCD-AVX512/blob/master/doc/su3_fcts.pdf) or by downloading a release.
864/// Note that the documentation above explain the algorithm for exp(X) here it is a modified version for
865/// exp(i X).
866///
867/// # Panic
868/// The input matrix must be an su(3) (Lie algebra of SU(3)) matrix or approximately su(3),
869/// otherwise the function will panic in debug mod, in release the output gives unexpected values.
870/// ```should_panic
871/// use lattice_qcd_rs::{
872///     assert_eq_matrix,
873///     su3::{matrix_su3_exp_i, MatrixExp},
874/// };
875/// use nalgebra::{Complex, Matrix3};
876/// let i = Complex::new(0_f64, 1_f64);
877/// let matrix = Matrix3::identity(); // this is NOT an su(3) matrix
878/// let output = matrix_su3_exp_i(matrix);
879/// // We panic in debug. In release the following asset will fail.
880/// // assert_eq_matrix!(output, (matrix* i).exp(), f64::EPSILON * 100_000_f64);
881/// ```
882#[inline]
883#[allow(clippy::as_conversions)] // no try into for f64
884pub fn matrix_su3_exp_i(matrix: CMatrix3) -> CMatrix3 {
885    debug_assert!(is_matrix_su3_lie(&matrix, f64::EPSILON * 100_f64));
886    const N_LOOP: usize = N - 1;
887    let mut q0: Complex =
888        Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(N_LOOP).unwrap() as f64);
889    let mut q1: Complex = Complex::from(0_f64);
890    let mut q2: Complex = Complex::from(0_f64);
891    let d: Complex = matrix.determinant() * I;
892    let t: Complex = -Complex::from(0.5_f64) * (matrix * matrix).trace();
893    for i in (0..N_LOOP).rev() {
894        let q0_n =
895            Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(i).unwrap() as f64)
896                + d * q2;
897        let q1_n = I * (q0 - t * q2);
898        let q2_n = I * q1;
899
900        q0 = q0_n;
901        q1 = q1_n;
902        q2 = q2_n;
903    }
904
905    CMatrix3::from_diagonal_element(q0) + matrix * q1 + matrix * matrix * q2
906}
907
908/// gives the value `exp(v^a T^a )`
909///
910/// The algorithm use is much more efficient the diagonalization method.
911/// It use the Cayley–Hamilton theorem. If you wish to find more about it you can read the
912/// [OpenQCD](https://luscher.web.cern.ch/luscher/openQCD/) documentation that can be found
913/// [here](https://github.com/sa2c/OpenQCD-AVX512/blob/master/doc/su3_fcts.pdf) or by downloading a release.
914#[inline]
915#[allow(clippy::as_conversions)] // no try into for f64
916pub fn su3_exp_r(su3_adj: Su3Adjoint) -> CMatrix3 {
917    const N_LOOP: usize = N - 1;
918    let mut q0: Complex =
919        Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(N_LOOP).unwrap() as f64);
920    let mut q1: Complex = Complex::from(0_f64);
921    let mut q2: Complex = Complex::from(0_f64);
922    let d: Complex = su3_adj.d();
923    let t: Complex = su3_adj.t().into();
924    for i in (0..N_LOOP).rev() {
925        let q0_n =
926            Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(i).unwrap() as f64)
927                - I * d * q2;
928        let q1_n = q0 - t * q2;
929        let q2_n = q1;
930
931        q0 = q0_n;
932        q1 = q1_n;
933        q2 = q2_n;
934    }
935
936    let m = su3_adj.to_matrix();
937    CMatrix3::from_diagonal_element(q0) + m * q1 + m * m * q2
938}
939
940/// the input must be a su(3) matrix (generator of SU(3)), gives the value `exp(v^a T^a )`
941///
942/// The algorithm use is much more efficient the diagonalization method.
943/// It use the Cayley–Hamilton theorem. If you wish to find more about it you can read the
944/// [OpenQCD](https://luscher.web.cern.ch/luscher/openQCD/) documentation that can be found
945/// [here](https://github.com/sa2c/OpenQCD-AVX512/blob/master/doc/su3_fcts.pdf) or by downloading a release.
946///
947/// # Panic
948/// The input matrix must be an su(3) (Lie algebra of SU(3)) matrix or approximately su(3),
949/// otherwise the function will panic in debug mod, in release the output gives unexpected values.
950/// ```should_panic
951/// use lattice_qcd_rs::{
952///     assert_eq_matrix,
953///     su3::{matrix_su3_exp_r, MatrixExp},
954/// };
955/// use nalgebra::{Complex, Matrix3};
956/// let i = Complex::new(0_f64, 1_f64);
957/// let matrix = Matrix3::identity(); // this is NOT an su(3)
958/// let output = matrix_su3_exp_r(matrix);
959/// // We panic in debug. In release the following asset will fail.
960/// // assert_eq_matrix!(output, matrix.exp(), f64::EPSILON * 100_000_f64);
961/// ```
962#[inline]
963#[allow(clippy::as_conversions)] // no try into for f64
964pub fn matrix_su3_exp_r(matrix: CMatrix3) -> CMatrix3 {
965    debug_assert!(is_matrix_su3_lie(&matrix, f64::EPSILON * 100_f64));
966    const N_LOOP: usize = N - 1;
967    let mut q0: Complex =
968        Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(N_LOOP).unwrap() as f64);
969    let mut q1: Complex = Complex::from(0_f64);
970    let mut q2: Complex = Complex::from(0_f64);
971    let d: Complex = matrix.determinant() * I;
972    let t: Complex = -Complex::from(0.5_f64) * (matrix * matrix).trace();
973    for i in (0..N_LOOP).rev() {
974        let q0_n =
975            Complex::from(1_f64 / *FACTORIAL_STORAGE_STAT.try_get_factorial(i).unwrap() as f64)
976                - I * d * q2;
977        let q1_n = q0 - t * q2;
978        let q2_n = q1;
979
980        q0 = q0_n;
981        q1 = q1_n;
982        q2 = q2_n;
983    }
984
985    CMatrix3::from_diagonal_element(q0) + matrix * q1 + matrix * matrix * q2
986}
987
988/// Return wether the input matrix is SU(3) up to epsilon.
989// TODO example
990pub fn is_matrix_su3(m: &CMatrix3, epsilon: f64) -> bool {
991    ((m.determinant() - Complex::from(1_f64)).modulus_squared() < epsilon)
992        && ((m * m.adjoint() - CMatrix3::identity()).norm() < epsilon)
993}
994
995/// Returns wether the given matrix is in the lie algebra su(3) that generates SU(3) up to epsilon.
996pub fn is_matrix_su3_lie(matrix: &CMatrix3, epsilon: Real) -> bool {
997    matrix.trace().modulus() < epsilon && (matrix - matrix.adjoint()).norm() < epsilon
998}
999
1000#[doc(hidden)]
1001/// Crate a random 3x3 Matrix
1002pub fn random_matrix_3<R: Rng + ?Sized>(rng: &mut R) -> CMatrix3 {
1003    let d = Uniform::from(-10_f64..10_f64);
1004    CMatrix3::from_fn(|_, _| Complex::new(d.sample(rng), d.sample(rng)))
1005}
1006
1007#[cfg(test)]
1008mod test {
1009    use rand::SeedableRng;
1010
1011    use super::*;
1012    use crate::{assert_eq_matrix, I};
1013
1014    const EPSILON: f64 = 0.000_000_001_f64;
1015    const SEED_RNG: u64 = 0x45_78_93_f4_4a_b0_67_f0;
1016
1017    #[test]
1018    /// test that [`N`] is indeed what we need
1019    fn test_constant() {
1020        assert_eq!(N, factorial_size_for_exp() + 1);
1021    }
1022
1023    #[test]
1024    fn test_factorial() {
1025        for i in 0..FACTORIAL_STORAGE_STAT_SIZE {
1026            assert_eq!(
1027                *FACTORIAL_STORAGE_STAT.try_get_factorial(i).unwrap(),
1028                utils::factorial(i)
1029            );
1030        }
1031    }
1032
1033    #[test]
1034    fn sub_block() {
1035        let m = CMatrix3::new(
1036            Complex::from(1_f64),
1037            Complex::from(2_f64),
1038            Complex::from(3_f64),
1039            // ---
1040            Complex::from(4_f64),
1041            Complex::from(5_f64),
1042            Complex::from(6_f64),
1043            // ---
1044            Complex::from(7_f64),
1045            Complex::from(8_f64),
1046            Complex::from(9_f64),
1047        );
1048        let r = CMatrix2::new(
1049            Complex::from(1_f64),
1050            Complex::from(2_f64),
1051            // ---
1052            Complex::from(4_f64),
1053            Complex::from(5_f64),
1054        );
1055        assert_eq_matrix!(get_sub_block_r(m), r, EPSILON);
1056
1057        let m_r = CMatrix3::new(
1058            Complex::from(1_f64),
1059            Complex::from(2_f64),
1060            Complex::from(0_f64),
1061            // ---
1062            Complex::from(4_f64),
1063            Complex::from(5_f64),
1064            Complex::from(0_f64),
1065            // ---
1066            Complex::from(0_f64),
1067            Complex::from(0_f64),
1068            Complex::from(1_f64),
1069        );
1070        assert_eq_matrix!(get_r(r), m_r, EPSILON);
1071
1072        let s = CMatrix2::new(
1073            Complex::from(1_f64),
1074            Complex::from(3_f64),
1075            // ---
1076            Complex::from(7_f64),
1077            Complex::from(9_f64),
1078        );
1079        assert_eq_matrix!(get_sub_block_s(m), s, EPSILON);
1080
1081        let m_s = CMatrix3::new(
1082            Complex::from(1_f64),
1083            Complex::from(0_f64),
1084            Complex::from(3_f64),
1085            // ---
1086            Complex::from(0_f64),
1087            Complex::from(1_f64),
1088            Complex::from(0_f64),
1089            // ---
1090            Complex::from(7_f64),
1091            Complex::from(0_f64),
1092            Complex::from(9_f64),
1093        );
1094        assert_eq_matrix!(get_s(s), m_s, EPSILON);
1095
1096        let t = CMatrix2::new(
1097            Complex::from(5_f64),
1098            Complex::from(6_f64),
1099            // ---
1100            Complex::from(8_f64),
1101            Complex::from(9_f64),
1102        );
1103        assert_eq_matrix!(get_sub_block_t(m), t, EPSILON);
1104
1105        let m_t = CMatrix3::new(
1106            Complex::from(1_f64),
1107            Complex::from(0_f64),
1108            Complex::from(0_f64),
1109            // ---
1110            Complex::from(0_f64),
1111            Complex::from(5_f64),
1112            Complex::from(6_f64),
1113            // ---
1114            Complex::from(0_f64),
1115            Complex::from(8_f64),
1116            Complex::from(9_f64),
1117        );
1118        assert_eq_matrix!(get_t(t), m_t, EPSILON);
1119
1120        for p in &extract_su2_unorm(m) {
1121            assert_matrix_is_su_2!(p / p.determinant().sqrt(), EPSILON);
1122        }
1123    }
1124
1125    #[test]
1126    #[allow(deprecated)]
1127    fn exp_matrix() {
1128        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
1129        let d = rand::distributions::Uniform::from(-1_f64..1_f64);
1130        for m in &GENERATORS {
1131            let output_r = matrix_su3_exp_r(**m);
1132            assert_eq_matrix!(output_r, m.exp(), EPSILON);
1133            let output_i = matrix_su3_exp_i(**m);
1134            assert_eq_matrix!(output_i, (*m * I).exp(), EPSILON);
1135        }
1136        for _ in 0_u32..100_u32 {
1137            let v = Su3Adjoint::random(&mut rng, &d);
1138            let m = v.to_matrix();
1139            let output_r = matrix_su3_exp_r(m);
1140            assert_eq_matrix!(output_r, m.exp(), EPSILON);
1141            let output_i = matrix_su3_exp_i(m);
1142            assert_eq_matrix!(output_i, (m * I).exp(), EPSILON);
1143        }
1144    }
1145    #[test]
1146    fn test_is_matrix_su3_lie() {
1147        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
1148        let d = rand::distributions::Uniform::from(-1_f64..1_f64);
1149        for m in &GENERATORS {
1150            assert!(is_matrix_su3_lie(*m, f64::EPSILON * 100_f64));
1151        }
1152        for _ in 0_u32..100_u32 {
1153            let v = Su3Adjoint::random(&mut rng, &d);
1154            let m = v.to_matrix();
1155            assert!(is_matrix_su3_lie(&m, f64::EPSILON * 100_f64));
1156        }
1157    }
1158
1159    #[test]
1160    fn test_gen() {
1161        assert_eq!(
1162            CMatrix3::new(
1163                Complex::new(1_f64, 0_f64),
1164                ZERO,
1165                ZERO,
1166                // ---
1167                ZERO,
1168                Complex::new(1_f64, 0_f64),
1169                ZERO,
1170                // ---
1171                ZERO,
1172                ZERO,
1173                Complex::new(-2_f64, 0_f64),
1174            ) * Complex::from(0.5_f64 / 3_f64.sqrt()),
1175            GENERATOR_8
1176        );
1177    }
1178
1179    #[test]
1180    fn test_is_matrix_su3() {
1181        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
1182        let m = CMatrix3::new(
1183            Complex::from(0_f64),
1184            Complex::from(0_f64),
1185            Complex::from(0_f64),
1186            Complex::from(0_f64),
1187            Complex::from(0_f64),
1188            Complex::from(0_f64),
1189            Complex::from(0_f64),
1190            Complex::from(0_f64),
1191            Complex::from(0_f64),
1192        );
1193        assert!(!is_matrix_su3(&m, f64::EPSILON * 100_f64));
1194        let m = CMatrix3::new(
1195            Complex::from(1_f64),
1196            Complex::from(0_f64),
1197            Complex::from(0_f64),
1198            Complex::from(0_f64),
1199            Complex::from(0_f64),
1200            Complex::from(0_f64),
1201            Complex::from(0_f64),
1202            Complex::from(0_f64),
1203            Complex::from(0_f64),
1204        );
1205        assert!(!is_matrix_su3(&m, f64::EPSILON * 100_f64),);
1206        let m = CMatrix3::new(
1207            Complex::from(1_f64),
1208            Complex::from(0_f64),
1209            Complex::from(1_f64),
1210            Complex::from(0_f64),
1211            Complex::from(1_f64),
1212            Complex::from(0_f64),
1213            Complex::from(1_f64),
1214            Complex::from(0_f64),
1215            Complex::from(1_f64),
1216        );
1217        assert!(!is_matrix_su3(&m, f64::EPSILON * 100_f64));
1218        for _ in 0_u32..10_u32 {
1219            let m = random_su3(&mut rng);
1220            assert!(is_matrix_su3(&m, f64::EPSILON * 100_f64));
1221        }
1222        for _ in 0_u32..10_u32 {
1223            let m = random_su3(&mut rng) * Complex::from(1.1_f64);
1224            assert!(!is_matrix_su3(&m, f64::EPSILON * 100_f64));
1225        }
1226    }
1227
1228    #[test]
1229    fn factorial_storage() {
1230        let storage = FactorialStorageStatic::new();
1231        for (i, f) in storage.data.iter().enumerate() {
1232            assert_eq!(utils::factorial(i), *f);
1233        }
1234        assert!(storage.try_get_factorial(35).is_none());
1235        assert_eq!(storage.try_get_factorial(34), Some(&utils::factorial(34)));
1236    }
1237}
lattice_qcd_rs/su3.rs

lattice_qcd_rs/
su3.rs
