1use 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
35pub trait CharacterProperties: ClassProperties
42where
43 Self::RowSymbol: LinearSpaceSymbol,
44 Self::CharTab: CharacterTable<RowSymbol = Self::RowSymbol, ColSymbol = Self::ClassSymbol>,
45{
46 type RowSymbol;
48
49 type CharTab;
52
53 fn character_table(&self) -> &Self::CharTab;
55
56 fn unitary_represented(&self) -> bool;
62}
63
64pub 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 fn set_irrep_character_table(&mut self, chartab: Self::CharTab);
88
89 #[allow(clippy::too_many_lines)]
103 fn construct_irrep_character_table(&mut self) {
104 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 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 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 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 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 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
616pub 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 fn set_ircorep_character_table(&mut self, chartab: Self::CharTab);
637
638 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 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 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 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 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 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
924impl<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
983impl<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}