sva_formula/spectral_sum/
kink.rs

1// Concern: splits a form at the crossing of a min or max of two lines | Non-concern: lowering either half (build.rs) | IO: (&Body) -> the same value as two crops, or nothing
2
3use crate::affine::exact_affine_read;
4use crate::closed_form::{Body, Edge, Fold, Part, map_children};
5use crate::complex::C64;
6use crate::through::Reads;
7
8pub fn split(f: &Body, reads: &dyn Reads) -> Option<Body> {
9    let (kink, (at, below, above)) = first(f, reads)?;
10    let half = |with: &Body, l: Edge, r: Edge| {
11        let scope = reads.scope();
12        Part::bare(Body::Crop {
13            of: Part::bare(replaced(f, (scope, &kink, with), reads)),
14            l,
15            r,
16            rise: 0.0,
17            fall: 0.0,
18        })
19    };
20    Some(Body::Add(vec![
21        half(&below, Edge::NegInf, Edge::at(at)),
22        half(&above, Edge::at(at), Edge::PosInf),
23    ]))
24}
25
26/// A constant, or a slope and the instant it is zero at: `k*(t - r)` keeps `r` itself. Each
27/// fold of a written constant rounds once, at most 2^-53 of its result, as `exact_affine`.
28#[derive(Clone, Copy, Debug)]
29pub enum Line {
30    Flat(f64),
31    Sloped { slope: f64, root: f64 },
32}
33
34impl Line {
35    fn slope(self) -> f64 {
36        match self {
37            Line::Flat(_) => 0.0,
38            Line::Sloped { slope, .. } => slope,
39        }
40    }
41
42    fn scaled(self, k: f64, divide: bool) -> Line {
43        let by = |v: f64| if divide { v / k } else { v * k };
44        match self {
45            Line::Flat(v) => Line::Flat(by(v)),
46            Line::Sloped { slope, root } => Line::Sloped {
47                slope: by(slope),
48                root,
49            },
50        }
51    }
52
53    fn plus(self, other: Line) -> Line {
54        match (self, other) {
55            (Line::Flat(a), Line::Flat(b)) => Line::Flat(a + b),
56            (Line::Flat(c), Line::Sloped { slope, root })
57            | (Line::Sloped { slope, root }, Line::Flat(c)) => match c == 0.0 {
58                true => Line::Sloped { slope, root },
59                false => Line::Sloped {
60                    slope,
61                    root: root - c / slope,
62                },
63            },
64            (Line::Sloped { slope: s, root: r }, Line::Sloped { slope: u, root: q }) => {
65                match (s + u, r == q) {
66                    (0.0, _) => Line::Flat(-(s * r + u * q)),
67                    (sum, true) => Line::Sloped {
68                        slope: sum,
69                        root: r,
70                    },
71                    (sum, false) => Line::Sloped {
72                        slope: sum,
73                        root: (s * r + u * q) / sum,
74                    },
75                }
76            }
77        }
78    }
79}
80
81pub(crate) fn line(f: &Body, reads: &dyn Reads) -> Option<Line> {
82    if let Body::Node(id) = f {
83        return reads.line(*id);
84    }
85    let line = |g: &Body| line(g, reads);
86    let real = |c: C64| c.is_real().then_some(c.re);
87    let read = match f {
88        Body::Const(c) => Some(Line::Flat(real(*c)?)),
89        Body::Line => Some(Line::Sloped {
90            slope: 1.0,
91            root: 0.0,
92        }),
93        Body::Shift { by, of } => match line(&of.body)? {
94            Line::Sloped { slope, root } => Some(Line::Sloped {
95                slope,
96                root: root + by,
97            }),
98            flat => Some(flat),
99        },
100        Body::Add(parts) => parts
101            .iter()
102            .try_fold(Line::Flat(0.0), |acc, p| Some(acc.plus(line(&p.body)?))),
103        Body::Mul(parts) => {
104            let mut sloped: Option<Line> = None;
105            let mut k = 1.0;
106            for p in parts {
107                match line(&p.body)? {
108                    Line::Flat(v) => k *= v,
109                    held if sloped.is_none() => sloped = Some(held),
110                    _ => return None,
111                }
112            }
113            Some(sloped.map_or(Line::Flat(k), |l| l.scaled(k, false)))
114        }
115        Body::Div(num, den) => match line(&den.body)? {
116            Line::Flat(k) if k != 0.0 => Some(line(&num.body)?.scaled(k, true)),
117            _ => None,
118        },
119        _ => None,
120    };
121    let read = read.or_else(|| {
122        let (a, b) = exact_affine_read(f, reads)?;
123        match (real(a)?, real(b)?) {
124            (0.0, b) => Some(Line::Flat(b)),
125            (a, b) => Some(Line::Sloped {
126                slope: a,
127                root: -b / a,
128            }),
129        }
130    })?;
131    match read {
132        Line::Sloped { slope: 0.0, .. } => Some(Line::Flat(0.0)),
133        held => Some(held),
134    }
135    .filter(|l| match *l {
136        Line::Flat(v) => v.is_finite(),
137        Line::Sloped { slope, root } => slope.is_finite() && root.is_finite(),
138    })
139}
140
141/// Each half is its own line, read from its root. Against a constant, the crossing is the
142/// first double where the line as evaluated passes it, so a `max` never falls under nor a
143/// `min` climbs over it, and a line meeting zero crosses at its root.
144fn crossing(f: &Body, reads: &dyn Reads) -> Option<Crossing> {
145    let Body::Fold(op @ (Fold::Min | Fold::Max), args) = f else {
146        return None;
147    };
148    let [p, q] = args.as_slice() else {
149        return None;
150    };
151    let (l1, l2) = (line(&p.body, reads)?, line(&q.body, reads)?);
152    if l1.slope() == l2.slope() {
153        return None;
154    }
155    let max = *op == Fold::Max;
156    match (l1, l2) {
157        (Line::Flat(c), sloped @ Line::Sloped { slope, root })
158        | (sloped @ Line::Sloped { slope, root }, Line::Flat(c)) => {
159            let taken = |t: f64| {
160                let v = slope * (t - root);
161                if max { v >= c } else { v <= c }
162            };
163            let flat = Body::Const(C64::real(c));
164            let estimate = root + c / slope;
165            match (slope > 0.0) == max {
166                true => Some((first_double(taken, estimate)?, flat, written(sloped))),
167                false => Some((
168                    first_double(|t| !taken(t), estimate)?,
169                    written(sloped),
170                    flat,
171                )),
172            }
173        }
174        (Line::Sloped { slope: s, root: r }, Line::Sloped { slope: u, root: v }) => {
175            let at = match r == v {
176                true => r,
177                false => (s * r - u * v) / (s - u),
178            };
179            let p_below = (s > u) != max;
180            let (below, above) = if p_below { (l1, l2) } else { (l2, l1) };
181            at.is_finite().then(|| (at, written(below), written(above)))
182        }
183        (Line::Flat(_), Line::Flat(_)) => None,
184    }
185}
186
187fn written(l: Line) -> Body {
188    match l {
189        Line::Flat(c) => Body::Const(C64::real(c)),
190        Line::Sloped { slope, root } => Body::Mul(vec![
191            Part::bare(Body::Const(C64::real(slope))),
192            Part::bare(Body::Shift {
193                by: root,
194                of: Part::bare(Body::Line),
195            }),
196        ]),
197    }
198}
199
200/// Doubles in numeric order as integers.
201fn ordered(x: i64) -> i64 {
202    x ^ (((x >> 63) as u64) >> 1) as i64
203}
204
205fn key(x: f64) -> i64 {
206    ordered(x.to_bits() as i64)
207}
208
209fn double(k: i64) -> f64 {
210    f64::from_bits(ordered(k) as u64)
211}
212
213/// The least finite double where `holds`, monotone over the doubles, turns true, galloping
214/// out from `near` and bisecting; nothing where it never turns within the finite doubles.
215fn first_double(holds: impl Fn(f64) -> bool, near: f64) -> Option<f64> {
216    let finite = |k: i64| double(k).is_finite().then_some(k);
217    let near = finite(key(near))?;
218    let mut step = 1i64;
219    let (mut no, mut yes) = (near, near);
220    match holds(double(near)) {
221        true => loop {
222            no = finite(yes.checked_sub(step)?)?;
223            if !holds(double(no)) {
224                break;
225            }
226            yes = no;
227            step = step.checked_mul(2)?;
228        },
229        false => loop {
230            yes = finite(no.checked_add(step)?)?;
231            if holds(double(yes)) {
232                break;
233            }
234            no = yes;
235            step = step.checked_mul(2)?;
236        },
237    }
238    while yes - no > 1 {
239        let mid = no + (yes - no) / 2;
240        match holds(double(mid)) {
241            true => yes = mid,
242            false => no = mid,
243        }
244    }
245    Some(double(yes))
246}
247
248/// A shift, warp or series moves the kink; a derivative would add a delta at the seam.
249pub(crate) fn in_time(f: &Body) -> bool {
250    matches!(
251        f,
252        Body::Add(_)
253            | Body::Mul(_)
254            | Body::Div(..)
255            | Body::Pow(..)
256            | Body::Apply(..)
257            | Body::Fold(..)
258            | Body::Crop { .. }
259            | Body::Join(_)
260            | Body::Channel(..)
261    )
262}
263
264pub type Crossing = (f64, Body, Body);
265
266/// The first min or max of two lines, searched through each ref as through its own form.
267pub(crate) fn first(f: &Body, reads: &dyn Reads) -> Option<(Body, Crossing)> {
268    if let Body::Node(id) = f {
269        return reads.kink(*id);
270    }
271    if let Some(found) = crossing(f, reads) {
272        return Some((f.clone(), found));
273    }
274    if !in_time(f) {
275        return None;
276    }
277    crate::closed_form::children(f)
278        .into_iter()
279        .find_map(|p| first(&p.body, reads))
280}
281
282/// Which split it is for, the kink, and the line read in its place.
283type Replacing<'a> = (usize, &'a Body, &'a Body);
284
285fn replaced(f: &Body, at: Replacing, reads: &dyn Reads) -> Body {
286    let (scope, kink, with) = at;
287    if f == kink {
288        return with.clone();
289    }
290    if let Body::Node(id) = f {
291        return reads
292            .replaced(*id, scope, kink, with)
293            .map_or_else(|| f.clone(), Body::Node);
294    }
295    if !in_time(f) {
296        return f.clone();
297    }
298    map_children(f, |p| replaced_part(p, at, reads))
299}
300
301pub(crate) fn replaced_part(p: &Part, at: Replacing, reads: &dyn Reads) -> Part {
302    Part::new(p.origin, replaced(&p.body, at, reads))
303}
sva_formula/spectral_sum/kink.rs

sva_formula/spectral_sum/
kink.rs
