outram_foam_basic_lib/polynomial/
roots.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3// Derived from OpenFOAM (www.openfoam.com)
4// Copyright (C) 2004-2023 OpenFOAM Foundation
5// Copyright (C) 2016-2023 OpenCFD Ltd.
6//
7// This file is part of OUTRAM PARK.
8//
9// OUTRAM PARK is free software: you can redistribute it and/or modify it
10// under the terms of the GNU General Public License as published by the
11// Free Software Foundation, either version 3 of the License, or (at your
12// option) any later version.
13//
14// OUTRAM PARK is distributed in the hope that it will be useful, but
15// WITHOUT ANY WARRANTY; without even the implied warranty of
16// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17// General Public License for more details.
18//
19// You should have received a copy of the GNU General Public License along
20// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
21
22/// Root classification tag, matching `Foam::roots::type`.
23#[derive(Debug, Clone, Copy, PartialEq, Eq)]
24#[repr(u64)]
25pub enum RootType {
26    Real = 0,
27    Complex = 1,
28    PosInf = 2,
29    NegInf = 3,
30    Nan = 4,
31}
32
33/// Tagged root container for N roots.
34/// Types are packed 3 bits per root into a u64, identical to C++ `Roots<N>`.
35/// For complex conjugate pairs, slot i holds the real part and slot i+1 holds
36/// the imaginary part; both slots are tagged `Complex`.
37#[derive(Debug, Clone, Copy)]
38pub struct Roots<const N: usize> {
39    pub(crate) values: [f64; N],
40    types: u64,
41}
42
43impl<const N: usize> Roots<N> {
44    /// Value stored at slot `i`.
45    #[inline]
46    pub fn get(&self, i: usize) -> f64 {
47        self.values[i]
48    }
49
50    /// Root type at slot `i`.
51    #[inline]
52    pub fn root_type(&self, i: usize) -> RootType {
53        match (self.types >> (3 * i)) & 7 {
54            0 => RootType::Real,
55            1 => RootType::Complex,
56            2 => RootType::PosInf,
57            3 => RootType::NegInf,
58            _ => RootType::Nan,
59        }
60    }
61
62    /// Overwrite the type at slot `i`.
63    #[inline]
64    pub fn set_type(&mut self, i: usize, t: RootType) {
65        let shift = 3 * i;
66        self.types = (self.types & !(7u64 << shift)) | ((t as u64) << shift);
67    }
68}
69
70impl<const N: usize> std::ops::Index<usize> for Roots<N> {
71    type Output = f64;
72    #[inline]
73    fn index(&self, i: usize) -> &f64 {
74        &self.values[i]
75    }
76}
77
78// ---------------------------------------------------------------------------
79// Constructors for each concrete arity.
80// The C++ template constructors (concatenation by prepend/append) are exposed
81// as named static methods on the concrete impl blocks.
82// ---------------------------------------------------------------------------
83
84impl Roots<1> {
85    /// Single root with the given type and value.
86    #[inline]
87    pub fn new(t: RootType, x: f64) -> Self {
88        Self {
89            values: [x],
90            types: t as u64,
91        }
92    }
93}
94
95impl Roots<2> {
96    /// Concatenate two single roots.  C++ `Roots<2>(Roots<1>, Roots<1>)`.
97    #[inline]
98    pub fn from_pair(a: Roots<1>, b: Roots<1>) -> Self {
99        Self {
100            values: [a.values[0], b.values[0]],
101            types: a.types | (b.types << 3),
102        }
103    }
104
105    /// `Roots<1>` followed by one additional root.
106    /// C++ `Roots<2>(Roots<1>, type, x)`.
107    #[inline]
108    pub fn with_tail(head: Roots<1>, t: RootType, x: f64) -> Self {
109        Self {
110            values: [head.values[0], x],
111            types: head.types | ((t as u64) << 3),
112        }
113    }
114
115    /// Duplicate a single root into both slots.
116    /// C++ `Roots<2>(r, r)`.
117    #[inline]
118    pub fn both(r: Roots<1>) -> Self {
119        Self::from_pair(r, r)
120    }
121}
122
123impl Roots<3> {
124    /// All three slots get the same type and value.
125    /// C++ `Roots<3>(type, x)`.
126    #[inline]
127    pub fn uniform(t: RootType, x: f64) -> Self {
128        let ti = t as u64;
129        Self {
130            values: [x, x, x],
131            types: ti | (ti << 3) | (ti << 6),
132        }
133    }
134
135    /// Concatenate `Roots<1>` then `Roots<2>`.
136    /// C++ `Roots<3>(Roots<1>, Roots<2>)`.
137    #[inline]
138    pub fn concat_1_2(a: Roots<1>, b: Roots<2>) -> Self {
139        Self {
140            values: [a.values[0], b.values[0], b.values[1]],
141            types: a.types | (b.types << 3),
142        }
143    }
144
145    /// Concatenate `Roots<2>` then `Roots<1>`.
146    /// C++ `Roots<3>(Roots<2>, Roots<1>)`.
147    #[inline]
148    pub fn concat_2_1(a: Roots<2>, b: Roots<1>) -> Self {
149        Self {
150            values: [a.values[0], a.values[1], b.values[0]],
151            types: a.types | (b.types << 6),
152        }
153    }
154
155    /// `Roots<2>` followed by one additional root.
156    /// C++ `Roots<3>(Roots<2>, type, x)`.
157    #[inline]
158    pub fn with_tail(head: Roots<2>, t: RootType, x: f64) -> Self {
159        Self::concat_2_1(head, Roots::<1>::new(t, x))
160    }
161}
162
163#[cfg(test)]
164mod tests {
165    use super::*;
166
167    #[test]
168    fn encoding_round_trip() {
169        let r = Roots::<1>::new(RootType::Real, 3.14);
170        assert_eq!(r.root_type(0), RootType::Real);
171        assert_eq!(r.get(0), 3.14);
172    }
173
174    #[test]
175    fn concat_2() {
176        let a = Roots::<1>::new(RootType::Real, 1.0);
177        let b = Roots::<1>::new(RootType::Nan, 0.0);
178        let r = Roots::<2>::from_pair(a, b);
179        assert_eq!(r.root_type(0), RootType::Real);
180        assert_eq!(r.root_type(1), RootType::Nan);
181        assert_eq!(r[0], 1.0);
182        assert_eq!(r[1], 0.0);
183    }
184
185    #[test]
186    fn concat_3() {
187        let a = Roots::<1>::new(RootType::Real, 1.0);
188        let b = Roots::<2>::from_pair(
189            Roots::<1>::new(RootType::Real, 2.0),
190            Roots::<1>::new(RootType::Nan, 0.0),
191        );
192        let r = Roots::<3>::concat_1_2(a, b);
193        assert_eq!(r.root_type(0), RootType::Real);
194        assert_eq!(r.root_type(1), RootType::Real);
195        assert_eq!(r.root_type(2), RootType::Nan);
196        assert_eq!(r[0], 1.0);
197        assert_eq!(r[1], 2.0);
198    }
199
200    #[test]
201    fn uniform_3() {
202        let r = Roots::<3>::uniform(RootType::Complex, -5.0);
203        for i in 0..3 {
204            assert_eq!(r.root_type(i), RootType::Complex);
205            assert_eq!(r[i], -5.0);
206        }
207    }
208}
outram_foam_basic_lib/polynomial/roots.rs

outram_foam_basic_lib/polynomial/
roots.rs
