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

sva_formula/
series.rs
