haar_lib/math/
stern_brocot.rs1use crate::math::continued_fraction::*;
11
12#[derive(Clone, Copy, Debug)]
14pub struct Frac(pub u64, pub u64);
15
16impl PartialEq for Frac {
17 fn eq(&self, Frac(c, d): &Self) -> bool {
18 let Frac(a, b) = self;
19 a * d == b * c
20 }
21}
22impl Eq for Frac {}
23impl PartialOrd for Frac {
24 fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
25 let Frac(a, b) = self;
26 let Frac(c, d) = other;
27 (a * d).partial_cmp(&(b * c))
28 }
29}
30
31#[derive(Clone, Copy, Debug, Eq, PartialEq, Hash)]
33pub struct SBNode(u64, u64, u64, u64);
34
35impl SBNode {
36 pub fn root() -> Self {
38 Self(0, 1, 1, 0)
39 }
40
41 pub fn quadruplet(self) -> (u64, u64, u64, u64) {
43 (self.0, self.1, self.2, self.3)
44 }
45
46 pub fn left_child(self, n: u64) -> Self {
48 let Self(p, q, r, s) = self;
49 Self(p, q, p * n + r, q * n + s)
50 }
51
52 pub fn right_child(self, n: u64) -> Self {
54 let Self(p, q, r, s) = self;
55 Self(p + r * n, q + s * n, r, s)
56 }
57
58 pub fn as_frac(self) -> Frac {
60 let Self(p, q, r, s) = self;
61 Frac(p + r, q + s)
62 }
63
64 pub fn as_f64(self) -> f64 {
66 let Self(p, q, r, s) = self;
67 ((p + r) as f64) / ((q + s) as f64)
68 }
69
70 pub fn range(self) -> (Frac, Frac) {
72 let Self(p, q, r, s) = self;
73 (Frac(p, q), Frac(r, s))
74 }
75
76 pub fn from_frac(a: Frac) -> Option<Self> {
78 Some(SBPath::decode(SBPath::encode(a)?))
79 }
80
81 pub fn lca(a: impl Into<Frac>, b: impl Into<Frac>) -> Option<Self> {
83 let pa = SBPath::encode(a)?;
84 let pb = SBPath::encode(b)?;
85
86 let mut pc = vec![];
87 for (a, b) in std::iter::zip(pa.0, pb.0) {
88 if a == b {
89 pc.push(a);
90 } else {
91 match (a, b) {
92 (SBMove::L(a), SBMove::L(b)) => pc.push(SBMove::L(a.min(b))),
93 (SBMove::R(a), SBMove::R(b)) => pc.push(SBMove::R(a.min(b))),
94 _ => {}
95 }
96 break;
97 }
98 }
99
100 Some(SBPath::decode(SBPath(pc)))
101 }
102
103 pub fn ancestor(a: impl Into<Frac>, mut d: u64) -> Option<Self> {
105 let path = SBPath::encode(a)?;
106
107 let mut path_d = vec![];
108 for m in path.0 {
109 match m {
110 SBMove::L(n) => {
111 path_d.push(SBMove::L(n.min(d)));
112 d = d.saturating_sub(n);
113 }
114 SBMove::R(n) => {
115 path_d.push(SBMove::R(n.min(d)));
116 d = d.saturating_sub(n);
117 }
118 }
119 if d == 0 {
120 break;
121 }
122 }
123 if d > 0 {
124 return None;
125 }
126
127 Some(SBPath::decode(SBPath(path_d)))
128 }
129}
130
131impl From<SBNode> for Frac {
132 fn from(value: SBNode) -> Self {
133 value.as_frac()
134 }
135}
136
137#[derive(Clone, Copy, Debug, Eq, PartialEq, Hash)]
139pub enum SBMove {
140 L(u64),
142 R(u64),
144}
145#[derive(Clone, Debug, Eq, PartialEq, Hash)]
147pub struct SBPath(pub Vec<SBMove>);
148
149impl SBPath {
150 pub fn encode(f: impl Into<Frac>) -> Option<Self> {
154 let Frac(p, q) = f.into();
155 if p == 0 || q == 0 {
156 return None;
157 }
158
159 let mut cfe = continued_fraction(p, q)?;
160 *cfe.last_mut().unwrap() -= 1;
161
162 let mut ret = vec![];
163 if let Some(&m) = cfe.first() {
164 if m != 0 {
165 ret.push(SBMove::R(m));
166 }
167 }
168
169 for (i, m) in cfe.into_iter().skip(1).enumerate() {
170 if i % 2 == 0 {
171 ret.push(SBMove::L(m));
172 } else {
173 ret.push(SBMove::R(m));
174 }
175 }
176
177 Some(Self(ret))
178 }
179
180 pub fn decode(self) -> SBNode {
182 let mut ret = SBNode::root();
183 for m in self.0 {
184 match m {
185 SBMove::L(n) => ret = ret.left_child(n),
186 SBMove::R(n) => ret = ret.right_child(n),
187 }
188 }
189 ret
190 }
191}