qualia_core_db/solvers/calculus/
quadrature.rs

1//! Adaptive scalar quadrature with bounded, non-recursive work stacks.
2
3#[derive(Debug, Clone, Copy, PartialEq)]
4pub enum QuadratureError {
5    InvalidDomain,
6    NonFiniteIntegrand { x: f64 },
7    EvaluationBudgetExceeded,
8    WorkspaceExceeded,
9}
10
11#[repr(C)]
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct QuadratureResult {
14    pub value: f64,
15    pub absolute_error: f64,
16    pub evaluations: u32,
17    pub intervals: u32,
18}
19
20#[derive(Clone, Copy, Default)]
21struct SimpsonPanel {
22    a: f64,
23    b: f64,
24    fa: f64,
25    fm: f64,
26    fb: f64,
27    whole: f64,
28    tolerance: f64,
29}
30
31const MAX_ADAPTIVE_PANELS: usize = 512;
32
33fn checked_eval<F>(function: &F, x: f64) -> Result<f64, QuadratureError>
34where
35    F: Fn(f64) -> f64,
36{
37    let value = function(x);
38    if value.is_finite() {
39        Ok(value)
40    } else {
41        Err(QuadratureError::NonFiniteIntegrand { x })
42    }
43}
44
45pub fn adaptive_simpson<F>(
46    function: F,
47    a: f64,
48    b: f64,
49    absolute_tolerance: f64,
50    max_evaluations: u32,
51) -> Result<QuadratureResult, QuadratureError>
52where
53    F: Fn(f64) -> f64,
54{
55    if !a.is_finite()
56        || !b.is_finite()
57        || b <= a
58        || !absolute_tolerance.is_finite()
59        || absolute_tolerance <= 0.0
60        || max_evaluations < 3
61    {
62        return Err(QuadratureError::InvalidDomain);
63    }
64
65    let midpoint = 0.5 * (a + b);
66    let fa = checked_eval(&function, a)?;
67    let fm = checked_eval(&function, midpoint)?;
68    let fb = checked_eval(&function, b)?;
69    let whole = (b - a) * (fa + 4.0 * fm + fb) / 6.0;
70    let mut evaluations = 3_u32;
71    let mut stack = [SimpsonPanel::default(); MAX_ADAPTIVE_PANELS];
72    stack[0] = SimpsonPanel {
73        a,
74        b,
75        fa,
76        fm,
77        fb,
78        whole,
79        tolerance: absolute_tolerance,
80    };
81    let mut stack_len = 1;
82    let mut value = 0.0;
83    let mut error = 0.0;
84    let mut intervals = 0;
85
86    while stack_len > 0 {
87        stack_len -= 1;
88        let panel = stack[stack_len];
89        if evaluations + 2 > max_evaluations {
90            return Err(QuadratureError::EvaluationBudgetExceeded);
91        }
92        let midpoint = 0.5 * (panel.a + panel.b);
93        let left_midpoint = 0.5 * (panel.a + midpoint);
94        let right_midpoint = 0.5 * (midpoint + panel.b);
95        let flm = checked_eval(&function, left_midpoint)?;
96        let frm = checked_eval(&function, right_midpoint)?;
97        evaluations += 2;
98
99        let left = (midpoint - panel.a) * (panel.fa + 4.0 * flm + panel.fm) / 6.0;
100        let right = (panel.b - midpoint) * (panel.fm + 4.0 * frm + panel.fb) / 6.0;
101        let delta = left + right - panel.whole;
102        let local_error = delta.abs() / 15.0;
103        if local_error <= panel.tolerance {
104            value += left + right + delta / 15.0;
105            error += local_error;
106            intervals += 2;
107            continue;
108        }
109
110        if stack_len + 2 > stack.len() {
111            return Err(QuadratureError::WorkspaceExceeded);
112        }
113        let child_tolerance = panel.tolerance * 0.5;
114        stack[stack_len] = SimpsonPanel {
115            a: midpoint,
116            b: panel.b,
117            fa: panel.fm,
118            fm: frm,
119            fb: panel.fb,
120            whole: right,
121            tolerance: child_tolerance,
122        };
123        stack[stack_len + 1] = SimpsonPanel {
124            a: panel.a,
125            b: midpoint,
126            fa: panel.fa,
127            fm: flm,
128            fb: panel.fm,
129            whole: left,
130            tolerance: child_tolerance,
131        };
132        stack_len += 2;
133    }
134
135    Ok(QuadratureResult {
136        value,
137        absolute_error: error,
138        evaluations,
139        intervals,
140    })
141}
142
143const GK15_ABSCISSAE: [f64; 8] = [
144    0.991_455_371_120_812_6,
145    0.949_107_912_342_758_5,
146    0.864_864_423_359_769_1,
147    0.741_531_185_599_394_5,
148    0.586_087_235_467_691_1,
149    0.405_845_151_377_397_2,
150    0.207_784_955_007_898_48,
151    0.0,
152];
153const GK15_WEIGHTS: [f64; 8] = [
154    0.022_935_322_010_529_224,
155    0.063_092_092_629_978_55,
156    0.104_790_010_322_250_19,
157    0.140_653_259_715_525_92,
158    0.169_004_726_639_267_9,
159    0.190_350_578_064_785_42,
160    0.204_432_940_075_298_89,
161    0.209_482_141_084_727_82,
162];
163const G7_WEIGHTS: [f64; 4] = [
164    0.129_484_966_168_869_7,
165    0.279_705_391_489_276_64,
166    0.381_830_050_505_118_9,
167    0.417_959_183_673_469_4,
168];
169
170#[derive(Clone, Copy, Default)]
171struct IntervalPanel {
172    a: f64,
173    b: f64,
174    tolerance: f64,
175}
176
177fn gauss_kronrod_15_panel<F>(
178    function: &F,
179    a: f64,
180    b: f64,
181) -> Result<(f64, f64, u32), QuadratureError>
182where
183    F: Fn(f64) -> f64,
184{
185    let center = 0.5 * (a + b);
186    let half = 0.5 * (b - a);
187    let center_value = checked_eval(function, center)?;
188    let mut kronrod = GK15_WEIGHTS[7] * center_value;
189    let mut gauss = G7_WEIGHTS[3] * center_value;
190
191    for index in 0..7 {
192        let offset = half * GK15_ABSCISSAE[index];
193        let pair =
194            checked_eval(function, center - offset)? + checked_eval(function, center + offset)?;
195        kronrod += GK15_WEIGHTS[index] * pair;
196        match index {
197            1 => gauss += G7_WEIGHTS[0] * pair,
198            3 => gauss += G7_WEIGHTS[1] * pair,
199            5 => gauss += G7_WEIGHTS[2] * pair,
200            _ => {}
201        }
202    }
203    let kronrod = kronrod * half;
204    let gauss = gauss * half;
205    Ok((kronrod, (kronrod - gauss).abs(), 15))
206}
207
208pub fn adaptive_gauss_kronrod_15<F>(
209    function: F,
210    a: f64,
211    b: f64,
212    absolute_tolerance: f64,
213    max_evaluations: u32,
214) -> Result<QuadratureResult, QuadratureError>
215where
216    F: Fn(f64) -> f64,
217{
218    if !a.is_finite()
219        || !b.is_finite()
220        || b <= a
221        || !absolute_tolerance.is_finite()
222        || absolute_tolerance <= 0.0
223        || max_evaluations < 15
224    {
225        return Err(QuadratureError::InvalidDomain);
226    }
227
228    let mut stack = [IntervalPanel::default(); MAX_ADAPTIVE_PANELS];
229    stack[0] = IntervalPanel {
230        a,
231        b,
232        tolerance: absolute_tolerance,
233    };
234    let mut stack_len = 1;
235    let mut evaluations = 0;
236    let mut intervals = 0;
237    let mut value = 0.0;
238    let mut error = 0.0;
239
240    while stack_len > 0 {
241        stack_len -= 1;
242        let panel = stack[stack_len];
243        if evaluations + 15 > max_evaluations {
244            return Err(QuadratureError::EvaluationBudgetExceeded);
245        }
246        let (estimate, local_error, used) = gauss_kronrod_15_panel(&function, panel.a, panel.b)?;
247        evaluations += used;
248        if local_error <= panel.tolerance {
249            value += estimate;
250            error += local_error;
251            intervals += 1;
252            continue;
253        }
254        if stack_len + 2 > stack.len() {
255            return Err(QuadratureError::WorkspaceExceeded);
256        }
257        let midpoint = 0.5 * (panel.a + panel.b);
258        let tolerance = panel.tolerance * 0.5;
259        stack[stack_len] = IntervalPanel {
260            a: midpoint,
261            b: panel.b,
262            tolerance,
263        };
264        stack[stack_len + 1] = IntervalPanel {
265            a: panel.a,
266            b: midpoint,
267            tolerance,
268        };
269        stack_len += 2;
270    }
271
272    Ok(QuadratureResult {
273        value,
274        absolute_error: error,
275        evaluations,
276        intervals,
277    })
278}
279
280#[cfg(test)]
281mod tests {
282    use super::*;
283
284    #[test]
285    fn adaptive_simpson_meets_tolerance_and_reports_budget_failure() {
286        let result =
287            adaptive_simpson(|x| x.sin(), 0.0, core::f64::consts::PI, 1e-12, 10_000).unwrap();
288        assert!((result.value - 2.0).abs() < 1e-12);
289        assert!(result.absolute_error < 1e-12);
290        assert_eq!(
291            adaptive_simpson(|x| x.sin(), 0.0, 1.0, 1e-14, 3),
292            Err(QuadratureError::EvaluationBudgetExceeded)
293        );
294    }
295
296    #[test]
297    fn gauss_kronrod_integrates_polynomial_and_oscillation() {
298        let polynomial =
299            adaptive_gauss_kronrod_15(|x| x.powi(12), 0.0, 1.0, 1e-13, 10_000).unwrap();
300        assert!((polynomial.value - 1.0 / 13.0).abs() < 1e-13);
301
302        let oscillatory =
303            adaptive_gauss_kronrod_15(|x| (50.0 * x).sin(), 0.0, 1.0, 1e-10, 50_000).unwrap();
304        let expected = (1.0 - 50.0_f64.cos()) / 50.0;
305        assert!((oscillatory.value - expected).abs() < 1e-10);
306    }
307
308    #[test]
309    fn quadrature_rejects_reversed_and_non_finite_domains() {
310        assert_eq!(
311            adaptive_simpson(|x| x, 1.0, 0.0, 1e-6, 100),
312            Err(QuadratureError::InvalidDomain)
313        );
314        assert!(matches!(
315            adaptive_simpson(|x| if x > 0.4 { f64::NAN } else { x }, 0.0, 1.0, 1e-6, 100),
316            Err(QuadratureError::NonFiniteIntegrand { .. })
317        ));
318    }
319}
qualia_core_db/solvers/calculus/quadrature.rs

qualia_core_db/solvers/calculus/
quadrature.rs
