qsym2/chartab/
chartab_group.rs

1//! Traits enabling groups to construct and possess character tables.
2
3use std::fmt;
4use std::hash::Hash;
5use std::ops::Mul;
6use std::str::FromStr;
7
8use approx;
9use log;
10use ndarray::{array, s, Array1, Array2, Zip};
11use num::integer::lcm;
12use num_modular::{ModularInteger, MontgomeryInt};
13use num_ord::NumOrd;
14use num_traits::{Inv, One, Pow, ToPrimitive, Zero};
15use primes::is_prime;
16use rayon::prelude::*;
17use serde::{de::DeserializeOwned, Serialize};
18
19use crate::chartab::character::Character;
20use crate::chartab::chartab_symbols::{
21    CollectionSymbol, LinearSpaceSymbol, ReducibleLinearSpaceSymbol,
22};
23use crate::chartab::modular_linalg::{
24    modular_eig, split_2d_space, split_space, weighted_hermitian_inprod,
25};
26use crate::chartab::reducedint::{IntoLinAlgReducedInt, LinAlgMontgomeryInt};
27use crate::chartab::unityroot::UnityRoot;
28use crate::chartab::{CharacterTable, CorepCharacterTable, RepCharacterTable};
29use crate::group::class::ClassProperties;
30use crate::group::{
31    FiniteOrder, GroupProperties, HasUnitarySubgroup, MagneticRepresentedGroup,
32    UnitaryRepresentedGroup,
33};
34
35// =================
36// Trait definitions
37// =================
38
39/// Trait to indicate the presence of character properties in a group and enable access to the
40/// character table of the group.
41pub trait CharacterProperties: ClassProperties
42where
43    Self::RowSymbol: LinearSpaceSymbol,
44    Self::CharTab: CharacterTable<RowSymbol = Self::RowSymbol, ColSymbol = Self::ClassSymbol>,
45{
46    /// Type of the row-labelling symbols in the associated character table.
47    type RowSymbol;
48
49    /// Type of the associated character table whose row-labelling symbol type is constrained to be
50    /// the same as [`Self::RowSymbol`].
51    type CharTab;
52
53    /// Returns a shared reference to the character table of this group.
54    fn character_table(&self) -> &Self::CharTab;
55
56    /// Returns a boolean indicating if this character table contains irreducible representations
57    /// of a unitary-represented group.
58    ///
59    /// If `false`, then some elements in the group are not unitary-represented and one has
60    /// corepresentations instead of representations.
61    fn unitary_represented(&self) -> bool;
62}
63
64/// Trait for the ability to construct an irrep character table for the group.
65///
66/// This trait comes with a default implementation of character table calculation based on the
67/// Burnside--Dixon algorithm with Schneider optimisation.
68pub trait IrrepCharTabConstruction:
69    CharacterProperties<
70    CharTab = RepCharacterTable<
71        <Self as CharacterProperties>::RowSymbol,
72        <Self as ClassProperties>::ClassSymbol,
73    >,
74>
75where
76    Self: Sync,
77    Self::GroupElement: Mul<Output = Self::GroupElement>
78        + Hash
79        + Eq
80        + Clone
81        + Sync
82        + fmt::Debug
83        + FiniteOrder<Int = u32>
84        + Pow<i32, Output = Self::GroupElement>,
85{
86    /// Sets the irrep character table internally.
87    fn set_irrep_character_table(&mut self, chartab: Self::CharTab);
88
89    /// Constructs the irrep character table for this group using the Burnside--Dixon algorithm
90    /// with Schneider optimisation.
91    ///
92    /// # References
93    ///
94    /// * Dixon, J. D. High speed computation of group characters. *Numerische Mathematik* **10**, 446–450 (1967).
95    /// * Schneider, G. J. A. Dixon’s character table algorithm revisited. *Journal of Symbolic Computation* **9**, 601–606 (1990).
96    /// * Grove, L. C. Groups and Characters. (John Wiley & Sons, Inc., 1997).
97    ///
98    /// # Panics
99    ///
100    /// Panics if the Frobenius--Schur indicator for any resulted irrep takes on unexpected values
101    /// and thus implies that the computed irrep is invalid.
102    #[allow(clippy::too_many_lines)]
103    fn construct_irrep_character_table(&mut self) {
104        // Variable definitions
105        // --------------------
106        // m: LCM of the orders of the elements in the group (i.e. the group
107        //    exponent)
108        // p: A prime greater than 2*sqrt(|G|) and m | (p - 1), which is
109        //    guaranteed to exist by Dirichlet's theorem. p is guaranteed to be
110        //    odd.
111        // z: An integer having multiplicative order m when viewed as an
112        //    element of Z*p, i.e. z^m ≡ 1 (mod p) but z^n !≡ 1 (mod p) for all
113        //    0 <= n < m.
114
115        log::debug!("=============================================");
116        log::debug!("Construction of irrep character table begins.");
117        log::debug!("      *** Burnside--Dixon algorithm ***      ");
118        log::debug!("      ** with Schneider optimisation **      ");
119        log::debug!("=============================================");
120
121        // Identify a suitable finite field
122        let m = (0..self.class_number())
123            .map(|i| {
124                self.get_cc_transversal(i)
125                    .unwrap_or_else(|| panic!("No representative of class index `{i}` found."))
126                    .order()
127            })
128            .reduce(lcm)
129            .expect("Unable to find the LCM for the orders of the elements in this group.");
130        let zeta = UnityRoot::new(1, m);
131        log::debug!("Found group exponent m = {}.", m);
132        log::debug!("Chosen primitive unity root ζ = {}.", zeta);
133
134        let rf64 = (2.0
135            * self
136                .order()
137                .to_f64()
138                .unwrap_or_else(|| panic!("Unable to convert `{}` to `f64`.", self.order()))
139                .sqrt()
140            / (f64::from(m)))
141        .ceil();
142        assert!(rf64.is_sign_positive());
143        assert!(rf64 <= f64::from(u32::MAX));
144        #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
145        let mut r = rf64 as u32;
146        if r == 0 {
147            r = 1;
148        };
149        let mut p = r * m + 1;
150        while !is_prime(u64::from(p)) {
151            log::debug!("Trying {}: not prime.", p);
152            r += 1;
153            p = r * m + 1;
154        }
155        log::debug!("Found r = {}.", r);
156        log::debug!("Found prime number p = r * m + 1 = {}.", p);
157        log::debug!("All arithmetic will now be carried out in GF({}).", p);
158
159        let modp = MontgomeryInt::<u32>::new(1, &p).linalg();
160        // p is prime, so there is guaranteed a z < p such that z^m ≡ 1 (mod p).
161        let mut i = 1u32;
162        while modp.convert(i).multiplicative_order().unwrap_or_else(|| {
163            panic!(
164                "Unable to find multiplicative order for `{}`",
165                modp.convert(i)
166            )
167        }) != m
168            && i < p
169        {
170            i += 1;
171        }
172        let z = modp.convert(i);
173        assert_eq!(
174            z.multiplicative_order()
175                .unwrap_or_else(|| panic!("Unable to find multiplicative order for `{z}`.")),
176            m
177        );
178        log::debug!("Found integer z = {} with multiplicative order {}.", z, m);
179
180        // Diagonalise class matrices
181        let class_sizes: Vec<_> = (0..self.class_number())
182            .map(|i| {
183                self.class_size(i)
184                    .expect("Not all class sizes can be found.")
185            })
186            .collect();
187        log::debug!("Class sizes: {:?}", class_sizes);
188        let sq_indices: Vec<usize> = (0..self.class_number())
189            .map(|i| {
190                let ele = self
191                    .get_cc_transversal(i)
192                    .unwrap_or_else(|| panic!("No representative of class index `{i}` found."));
193                let elep2 = ele.pow(2);
194                let elep2_i = self
195                    .get_index_of(&elep2)
196                    .unwrap_or_else(|| panic!("Element {elep2:?} not found."));
197                self.get_cc_of_element_index(elep2_i)
198                    .unwrap_or_else(|| panic!("Conjugacy class for element {elep2:?} not found."))
199            })
200            .collect();
201        let inverse_conjugacy_classes = Some(
202            (0..self.class_number())
203                .map(|i| {
204                    self.get_inverse_cc(i)
205                        .expect("Not all class inverses could be obtained.")
206                })
207                .collect::<Vec<_>>(),
208        );
209        log::debug!(
210            "Inverse conjugacy classes: {:?}",
211            inverse_conjugacy_classes.as_ref().unwrap()
212        );
213
214        let mut eigvecs_1d: Vec<Array1<LinAlgMontgomeryInt<u32>>> = vec![];
215        let mut success = false;
216        let mut rotate_count = 0;
217        let rs_original = (1..self.class_number()).collect::<Vec<_>>();
218
219        while !success && rotate_count < self.class_number() {
220            let mut rs = rs_original.clone();
221            rs.rotate_left(rotate_count);
222            log::debug!("Class matrix consideration order: {rs:?}");
223            rotate_count += 1;
224            let mut r_iter = rs.into_iter();
225            eigvecs_1d = vec![];
226
227            if self.class_number() == 1 {
228                eigvecs_1d.push(array![modp.convert(1)]);
229                success = true;
230            } else {
231                let mut degenerate_subspaces: Vec<Vec<Array1<LinAlgMontgomeryInt<u32>>>> = vec![];
232                let ctb = self.cayley_table();
233                let r = r_iter
234                    .next()
235                    .expect("Unable to obtain any `r` indices for class matrices.");
236                log::debug!("Considering class matrix N{r}...");
237                let nmat_1 = self.class_matrix(ctb, r).map(|&i| {
238                    modp.convert(
239                        u32::try_from(i)
240                            .unwrap_or_else(|_| panic!("Unable to convert `{i}` to `u32`.")),
241                    )
242                });
243                let eigs_1 = modular_eig(&nmat_1).unwrap_or_else(|err| {
244                    log::error!("{err}");
245                    panic!("{err}");
246                });
247                eigvecs_1d.extend(eigs_1.iter().filter_map(|(_, eigvecs)| {
248                    if eigvecs.len() == 1 {
249                        Some(eigvecs[0].clone())
250                    } else {
251                        None
252                    }
253                }));
254                degenerate_subspaces.extend(
255                    eigs_1
256                        .iter()
257                        .filter_map(|(_, eigvecs)| {
258                            if eigvecs.len() > 1 {
259                                Some(eigvecs)
260                            } else {
261                                None
262                            }
263                        })
264                        .cloned(),
265                );
266
267                while !degenerate_subspaces.is_empty() {
268                    log::debug!(
269                        "Number of 1-D eigenvectors found: {} / {}.",
270                        eigvecs_1d.len(),
271                        self.class_number()
272                    );
273                    log::debug!(
274                        "Number of degenerate subspaces found: {}.",
275                        degenerate_subspaces.len(),
276                    );
277
278                    let mut degenerate_2d_subspaces = degenerate_subspaces
279                        .iter()
280                        .filter(|subspace| subspace.len() == 2)
281                        .cloned()
282                        .collect::<Vec<_>>();
283                    if !degenerate_2d_subspaces.is_empty() {
284                        degenerate_subspaces.retain(|subspace| subspace.len() > 2);
285                        log::debug!(
286                            "Number of 2-D degenerate subspaces found: {}",
287                            degenerate_2d_subspaces.len()
288                        );
289                        log::debug!(
290                            "Schneider's greedy algorithm for splitting 2-D spaces will be attempted."
291                        );
292                        while let Some(subspace) = degenerate_2d_subspaces.pop() {
293                            if let Ok(subsubspaces) = split_2d_space(
294                                &subspace,
295                                &class_sizes,
296                                &sq_indices,
297                                inverse_conjugacy_classes.as_ref(),
298                            ) {
299                                eigvecs_1d.extend(subsubspaces.iter().filter_map(|subsubspace| {
300                                    if subsubspace.len() == 1 {
301                                        Some(subsubspace[0].clone())
302                                    } else {
303                                        None
304                                    }
305                                }));
306                                log::debug!(
307                                    "2-D subspace index {} successfully split.",
308                                    degenerate_2d_subspaces.len()
309                                );
310                            } else {
311                                log::debug!(
312                                    "2-D subspace index {} cannot be split greedily.",
313                                    degenerate_2d_subspaces.len()
314                                );
315                                log::debug!(
316                                    "Stashing this 2-D subspace for splitting with class matrices..."
317                                );
318                                degenerate_subspaces.push(subspace);
319                            }
320                        }
321                    }
322
323                    if !degenerate_subspaces.is_empty() {
324                        if let Some(r) = r_iter.next() {
325                            log::debug!("Considering class matrix N{}...", r);
326                            let nmat_r = self.class_matrix(ctb, r).map(|&i| {
327                                modp.convert(u32::try_from(i).unwrap_or_else(|_| {
328                                    panic!("Unable to convert `{i}` to `u32`.")
329                                }))
330                            });
331
332                            let mut remaining_degenerate_subspaces: Vec<
333                                Vec<Array1<LinAlgMontgomeryInt<u32>>>,
334                            > = vec![];
335                            while !degenerate_subspaces.is_empty() {
336                                let subspace = degenerate_subspaces
337                                    .pop()
338                                    .expect("Unexpected empty `degenerate_subspaces`.");
339
340                                if let Ok(subsubspaces) = split_space(
341                                    &nmat_r,
342                                    &subspace,
343                                    &class_sizes,
344                                    inverse_conjugacy_classes.as_ref(),
345                                ) {
346                                    eigvecs_1d.extend(subsubspaces.iter().filter_map(
347                                        |subsubspace| {
348                                            if subsubspace.len() == 1 {
349                                                Some(subsubspace[0].clone())
350                                            } else {
351                                                None
352                                            }
353                                        },
354                                    ));
355                                    remaining_degenerate_subspaces.extend(
356                                        subsubspaces
357                                            .iter()
358                                            .filter(|subsubspace| subsubspace.len() > 1)
359                                            .cloned(),
360                                    );
361                                } else {
362                                    log::warn!(
363                                        "Class matrix N{} failed to split degenerate subspace {}.",
364                                        r,
365                                        degenerate_subspaces.len()
366                                    );
367                                    log::warn!(
368                                        "Stashing this subspace for the next class matrices..."
369                                    );
370                                    remaining_degenerate_subspaces.push(subspace);
371                                }
372                            }
373                            degenerate_subspaces = remaining_degenerate_subspaces;
374                        } else {
375                            log::warn!(
376                                "Class matrices exhausted before degenerate subspaces are fully resolved."
377                            );
378                            log::warn!(
379                                "Restarting the solver using a different class matrix order..."
380                            );
381                            break;
382                        }
383                    }
384                }
385            }
386            if eigvecs_1d.len() == self.class_number() {
387                success = true;
388            }
389        }
390
391        if eigvecs_1d.len() != self.class_number() {
392            log::error!("Degenerate subspaces failed to be fully resolved for all cyclic permutations of class matrices N1, ..., N{}.", self.class_number());
393        }
394        assert_eq!(
395            eigvecs_1d.len(),
396            self.class_number(),
397            "Degenerate subspaces failed to be fully resolved."
398        );
399        log::debug!(
400            "Successfully found {} / {} one-dimensional eigenvectors for the class matrices.",
401            eigvecs_1d.len(),
402            self.class_number()
403        );
404        for (i, vec) in eigvecs_1d.iter().enumerate() {
405            log::debug!("Eigenvector {}: {}", i, vec);
406        }
407
408        // Lift characters back to the complex field
409        log::debug!("Lifting characters from GF({p}) back to the complex field...",);
410
411        let chars: Vec<_> = eigvecs_1d
412            .par_iter()
413            .flat_map(|vec_i| {
414                let vec_i_inprod = weighted_hermitian_inprod(
415                    (vec_i, vec_i),
416                    &class_sizes,
417                    inverse_conjugacy_classes.as_ref(),
418                );
419                let dim_i = (1..=p.div_euclid(2))
420                    .map(|d| vec_i_inprod.convert(d))
421                    .find(|d_modp| {
422                        Some(vec_i_inprod) == d_modp.square().inv()
423                    })
424                    .unwrap_or_else(|| {
425                        log::error!("Unable to deduce the irrep dimensionality from ⟨θvi, θvi⟩ = {vec_i_inprod} where vi = {vec_i}.");
426                        panic!("Unable to deduce the irrep dimensionality from ⟨θvi, θvi⟩ = {vec_i_inprod} where vi = {vec_i}.");
427                    });
428                log::debug!("⟨θvi, θvi⟩ = {vec_i_inprod} yields irrep dimensionality {}.", dim_i.residue());
429
430                let tchar_i =
431                    Zip::from(vec_i)
432                        .and(class_sizes.as_slice())
433                        .par_map_collect(|&v, &k| {
434                            v * dim_i
435                                / modp.convert(u32::try_from(k).unwrap_or_else(|_| {
436                                    panic!("Unable to convert `{k}` to `u32`.")
437                                }))
438                        });
439                let char_i: Vec<_> = (0..self.class_number())
440                    .into_par_iter()
441                    .map(|cc_i| {
442                        let x = self
443                            .get_cc_transversal(cc_i)
444                            .unwrap_or_else(|| {
445                                panic!("Representative of class index `{cc_i}` not found.")
446                            });
447                        let k = x.order();
448                        let xi = zeta.clone().pow(
449                            i32::try_from(m.checked_div_euclid(k).unwrap_or_else(|| {
450                                panic!("`{m}` cannot be Euclid-divided by `{k}`.")
451                            }))
452                            .unwrap_or_else(|_| {
453                                panic!(
454                                    "The Euclid division `{m} / {k}` cannot be converted to `i32`."
455                                )
456                            }),
457                        );
458                        let char_ij_terms: Vec<_> = (0..k)
459                            .into_par_iter()
460                            .map(|s| {
461                                let mu_s = (0..k).fold(modp.convert(0), |acc, l| {
462                                    let x_l =
463                                        x.clone().pow(i32::try_from(l).unwrap_or_else(|_| {
464                                            panic!("Unable to convert `{l}` to `i32`.")
465                                        }));
466                                    let x_l_idx = self
467                                        .get_index_of(&x_l)
468                                        .unwrap_or_else(|| panic!("Element {x_l:?} not found."));
469                                    let x_l_class_idx =
470                                        self.get_cc_of_element_index(x_l_idx).unwrap_or_else(|| {
471                                            panic!("Element `{x_l_idx}` does not have a conjugacy class.")
472                                        });
473                                    let tchar_i_x_l = tchar_i[x_l_class_idx];
474                                    acc + tchar_i_x_l
475                                        * z.pow(
476                                            s * l
477                                                * m.checked_div_euclid(k).unwrap_or_else(|| {
478                                                    panic!(
479                                                        "`{m}` cannot be Euclid-divided by `{k}`."
480                                                    )
481                                                }),
482                                        )
483                                        .inv()
484                                        .expect("Unable to find the modular inverse of z^(slm//k).")
485                                }) / k;
486                                (
487                                    xi.pow(i32::try_from(s).unwrap_or_else(|_| {
488                                        panic!("Unable to convert `{s}` to `i32`.")
489                                    })),
490                                    usize::try_from(mu_s.residue()).unwrap_or_else(|_| {
491                                        panic!("Unable to convert `{}` to `usize`.", mu_s.residue())
492                                    }),
493                                )
494                            })
495                            .collect();
496                        // We do not wish to simplify the character here, even if it can be
497                        // simplified (e.g. (E8)^2 + (E8)^6 can be simplified to 0). This is so
498                        // that the full unity-root-term-structure can be used in the ordering of
499                        // irreps.
500                        Character::new(&char_ij_terms)
501                    })
502                    .collect();
503                char_i
504            })
505            .collect();
506
507        let char_arr = Array2::from_shape_vec((self.class_number(), self.class_number()), chars)
508            .expect("Unable to construct the two-dimensional array of characters.");
509        log::debug!("Lifting characters from GF({p}) back to the complex field... Done.",);
510
511        let default_irrep_symbols = char_arr
512            .rows()
513            .into_iter()
514            .enumerate()
515            .map(|(irrep_i, _)| {
516                Self::RowSymbol::from_str(&format!("||Λ|_({irrep_i})|"))
517                    .ok()
518                    .expect("Unable to construct default irrep symbols.")
519            })
520            .collect::<Vec<_>>();
521        let default_principal_classes = vec![self
522            .get_cc_symbol_of_index(0)
523            .expect("No conjugacy class symbols found.")];
524
525        log::debug!("Computing the Frobenius--Schur indicators in GF({p})...");
526        let group_order = class_sizes.iter().sum::<usize>();
527        let group_order_u32 = u32::try_from(group_order).unwrap_or_else(|_| {
528            panic!("Unable to convert the group order {group_order} to `u32`.")
529        });
530        log::debug!("Indices of squared conjugacy classes: {sq_indices:?}");
531        let frobenius_schur_indicators: Vec<i8> = eigvecs_1d
532            .par_iter()
533            .map(|vec_i| {
534                let vec_i_inprod = weighted_hermitian_inprod(
535                    (vec_i, vec_i),
536                    &class_sizes,
537                    inverse_conjugacy_classes.as_ref(),
538                );
539                let dim_i = (1..=p.div_euclid(2))
540                    .map(|d| vec_i_inprod.convert(d))
541                    .find(|d_modp| {
542                        Some(vec_i_inprod) == d_modp.square().inv()
543                    })
544                    .unwrap_or_else(|| {
545                        log::error!("Unable to deduce the irrep dimensionality from ⟨θvi, θvi⟩ = {vec_i_inprod} where vi = {vec_i}.");
546                        panic!("Unable to deduce the irrep dimensionality from ⟨θvi, θvi⟩ = {vec_i_inprod} where vi = {vec_i}.");
547                    });
548
549                let tchar_i =
550                    Zip::from(vec_i)
551                        .and(class_sizes.as_slice())
552                        .par_map_collect(|&v, &k| {
553                            v * dim_i
554                                / modp.convert(u32::try_from(k).unwrap_or_else(|_| {
555                                    panic!("Unable to convert `{k}` to `u32`.")
556                                }))
557                        });
558                let fs_i = sq_indices
559                    .iter()
560                    .zip(class_sizes.iter())
561                    .fold(modp.convert(0), |acc, (&sq_idx, &k)| {
562                        let k_u32 = u32::try_from(k).unwrap_or_else(|_| {
563                            panic!("Unable to convert the class size {k} to `u32`.");
564                        });
565                        acc + modp.convert(k_u32) * tchar_i[sq_idx]
566                    }) / modp.convert(group_order_u32);
567
568                if fs_i.is_one() {
569                    1i8
570                } else if Zero::is_zero(&fs_i) {
571                    0i8
572                } else if fs_i == modp.convert(modp.modulus() - 1) {
573                    -1i8
574                } else {
575                    log::error!("Invalid Frobenius -- Schur indicator: `{fs_i}`.");
576                    panic!("Invalid Frobenius -- Schur indicator: `{fs_i}`.");
577                }
578            }).collect();
579        log::debug!("Computing the Frobenius--Schur indicators in GF({p})... Done.");
580
581        let chartab_name = if let Some(finite_name) = self.finite_subgroup_name().as_ref() {
582            format!("{} > {finite_name}", self.name())
583        } else {
584            self.name()
585        };
586        let ccsyms = (0..self.class_number())
587            .map(|i| {
588                self.get_cc_symbol_of_index(i)
589                    .unwrap_or_else(|| {
590                        let rep = self
591                            .get_cc_transversal(i)
592                            .unwrap_or_else(||
593                                panic!("Unable to obtain a representative for conjugacy class `{i}`.")
594                            );
595                        panic!("Class symbol for conjugacy class `{i}` with representative element `{rep:?}` cannot be found.")
596                    })
597            })
598            .collect::<Vec<_>>();
599        self.set_irrep_character_table(RepCharacterTable::new(
600            chartab_name.as_str(),
601            &default_irrep_symbols,
602            &ccsyms,
603            &default_principal_classes,
604            char_arr,
605            &frobenius_schur_indicators,
606        ));
607
608        log::debug!("===========================================");
609        log::debug!("     *** Burnside--Dixon algorithm ***     ");
610        log::debug!("     ** with Schneider optimisation **     ");
611        log::debug!("Construction of irrep character table ends.");
612        log::debug!("===========================================");
613    }
614}
615
616/// Trait for the ability to construct an ircorep character table for the group.
617///
618/// This trait comes with a default implementation of ircorep character table calculation based on
619/// the irreps of the unitary subgroup.
620pub trait IrcorepCharTabConstruction: HasUnitarySubgroup + CharacterProperties<
621    CharTab = CorepCharacterTable<
622        <Self as CharacterProperties>::RowSymbol,
623        <<Self as HasUnitarySubgroup>::UnitarySubgroup as CharacterProperties>::CharTab,
624    >
625>
626where
627    Self: ClassProperties<ClassSymbol = <<Self as HasUnitarySubgroup>::UnitarySubgroup as ClassProperties>::ClassSymbol>,
628    Self::RowSymbol: ReducibleLinearSpaceSymbol<
629        Subspace = <
630            <Self as HasUnitarySubgroup>::UnitarySubgroup as CharacterProperties
631        >::RowSymbol
632    >,
633    <<Self as HasUnitarySubgroup>::UnitarySubgroup as ClassProperties>::ClassSymbol: Serialize + DeserializeOwned,
634{
635    /// Sets the ircorep character table internally.
636    fn set_ircorep_character_table(&mut self, chartab: Self::CharTab);
637
638    /// Constructs the ircorep character table for this group.
639    ///
640    /// For each irrep in the unitary subgroup, the type of the ircorep it induces is determined
641    /// using the Dimmock--Wheeler character test, then the ircorep's characters in the
642    /// unitary-represented part of the full group are determined to give a square character table.
643    ///
644    /// # References
645    ///
646    /// * Bradley, C. J. & Davies, B. L. Magnetic Groups and Their Corepresentations. *Rev. Mod. Phys.* **40**, 359–379 (1968).
647    /// * Newmarch, J. D. & Golding, R. M. The character table for the corepresentations of magnetic groups. *J. Math. Phys.* **23**, 695–704 (1982).
648    /// * Newmarch, J. D. Some character theory for groups of linear and antilinear operators. *J. Math. Phys.* **24**, 742–756 (1983).
649    ///
650    /// # Panics
651    ///
652    /// Panics if any calculated ircoreps are found to be invalid.
653    fn construct_ircorep_character_table(&mut self) {
654        log::debug!("===============================================");
655        log::debug!("Construction of ircorep character table begins.");
656        log::debug!("===============================================");
657
658        if self.unitary_subgroup().order() == self.order() {
659            log::debug!(
660                "The unitary subgroup order and the full group order are both {}. This is not a magnetic group.", self.order()
661            );
662            // This is not a magnetic group. There is nothing to do.
663            return;
664        }
665
666        debug_assert_eq!(self.order() % 2, 0);
667        debug_assert_eq!(self.order().div_euclid(2), self.unitary_subgroup().order());
668        let unitary_order: i32 = self
669            .order()
670            .div_euclid(2)
671            .try_into()
672            .expect("Unable to convert the unitary group order to `i32`.");
673        let unitary_chartab = self.unitary_subgroup().character_table();
674
675        let mag_ctb = self.cayley_table().expect("Cayley table not found for the magnetic group.");
676
677        let mag = &self;
678        let uni = self.unitary_subgroup();
679        let mag_ccsyms = (0..mag.class_number()) .map(|i| {
680            mag.get_cc_symbol_of_index(i).expect("Unable to retrieve all class symbols of the magnetic group.")
681        })
682        .collect::<Vec<_>>();
683        let uni_ccsyms = (0..uni.class_number()).map(|i| {
684            uni.get_cc_symbol_of_index(i).expect("Unable to retrieve all class symbols of the unitary subgroup.")
685        }).collect::<Vec<_>>();
686        let (a0_mag_idx, _) = self
687            .elements()
688            .clone()
689            .into_iter()
690            .enumerate()
691            .find(|(_, op)| self.check_elem_antiunitary(op))
692            .expect("No antiunitary elements found in the magnetic group.");
693
694        let mut remaining_irreps = unitary_chartab.get_all_rows();
695        remaining_irreps.reverse();
696
697        let mut ircoreps_ins: Vec<(Self::RowSymbol, u8)> = Vec::new();
698        while let Some(irrep) = remaining_irreps.pop() {
699            log::debug!("Considering irrep {irrep} of the unitary subgroup...");
700            let char_sum = self
701                .elements()
702                .clone()
703                .into_iter()
704                .enumerate()
705                .filter(|(_, op)| self.check_elem_antiunitary(op))
706                .fold(Character::zero(), |acc, (a_mag_idx, _)| {
707                    let a2_mag_idx = mag_ctb[(a_mag_idx, a_mag_idx)];
708                    let a2 = self.get_index(a2_mag_idx).unwrap_or_else(|| {
709                        panic!("Element index `{a2_mag_idx}` not found in the magnetic group.")
710                    });
711                    let a2_uni_idx = self.unitary_subgroup().get_index_of(&a2).unwrap_or_else(|| {
712                        panic!("Element `{a2:?}` not found in the unitary subgroup.")
713                    });
714                    let a2_uni_class = &uni_ccsyms[
715                        uni.get_cc_of_element_index(a2_uni_idx).unwrap_or_else(|| {
716                            panic!("Conjugacy class for `{a2:?}` not found in the unitary subgroup.")
717                        })
718                    ];
719                    acc + unitary_chartab.get_character(&irrep, a2_uni_class)
720                })
721                .simplify();
722            log::debug!("  Dimmock--Wheeler indicator for {irrep}: {char_sum}");
723            let char_sum_c128 = char_sum.complex_value();
724            approx::assert_relative_eq!(
725                char_sum_c128.im,
726                0.0,
727                max_relative = char_sum.threshold()
728                    * unitary_order
729                        .to_f64()
730                        .expect("Unable to convert the unitary order to `f64`.")
731                        .sqrt(),
732                epsilon = char_sum.threshold()
733                    * unitary_order
734                        .to_f64()
735                        .expect("Unable to convert the unitary order to `f64`.")
736                        .sqrt(),
737            );
738            approx::assert_relative_eq!(
739                char_sum_c128.re,
740                char_sum_c128.re.round(),
741                max_relative = char_sum.threshold()
742                    * unitary_order
743                        .to_f64()
744                        .expect("Unable to convert the unitary order to `f64`.")
745                        .sqrt(),
746                epsilon = char_sum.threshold()
747                    * unitary_order
748                        .to_f64()
749                        .expect("Unable to convert the unitary order to `f64`.")
750                        .sqrt(),
751            );
752            let char_sum = char_sum_c128.re.round();
753
754            let (intertwining_number, ircorep) = if NumOrd(char_sum) == NumOrd(unitary_order) {
755                // Irreducible corepresentation type a
756                // Δ(u) is equivalent to Δ*[a^(-1)ua].
757                // Δ(u) is contained once in the induced irreducible corepresentation.
758                log::debug!(
759                    "  Ircorep induced by {irrep} is of type (a) with intertwining number 1."
760                );
761                (1u8, Self::RowSymbol::from_subspaces(&[(irrep, 1)]))
762            } else if NumOrd(char_sum) == NumOrd(-unitary_order) {
763                // Irreducible corepresentation type b
764                // Δ(u) is equivalent to Δ*[a^(-1)ua].
765                // Δ(u) is contained twice in the induced irreducible corepresentation.
766                log::debug!(
767                    "  Ircorep induced by {irrep} is of type (b) with intertwining number 4."
768                );
769                (4u8, Self::RowSymbol::from_subspaces(&[(irrep, 2)]))
770            } else if NumOrd(char_sum) == NumOrd(0i8) {
771                // Irreducible corepresentation type c
772                // Δ(u) is inequivalent to Δ*[a^(-1)ua].
773                // Δ(u) and Δ*[a^(-1)ua] are contained the induced irreducible corepresentation.
774                let irrep_conj_chars: Vec<Character> = unitary_chartab
775                    .get_all_cols()
776                    .iter()
777                    .enumerate()
778                    .map(|(cc_idx, cc)| {
779                        let u_cc = self
780                            .unitary_subgroup()
781                            .get_cc_index(cc_idx)
782                            .unwrap_or_else(|| panic!("Conjugacy class `{cc}` not found."));
783                        let u_unitary_idx = u_cc
784                            .iter()
785                            .next()
786                            .unwrap_or_else(|| panic!("No unitary elements found for conjugacy class `{cc}`."));
787                        let u = self.unitary_subgroup()
788                            .get_index(*u_unitary_idx)
789                            .unwrap_or_else(|| panic!("Unitary element with index `{u_unitary_idx}` cannot be retrieved."));
790                        let u_mag_idx = self
791                            .get_index_of(&u)
792                            .unwrap_or_else(|| panic!("Unable to retrieve the index of unitary element `{u:?}` in the magnetic group."));
793                        let ua0_mag_idx = mag_ctb[(u_mag_idx, a0_mag_idx)];
794                        let mag_ctb_a0x = mag_ctb.slice(s![a0_mag_idx, ..]);
795                        let a0invua0_mag_idx = mag_ctb_a0x.iter().position(|&x| x == ua0_mag_idx).unwrap_or_else(|| {
796                            panic!("No element `{ua0_mag_idx}` can be found in row `{a0_mag_idx}` of the magnetic group Cayley table.")
797                        });
798                        let a0invua0 = self
799                            .get_index(a0invua0_mag_idx)
800                            .unwrap_or_else(|| {
801                                panic!("Unable to retrieve element with index `{a0invua0_mag_idx}` in the magnetic group.")
802                            });
803                        let a0invua0_unitary_idx = self.unitary_subgroup()
804                            .get_index_of(&a0invua0)
805                            .unwrap_or_else(|| {
806                                panic!("Unable to retrieve the index of element `{a0invua0:?}` in the unitary subgroup.")
807                            });
808                        let a0invua0_unitary_class = &uni_ccsyms[
809                            uni.get_cc_of_element_index(a0invua0_unitary_idx).unwrap_or_else(|| {
810                                panic!("Unable to retrieve the class for `{a0invua0:?}` in the unitary subgroup.")
811                            })
812                        ];
813                        unitary_chartab.get_character(&irrep, a0invua0_unitary_class).complex_conjugate()
814                }).collect();
815                let all_irreps = unitary_chartab.get_all_rows();
816                let (_, conj_irrep) = all_irreps
817                    .iter()
818                    .enumerate()
819                    .find(|(irrep_idx, _)| {
820                        unitary_chartab.array().row(*irrep_idx).to_vec() == irrep_conj_chars
821                    })
822                    .unwrap_or_else(|| panic!("Conjugate irrep for {irrep} not found."));
823                assert!(remaining_irreps.shift_remove(conj_irrep));
824
825                log::debug!("  The Wigner-conjugate irrep of {irrep} is {conj_irrep}.");
826                log::debug!("  Ircorep induced by {irrep} and {conj_irrep} is of type (c) with intertwining number 2.");
827                (
828                    2u8,
829                    Self::RowSymbol::from_subspaces(&[(irrep, 1), (conj_irrep.to_owned(), 1)]),
830                )
831            } else {
832                log::error!(
833                    "Unexpected `char_sum`: {char_sum}. This can only be ±{unitary_order} or 0."
834                );
835                panic!("Unexpected `char_sum`: {char_sum}. This can only be ±{unitary_order} or 0.")
836            };
837            ircoreps_ins.push((ircorep, intertwining_number));
838        }
839
840        let mut char_arr: Array2<Character> =
841            Array2::zeros((ircoreps_ins.len(), self.class_number()));
842        for (i, (ircorep, intertwining_number)) in ircoreps_ins.iter().enumerate() {
843            for cc_idx in 0..mag.class_number() {
844                let mag_cc = &mag_ccsyms[cc_idx];
845                let mag_cc_rep = mag_cc.representative().unwrap_or_else(|| {
846                    panic!(
847                        "No representative element found for magnetic conjugacy class {mag_cc}."
848                    );
849                });
850                let mag_cc_uni_idx = self
851                    .unitary_subgroup()
852                    .get_index_of(mag_cc_rep)
853                    .unwrap_or_else(|| {
854                        panic!(
855                            "Index for element {mag_cc_rep:?} not found in the unitary subgroup."
856                        );
857                    });
858                let uni_cc = &uni_ccsyms[
859                    uni.get_cc_of_element_index(mag_cc_uni_idx).unwrap_or_else(|| {
860                        panic!("Unable to find the conjugacy class of element {mag_cc_rep:?} in the unitary subgroup.");
861                    })
862                ];
863
864                char_arr[(i, cc_idx)] = ircorep
865                    .subspaces()
866                    .iter()
867                    .fold(Character::zero(), |acc, (irrep, _)| {
868                        acc + unitary_chartab.get_character(irrep, uni_cc)
869                    });
870                if *intertwining_number == 4 {
871                    // Irreducible corepresentation type b
872                    // The inducing irrep appears twice.
873                    char_arr[(i, cc_idx)] *= 2;
874                }
875            }
876        }
877
878        let principal_classes = (0..mag.class_number())
879            .filter_map(|i| {
880                let mag_cc = &mag_ccsyms[i];
881                let mag_cc_rep = mag_cc.representative().unwrap_or_else(|| {
882                    panic!("No representative element found for magnetic conjugacy class {mag_cc}.");
883                });
884                let mag_cc_uni_idx = self.unitary_subgroup().get_index_of(mag_cc_rep).unwrap_or_else(|| {
885                    panic!("Index for element {mag_cc_rep:?} not found in the unitary subgroup.");
886                });
887                let uni_cc = uni.get_cc_symbol_of_index(
888                    uni.get_cc_of_element_index(mag_cc_uni_idx).unwrap_or_else(|| {
889                        panic!("Unable to find the conjugacy class of element {mag_cc_rep:?} in the unitary subgroup.");
890                    })
891                ).unwrap_or_else(|| {
892                    panic!("Unable to find the conjugacy class symbol of element {mag_cc_rep:?} in the unitary subgroup.");
893                });
894                if unitary_chartab.get_principal_cols().contains(&uni_cc) {
895                    Some(mag_cc.clone())
896                } else {
897                    None
898                }
899            }).collect::<Vec<_>>();
900
901        let (ircoreps, ins): (Vec<Self::RowSymbol>, Vec<u8>) = ircoreps_ins.into_iter().unzip();
902
903        let chartab_name = if let Some(finite_name) = self.finite_subgroup_name().as_ref() {
904            format!("{} > {finite_name}", self.name())
905        } else {
906            self.name()
907        };
908        self.set_ircorep_character_table(Self::CharTab::new(
909            chartab_name.as_str(),
910            unitary_chartab.clone(),
911            &ircoreps,
912            &mag_ccsyms,
913            &principal_classes,
914            char_arr,
915            &ins,
916        ));
917
918        log::debug!("=============================================");
919        log::debug!("Construction of ircorep character table ends.");
920        log::debug!("=============================================");
921    }
922}
923
924// =====================
925// Trait implementations
926// =====================
927
928// ---------------------------------------------
929// UnitaryRepresentedGroup trait implementations
930// ---------------------------------------------
931
932impl<T, RowSymbol, ColSymbol> CharacterProperties
933    for UnitaryRepresentedGroup<T, RowSymbol, ColSymbol>
934where
935    RowSymbol: LinearSpaceSymbol + Sync,
936    ColSymbol: CollectionSymbol<CollectionElement = T> + Sync,
937    T: Mul<Output = T>
938        + Inv<Output = T>
939        + Hash
940        + Eq
941        + Clone
942        + Sync
943        + fmt::Debug
944        + FiniteOrder<Int = u32>
945        + Pow<i32, Output = T>,
946    for<'a, 'b> &'b T: Mul<&'a T, Output = T>,
947{
948    type RowSymbol = RowSymbol;
949    type CharTab = RepCharacterTable<RowSymbol, ColSymbol>;
950
951    fn character_table(&self) -> &Self::CharTab {
952        self.irrep_character_table
953            .as_ref()
954            .expect("Irrep character table not found for this group.")
955    }
956
957    fn unitary_represented(&self) -> bool {
958        true
959    }
960}
961
962impl<T, RowSymbol, ColSymbol> IrrepCharTabConstruction
963    for UnitaryRepresentedGroup<T, RowSymbol, ColSymbol>
964where
965    RowSymbol: LinearSpaceSymbol + Sync,
966    ColSymbol: CollectionSymbol<CollectionElement = T> + Sync,
967    T: Mul<Output = T>
968        + Inv<Output = T>
969        + Hash
970        + Eq
971        + Clone
972        + Sync
973        + fmt::Debug
974        + FiniteOrder<Int = u32>
975        + Pow<i32, Output = T>,
976    for<'a, 'b> &'b T: Mul<&'a T, Output = T>,
977{
978    fn set_irrep_character_table(&mut self, chartab: Self::CharTab) {
979        self.irrep_character_table = Some(chartab)
980    }
981}
982
983// ----------------------------------------------
984// MagneticRepresentedGroup trait implementations
985// ----------------------------------------------
986
987impl<T, RowSymbol, UG> CharacterProperties for MagneticRepresentedGroup<T, UG, RowSymbol>
988where
989    RowSymbol: ReducibleLinearSpaceSymbol<Subspace = UG::RowSymbol> + Serialize + DeserializeOwned,
990    T: Mul<Output = T>
991        + Inv<Output = T>
992        + Hash
993        + Eq
994        + Clone
995        + Sync
996        + fmt::Debug
997        + FiniteOrder<Int = u32>
998        + Pow<i32, Output = T>,
999    for<'a, 'b> &'b T: Mul<&'a T, Output = T>,
1000    UG: Clone
1001        + GroupProperties<GroupElement = T>
1002        + ClassProperties<GroupElement = T>
1003        + CharacterProperties,
1004    <UG as ClassProperties>::ClassSymbol: Serialize + DeserializeOwned,
1005    <UG as CharacterProperties>::CharTab: Serialize + DeserializeOwned,
1006    CorepCharacterTable<RowSymbol, <UG as CharacterProperties>::CharTab>:
1007        Serialize + DeserializeOwned,
1008{
1009    type RowSymbol = RowSymbol;
1010    type CharTab = CorepCharacterTable<Self::RowSymbol, UG::CharTab>;
1011
1012    fn character_table(&self) -> &Self::CharTab {
1013        self.ircorep_character_table
1014            .as_ref()
1015            .expect("Ircorep character table not found for this group.")
1016    }
1017
1018    fn unitary_represented(&self) -> bool {
1019        false
1020    }
1021}
1022
1023impl<T, RowSymbol, UG> IrcorepCharTabConstruction for MagneticRepresentedGroup<T, UG, RowSymbol>
1024where
1025    RowSymbol: ReducibleLinearSpaceSymbol<Subspace = UG::RowSymbol> + Serialize + DeserializeOwned,
1026    T: Mul<Output = T>
1027        + Inv<Output = T>
1028        + Hash
1029        + Eq
1030        + Clone
1031        + Sync
1032        + fmt::Debug
1033        + FiniteOrder<Int = u32>
1034        + Pow<i32, Output = T>,
1035    for<'a, 'b> &'b T: Mul<&'a T, Output = T>,
1036    UG: Clone
1037        + GroupProperties<GroupElement = T>
1038        + ClassProperties<GroupElement = T>
1039        + CharacterProperties,
1040    <UG as ClassProperties>::ClassSymbol: Serialize + DeserializeOwned,
1041    <UG as CharacterProperties>::CharTab: Serialize + DeserializeOwned,
1042    CorepCharacterTable<RowSymbol, <UG as CharacterProperties>::CharTab>:
1043        Serialize + DeserializeOwned,
1044{
1045    fn set_ircorep_character_table(&mut self, chartab: Self::CharTab) {
1046        self.ircorep_character_table = Some(chartab);
1047    }
1048}