use std::collections::{BTreeMap, HashSet};
use std::fmt;
use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::One;
use crate::group::{CenterSubgroup, GlobalForm, GroupId, RootSystem};
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub enum Series {
B,
C,
D,
}
impl Series {
fn name(self) -> &'static str {
match self {
Series::B => "B (SO(2r+1))",
Series::C => "C (Sp(2r))",
Series::D => "D (SO(2r))",
}
}
fn min_rank(self) -> usize {
match self {
Series::B | Series::C => 2,
Series::D => 3,
}
}
pub fn root_system(self, r: usize) -> RootSystem {
match self {
Series::B => RootSystem::B(r),
Series::C => RootSystem::C(r),
Series::D => RootSystem::D(r),
}
}
fn published_group(self, r: usize) -> GroupId {
GroupId {
root_system: self.root_system(r),
form: match self {
Series::B => GlobalForm::Quotient(CenterSubgroup::Z2),
Series::C => GlobalForm::SimplyConnected,
Series::D => GlobalForm::Quotient(CenterSubgroup::DVector),
},
}
}
pub(crate) fn cover_group(self, r: usize) -> GroupId {
GroupId {
root_system: self.root_system(r),
form: GlobalForm::SimplyConnected,
}
}
fn low_rank_redirect(self, form: GlobalForm) -> &'static str {
match (self, form) {
(Series::B | Series::C, GlobalForm::SimplyConnected) => "use SU(2) instead",
(Series::B | Series::C, GlobalForm::Quotient(_)) => {
"use SU(2) with integer j only instead"
}
(Series::D, GlobalForm::SimplyConnected) => "use SU(2)×SU(2) instead",
(Series::D, GlobalForm::Quotient(_)) => "use SU(2)×SU(2) with j₁+j₂ integer instead",
}
}
}
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum BcdError {
EmptyLabel,
NegativeDynkin {
dynkin: Vec<i64>,
},
ExcludedRank {
series: Series,
rank: usize,
redirect: &'static str,
},
UnsupportedRootSystem {
root_system: RootSystem,
},
RankMismatch {
expected: usize,
got: usize,
},
NotAdmissible {
group: GroupId,
dynkin: Vec<i64>,
},
GroupMismatch {
a: (Series, usize),
b: (Series, usize),
},
CommutatorViolation {
series: Series,
relation: &'static str,
i: usize,
j: usize,
},
}
impl fmt::Display for BcdError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
BcdError::EmptyLabel => write!(f, "B/C/D irrep label must be non-empty"),
BcdError::NegativeDynkin { dynkin } => {
write!(f, "Dynkin label has a negative component: {dynkin:?}")
}
BcdError::ExcludedRank {
series,
rank,
redirect,
} => write!(
f,
"series {} rank {rank} is an excluded low-rank isomorphism: {redirect}",
series.name()
),
BcdError::UnsupportedRootSystem { root_system } => write!(
f,
"root system {root_system} is not one of the B/C/D families of this module"
),
BcdError::RankMismatch { expected, got } => write!(
f,
"Dynkin label of length {got} for a rank-{expected} group"
),
BcdError::NotAdmissible { group, dynkin } => write!(
f,
"Dynkin label {dynkin:?} is a dominant integral weight of the cover but not \
a representation of {group:?}: its central character is non-trivial on the \
quotiented-out subgroup of the center — construct it through \
Irrep::from_dynkin_in with the simply-connected form"
),
BcdError::GroupMismatch { a, b } => write!(
f,
"directproduct across distinct groups {:?} and {:?}",
a, b
),
BcdError::CommutatorViolation {
series,
relation,
i,
j,
} => write!(
f,
"series {} defining-seed commutator self-check failed: {relation} \
(generators i={i}, j={j})",
series.name()
),
}
}
}
impl std::error::Error for BcdError {}
#[derive(Clone, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub struct Irrep {
series: Series,
two_weight: Box<[i64]>,
}
impl crate::cache::CacheKeyCharge for Irrep {
fn key_bytes(&self) -> usize {
std::mem::size_of::<Self>().saturating_add(std::mem::size_of_val(self.two_weight.as_ref()))
}
}
impl Irrep {
pub fn from_dynkin(series: Series, dynkin: &[i64]) -> Result<Self, BcdError> {
if dynkin.is_empty() {
return Err(BcdError::EmptyLabel);
}
let group = series.published_group(dynkin.len());
Self::from_dynkin_in(&group, dynkin)
}
pub fn from_dynkin_in(group: &GroupId, dynkin: &[i64]) -> Result<Self, BcdError> {
if dynkin.is_empty() {
return Err(BcdError::EmptyLabel);
}
let (series, r) = match group.root_system {
RootSystem::B(r) => (Series::B, r),
RootSystem::C(r) => (Series::C, r),
RootSystem::D(r) => (Series::D, r),
other => return Err(BcdError::UnsupportedRootSystem { root_system: other }),
};
if dynkin.len() != r {
return Err(BcdError::RankMismatch {
expected: r,
got: dynkin.len(),
});
}
if r < series.min_rank() {
return Err(BcdError::ExcludedRank {
series,
rank: r,
redirect: series.low_rank_redirect(group.form),
});
}
if dynkin.iter().any(|&a| a < 0) {
return Err(BcdError::NegativeDynkin {
dynkin: dynkin.to_vec(),
});
}
if !group.admits(dynkin) {
return Err(BcdError::NotAdmissible {
group: *group,
dynkin: dynkin.to_vec(),
});
}
Ok(Irrep {
series,
two_weight: dynkin_to_two_partition(series, dynkin).into_boxed_slice(),
})
}
pub fn trivial(series: Series, r: usize) -> Result<Self, BcdError> {
Self::from_dynkin(series, &vec![0i64; r])
}
pub(crate) fn from_two_weight(series: Series, two_weight: Vec<i64>) -> Irrep {
Irrep {
series,
two_weight: two_weight.into_boxed_slice(),
}
}
pub fn series(&self) -> Series {
self.series
}
pub fn rank(&self) -> usize {
self.two_weight.len()
}
pub fn two_partition(&self) -> &[i64] {
&self.two_weight
}
pub fn partition(&self) -> Option<Vec<i64>> {
if self.two_weight.iter().any(|x| x % 2 != 0) {
return None;
}
Some(self.two_weight.iter().map(|x| x / 2).collect())
}
pub fn is_spinor(&self) -> bool {
self.two_weight.iter().any(|x| x % 2 != 0)
}
pub fn dynkin(&self) -> Vec<i64> {
two_partition_to_dynkin(self.series, &self.two_weight)
}
pub fn dim(&self) -> BigInt {
let r = self.rank();
let two_rho = two_rho(self.series, r);
let mut acc = Ratio::<BigInt>::one();
for alpha in positive_roots(self.series, r) {
let two_lam = dot(&self.two_weight, &alpha);
let two_rho_a = dot(&two_rho, &alpha);
acc *= Ratio::new(BigInt::from(two_lam + two_rho_a), BigInt::from(two_rho_a));
}
acc.to_integer()
}
pub fn dual(&self) -> Irrep {
let mut w = self.two_weight.to_vec();
if self.series == Series::D && self.rank() % 2 == 1 {
let last = self.rank() - 1;
w[last] = -w[last];
}
Irrep {
series: self.series,
two_weight: w.into_boxed_slice(),
}
}
pub fn frobenius_schur(&self) -> i32 {
if *self != self.dual() {
return 0;
}
let twice = dot(&self.two_weight, &two_rho_vee(self.series, self.rank()));
debug_assert_eq!(twice % 2, 0, "⟨λ,2ρ^∨⟩ must be an integer");
if (twice / 2).rem_euclid(2) == 0 {
1
} else {
-1
}
}
pub fn two_weight_multiplicities(&self) -> BTreeMap<Vec<i64>, u64> {
freudenthal(self.series, &self.two_weight)
}
pub fn weight_multiplicities(&self) -> Option<BTreeMap<Vec<i64>, u64>> {
if self.is_spinor() {
return None;
}
Some(
self.two_weight_multiplicities()
.into_iter()
.map(|(k, v)| (k.iter().map(|x| x / 2).collect(), v))
.collect(),
)
}
}
fn dynkin_to_two_partition(series: Series, a: &[i64]) -> Vec<i64> {
let r = a.len();
let mut lam = vec![0i64; r];
match series {
Series::B => {
lam[r - 1] = a[r - 1];
for i in (0..r - 1).rev() {
lam[i] = lam[i + 1] + 2 * a[i];
}
}
Series::C => {
lam[r - 1] = 2 * a[r - 1];
for i in (0..r - 1).rev() {
lam[i] = lam[i + 1] + 2 * a[i];
}
}
Series::D => {
lam[r - 1] = a[r - 1] - a[r - 2];
lam[r - 2] = a[r - 1] + a[r - 2];
for i in (0..r - 2).rev() {
lam[i] = lam[i + 1] + 2 * a[i];
}
}
}
lam
}
fn two_partition_to_dynkin(series: Series, two_lam: &[i64]) -> Vec<i64> {
let r = two_lam.len();
let mut a = vec![0i64; r];
for i in 0..r - 1 {
a[i] = (two_lam[i] - two_lam[i + 1]) / 2;
}
match series {
Series::B => a[r - 1] = two_lam[r - 1],
Series::C => a[r - 1] = two_lam[r - 1] / 2,
Series::D => {
a[r - 2] = (two_lam[r - 2] - two_lam[r - 1]) / 2;
a[r - 1] = (two_lam[r - 2] + two_lam[r - 1]) / 2;
}
}
a
}
fn dot(u: &[i64], v: &[i64]) -> i64 {
u.iter().zip(v).map(|(a, b)| a * b).sum()
}
fn two_rho(series: Series, r: usize) -> Vec<i64> {
(0..r)
.map(|i0| {
let i = i0 as i64 + 1;
let rr = r as i64;
match series {
Series::B => 2 * rr - 2 * i + 1,
Series::C => 2 * rr - 2 * i + 2,
Series::D => 2 * rr - 2 * i,
}
})
.collect()
}
fn two_rho_vee(series: Series, r: usize) -> Vec<i64> {
(0..r)
.map(|i0| {
let i = i0 as i64 + 1;
let rr = r as i64;
match series {
Series::B => 2 * (rr - i) + 2,
Series::C => 2 * (rr - i) + 1,
Series::D => 2 * (rr - i),
}
})
.collect()
}
fn positive_roots(series: Series, r: usize) -> Vec<Vec<i64>> {
let mut roots = Vec::new();
for i in 0..r {
for j in i + 1..r {
let mut minus = vec![0i64; r];
minus[i] = 1;
minus[j] = -1;
roots.push(minus);
let mut plus = vec![0i64; r];
plus[i] = 1;
plus[j] = 1;
roots.push(plus);
}
}
match series {
Series::B => {
for i in 0..r {
let mut e = vec![0i64; r];
e[i] = 1;
roots.push(e);
}
}
Series::C => {
for i in 0..r {
let mut e = vec![0i64; r];
e[i] = 2;
roots.push(e);
}
}
Series::D => {}
}
roots
}
fn sort_desc_with_parity(v: &[i64]) -> (Vec<i64>, i32) {
let n = v.len();
let mut inv = 0usize;
for i in 0..n {
for j in i + 1..n {
if v[i] < v[j] {
inv += 1;
}
}
}
let mut out = v.to_vec();
out.sort_unstable_by(|a, b| b.cmp(a));
(out, if inv.is_multiple_of(2) { 1 } else { -1 })
}
fn weyl_dominant(series: Series, v: &[i64]) -> Vec<i64> {
let negcount = v.iter().filter(|&&x| x < 0).count();
let absv: Vec<i64> = v.iter().map(|x| x.abs()).collect();
let (mut sorted, _) = sort_desc_with_parity(&absv);
if series == Series::D && !negcount.is_multiple_of(2) {
let last = sorted.len() - 1;
if sorted[last] != 0 {
sorted[last] = -sorted[last];
}
}
sorted
}
fn dominant_conjugate_signed(series: Series, two_v: &[i64]) -> Option<(Vec<i64>, i32)> {
let negcount = two_v.iter().filter(|&&x| x < 0).count();
let absv: Vec<i64> = two_v.iter().map(|x| x.abs()).collect();
for i in 0..absv.len() {
for j in i + 1..absv.len() {
if absv[i] == absv[j] {
return None;
}
}
}
let (mut sorted, perm_sign) = sort_desc_with_parity(&absv);
match series {
Series::B | Series::C => {
if absv.contains(&0) {
return None;
}
let sign = perm_sign * if negcount.is_multiple_of(2) { 1 } else { -1 };
Some((sorted, sign))
}
Series::D => {
let last = sorted.len() - 1;
if !negcount.is_multiple_of(2) && sorted[last] != 0 {
sorted[last] = -sorted[last];
}
Some((sorted, perm_sign))
}
}
}
fn freudenthal(series: Series, lambda: &[i64]) -> BTreeMap<Vec<i64>, u64> {
let r = lambda.len();
let two_rho: Vec<i64> = two_rho(series, r).into_iter().map(|x| 2 * x).collect();
let roots: Vec<Vec<i64>> = positive_roots(series, r)
.into_iter()
.map(|a| a.into_iter().map(|x| 2 * x).collect())
.collect();
let mut doms: Vec<(i64, Vec<i64>)> = enumerate_dominant_below(series, lambda)
.into_iter()
.map(|mu| (depth(series, lambda, &mu), mu))
.collect();
doms.sort();
let casimir = |w: &[i64]| -> i128 { (dot(w, w) + dot(w, &two_rho)) as i128 };
let cas_lambda = casimir(lambda);
let mut mult: BTreeMap<Vec<i64>, u64> = BTreeMap::new();
for (_, mu) in &doms {
if mu == lambda {
mult.insert(mu.clone(), 1);
continue;
}
let denom = cas_lambda - casimir(mu);
debug_assert!(denom > 0, "Freudenthal denominator must be positive");
let mut num: i128 = 0;
for alpha in &roots {
let aa = dot(alpha, alpha) as i128;
let mu_a = dot(mu, alpha) as i128;
let mut k: i128 = 1;
loop {
let shifted: Vec<i64> = mu
.iter()
.zip(alpha)
.map(|(&m, &al)| m + (k as i64) * al)
.collect();
let dom = weyl_dominant(series, &shifted);
match mult.get(&dom) {
Some(&m) if m > 0 => {
num += 2 * (mu_a + k * aa) * (m as i128);
k += 1;
}
_ => break,
}
}
}
debug_assert_eq!(num % denom, 0, "Freudenthal must divide exactly");
let m = num / denom;
if m > 0 {
mult.insert(mu.clone(), m as u64);
}
}
mult
}
fn depth(series: Series, lambda: &[i64], mu: &[i64]) -> i64 {
let d: Vec<i64> = lambda.iter().zip(mu).map(|(&l, &m)| l - m).collect();
simple_root_coeffs(series, &d)
.map(|c| c.iter().sum())
.unwrap_or(-1)
}
fn simple_root_coeffs(series: Series, dd: &[i64]) -> Option<Vec<i64>> {
let r = dd.len();
let mut c = vec![0i64; r];
let mut acc = 0i64;
for i in 0..r {
acc += dd[i];
if acc % 2 != 0 {
return None;
}
c[i] = acc / 2; }
let total = acc; match series {
Series::B => {
}
Series::C => {
if total % 4 != 0 {
return None;
}
c[r - 1] = total / 4;
}
Series::D => {
if total % 4 != 0 {
return None;
}
c[r - 1] = total / 4;
c[r - 2] = total / 4 - dd[r - 1] / 2;
}
}
if c.iter().all(|&x| x >= 0) {
Some(c)
} else {
None
}
}
fn enumerate_dominant_below(series: Series, lambda: &[i64]) -> Vec<Vec<i64>> {
let r = lambda.len();
let hi = lambda[0]; let mut out = Vec::new();
let mut cur = vec![0i64; r];
enum_dom_rec(series, lambda, hi, 0, &mut cur, &mut out);
out
}
fn enum_dom_rec(
series: Series,
lambda: &[i64],
hi: i64,
pos: usize,
cur: &mut Vec<i64>,
out: &mut Vec<Vec<i64>>,
) {
let r = cur.len();
if pos == r {
if simple_root_coeffs(series, &sub(lambda, cur)).is_some() {
out.push(cur.clone());
}
return;
}
let upper = if pos == 0 { hi } else { cur[pos - 1] };
let lower = if series == Series::D && pos == r - 1 {
-cur[pos - 1]
} else {
0
};
let mut v = upper;
while v >= lower {
cur[pos] = v;
enum_dom_rec(series, lambda, hi, pos + 1, cur, out);
v -= 2;
}
cur[pos] = 0;
}
fn sub(a: &[i64], b: &[i64]) -> Vec<i64> {
a.iter().zip(b).map(|(&x, &y)| x - y).collect()
}
fn weyl_orbit(series: Series, mu: &[i64]) -> Vec<Vec<i64>> {
let r = mu.len();
let mut set: HashSet<Vec<i64>> = HashSet::new();
for signs in 0u32..(1u32 << r) {
let flips = signs.count_ones() as usize;
if series == Series::D && !flips.is_multiple_of(2) {
continue;
}
let signed: Vec<i64> = (0..r)
.map(|i| if signs & (1 << i) != 0 { -mu[i] } else { mu[i] })
.collect();
permute_into(&signed, &mut set);
}
set.into_iter().collect()
}
fn permute_into(v: &[i64], set: &mut HashSet<Vec<i64>>) {
let mut idx: Vec<usize> = (0..v.len()).collect();
permute_rec(v, &mut idx, 0, set);
}
fn permute_rec(v: &[i64], idx: &mut Vec<usize>, k: usize, set: &mut HashSet<Vec<i64>>) {
let n = idx.len();
if k == n {
set.insert(idx.iter().map(|&i| v[i]).collect());
return;
}
for i in k..n {
idx.swap(k, i);
permute_rec(v, idx, k + 1, set);
idx.swap(k, i);
}
}
pub fn directproduct(a: &Irrep, b: &Irrep) -> Result<BTreeMap<Irrep, u32>, BcdError> {
if a.series != b.series || a.rank() != b.rank() {
return Err(BcdError::GroupMismatch {
a: (a.series, a.rank()),
b: (b.series, b.rank()),
});
}
let series = a.series;
let r = a.rank();
let two_rho = two_rho(series, r);
let two_a_rho: Vec<i64> = a
.two_weight
.iter()
.zip(&two_rho)
.map(|(&av, &tr)| av + tr)
.collect();
let mut acc: BTreeMap<Vec<i64>, i64> = BTreeMap::new();
for (mu, &m) in &b.two_weight_multiplicities() {
let m = m as i64;
for omega in weyl_orbit(series, mu) {
let two_xi: Vec<i64> = two_a_rho
.iter()
.zip(&omega)
.map(|(&ar, &w)| ar + w)
.collect();
if let Some((dom, sign)) = dominant_conjugate_signed(series, &two_xi) {
let c: Vec<i64> = dom.iter().zip(&two_rho).map(|(&d, &tr)| d - tr).collect();
*acc.entry(c).or_insert(0) += sign as i64 * m;
}
}
}
let mut result: BTreeMap<Irrep, u32> = BTreeMap::new();
for (c, n) in acc {
debug_assert!(n >= 0, "Racah–Speiser multiplicity must be non-negative");
if n > 0 {
result.insert(
Irrep {
series,
two_weight: c.into_boxed_slice(),
},
n as u32,
);
}
}
Ok(result)
}
mod seeds;
pub use seeds::{check_commutators, defining_seed, spinor_seeds, CommReport, Seed};
#[cfg(feature = "cgc-gen")]
mod linalg;
#[cfg(feature = "cgc-gen")]
mod sweep;
#[cfg(feature = "cgc-gen")]
pub use sweep::{
decompose, decompose_defining_product, Block, Decomposition, Generators, SweepError,
};
#[cfg(feature = "cgc-gen")]
mod catalog;
#[cfg(feature = "cgc-gen")]
pub use catalog::{CanonicalCatalog, CatalogCgc, CatalogError};
#[cfg(feature = "cgc-gen")]
mod fr;
#[cfg(feature = "cgc-gen")]
pub use fr::{
cgc_sweeps, check_f_unitarity, check_hexagon, check_pentagon, f_symbol, r_symbol, FBlock,
FrError, RBlock,
};
#[cfg(feature = "cgc-gen")]
pub fn bcd_authority_fingerprint() -> &'static [u8] {
b"racah:bcd-bootstrap:ref=qspace-v4-dd2cc7e:kron=a-fast:parent=canonical-parent:sweep=gs2-qrpos-posdiag:sort=maxweight-desc:sign=first-significant-positive:align=procrustes-canonical:tol=cg-eps-tier:epoch=2"
}
#[cfg(test)]
mod tests;
#[cfg(all(test, feature = "cgc-gen"))]
mod qspace_oracle_tests;