sva_formula/
series.rs

1// Concern: decides a series' convergence in A and enumerates the lines it yields | Non-concern: placing them on a grid (sva-samples) | IO: (&Series, ceiling) -> bool, Lines, Enumerated, a spacing
2
3use std::f64::consts::TAU;
4
5use crate::affine::{Axis, Reading, affine_in, axis, exact_constant};
6use crate::closed_form::{Body, Bound, IndexId, Part, Series, Unary, children, map_children};
7use crate::complex::C64;
8use crate::env::Env;
9use crate::spectral_sum::atom::{Exp, Factors, Singular, SpectralAtom};
10use crate::table::series::{Shape, read};
11
12/// A tempered limit needs the coefficient polynomially bounded: `1/k` is, `2^k` is not.
13#[derive(Clone, Copy, Debug, PartialEq, Eq)]
14pub enum IndexGrowth {
15    Polynomial,
16    Unbounded,
17}
18
19pub fn summable(s: &Series, env: &dyn Env) -> bool {
20    match s.hi {
21        Bound::Finite(_) => true,
22        Bound::Infinite => growth(&s.term.body, s.index, env) == IndexGrowth::Polynomial,
23    }
24}
25
26/// Decided from where the index sits, never by evaluating a term.
27pub fn growth(f: &Body, k: IndexId, env: &dyn Env) -> IndexGrowth {
28    if !mentions(f, k) {
29        return IndexGrowth::Polynomial;
30    }
31    match f {
32        Body::Index(_) | Body::Line | Body::Const(_) => IndexGrowth::Polynomial,
33        Body::Add(parts) | Body::Mul(parts) | Body::Join(parts) => join(parts, k, env),
34        Body::Div(a, b) => {
35            join(std::slice::from_ref(a), k, env).and(join(std::slice::from_ref(b), k, env))
36        }
37        Body::Pow(base, _) => growth(&base.body, k, env),
38        Body::Keyed { .. } => IndexGrowth::Polynomial,
39        Body::Apply(Unary::Sin | Unary::Cos, arg) => bounded_along(arg, k, Axis::Real, env),
40        Body::Apply(Unary::Exp, arg) if logarithmic(&arg.body, k) => IndexGrowth::Polynomial,
41        Body::Apply(Unary::Exp, arg) => exponential_in(arg, k, env),
42        Body::Channel(of, _) | Body::Crop { of, .. } => growth(&of.body, k, env),
43        Body::Shift { of, .. } | Body::Deriv { of, .. } => growth(&of.body, k, env),
44        Body::Pv(at) | Body::Delta { at, .. } => growth(&at.body, k, env),
45        Body::Series(inner) => growth(&inner.term.body, k, env),
46        _ => IndexGrowth::Unbounded,
47    }
48}
49
50/// A turning exponent is bounded and a falling one is a geometric decay; only a rising real
51/// part outgrows every polynomial.
52fn exponential_in(arg: &Part, k: IndexId, env: &dyn Env) -> IndexGrowth {
53    if let bounded @ IndexGrowth::Polynomial = bounded_along(arg, k, Axis::Imaginary, env) {
54        return bounded;
55    }
56    match affine_in(&arg.body, Reading::Index(k)).and_then(|(slope, _)| slope.exact()) {
57        Some(slope) if slope.re < 0.0 => IndexGrowth::Polynomial,
58        _ => IndexGrowth::Unbounded,
59    }
60}
61
62fn bounded_along(arg: &Part, k: IndexId, wanted: Axis, env: &dyn Env) -> IndexGrowth {
63    if !mentions(&arg.body, k) || axis(&arg.body, env) == wanted {
64        IndexGrowth::Polynomial
65    } else {
66        IndexGrowth::Unbounded
67    }
68}
69
70/// An exponent reaching the index only through a logarithm is a power of it.
71fn logarithmic(f: &Body, k: IndexId) -> bool {
72    match f {
73        _ if !mentions(f, k) => true,
74        Body::Apply(Unary::Log, _) => true,
75        Body::Add(parts) | Body::Mul(parts) => parts.iter().all(|p| logarithmic(&p.body, k)),
76        Body::Div(a, b) => logarithmic(&a.body, k) && logarithmic(&b.body, k),
77        _ => false,
78    }
79}
80
81fn join(parts: &[Part], k: IndexId, env: &dyn Env) -> IndexGrowth {
82    parts.iter().fold(IndexGrowth::Polynomial, |acc, p| {
83        acc.and(growth(&p.body, k, env))
84    })
85}
86
87impl IndexGrowth {
88    fn and(self, other: IndexGrowth) -> IndexGrowth {
89        match (self, other) {
90            (IndexGrowth::Polynomial, IndexGrowth::Polynomial) => IndexGrowth::Polynomial,
91            _ => IndexGrowth::Unbounded,
92        }
93    }
94}
95
96/// What separates a series term's coefficient from its wave.
97pub fn mentions_line(f: &Body) -> bool {
98    reaches(f, &|x| matches!(x, Body::Line))
99}
100
101pub fn mentions(f: &Body, k: IndexId) -> bool {
102    reaches(f, &|x| matches!(x, Body::Index(i) if *i == k))
103}
104
105fn reaches(f: &Body, leaf: &dyn Fn(&Body) -> bool) -> bool {
106    leaf(f) || children(f).iter().any(|p| reaches(&p.body, leaf))
107}
108
109#[derive(Clone, Copy, Debug, PartialEq)]
110pub struct Line {
111    pub hz: f64,
112    pub amp: C64,
113}
114
115/// `tail_db` is the loudest dropped line against the loudest taken one.
116#[derive(Clone, Debug, PartialEq)]
117pub struct Lines {
118    pub taken: Vec<Line>,
119    pub dropped: Vec<Line>,
120    pub tail_db: f64,
121}
122
123pub const AUDIBLE_CEILING_HZ: f64 = 20_000.0;
124
125/// Stops at the ceiling where the frequency closed form leaves the band. Where it never does, every
126/// term piles onto one line, and the coefficient falling below the floor stops it.
127pub fn lines(s: &Series, ceiling: f64, floor_db: f64) -> Lines {
128    let ceiling = ceiling.min(AUDIBLE_CEILING_HZ);
129    let Some(shape) = read(&s.term.body) else {
130        return Lines {
131            taken: Vec::new(),
132            dropped: Vec::new(),
133            tail_db: f64::NEG_INFINITY,
134        };
135    };
136    let voices = places(&shape);
137    let band = voices
138        .iter()
139        .filter_map(|(place, _)| leaves_band(place, s.index, ceiling))
140        .fold(None, |held: Option<i64>, next| {
141            Some(held.map_or(next, |held| held.max(next)))
142        });
143    let hi = match (s.hi, band) {
144        (Bound::Finite(n), Some(last)) => n.min(last),
145        (Bound::Finite(n), None) => n,
146        (Bound::Infinite, Some(last)) => last,
147        (Bound::Infinite, None) => s.lo.saturating_add(MAX_TERMS),
148    };
149    let floor = band.is_none().then(|| 10f64.powf(floor_db / 20.0));
150
151    let mut taken = Vec::new();
152    let mut dropped = Vec::new();
153    let mut first = 0.0f64;
154    for k in s.lo..=hi {
155        let mut here = Vec::new();
156        for (place, weight) in &voices {
157            let (Some(hz), Some(amp)) = (
158                at_index(place, s.index, k).map(|c| c.re),
159                at_index(weight, s.index, k),
160            ) else {
161                continue;
162            };
163            here.push(Line { hz, amp });
164        }
165        let loudest = here.iter().map(|l| l.amp.abs()).fold(0.0f64, f64::max);
166        if k == s.lo {
167            first = loudest;
168        }
169        if let Some(floor) = floor
170            && k > s.lo
171            && first > 0.0
172            && loudest < first * floor
173        {
174            dropped.extend(here);
175            break;
176        }
177        for line in here {
178            if line.hz.abs() <= ceiling {
179                taken.push(line);
180            } else {
181                dropped.push(line);
182            }
183        }
184    }
185    let loudest = |set: &[Line]| set.iter().map(|l| l.amp.abs()).fold(0.0f64, f64::max);
186    let (kept, gone) = (loudest(&taken), loudest(&dropped));
187    Lines {
188        tail_db: if gone > 0.0 && kept > 0.0 {
189            20.0 * (gone / kept).log10()
190        } else {
191            f64::NEG_INFINITY
192        },
193        taken,
194        dropped,
195    }
196}
197
198/// The step a series' own frequency walks: every term lands on a multiple of it. A
199/// delta train's places are instants, not frequencies, and name no such step.
200pub fn spacing(s: &Series) -> Option<f64> {
201    let Some(shape @ Shape::Lines(_)) = read(&s.term.body) else {
202        return None;
203    };
204    let mut held: Option<f64> = None;
205    for (place, _) in places(&shape) {
206        let (slope, offset) = affine_in(&place, Reading::Index(s.index))?;
207        let (slope, offset) = (slope.exact()?.re, offset.exact()?.re);
208        if slope == 0.0 || !slope.is_finite() || !offset.is_finite() {
209            return None;
210        }
211        let steps = offset / slope;
212        if (steps.round() - steps).abs() > TURN_EPSILON * steps.abs().max(1.0) {
213            return None;
214        }
215        match held {
216            Some(step) if step != slope.abs() => return None,
217            _ => held = Some(slope.abs()),
218        }
219    }
220    held
221}
222
223#[derive(Clone, Debug, PartialEq)]
224pub struct Enumerated {
225    pub atoms: Vec<SpectralAtom>,
226    pub dropped: Vec<Line>,
227}
228
229/// A crop of a series is the series of cropped terms: the window lifts off, goes back on
230/// each. `None` where no line closed form reads under it. A delta's `hz` is an instant, not a pitch.
231pub fn line_atoms(s: &Series, ceiling: f64, floor_db: f64) -> Option<Enumerated> {
232    let (body, window) = crate::spectral_sum::image::crop_peeled(&s.term.body);
233    let bare = Series {
234        term: Part::new(s.term.origin, body),
235        ..s.clone()
236    };
237    let singular = match read(&bare.term.body)? {
238        Shape::Deltas(_) => true,
239        Shape::Lines(_) => false,
240    };
241    let found = lines(&bare, ceiling, floor_db);
242    let atoms = found
243        .taken
244        .into_iter()
245        .filter_map(|l| match singular {
246            true => window.is_none_or(|w| w.contains(l.hz)).then(|| {
247                SpectralAtom::new(
248                    l.amp,
249                    Factors::NONE,
250                    Singular::Delta { at: l.hz, order: 0 },
251                    s.term.origin,
252                )
253            }),
254            false => Some(SpectralAtom::new(
255                l.amp,
256                Factors {
257                    exp: Some(Exp::at(0.0, TAU * l.hz)),
258                    ind: window,
259                    ..Factors::NONE
260                },
261                Singular::Regular,
262                s.term.origin,
263            )),
264        })
265        .collect();
266    Some(Enumerated {
267        atoms,
268        dropped: found.dropped,
269    })
270}
271
272/// A whole turn count to floating precision: a tolerance would put a line on a neighbouring
273/// bin, which `exact` cannot carry.
274pub fn commensurate(hz: f64, horizon: f64) -> bool {
275    let turns = hz * horizon;
276    (turns.round() - turns).abs() <= TURN_EPSILON * turns.abs().max(1.0)
277}
278
279const TURN_EPSILON: f64 = 1e-9;
280
281fn places(shape: &Shape) -> Vec<(Body, Body)> {
282    match shape {
283        Shape::Lines(lines) => lines
284            .iter()
285            .map(|l| (l.freq.clone(), l.amp.clone()))
286            .collect(),
287        Shape::Deltas(deltas) => deltas
288            .iter()
289            .map(|d| (d.at.clone(), d.weight.clone()))
290            .collect(),
291    }
292}
293
294/// The last index whose frequency still fits the band, solving `|slope*k + offset| <= ceiling`
295/// at both signs: an offset opposing the slope carries the line back in before it leaves.
296fn leaves_band(place: &Body, k: IndexId, ceiling: f64) -> Option<i64> {
297    let (slope, offset) = affine_in(place, Reading::Index(k))?;
298    let (slope, offset) = (slope.exact()?.re, offset.exact()?.re);
299    if slope == 0.0 {
300        return None;
301    }
302    let ends = [(ceiling - offset) / slope, (-ceiling - offset) / slope];
303    let last = ends[0].max(ends[1]).floor();
304    Some(last.clamp(0.0, MAX_TERMS as f64) as i64 + 1)
305}
306
307const MAX_TERMS: i64 = 1 << 20;
308
309fn at_index(f: &Body, k: IndexId, value: i64) -> Option<C64> {
310    exact_constant(&substitute(f, k, value as f64))
311}
312
313pub fn substitute(f: &Body, k: IndexId, value: f64) -> Body {
314    match f {
315        Body::Index(i) if *i == k => Body::Const(C64::real(value)),
316        other => map_children(other, |p| {
317            Part::new(p.origin, substitute(&p.body, k, value))
318        }),
319    }
320}
sva_formula/series.rs

sva_formula/
series.rs
