peroxide/numerical/
integral.rs

1use crate::structure::polynomial::{lagrange_polynomial, Calculus, Polynomial};
2use crate::traits::fp::FPVector;
3use crate::util::non_macro::seq;
4
5#[derive(Debug, Copy, Clone, PartialEq)]
6pub enum Integral {
7    GaussLegendre(usize),
8    NewtonCotes(usize),
9    G7K15(f64, u32),
10    G10K21(f64, u32),
11    G15K31(f64, u32),
12    G20K41(f64, u32),
13    G25K51(f64, u32),
14    G30K61(f64, u32),
15    G7K15R(f64, u32),
16    G10K21R(f64, u32),
17    G15K31R(f64, u32),
18    G20K41R(f64, u32),
19    G25K51R(f64, u32),
20    G30K61R(f64, u32),
21}
22
23impl Integral {
24    pub fn get_num_node(&self) -> usize {
25        match self {
26            Integral::GaussLegendre(n) => *n,
27            Integral::NewtonCotes(n) => *n,
28            _ => panic!("This method does not have a fixed number of nodes."),
29        }
30    }
31
32    pub fn get_tol(&self) -> f64 {
33        match self {
34            Integral::G7K15(tol, _) => *tol,
35            Integral::G10K21(tol, _) => *tol,
36            Integral::G15K31(tol, _) => *tol,
37            Integral::G20K41(tol, _) => *tol,
38            Integral::G25K51(tol, _) => *tol,
39            Integral::G30K61(tol, _) => *tol,
40            Integral::G7K15R(tol, _) => *tol,
41            Integral::G10K21R(tol, _) => *tol,
42            Integral::G15K31R(tol, _) => *tol,
43            Integral::G20K41R(tol, _) => *tol,
44            Integral::G25K51R(tol, _) => *tol,
45            Integral::G30K61R(tol, _) => *tol,
46            _ => panic!("This method does not have a tolerance."),
47        }
48    }
49
50    pub fn get_max_iter(&self) -> u32 {
51        match self {
52            Integral::G7K15(_, max_iter) => *max_iter,
53            Integral::G10K21(_, max_iter) => *max_iter,
54            Integral::G15K31(_, max_iter) => *max_iter,
55            Integral::G20K41(_, max_iter) => *max_iter,
56            Integral::G25K51(_, max_iter) => *max_iter,
57            Integral::G30K61(_, max_iter) => *max_iter,
58            Integral::G7K15R(_, max_iter) => *max_iter,
59            Integral::G10K21R(_, max_iter) => *max_iter,
60            Integral::G15K31R(_, max_iter) => *max_iter,
61            Integral::G20K41R(_, max_iter) => *max_iter,
62            Integral::G25K51R(_, max_iter) => *max_iter,
63            Integral::G30K61R(_, max_iter) => *max_iter,
64            _ => panic!("This method does not have a maximum number of iterations."),
65        }
66    }
67
68    pub fn get_gauss_kronrod_order(&self) -> (u8, u8) {
69        match self {
70            Integral::G7K15(_, _) => (7, 15),
71            Integral::G10K21(_, _) => (10, 21),
72            Integral::G15K31(_, _) => (15, 31),
73            Integral::G20K41(_, _) => (20, 41),
74            Integral::G25K51(_, _) => (25, 51),
75            Integral::G30K61(_, _) => (30, 61),
76            Integral::G7K15R(_, _) => (7, 15),
77            Integral::G10K21R(_, _) => (10, 21),
78            Integral::G15K31R(_, _) => (15, 31),
79            Integral::G20K41R(_, _) => (20, 41),
80            Integral::G25K51R(_, _) => (25, 51),
81            Integral::G30K61R(_, _) => (30, 61),
82            _ => panic!("This method does not have a Gauss-Kronrod order."),
83        }
84    }
85
86    pub fn is_relative(&self) -> bool {
87        match self {
88            Integral::G7K15R(_, _) => true,
89            Integral::G10K21R(_, _) => true,
90            Integral::G15K31R(_, _) => true,
91            Integral::G20K41R(_, _) => true,
92            Integral::G25K51R(_, _) => true,
93            Integral::G30K61R(_, _) => true,
94            _ => false,
95        }
96    }
97
98    pub fn change_tol(&self, tol: f64) -> Self {
99        match self {
100            Integral::G7K15(_, max_iter) => Integral::G7K15(tol, *max_iter),
101            Integral::G10K21(_, max_iter) => Integral::G10K21(tol, *max_iter),
102            Integral::G15K31(_, max_iter) => Integral::G15K31(tol, *max_iter),
103            Integral::G20K41(_, max_iter) => Integral::G20K41(tol, *max_iter),
104            Integral::G25K51(_, max_iter) => Integral::G25K51(tol, *max_iter),
105            Integral::G30K61(_, max_iter) => Integral::G30K61(tol, *max_iter),
106            Integral::G7K15R(_, max_iter) => Integral::G7K15R(tol, *max_iter),
107            Integral::G10K21R(_, max_iter) => Integral::G10K21R(tol, *max_iter),
108            Integral::G15K31R(_, max_iter) => Integral::G15K31R(tol, *max_iter),
109            Integral::G20K41R(_, max_iter) => Integral::G20K41R(tol, *max_iter),
110            Integral::G25K51R(_, max_iter) => Integral::G25K51R(tol, *max_iter),
111            Integral::G30K61R(_, max_iter) => Integral::G30K61R(tol, *max_iter),
112            _ => panic!("This method does not have a tolerance."),
113        }
114    }
115
116    pub fn change_max_iter(&self, max_iter: u32) -> Self {
117        match self {
118            Integral::G7K15(tol, _) => Integral::G7K15(*tol, max_iter),
119            Integral::G10K21(tol, _) => Integral::G10K21(*tol, max_iter),
120            Integral::G15K31(tol, _) => Integral::G15K31(*tol, max_iter),
121            Integral::G20K41(tol, _) => Integral::G20K41(*tol, max_iter),
122            Integral::G25K51(tol, _) => Integral::G25K51(*tol, max_iter),
123            Integral::G30K61(tol, _) => Integral::G30K61(*tol, max_iter),
124            Integral::G7K15R(tol, _) => Integral::G7K15R(*tol, max_iter),
125            Integral::G10K21R(tol, _) => Integral::G10K21R(*tol, max_iter),
126            Integral::G15K31R(tol, _) => Integral::G15K31R(*tol, max_iter),
127            Integral::G20K41R(tol, _) => Integral::G20K41R(*tol, max_iter),
128            Integral::G25K51R(tol, _) => Integral::G25K51R(*tol, max_iter),
129            Integral::G30K61R(tol, _) => Integral::G30K61R(*tol, max_iter),
130            _ => panic!("This method does not have a maximum number of iterations."),
131        }
132    }
133}
134
135/// Type that can be integrated using Newton Cotes Quadrature
136pub trait NCIntegrable: std::ops::Sub<Self, Output = Self> + Sized {
137    type NodeY;
138    type NCPolynomial;
139
140    /// Returns the image of `node_x` under function `f`, in a representation
141    /// (`Self::NodeY`) suitable for separately computing one Lagrange polynomial
142    /// for each of the degrees of freedom of `Self` (e.g. a `Vec` with as many
143    /// entries as the number of degrees of freedom).
144    fn compute_node_y<F>(f: F, node_x: &[f64]) -> Self::NodeY
145    where
146        F: Fn(f64) -> Self;
147
148    /// Computes one [`Lagrange polynomial`](lagrange_polynomial) for each one
149    /// of the degrees of freedom of `Self`, returning them in a representation
150    /// (`Self::NCPolynomial`) suitable to separately performing operation on
151    /// them (e.g. [`Vec<Polynomial>`](Polynomial)).
152    fn compute_polynomial(node_x: &[f64], node_y: &Self::NodeY) -> Self::NCPolynomial;
153
154    /// Separately integrates each of the polynomials obtained using
155    /// [`compute_polynomial`](NCIntegrable::compute_polynomial).
156    fn integrate_polynomial(p: &Self::NCPolynomial) -> Self::NCPolynomial;
157
158    /// Separately evaluates each of the polynomial integrated using
159    /// [`integrate_polynomial`](NCIntegrable::integrate_polynomial) at `x`,
160    /// then recombines the result into a `Self`.
161    fn evaluate_polynomial(p: &Self::NCPolynomial, x: f64) -> Self;
162}
163
164impl NCIntegrable for f64 {
165    type NodeY = Vec<f64>;
166    type NCPolynomial = Polynomial;
167
168    fn compute_node_y<F>(f: F, node_x: &[f64]) -> Vec<f64>
169    where
170        F: Fn(f64) -> Self,
171    {
172        node_x.to_vec().fmap(|x| f(x))
173    }
174
175    fn compute_polynomial(node_x: &[f64], node_y: &Vec<f64>) -> Polynomial {
176        lagrange_polynomial(node_x.to_vec(), node_y.to_vec())
177    }
178
179    fn integrate_polynomial(p: &Polynomial) -> Polynomial {
180        p.integral()
181    }
182
183    fn evaluate_polynomial(p: &Polynomial, x: f64) -> Self {
184        p.eval(x)
185    }
186}
187
188/// Type that can be integrated using Gauss Legendre or Kronrod Quadrature
189pub trait GLKIntegrable:
190    std::ops::Add<Self, Output = Self> + std::ops::Mul<f64, Output = Self> + Sized
191{
192    /// The neutral element of addition.
193    const ZERO: Self;
194}
195
196impl GLKIntegrable for f64 {
197    const ZERO: Self = 0f64;
198}
199
200/// Type that can be integrated using Gauss Kronrod Quadrature
201pub trait GKIntegrable: GLKIntegrable + std::ops::Sub<Self, Output = Self> + Sized + Clone {
202    /// Returns `true` if this object is neither infinite nor NaN.
203    fn is_finite(&self) -> bool;
204
205    /// Returns the value of the norm used by Gauss Kronrod Quadrature to
206    /// compute relative tolerances.
207    fn gk_norm(&self) -> f64;
208
209    /// Returns true if this object is less than the tolerance passed as
210    /// the `tol` argument. By default, returns true if its [`gk_norm`](GKIntegrable::gk_norm)
211    /// is less than `tol`.
212    fn is_within_tol(&self, tol: f64) -> bool {
213        self.gk_norm() < tol
214    }
215}
216
217impl GKIntegrable for f64 {
218    fn is_finite(&self) -> bool {
219        f64::is_finite(*self)
220    }
221
222    fn gk_norm(&self) -> f64 {
223        self.abs()
224    }
225}
226
227/// Numerical Integration
228///
229/// # Description
230/// `fn integrate(f, (a,b), method) -> Y`
231///
232/// `Y` must implement [`GKIntegrable`] and [`NCIntegrable`], like `f64` or
233/// `C64` (`complex` feature only).
234///
235/// * `f`: Target function (`Fn(f64) -> Y`)
236/// * `(a,b)` : Target interval
237/// * `method` : Numerical integration method
238///
239/// # Method
240///
241/// * Gauss-Legendre Quadrature (up to order 30) : `GaussLegendre(usize)`
242/// * Newton-Cotes Quadrature: `NewtonCotes(usize)`
243/// * Gauss-Kronrod Quadrature
244///     * `G7K15(tol, max_iter)`
245///     * `G10K21`
246///     * `G15K31`
247///     * `G20K41`
248///     * `G25K51`
249///     * `G30K61`
250/// * Gauss-Kronrod Quarature with Relative Error
251///     * `G7K15R(rtol, max_iter)`
252///     * `G10K21R`
253///     * `G15K31R`
254///     * `G20K41R`
255///     * `G25K51R`
256///     * `G30K61R`
257pub fn integrate<F, Y>(f: F, (a, b): (f64, f64), method: Integral) -> Y
258where
259    F: Fn(f64) -> Y + Copy,
260    Y: GKIntegrable + NCIntegrable,
261{
262    match method {
263        Integral::GaussLegendre(n) => gauss_legendre_quadrature(f, n, (a, b)),
264        Integral::NewtonCotes(n) => newton_cotes_quadrature(f, n, (a, b)),
265        method => gauss_kronrod_quadrature_optimized(f, (a, b), method),
266    }
267}
268
269/// Newton Cotes Quadrature
270pub fn newton_cotes_quadrature<F, Y>(f: F, n: usize, (a, b): (f64, f64)) -> Y
271where
272    F: Fn(f64) -> Y,
273    Y: NCIntegrable,
274{
275    let h = (b - a) / (n as f64);
276    let node_x = seq(a, b, h);
277    let node_y = Y::compute_node_y(f, &node_x);
278    let p = Y::compute_polynomial(&node_x, &node_y);
279    let q = Y::integrate_polynomial(&p);
280
281    Y::evaluate_polynomial(&q, b) - Y::evaluate_polynomial(&q, a)
282}
283
284/// Gauss Legendre Quadrature
285///
286/// # Type
287/// * `f, n, (a,b) -> Y`
288///     * `Y` implements [`GLKIntegrable`], like `f64` or `C64` (`complex` feature
289///       only)
290///     * `f`: Numerical function (`Fn(f64) -> Y`)
291///     * `n`: Order of Legendre polynomial (up to 16)
292///     * `(a,b)`: Interval of integration
293///
294/// # Reference
295/// * A. N. Lowan et al. (1942)
296pub fn gauss_legendre_quadrature<F, Y>(f: F, n: usize, (a, b): (f64, f64)) -> Y
297where
298    F: Fn(f64) -> Y,
299    Y: GLKIntegrable,
300{
301    unit_gauss_legendre_quadrature(|x| f(x * (b - a) / 2f64 + (a + b) / 2f64), n) * ((b - a) / 2f64)
302}
303
304/// Gauss Kronrod Quadrature
305///
306/// # Type
307/// * `f, (a,b), method -> Y`
308///     * `Y` implements [`GKIntegrable`], like `f64` or `C64` (`complex` feature
309///       only)
310///     * `f`: Numerical function (`Fn(f64) -> Y + Copy`)
311///     * `(a,b)`: Interval of integration
312///     * `method`: Integration method
313///         * `G7K15(tol)`
314///
315/// # Reference
316/// * arXiv: [1003.4629](https://arxiv.org/abs/1003.4629)
317#[allow(non_snake_case)]
318pub fn gauss_kronrod_quadrature<F, T, S, Y>(f: F, (a, b): (T, S), method: Integral) -> Y
319where
320    F: Fn(f64) -> Y + Copy,
321    T: Into<f64>,
322    S: Into<f64>,
323    Y: GKIntegrable,
324{
325    let (g, k) = method.get_gauss_kronrod_order();
326    let tol = method.get_tol();
327    let max_iter = method.get_max_iter();
328    let mut I = Y::ZERO;
329    let mut S: Vec<(f64, f64, f64, u32)> = vec![];
330    S.push((a.into(), b.into(), tol, max_iter));
331
332    loop {
333        match S.pop() {
334            Some((a, b, tol, max_iter)) => {
335                let G = gauss_legendre_quadrature(f, g as usize, (a, b));
336                let K = kronrod_quadrature(f, k as usize, (a, b));
337                let c = (a + b) / 2f64;
338                let tol_curr = if method.is_relative() {
339                    tol * G.gk_norm()
340                } else {
341                    tol
342                };
343                if (G.clone() - K).is_within_tol(tol_curr) || a == b || max_iter == 0 {
344                    if !G.is_finite() {
345                        return G;
346                    }
347                    I = I + G;
348                } else {
349                    S.push((a, c, tol / 2f64, max_iter - 1));
350                    S.push((c, b, tol / 2f64, max_iter - 1));
351                }
352            }
353            None => break,
354        }
355    }
356    I
357}
358
359pub fn kronrod_quadrature<F, Y>(f: F, n: usize, (a, b): (f64, f64)) -> Y
360where
361    F: Fn(f64) -> Y,
362    Y: GLKIntegrable,
363{
364    unit_kronrod_quadrature(|x| f(x * (b - a) / 2f64 + (a + b) / 2f64), n) * ((b - a) / 2f64)
365}
366
367#[allow(non_snake_case)]
368pub fn gauss_kronrod_quadrature_optimized<F, T, S, Y>(f: F, (a, b): (T, S), method: Integral) -> Y
369where
370    F: Fn(f64) -> Y,
371    T: Into<f64>,
372    S: Into<f64>,
373    Y: GKIntegrable, // Requires Clone, Add, Sub, Mul<f64>, ZERO, gk_norm, is_finite
374{
375    // Extract parameters from the method enum
376    let (g, k) = method.get_gauss_kronrod_order();
377    let tol = method.get_tol(); // Base tolerance
378    let max_iter = method.get_max_iter(); // Max recursion depth
379    let is_relative = method.is_relative(); // Use relative tolerance?
380
381    let mut I = Y::ZERO; // Accumulator for the integral result
382    let mut S: Vec<(f64, f64, f64, u32)> = Vec::with_capacity(max_iter as usize + 10); // Stack for subintervals: (a, b, tolerance, depth)
383
384    let a_f64 = a.into();
385    let b_f64 = b.into();
386
387    // Handle zero-width interval initially
388    if a_f64 == b_f64 {
389        return Y::ZERO;
390    }
391
392    // Push the initial interval onto the stack
393    S.push((a_f64, b_f64, tol, max_iter));
394
395    // Adaptive quadrature loop
396    while let Some((a_curr, b_curr, tol_curr, iter_left)) = S.pop() {
397        // Call the helper function to get G, K, and error estimate for the current interval
398        // This performs all necessary function evaluations efficiently
399        let (G, K, error_estimate_raw) =
400            compute_gauss_kronrod_sums_stored(&f, (a_curr, b_curr), g as usize, k as usize);
401
402        // Calculate the norm of the error estimate |G - K|
403        let error_norm = error_estimate_raw.gk_norm();
404
405        // Determine the tolerance threshold for this interval
406        let current_tolerance_threshold = if is_relative {
407            // Scale relative tolerance by the norm of the higher-order estimate (K)
408            tol_curr * K.clone().gk_norm()
409        } else {
410            // Use absolute tolerance directly
411            tol_curr
412        };
413
414        // Check termination conditions for this interval:
415        // 1. Error is within tolerance
416        // 2. Interval width is zero (or numerically close)
417        // 3. Maximum iteration depth reached
418        if error_norm <= current_tolerance_threshold || a_curr == b_curr || iter_left == 0 {
419            // Interval meets criteria, add its contribution to the total integral I
420            // Prefer the higher-order Kronrod estimate (K) if it's finite.
421            if !K.is_finite() {
422                // If K is not finite, check G
423                if G.is_finite() {
424                    // If K is not finite but G is, add G
425                    I = I + G.clone();
426                }
427            } else {
428                // K is finite, add it to the total integral
429                I = I + K.clone();
430            }
431        } else {
432            // Interval does not meet criteria, subdivide it
433            let c = (a_curr + b_curr) * 0.5; // Midpoint
434
435            // Avoid infinite loops if midpoint calculation doesn't change due to precision
436            if c <= a_curr || c >= b_curr {
437                // Cannot subdivide further, precision limit reached.
438                // Add the best estimate (K if finite, else G if finite) for this interval and stop subdividing.
439                if !K.is_finite() {
440                    if G.is_finite() {
441                        I = I + G.clone();
442                    }
443                } else {
444                    I = I + K.clone();
445                }
446                continue; // Skip pushing to stack
447            }
448
449            // Tolerance propagation for subintervals: divide by sqrt(2)
450            // This is often preferred over dividing by 2 for better error distribution.
451            let next_tol = tol_curr / 1.4142135623730951; // tol / sqrt(2)
452            let next_iter = iter_left - 1; // Decrement remaining iterations
453
454            // Push the two subintervals onto the stack for further processing
455            S.push((a_curr, c, next_tol, next_iter)); // Left subinterval [a, c]
456            S.push((c, b_curr, next_tol, next_iter)); // Right subinterval [c, b]
457        }
458    }
459
460    I
461}
462
463/// Helper function to compute G and K integrals reusing function evaluations.
464/// Optimized to store all function evaluations.
465/// Returns (Gauss Integral G, Kronrod Integral K, Error Estimate G - K).
466#[allow(non_snake_case)]
467fn compute_gauss_kronrod_sums_stored<F, Y>(
468    f: F,
469    (a, b): (f64, f64),
470    g: usize, // Gauss order
471    k: usize, // Kronrod order
472) -> (Y, Y, Y)
473where
474    F: Fn(f64) -> Y,
475    Y: GKIntegrable,
476{
477    // 1. Get Kronrod nodes and weights
478    // Assumes nodes ordered [-xn,...,0,...,+xn] and weights correspond.
479    let (kronrod_weights, kronrod_nodes) = kronrod_table(k);
480    // 2. Get Gauss weights
481    // Assumes weights ordered [-wgn,...,wg0,...,+wgn] and symmetric w(+x) == w(-x).
482    let (gauss_weights, _gauss_nodes_unused) = gauss_legendre_table(g);
483
484    // 3. Calculate interval midpoint and half-length
485    let xm = (a + b) * 0.5;
486    let xh = (b - a) * 0.5;
487
488    // Handle zero interval length
489    if xh == 0.0 {
490        let center_f = f(xm);
491        return (center_f.clone(), center_f, Y::ZERO);
492    }
493
494    // 4. Evaluate f at all Kronrod nodes and store results (matching node order)
495    let mut f_evals = Vec::with_capacity(k);
496    for i in 0..k {
497        let node = kronrod_nodes[i];
498        f_evals.push(f(xm + xh * node));
499    }
500
501    // 5. Calculate Kronrod sum (K)
502    let mut kronrod_sum = Y::ZERO;
503    for i in 0..k {
504        kronrod_sum = kronrod_sum + f_evals[i].clone() * kronrod_weights[i];
505    }
506    let K = kronrod_sum * xh;
507
508    // 6. Calculate Gauss sum (G) using stored evaluations and Gauss weights
509    // Assumes G(g) nodes are subset of K(k), g is odd, standard mapping: +/- xg_j maps to +/- xk_{2j}.
510    let kronrod_center_idx = k / 2; // Index of 0.0 node in Kronrod/f_evals
511    let gauss_center_idx = g / 2; // Index of 0.0 weight in Gauss weights
512    let num_gauss_pairs = (g - 1) / 2;
513
514    // Contribution from center node (0.0)
515    let mut gauss_sum = f_evals[kronrod_center_idx].clone() * gauss_weights[gauss_center_idx];
516
517    // Contribution from paired Gauss nodes (+/- xg_j)
518    for j in 1..=num_gauss_pairs {
519        if 2 * j > kronrod_center_idx {
520            panic!(
521                "Kronrod node index underflow. k={}, center_idx={}, j={}",
522                k, kronrod_center_idx, j
523            );
524        }
525
526        // Indices in f_evals corresponding to Kronrod nodes +/- xk_{2j}
527        let kronrod_node_pos_idx = kronrod_center_idx + 2 * j;
528        let kronrod_node_neg_idx = kronrod_center_idx - 2 * j;
529
530        // Index in gauss_weights for the symmetric weight of the pair +/- xg_j
531        let gauss_weight_pair_idx = gauss_center_idx + j;
532
533        // Safety checks
534        if kronrod_node_pos_idx >= k || kronrod_node_neg_idx >= k {
535            panic!("Kronrod node index for Gauss sum out of bounds. Check mapping/indices. k={}, pos_idx={}, neg_idx={}", k, kronrod_node_pos_idx, kronrod_node_neg_idx);
536        }
537        if gauss_weight_pair_idx >= g {
538            panic!("Gauss weight index out of bounds. Check gauss_weights structure. g={}, weight_idx={}", g, gauss_weight_pair_idx);
539        }
540
541        // Retrieve stored function evaluations
542        let f_pos = f_evals[kronrod_node_pos_idx].clone();
543        let f_neg = f_evals[kronrod_node_neg_idx].clone();
544
545        // Retrieve the Gauss weight for the pair
546        let weight_g_pair = gauss_weights[gauss_weight_pair_idx];
547
548        // Add contribution: (f(+x) + f(-x)) * w_pair
549        let f_sum_pair = f_pos + f_neg;
550        gauss_sum = gauss_sum + f_sum_pair * weight_g_pair;
551    }
552
553    let G = gauss_sum * xh;
554
555    // 7. Calculate Error Estimate (G - K)
556    let error_estimate = G.clone() - K.clone();
557
558    (G, K, error_estimate)
559}
560
561// =============================================================================
562// Gauss Legendre Backends
563// =============================================================================
564fn unit_gauss_legendre_quadrature<F, Y>(f: F, n: usize) -> Y
565where
566    F: Fn(f64) -> Y,
567    Y: GLKIntegrable,
568{
569    let (a, x) = gauss_legendre_table(n);
570    let mut s = Y::ZERO;
571    for i in 0..a.len() {
572        s = s + f(x[i]) * a[i];
573    }
574    s
575}
576
577fn gauss_legendre_table(n: usize) -> (&'static [f64], &'static [f64]) {
578    let root: &'static [f64] = match n {
579        2 => &LEGENDRE_ROOT_2[..],
580        3 => &LEGENDRE_ROOT_3[..],
581        4 => &LEGENDRE_ROOT_4[..],
582        5 => &LEGENDRE_ROOT_5[..],
583        6 => &LEGENDRE_ROOT_6[..],
584        7 => &LEGENDRE_ROOT_7[..],
585        8 => &LEGENDRE_ROOT_8[..],
586        9 => &LEGENDRE_ROOT_9[..],
587        10 => &LEGENDRE_ROOT_10[..],
588        11 => &LEGENDRE_ROOT_11[..],
589        12 => &LEGENDRE_ROOT_12[..],
590        13 => &LEGENDRE_ROOT_13[..],
591        14 => &LEGENDRE_ROOT_14[..],
592        15 => &LEGENDRE_ROOT_15[..],
593        16 => &LEGENDRE_ROOT_16[..],
594        17 => &LEGENDRE_ROOT_17[..],
595        18 => &LEGENDRE_ROOT_18[..],
596        19 => &LEGENDRE_ROOT_19[..],
597        20 => &LEGENDRE_ROOT_20[..],
598        21 => &LEGENDRE_ROOT_21[..],
599        22 => &LEGENDRE_ROOT_22[..],
600        23 => &LEGENDRE_ROOT_23[..],
601        24 => &LEGENDRE_ROOT_24[..],
602        25 => &LEGENDRE_ROOT_25[..],
603        26 => &LEGENDRE_ROOT_26[..],
604        27 => &LEGENDRE_ROOT_27[..],
605        28 => &LEGENDRE_ROOT_28[..],
606        29 => &LEGENDRE_ROOT_29[..],
607        30 => &LEGENDRE_ROOT_30[..],
608        _ => panic!("Legendre quadrature is limited up to n = 16"),
609    };
610
611    let weight: &'static [f64] = match n {
612        2 => &LEGENDRE_WEIGHT_2[..],
613        3 => &LEGENDRE_WEIGHT_3[..],
614        4 => &LEGENDRE_WEIGHT_4[..],
615        5 => &LEGENDRE_WEIGHT_5[..],
616        6 => &LEGENDRE_WEIGHT_6[..],
617        7 => &LEGENDRE_WEIGHT_7[..],
618        8 => &LEGENDRE_WEIGHT_8[..],
619        9 => &LEGENDRE_WEIGHT_9[..],
620        10 => &LEGENDRE_WEIGHT_10[..],
621        11 => &LEGENDRE_WEIGHT_11[..],
622        12 => &LEGENDRE_WEIGHT_12[..],
623        13 => &LEGENDRE_WEIGHT_13[..],
624        14 => &LEGENDRE_WEIGHT_14[..],
625        15 => &LEGENDRE_WEIGHT_15[..],
626        16 => &LEGENDRE_WEIGHT_16[..],
627        17 => &LEGENDRE_WEIGHT_17[..],
628        18 => &LEGENDRE_WEIGHT_18[..],
629        19 => &LEGENDRE_WEIGHT_19[..],
630        20 => &LEGENDRE_WEIGHT_20[..],
631        21 => &LEGENDRE_WEIGHT_21[..],
632        22 => &LEGENDRE_WEIGHT_22[..],
633        23 => &LEGENDRE_WEIGHT_23[..],
634        24 => &LEGENDRE_WEIGHT_24[..],
635        25 => &LEGENDRE_WEIGHT_25[..],
636        26 => &LEGENDRE_WEIGHT_26[..],
637        27 => &LEGENDRE_WEIGHT_27[..],
638        28 => &LEGENDRE_WEIGHT_28[..],
639        29 => &LEGENDRE_WEIGHT_29[..],
640        30 => &LEGENDRE_WEIGHT_30[..],
641        _ => panic!("Legendre quadrature is limited up to n = 16"),
642    };
643
644    (weight, root)
645}
646
647// =============================================================================
648// Gauss Kronrod Backends
649// =============================================================================
650
651fn unit_kronrod_quadrature<F, Y>(f: F, n: usize) -> Y
652where
653    F: Fn(f64) -> Y,
654    Y: GLKIntegrable,
655{
656    let (a, x) = kronrod_table(n);
657    let mut s = Y::ZERO;
658    for i in 0..a.len() {
659        s = s + f(x[i]) * a[i];
660    }
661    s
662}
663
664fn kronrod_table(n: usize) -> (&'static [f64], &'static [f64]) {
665    let node: &'static [f64] = match n {
666        15 => &KRONROD_NODES_15[..],
667        21 => &KRONROD_NODES_21[..],
668        31 => &KRONROD_NODES_31[..],
669        41 => &KRONROD_NODES_41[..],
670        51 => &KRONROD_NODES_51[..],
671        61 => &KRONROD_NODES_61[..],
672        _ => panic!("Not yet implemented"),
673    };
674    let weight: &'static [f64] = match n {
675        15 => &KRONROD_WEIGHTS_15[..],
676        21 => &KRONROD_WEIGHTS_21[..],
677        31 => &KRONROD_WEIGHTS_31[..],
678        41 => &KRONROD_WEIGHTS_41[..],
679        51 => &KRONROD_WEIGHTS_51[..],
680        61 => &KRONROD_WEIGHTS_61[..],
681        _ => panic!("Not yet implemented"),
682    };
683
684    (weight, node)
685}
686
687// =============================================================================
688// Table for Gauss-Legendre Quadrature (ref. A. N. Lowan et al. (1942))
689// =============================================================================
690const LEGENDRE_ROOT_2: [f64; 2] = [-0.577350269189626, 0.577350269189626];
691const LEGENDRE_ROOT_3: [f64; 3] = [-0.774596669241483, 0.0, 0.774596669241483];
692const LEGENDRE_ROOT_4: [f64; 4] = [
693    -0.861136311594053,
694    -0.339981043584856,
695    0.339981043584856,
696    0.861136311594053,
697];
698const LEGENDRE_ROOT_5: [f64; 5] = [
699    -0.906179845938664,
700    -0.538469310105683,
701    0.0,
702    0.538469310105683,
703    0.906179845938664,
704];
705const LEGENDRE_ROOT_6: [f64; 6] = [
706    -0.932469514203152,
707    -0.661209386466265,
708    -0.238619186083197,
709    0.238619186083197,
710    0.661209386466265,
711    0.932469514203152,
712];
713const LEGENDRE_ROOT_7: [f64; 7] = [
714    -0.949107912342759,
715    -0.741531185599394,
716    -0.405845151377397,
717    0.0,
718    0.405845151377397,
719    0.741531185599394,
720    0.949107912342759,
721];
722
723const LEGENDRE_ROOT_8: [f64; 8] = [
724    -0.960289856497536,
725    -0.796666477413627,
726    -0.525532409916329,
727    -0.18343464249565,
728    0.18343464249565,
729    0.525532409916329,
730    0.796666477413627,
731    0.960289856497536,
732];
733const LEGENDRE_ROOT_9: [f64; 9] = [
734    -0.968160239507626,
735    -0.836031107326636,
736    -0.61337143270059,
737    -0.324253423403809,
738    0.0,
739    0.324253423403809,
740    0.61337143270059,
741    0.836031107326636,
742    0.968160239507626,
743];
744const LEGENDRE_ROOT_10: [f64; 10] = [
745    -0.973906528517172,
746    -0.865063366688985,
747    -0.679409568299024,
748    -0.433395394129247,
749    -0.148874338981631,
750    0.148874338981631,
751    0.433395394129247,
752    0.679409568299024,
753    0.865063366688985,
754    0.973906528517172,
755];
756const LEGENDRE_ROOT_11: [f64; 11] = [
757    -0.978228658146057,
758    -0.887062599768095,
759    -0.730152005574049,
760    -0.519096129110681,
761    -0.269543155952345,
762    0.0,
763    0.269543155952345,
764    0.519096129110681,
765    0.730152005574049,
766    0.887062599768095,
767    0.978228658146057,
768];
769const LEGENDRE_ROOT_12: [f64; 12] = [
770    -0.981560634246719,
771    -0.904117256370475,
772    -0.769902674194305,
773    -0.587317954286617,
774    -0.36783149891818,
775    -0.125333408511469,
776    0.125333408511469,
777    0.36783149891818,
778    0.587317954286617,
779    0.769902674194305,
780    0.904117256370475,
781    0.981560634246719,
782];
783const LEGENDRE_ROOT_13: [f64; 13] = [
784    -0.984183054718588,
785    -0.917598399222978,
786    -0.80157809073331,
787    -0.64234933944034,
788    -0.448492751036447,
789    -0.230458315955135,
790    0.0,
791    0.230458315955135,
792    0.448492751036447,
793    0.64234933944034,
794    0.80157809073331,
795    0.917598399222978,
796    0.984183054718588,
797];
798const LEGENDRE_ROOT_14: [f64; 14] = [
799    -0.986283808696812,
800    -0.928434883663574,
801    -0.827201315069765,
802    -0.687292904811685,
803    -0.515248636358154,
804    -0.31911236892789,
805    -0.108054948707344,
806    0.108054948707344,
807    0.31911236892789,
808    0.515248636358154,
809    0.687292904811685,
810    0.827201315069765,
811    0.928434883663574,
812    0.986283808696812,
813];
814const LEGENDRE_ROOT_15: [f64; 15] = [
815    -0.987992518020485,
816    -0.937273392400706,
817    -0.848206583410427,
818    -0.72441773136017,
819    -0.570972172608539,
820    -0.394151347077563,
821    -0.201194093997435,
822    0.0,
823    0.201194093997435,
824    0.394151347077563,
825    0.570972172608539,
826    0.72441773136017,
827    0.848206583410427,
828    0.937273392400706,
829    0.987992518020485,
830];
831const LEGENDRE_ROOT_16: [f64; 16] = [
832    -0.98940093499165,
833    -0.944575023073233,
834    -0.865631202387832,
835    -0.755404408355003,
836    -0.617876244402644,
837    -0.458016777657227,
838    -0.281603550779259,
839    -0.095012509837637,
840    0.095012509837637,
841    0.281603550779259,
842    0.458016777657227,
843    0.617876244402644,
844    0.755404408355003,
845    0.865631202387832,
846    0.944575023073233,
847    0.98940093499165,
848];
849const LEGENDRE_ROOT_17: [f64; 17] = [
850    -0.9905754753144174,
851    -0.9506755217687678,
852    -0.8802391537269859,
853    -0.7815140038968014,
854    -0.6576711592166907,
855    -0.5126905370864769,
856    -0.3512317634538763,
857    -0.17848418149584785,
858    0.0,
859    0.17848418149584785,
860    0.3512317634538763,
861    0.5126905370864769,
862    0.6576711592166907,
863    0.7815140038968014,
864    0.8802391537269859,
865    0.9506755217687678,
866    0.9905754753144174,
867];
868const LEGENDRE_ROOT_18: [f64; 18] = [
869    -0.9915651684209309,
870    -0.9558239495713977,
871    -0.8926024664975557,
872    -0.8037049589725231,
873    -0.6916870430603532,
874    -0.5597708310739475,
875    -0.41175116146284263,
876    -0.2518862256915055,
877    -0.0847750130417353,
878    0.0847750130417353,
879    0.2518862256915055,
880    0.41175116146284263,
881    0.5597708310739475,
882    0.6916870430603532,
883    0.8037049589725231,
884    0.8926024664975557,
885    0.9558239495713977,
886    0.9915651684209309,
887];
888const LEGENDRE_ROOT_19: [f64; 19] = [
889    -0.9924068438435844,
890    -0.96020815213483,
891    -0.9031559036148179,
892    -0.8227146565371428,
893    -0.7209661773352294,
894    -0.600545304661681,
895    -0.46457074137596094,
896    -0.31656409996362983,
897    -0.16035864564022537,
898    0.0,
899    0.16035864564022537,
900    0.31656409996362983,
901    0.46457074137596094,
902    0.600545304661681,
903    0.7209661773352294,
904    0.8227146565371428,
905    0.9031559036148179,
906    0.96020815213483,
907    0.9924068438435844,
908];
909const LEGENDRE_ROOT_20: [f64; 20] = [
910    -0.9931285991850949,
911    -0.9639719272779138,
912    -0.912234428251326,
913    -0.8391169718222188,
914    -0.7463319064601508,
915    -0.636053680726515,
916    -0.5108670019508271,
917    -0.37370608871541955,
918    -0.22778585114164507,
919    -0.07652652113349734,
920    0.07652652113349734,
921    0.22778585114164507,
922    0.37370608871541955,
923    0.5108670019508271,
924    0.636053680726515,
925    0.7463319064601508,
926    0.8391169718222188,
927    0.912234428251326,
928    0.9639719272779138,
929    0.9931285991850949,
930];
931const LEGENDRE_ROOT_21: [f64; 21] = [
932    -0.9937521706203895,
933    -0.9672268385663063,
934    -0.9200993341504008,
935    -0.8533633645833173,
936    -0.7684399634756779,
937    -0.6671388041974123,
938    -0.5516188358872198,
939    -0.4243421202074388,
940    -0.2880213168024011,
941    -0.1455618541608951,
942    0.0,
943    0.1455618541608951,
944    0.2880213168024011,
945    0.4243421202074388,
946    0.5516188358872198,
947    0.6671388041974123,
948    0.7684399634756779,
949    0.8533633645833173,
950    0.9200993341504008,
951    0.9672268385663063,
952    0.9937521706203895,
953];
954const LEGENDRE_ROOT_22: [f64; 22] = [
955    -0.9942945854823992,
956    -0.9700604978354287,
957    -0.926956772187174,
958    -0.8658125777203002,
959    -0.7878168059792081,
960    -0.6944872631866827,
961    -0.5876404035069116,
962    -0.469355837986757,
963    -0.34193582089208424,
964    -0.20786042668822127,
965    -0.06973927331972223,
966    0.06973927331972223,
967    0.20786042668822127,
968    0.34193582089208424,
969    0.469355837986757,
970    0.5876404035069116,
971    0.6944872631866827,
972    0.7878168059792081,
973    0.8658125777203002,
974    0.926956772187174,
975    0.9700604978354287,
976    0.9942945854823992,
977];
978const LEGENDRE_ROOT_23: [f64; 23] = [
979    -0.9947693349975522,
980    -0.9725424712181152,
981    -0.9329710868260161,
982    -0.8767523582704416,
983    -0.8048884016188399,
984    -0.7186613631319502,
985    -0.6196098757636461,
986    -0.5095014778460075,
987    -0.3903010380302908,
988    -0.26413568097034495,
989    -0.1332568242984661,
990    0.0,
991    0.1332568242984661,
992    0.26413568097034495,
993    0.3903010380302908,
994    0.5095014778460075,
995    0.6196098757636461,
996    0.7186613631319502,
997    0.8048884016188399,
998    0.8767523582704416,
999    0.9329710868260161,
1000    0.9725424712181152,
1001    0.9947693349975522,
1002];
1003const LEGENDRE_ROOT_24: [f64; 24] = [
1004    -0.9951872199970213,
1005    -0.9747285559713095,
1006    -0.9382745520027328,
1007    -0.8864155270044011,
1008    -0.820001985973903,
1009    -0.7401241915785544,
1010    -0.6480936519369755,
1011    -0.5454214713888396,
1012    -0.4337935076260451,
1013    -0.3150426796961634,
1014    -0.1911188674736163,
1015    -0.06405689286260563,
1016    0.06405689286260563,
1017    0.1911188674736163,
1018    0.3150426796961634,
1019    0.4337935076260451,
1020    0.5454214713888396,
1021    0.6480936519369755,
1022    0.7401241915785544,
1023    0.820001985973903,
1024    0.8864155270044011,
1025    0.9382745520027328,
1026    0.9747285559713095,
1027    0.9951872199970213,
1028];
1029const LEGENDRE_ROOT_25: [f64; 25] = [
1030    -0.9955569697904981,
1031    -0.9766639214595175,
1032    -0.9429745712289743,
1033    -0.8949919978782753,
1034    -0.833442628760834,
1035    -0.7592592630373576,
1036    -0.6735663684734684,
1037    -0.577662930241223,
1038    -0.473002731445715,
1039    -0.36117230580938786,
1040    -0.24386688372098844,
1041    -0.1228646926107104,
1042    0.0,
1043    0.1228646926107104,
1044    0.24386688372098844,
1045    0.36117230580938786,
1046    0.473002731445715,
1047    0.577662930241223,
1048    0.6735663684734684,
1049    0.7592592630373576,
1050    0.833442628760834,
1051    0.8949919978782753,
1052    0.9429745712289743,
1053    0.9766639214595175,
1054    0.9955569697904981,
1055];
1056const LEGENDRE_ROOT_26: [f64; 26] = [
1057    -0.9958857011456169,
1058    -0.978385445956471,
1059    -0.9471590666617142,
1060    -0.9026378619843071,
1061    -0.845445942788498,
1062    -0.7763859488206789,
1063    -0.6964272604199573,
1064    -0.6066922930176181,
1065    -0.5084407148245057,
1066    -0.4030517551234863,
1067    -0.2920048394859569,
1068    -0.17685882035689018,
1069    -0.05923009342931321,
1070    0.05923009342931321,
1071    0.17685882035689018,
1072    0.2920048394859569,
1073    0.4030517551234863,
1074    0.5084407148245057,
1075    0.6066922930176181,
1076    0.6964272604199573,
1077    0.7763859488206789,
1078    0.845445942788498,
1079    0.9026378619843071,
1080    0.9471590666617142,
1081    0.978385445956471,
1082    0.9958857011456169,
1083];
1084const LEGENDRE_ROOT_27: [f64; 27] = [
1085    -0.9961792628889886,
1086    -0.9799234759615012,
1087    -0.9509005578147051,
1088    -0.9094823206774911,
1089    -0.8562079080182945,
1090    -0.7917716390705082,
1091    -0.7170134737394237,
1092    -0.6329079719464952,
1093    -0.5405515645794569,
1094    -0.44114825175002687,
1095    -0.3359939036385089,
1096    -0.22645936543953685,
1097    -0.11397258560952997,
1098    0.0,
1099    0.11397258560952997,
1100    0.22645936543953685,
1101    0.3359939036385089,
1102    0.44114825175002687,
1103    0.5405515645794569,
1104    0.6329079719464952,
1105    0.7170134737394237,
1106    0.7917716390705082,
1107    0.8562079080182945,
1108    0.9094823206774911,
1109    0.9509005578147051,
1110    0.9799234759615012,
1111    0.9961792628889886,
1112];
1113const LEGENDRE_ROOT_28: [f64; 28] = [
1114    -0.9964424975739544,
1115    -0.9813031653708727,
1116    -0.9542592806289382,
1117    -0.9156330263921321,
1118    -0.8658925225743951,
1119    -0.8056413709171791,
1120    -0.7356108780136318,
1121    -0.656651094038865,
1122    -0.5697204718114017,
1123    -0.4758742249551183,
1124    -0.3762515160890787,
1125    -0.2720616276351781,
1126    -0.16456928213338076,
1127    -0.05507928988403427,
1128    0.05507928988403427,
1129    0.16456928213338076,
1130    0.2720616276351781,
1131    0.3762515160890787,
1132    0.4758742249551183,
1133    0.5697204718114017,
1134    0.656651094038865,
1135    0.7356108780136318,
1136    0.8056413709171791,
1137    0.8658925225743951,
1138    0.9156330263921321,
1139    0.9542592806289382,
1140    0.9813031653708727,
1141    0.9964424975739544,
1142];
1143const LEGENDRE_ROOT_29: [f64; 29] = [
1144    -0.9966794422605966,
1145    -0.9825455052614132,
1146    -0.9572855957780877,
1147    -0.9211802329530587,
1148    -0.8746378049201028,
1149    -0.8181854876152524,
1150    -0.7524628517344771,
1151    -0.6782145376026865,
1152    -0.5962817971382278,
1153    -0.5075929551242276,
1154    -0.41315288817400864,
1155    -0.31403163786763993,
1156    -0.21135228616600107,
1157    -0.10627823013267923,
1158    0.0,
1159    0.10627823013267923,
1160    0.21135228616600107,
1161    0.31403163786763993,
1162    0.41315288817400864,
1163    0.5075929551242276,
1164    0.5962817971382278,
1165    0.6782145376026865,
1166    0.7524628517344771,
1167    0.8181854876152524,
1168    0.8746378049201028,
1169    0.9211802329530587,
1170    0.9572855957780877,
1171    0.9825455052614132,
1172    0.9966794422605966,
1173];
1174const LEGENDRE_ROOT_30: [f64; 30] = [
1175    -0.9968934840746495,
1176    -0.9836681232797472,
1177    -0.9600218649683075,
1178    -0.9262000474292743,
1179    -0.8825605357920527,
1180    -0.8295657623827684,
1181    -0.7677774321048262,
1182    -0.6978504947933158,
1183    -0.6205261829892429,
1184    -0.5366241481420199,
1185    -0.44703376953808915,
1186    -0.3527047255308781,
1187    -0.25463692616788985,
1188    -0.15386991360858354,
1189    -0.0514718425553177,
1190    0.0514718425553177,
1191    0.15386991360858354,
1192    0.25463692616788985,
1193    0.3527047255308781,
1194    0.44703376953808915,
1195    0.5366241481420199,
1196    0.6205261829892429,
1197    0.6978504947933158,
1198    0.7677774321048262,
1199    0.8295657623827684,
1200    0.8825605357920527,
1201    0.9262000474292743,
1202    0.9600218649683075,
1203    0.9836681232797472,
1204    0.9968934840746495,
1205];
1206
1207const LEGENDRE_WEIGHT_2: [f64; 2] = [1.0, 1.0];
1208const LEGENDRE_WEIGHT_3: [f64; 3] = [0.555555555555556, 0.888888888888889, 0.555555555555556];
1209const LEGENDRE_WEIGHT_4: [f64; 4] = [
1210    0.347854845137454,
1211    0.652145154862546,
1212    0.652145154862546,
1213    0.347854845137454,
1214];
1215const LEGENDRE_WEIGHT_5: [f64; 5] = [
1216    0.236926885056189,
1217    0.478628670499366,
1218    0.568888888888889,
1219    0.478628670499366,
1220    0.236926885056189,
1221];
1222const LEGENDRE_WEIGHT_6: [f64; 6] = [
1223    0.17132449237917,
1224    0.360761573048139,
1225    0.467913934572691,
1226    0.467913934572691,
1227    0.360761573048139,
1228    0.17132449237917,
1229];
1230const LEGENDRE_WEIGHT_7: [f64; 7] = [
1231    0.12948496616887,
1232    0.279705391489277,
1233    0.381830050505119,
1234    0.417959183673469,
1235    0.381830050505119,
1236    0.279705391489277,
1237    0.12948496616887,
1238];
1239const LEGENDRE_WEIGHT_8: [f64; 8] = [
1240    0.101228536290376,
1241    0.222381034453374,
1242    0.313706645877887,
1243    0.362683783378362,
1244    0.362683783378362,
1245    0.313706645877887,
1246    0.222381034453374,
1247    0.101228536290376,
1248];
1249const LEGENDRE_WEIGHT_9: [f64; 9] = [
1250    0.081274388361574,
1251    0.180648160694857,
1252    0.260610696402935,
1253    0.312347077040003,
1254    0.33023935500126,
1255    0.312347077040003,
1256    0.260610696402935,
1257    0.180648160694857,
1258    0.081274388361574,
1259];
1260const LEGENDRE_WEIGHT_10: [f64; 10] = [
1261    0.066671344308688,
1262    0.149451349150581,
1263    0.219086362515982,
1264    0.269266719309996,
1265    0.295524224714753,
1266    0.295524224714753,
1267    0.269266719309996,
1268    0.219086362515982,
1269    0.149451349150581,
1270    0.066671344308688,
1271];
1272const LEGENDRE_WEIGHT_11: [f64; 11] = [
1273    0.055668567116174,
1274    0.125580369464905,
1275    0.186290210927734,
1276    0.23319376459199,
1277    0.262804544510247,
1278    0.272925086777901,
1279    0.262804544510247,
1280    0.23319376459199,
1281    0.186290210927734,
1282    0.125580369464905,
1283    0.055668567116174,
1284];
1285const LEGENDRE_WEIGHT_12: [f64; 12] = [
1286    0.047175336386512,
1287    0.106939325995318,
1288    0.160078328543346,
1289    0.203167426723066,
1290    0.233492536538355,
1291    0.249147045813403,
1292    0.249147045813403,
1293    0.233492536538355,
1294    0.203167426723066,
1295    0.160078328543346,
1296    0.106939325995318,
1297    0.047175336386512,
1298];
1299const LEGENDRE_WEIGHT_13: [f64; 13] = [
1300    0.040484004765316,
1301    0.092121499837728,
1302    0.138873510219787,
1303    0.178145980761946,
1304    0.207816047536889,
1305    0.226283180262897,
1306    0.232551553230874,
1307    0.226283180262897,
1308    0.207816047536889,
1309    0.178145980761946,
1310    0.138873510219787,
1311    0.092121499837728,
1312    0.040484004765316,
1313];
1314const LEGENDRE_WEIGHT_14: [f64; 14] = [
1315    0.035119460331752,
1316    0.08015808715976,
1317    0.121518570687903,
1318    0.157203167158194,
1319    0.185538397477938,
1320    0.20519846372129,
1321    0.215263853463158,
1322    0.215263853463158,
1323    0.20519846372129,
1324    0.185538397477938,
1325    0.157203167158194,
1326    0.121518570687903,
1327    0.08015808715976,
1328    0.035119460331752,
1329];
1330const LEGENDRE_WEIGHT_15: [f64; 15] = [
1331    0.030753241996117,
1332    0.070366047488108,
1333    0.107159220467172,
1334    0.139570677926154,
1335    0.166269205816994,
1336    0.186161000015562,
1337    0.198431485327111,
1338    0.202578241925561,
1339    0.198431485327111,
1340    0.186161000015562,
1341    0.166269205816994,
1342    0.139570677926154,
1343    0.107159220467172,
1344    0.070366047488108,
1345    0.030753241996117,
1346];
1347const LEGENDRE_WEIGHT_16: [f64; 16] = [
1348    0.027152459411754,
1349    0.062253523938648,
1350    0.095158511682493,
1351    0.124628971255534,
1352    0.149595988816577,
1353    0.169156519395003,
1354    0.182603415044924,
1355    0.189450610455069,
1356    0.189450610455069,
1357    0.182603415044924,
1358    0.169156519395003,
1359    0.149595988816577,
1360    0.124628971255534,
1361    0.095158511682493,
1362    0.062253523938648,
1363    0.027152459411754,
1364];
1365const LEGENDRE_WEIGHT_17: [f64; 17] = [
1366    0.02414830286854793,
1367    0.0554595293739872,
1368    0.08503614831717918,
1369    0.11188384719340397,
1370    0.13513636846852548,
1371    0.15404576107681028,
1372    0.16800410215645004,
1373    0.17656270536699264,
1374    0.17944647035620653,
1375    0.17656270536699264,
1376    0.16800410215645004,
1377    0.15404576107681028,
1378    0.13513636846852548,
1379    0.11188384719340397,
1380    0.08503614831717918,
1381    0.0554595293739872,
1382    0.02414830286854793,
1383];
1384const LEGENDRE_WEIGHT_18: [f64; 18] = [
1385    0.02161601352648331,
1386    0.0497145488949698,
1387    0.07642573025488907,
1388    0.10094204410628717,
1389    0.12255520671147846,
1390    0.14064291467065065,
1391    0.15468467512626524,
1392    0.16427648374583273,
1393    0.1691423829631436,
1394    0.1691423829631436,
1395    0.16427648374583273,
1396    0.15468467512626524,
1397    0.14064291467065065,
1398    0.12255520671147846,
1399    0.10094204410628717,
1400    0.07642573025488907,
1401    0.0497145488949698,
1402    0.02161601352648331,
1403];
1404const LEGENDRE_WEIGHT_19: [f64; 19] = [
1405    0.019461788229726478,
1406    0.0448142267656996,
1407    0.06904454273764123,
1408    0.09149002162245,
1409    0.11156664554733399,
1410    0.12875396253933621,
1411    0.1426067021736066,
1412    0.15276604206585967,
1413    0.15896884339395434,
1414    0.1610544498487837,
1415    0.15896884339395434,
1416    0.15276604206585967,
1417    0.1426067021736066,
1418    0.12875396253933621,
1419    0.11156664554733399,
1420    0.09149002162245,
1421    0.06904454273764123,
1422    0.0448142267656996,
1423    0.019461788229726478,
1424];
1425const LEGENDRE_WEIGHT_20: [f64; 20] = [
1426    0.017614007139152118,
1427    0.04060142980038694,
1428    0.06267204833410907,
1429    0.08327674157670475,
1430    0.10193011981724044,
1431    0.11819453196151841,
1432    0.13168863844917664,
1433    0.14209610931838204,
1434    0.14917298647260374,
1435    0.15275338713072584,
1436    0.15275338713072584,
1437    0.14917298647260374,
1438    0.14209610931838204,
1439    0.13168863844917664,
1440    0.11819453196151841,
1441    0.10193011981724044,
1442    0.08327674157670475,
1443    0.06267204833410907,
1444    0.04060142980038694,
1445    0.017614007139152118,
1446];
1447const LEGENDRE_WEIGHT_21: [f64; 21] = [
1448    0.016017228257774335,
1449    0.036953789770852494,
1450    0.05713442542685721,
1451    0.0761001136283793,
1452    0.09344442345603386,
1453    0.10879729916714838,
1454    0.12183141605372853,
1455    0.13226893863333747,
1456    0.13988739479107315,
1457    0.14452440398997005,
1458    0.14608113364969041,
1459    0.14452440398997005,
1460    0.13988739479107315,
1461    0.13226893863333747,
1462    0.12183141605372853,
1463    0.10879729916714838,
1464    0.09344442345603386,
1465    0.0761001136283793,
1466    0.05713442542685721,
1467    0.036953789770852494,
1468    0.016017228257774335,
1469];
1470const LEGENDRE_WEIGHT_22: [f64; 22] = [
1471    0.014627995298272202,
1472    0.03377490158481416,
1473    0.052293335152683286,
1474    0.06979646842452049,
1475    0.08594160621706773,
1476    0.10041414444288096,
1477    0.11293229608053922,
1478    0.12325237681051242,
1479    0.13117350478706238,
1480    0.13654149834601517,
1481    0.13925187285563198,
1482    0.13925187285563198,
1483    0.13654149834601517,
1484    0.13117350478706238,
1485    0.12325237681051242,
1486    0.11293229608053922,
1487    0.10041414444288096,
1488    0.08594160621706773,
1489    0.06979646842452049,
1490    0.052293335152683286,
1491    0.03377490158481416,
1492    0.014627995298272202,
1493];
1494const LEGENDRE_WEIGHT_23: [f64; 23] = [
1495    0.013411859487141771,
1496    0.030988005856979445,
1497    0.04803767173108467,
1498    0.06423242140852585,
1499    0.07928141177671895,
1500    0.09291576606003515,
1501    0.10489209146454141,
1502    0.11499664022241136,
1503    0.12304908430672953,
1504    0.12890572218808216,
1505    0.1324620394046966,
1506    0.13365457218610619,
1507    0.1324620394046966,
1508    0.12890572218808216,
1509    0.12304908430672953,
1510    0.11499664022241136,
1511    0.10489209146454141,
1512    0.09291576606003515,
1513    0.07928141177671895,
1514    0.06423242140852585,
1515    0.04803767173108467,
1516    0.030988005856979445,
1517    0.013411859487141771,
1518];
1519const LEGENDRE_WEIGHT_24: [f64; 24] = [
1520    0.0123412297999872,
1521    0.028531388628933663,
1522    0.04427743881741981,
1523    0.05929858491543678,
1524    0.07334648141108031,
1525    0.08619016153195327,
1526    0.09761865210411388,
1527    0.10744427011596563,
1528    0.1155056680537256,
1529    0.12167047292780339,
1530    0.1258374563468283,
1531    0.12793819534675216,
1532    0.12793819534675216,
1533    0.1258374563468283,
1534    0.12167047292780339,
1535    0.1155056680537256,
1536    0.10744427011596563,
1537    0.09761865210411388,
1538    0.08619016153195327,
1539    0.07334648141108031,
1540    0.05929858491543678,
1541    0.04427743881741981,
1542    0.028531388628933663,
1543    0.0123412297999872,
1544];
1545const LEGENDRE_WEIGHT_25: [f64; 25] = [
1546    0.011393798501026288,
1547    0.026354986615032137,
1548    0.040939156701306316,
1549    0.054904695975835194,
1550    0.06803833381235691,
1551    0.08014070033500102,
1552    0.09102826198296365,
1553    0.10053594906705064,
1554    0.10851962447426365,
1555    0.11485825914571164,
1556    0.11945576353578477,
1557    0.12224244299031004,
1558    0.12317605372671545,
1559    0.12224244299031004,
1560    0.11945576353578477,
1561    0.11485825914571164,
1562    0.10851962447426365,
1563    0.10053594906705064,
1564    0.09102826198296365,
1565    0.08014070033500102,
1566    0.06803833381235691,
1567    0.054904695975835194,
1568    0.040939156701306316,
1569    0.026354986615032137,
1570    0.011393798501026288,
1571];
1572const LEGENDRE_WEIGHT_26: [f64; 26] = [
1573    0.01055137261734301,
1574    0.02441785109263191,
1575    0.03796238329436276,
1576    0.05097582529714781,
1577    0.06327404632957484,
1578    0.07468414976565975,
1579    0.08504589431348523,
1580    0.09421380035591415,
1581    0.10205916109442542,
1582    0.10847184052857659,
1583    0.11336181654631967,
1584    0.11666044348529658,
1585    0.11832141527926228,
1586    0.11832141527926228,
1587    0.11666044348529658,
1588    0.11336181654631967,
1589    0.10847184052857659,
1590    0.10205916109442542,
1591    0.09421380035591415,
1592    0.08504589431348523,
1593    0.07468414976565975,
1594    0.06327404632957484,
1595    0.05097582529714781,
1596    0.03796238329436276,
1597    0.02441785109263191,
1598    0.01055137261734301,
1599];
1600const LEGENDRE_WEIGHT_27: [f64; 27] = [
1601    0.00979899605129436,
1602    0.02268623159618062,
1603    0.03529705375741971,
1604    0.04744941252061506,
1605    0.0589835368598336,
1606    0.0697488237662456,
1607    0.07960486777305777,
1608    0.08842315854375694,
1609    0.0960887273700285,
1610    0.1025016378177458,
1611    0.10757828578853319,
1612    0.11125248835684519,
1613    0.11347634610896515,
1614    0.114220867378957,
1615    0.11347634610896515,
1616    0.11125248835684519,
1617    0.10757828578853319,
1618    0.1025016378177458,
1619    0.0960887273700285,
1620    0.08842315854375694,
1621    0.07960486777305777,
1622    0.0697488237662456,
1623    0.0589835368598336,
1624    0.04744941252061506,
1625    0.03529705375741971,
1626    0.02268623159618062,
1627    0.00979899605129436,
1628];
1629const LEGENDRE_WEIGHT_28: [f64; 28] = [
1630    0.00912428259309452,
1631    0.02113211259277126,
1632    0.03290142778230438,
1633    0.04427293475900423,
1634    0.05510734567571675,
1635    0.0652729239669996,
1636    0.07464621423456878,
1637    0.08311341722890121,
1638    0.09057174439303284,
1639    0.09693065799792992,
1640    0.10211296757806076,
1641    0.10605576592284642,
1642    0.10871119225829413,
1643    0.1100470130164752,
1644    0.1100470130164752,
1645    0.10871119225829413,
1646    0.10605576592284642,
1647    0.10211296757806076,
1648    0.09693065799792992,
1649    0.09057174439303284,
1650    0.08311341722890121,
1651    0.07464621423456878,
1652    0.0652729239669996,
1653    0.05510734567571675,
1654    0.04427293475900423,
1655    0.03290142778230438,
1656    0.02113211259277126,
1657    0.00912428259309452,
1658];
1659const LEGENDRE_WEIGHT_29: [f64; 29] = [
1660    0.00851690387874641,
1661    0.019732085056122707,
1662    0.030740492202093624,
1663    0.041402062518682836,
1664    0.05159482690249792,
1665    0.061203090657079136,
1666    0.07011793325505128,
1667    0.07823832713576379,
1668    0.08547225736617253,
1669    0.09173775713925876,
1670    0.0969638340944086,
1671    0.10109127375991497,
1672    0.10407331007772938,
1673    0.10587615509732094,
1674    0.10647938171831424,
1675    0.10587615509732094,
1676    0.10407331007772938,
1677    0.10109127375991497,
1678    0.0969638340944086,
1679    0.09173775713925876,
1680    0.08547225736617253,
1681    0.07823832713576379,
1682    0.07011793325505128,
1683    0.061203090657079136,
1684    0.05159482690249792,
1685    0.041402062518682836,
1686    0.030740492202093624,
1687    0.019732085056122707,
1688    0.00851690387874641,
1689];
1690const LEGENDRE_WEIGHT_30: [f64; 30] = [
1691    0.007968192496166607,
1692    0.01846646831109096,
1693    0.02878470788332337,
1694    0.03879919256962705,
1695    0.04840267283059405,
1696    0.057493156217619065,
1697    0.0659742298821805,
1698    0.0737559747377052,
1699    0.08075589522942021,
1700    0.08689978720108298,
1701    0.09212252223778614,
1702    0.09636873717464425,
1703    0.09959342058679527,
1704    0.1017623897484055,
1705    0.10285265289355884,
1706    0.10285265289355884,
1707    0.1017623897484055,
1708    0.09959342058679527,
1709    0.09636873717464425,
1710    0.09212252223778614,
1711    0.08689978720108298,
1712    0.08075589522942021,
1713    0.0737559747377052,
1714    0.0659742298821805,
1715    0.057493156217619065,
1716    0.04840267283059405,
1717    0.03879919256962705,
1718    0.02878470788332337,
1719    0.01846646831109096,
1720    0.007968192496166607,
1721];
1722
1723// =============================================================================
1724// Table for Kronrod Quadrature
1725// =============================================================================
1726const KRONROD_NODES_15: [f64; 15] = [
1727    -0.9914553711208126,
1728    -0.9491079123427585,
1729    -0.8648644233597691,
1730    -0.7415311855993945,
1731    -0.5860872354676911,
1732    -0.4058451513773972,
1733    -0.20778495500789848,
1734    0.0,
1735    0.20778495500789848,
1736    0.4058451513773972,
1737    0.5860872354676911,
1738    0.7415311855993945,
1739    0.8648644233597691,
1740    0.9491079123427585,
1741    0.9914553711208126,
1742];
1743
1744const KRONROD_WEIGHTS_15: [f64; 15] = [
1745    0.022935322010529224,
1746    0.06309209262997854,
1747    0.10479001032225017,
1748    0.14065325971552592,
1749    0.1690047266392679,
1750    0.19035057806478542,
1751    0.20443294007529889,
1752    0.20948214108472782,
1753    0.20443294007529889,
1754    0.19035057806478542,
1755    0.1690047266392679,
1756    0.14065325971552592,
1757    0.10479001032225017,
1758    0.06309209262997854,
1759    0.022935322010529224,
1760];
1761
1762const KRONROD_NODES_21: [f64; 21] = [
1763    -0.9956571630258081,
1764    -0.9739065285171717,
1765    -0.9301574913557082,
1766    -0.8650633666889845,
1767    -0.7808177265864169,
1768    -0.6794095682990244,
1769    -0.5627571346686047,
1770    -0.4333953941292472,
1771    -0.2943928627014602,
1772    -0.14887433898163122,
1773    0.0,
1774    0.14887433898163122,
1775    0.2943928627014602,
1776    0.4333953941292472,
1777    0.5627571346686047,
1778    0.6794095682990244,
1779    0.7808177265864169,
1780    0.8650633666889845,
1781    0.9301574913557082,
1782    0.9739065285171717,
1783    0.9956571630258081,
1784];
1785
1786const KRONROD_WEIGHTS_21: [f64; 21] = [
1787    0.011694638867371874,
1788    0.032558162307964725,
1789    0.054755896574351995,
1790    0.07503967481091996,
1791    0.09312545458369761,
1792    0.10938715880229764,
1793    0.12349197626206584,
1794    0.13470921731147334,
1795    0.14277593857706009,
1796    0.14773910490133849,
1797    0.1494455540029169,
1798    0.14773910490133849,
1799    0.14277593857706009,
1800    0.13470921731147334,
1801    0.12349197626206584,
1802    0.10938715880229764,
1803    0.09312545458369761,
1804    0.07503967481091996,
1805    0.054755896574351995,
1806    0.032558162307964725,
1807    0.011694638867371874,
1808];
1809
1810const KRONROD_NODES_31: [f64; 31] = [
1811    -0.9980022986933971,
1812    -0.9879925180204854,
1813    -0.9677390756791391,
1814    -0.937273392400706,
1815    -0.8972645323440819,
1816    -0.8482065834104272,
1817    -0.790418501442466,
1818    -0.7244177313601701,
1819    -0.650996741297417,
1820    -0.5709721726085388,
1821    -0.4850818636402397,
1822    -0.3941513470775634,
1823    -0.29918000715316884,
1824    -0.20119409399743451,
1825    -0.1011420669187175,
1826    0.0,
1827    0.1011420669187175,
1828    0.20119409399743451,
1829    0.29918000715316884,
1830    0.3941513470775634,
1831    0.4850818636402397,
1832    0.5709721726085388,
1833    0.650996741297417,
1834    0.7244177313601701,
1835    0.790418501442466,
1836    0.8482065834104272,
1837    0.8972645323440819,
1838    0.937273392400706,
1839    0.9677390756791391,
1840    0.9879925180204854,
1841    0.9980022986933971,
1842];
1843
1844const KRONROD_WEIGHTS_31: [f64; 31] = [
1845    0.005377479872923349,
1846    0.015007947329316124,
1847    0.02546084732671532,
1848    0.03534636079137585,
1849    0.04458975132476488,
1850    0.05348152469092809,
1851    0.06200956780067064,
1852    0.06985412131872826,
1853    0.07684968075772038,
1854    0.08308050282313302,
1855    0.08856444305621176,
1856    0.09312659817082532,
1857    0.09664272698362368,
1858    0.09917359872179196,
1859    0.10076984552387559,
1860    0.10133000701479154,
1861    0.10076984552387559,
1862    0.09917359872179196,
1863    0.09664272698362368,
1864    0.09312659817082532,
1865    0.08856444305621176,
1866    0.08308050282313302,
1867    0.07684968075772038,
1868    0.06985412131872826,
1869    0.06200956780067064,
1870    0.05348152469092809,
1871    0.04458975132476488,
1872    0.03534636079137585,
1873    0.02546084732671532,
1874    0.015007947329316124,
1875    0.005377479872923349,
1876];
1877
1878const KRONROD_NODES_41: [f64; 41] = [
1879    -0.9988590315882777,
1880    -0.9931285991850949,
1881    -0.9815078774502503,
1882    -0.9639719272779138,
1883    -0.9408226338317548,
1884    -0.912234428251326,
1885    -0.878276811252282,
1886    -0.8391169718222188,
1887    -0.7950414288375512,
1888    -0.7463319064601508,
1889    -0.6932376563347514,
1890    -0.636053680726515,
1891    -0.5751404468197103,
1892    -0.5108670019508271,
1893    -0.4435931752387251,
1894    -0.37370608871541955,
1895    -0.301627868114913,
1896    -0.22778585114164507,
1897    -0.15260546524092267,
1898    -0.07652652113349734,
1899    0.0,
1900    0.07652652113349734,
1901    0.15260546524092267,
1902    0.22778585114164507,
1903    0.301627868114913,
1904    0.37370608871541955,
1905    0.4435931752387251,
1906    0.5108670019508271,
1907    0.5751404468197103,
1908    0.636053680726515,
1909    0.6932376563347514,
1910    0.7463319064601508,
1911    0.7950414288375512,
1912    0.8391169718222188,
1913    0.878276811252282,
1914    0.912234428251326,
1915    0.9408226338317548,
1916    0.9639719272779138,
1917    0.9815078774502503,
1918    0.9931285991850949,
1919    0.9988590315882777,
1920];
1921
1922const KRONROD_WEIGHTS_41: [f64; 41] = [
1923    0.00307358371852053,
1924    0.008600269855642943,
1925    0.014626169256971253,
1926    0.020388373461266523,
1927    0.02588213360495116,
1928    0.0312873067770328,
1929    0.036600169758200796,
1930    0.04166887332797369,
1931    0.04643482186749767,
1932    0.05094457392372869,
1933    0.05519510534828599,
1934    0.05911140088063957,
1935    0.06265323755478117,
1936    0.06583459713361842,
1937    0.06864867292852161,
1938    0.07105442355344407,
1939    0.07303069033278667,
1940    0.07458287540049918,
1941    0.07570449768455667,
1942    0.07637786767208074,
1943    0.07660071191799966,
1944    0.07637786767208074,
1945    0.07570449768455667,
1946    0.07458287540049918,
1947    0.07303069033278667,
1948    0.07105442355344407,
1949    0.06864867292852161,
1950    0.06583459713361842,
1951    0.06265323755478117,
1952    0.05911140088063957,
1953    0.05519510534828599,
1954    0.05094457392372869,
1955    0.04643482186749767,
1956    0.04166887332797369,
1957    0.036600169758200796,
1958    0.0312873067770328,
1959    0.02588213360495116,
1960    0.020388373461266523,
1961    0.014626169256971253,
1962    0.008600269855642943,
1963    0.00307358371852053,
1964];
1965
1966const KRONROD_NODES_51: [f64; 51] = [
1967    -0.9992621049926098,
1968    -0.9955569697904981,
1969    -0.9880357945340772,
1970    -0.9766639214595175,
1971    -0.9616149864258425,
1972    -0.9429745712289743,
1973    -0.9207471152817016,
1974    -0.8949919978782753,
1975    -0.8658470652932756,
1976    -0.833442628760834,
1977    -0.7978737979985001,
1978    -0.7592592630373576,
1979    -0.7177664068130843,
1980    -0.6735663684734684,
1981    -0.6268100990103174,
1982    -0.577662930241223,
1983    -0.5263252843347191,
1984    -0.473002731445715,
1985    -0.4178853821930377,
1986    -0.36117230580938786,
1987    -0.30308953893110785,
1988    -0.24386688372098844,
1989    -0.1837189394210489,
1990    -0.1228646926107104,
1991    -0.06154448300568508,
1992    0.0,
1993    0.06154448300568508,
1994    0.1228646926107104,
1995    0.1837189394210489,
1996    0.24386688372098844,
1997    0.30308953893110785,
1998    0.36117230580938786,
1999    0.4178853821930377,
2000    0.473002731445715,
2001    0.5263252843347191,
2002    0.577662930241223,
2003    0.6268100990103174,
2004    0.6735663684734684,
2005    0.7177664068130843,
2006    0.7592592630373576,
2007    0.7978737979985001,
2008    0.833442628760834,
2009    0.8658470652932756,
2010    0.8949919978782753,
2011    0.9207471152817016,
2012    0.9429745712289743,
2013    0.9616149864258425,
2014    0.9766639214595175,
2015    0.9880357945340772,
2016    0.9955569697904981,
2017    0.9992621049926098,
2018];
2019
2020const KRONROD_WEIGHTS_51: [f64; 51] = [
2021    0.001987383892330316,
2022    0.005561932135356714,
2023    0.009473973386174152,
2024    0.013236229195571676,
2025    0.0168478177091283,
2026    0.02043537114588284,
2027    0.02400994560695322,
2028    0.02747531758785174,
2029    0.030792300167387487,
2030    0.03400213027432934,
2031    0.03711627148341554,
2032    0.04008382550403238,
2033    0.04287284502017005,
2034    0.04550291304992179,
2035    0.04798253713883671,
2036    0.05027767908071567,
2037    0.05236288580640747,
2038    0.05425112988854549,
2039    0.055950811220412316,
2040    0.057437116361567835,
2041    0.058689680022394206,
2042    0.05972034032417406,
2043    0.06053945537604586,
2044    0.06112850971705305,
2045    0.061471189871425316,
2046    0.061580818067832936,
2047    0.061471189871425316,
2048    0.06112850971705305,
2049    0.06053945537604586,
2050    0.05972034032417406,
2051    0.058689680022394206,
2052    0.057437116361567835,
2053    0.055950811220412316,
2054    0.05425112988854549,
2055    0.05236288580640747,
2056    0.05027767908071567,
2057    0.04798253713883671,
2058    0.04550291304992179,
2059    0.04287284502017005,
2060    0.04008382550403238,
2061    0.03711627148341554,
2062    0.03400213027432934,
2063    0.030792300167387487,
2064    0.02747531758785174,
2065    0.02400994560695322,
2066    0.02043537114588284,
2067    0.0168478177091283,
2068    0.013236229195571676,
2069    0.009473973386174152,
2070    0.005561932135356714,
2071    0.001987383892330316,
2072];
2073
2074const KRONROD_NODES_61: [f64; 61] = [
2075    -0.9994844100504906,
2076    -0.9968934840746495,
2077    -0.9916309968704046,
2078    -0.9836681232797472,
2079    -0.9731163225011262,
2080    -0.9600218649683075,
2081    -0.94437444474856,
2082    -0.9262000474292743,
2083    -0.9055733076999078,
2084    -0.8825605357920527,
2085    -0.8572052335460612,
2086    -0.8295657623827684,
2087    -0.799727835821839,
2088    -0.7677774321048262,
2089    -0.7337900624532268,
2090    -0.6978504947933158,
2091    -0.6600610641266269,
2092    -0.6205261829892429,
2093    -0.5793452358263617,
2094    -0.5366241481420199,
2095    -0.49248046786177857,
2096    -0.44703376953808915,
2097    -0.4004012548303944,
2098    -0.3527047255308781,
2099    -0.30407320227362505,
2100    -0.25463692616788985,
2101    -0.20452511668230988,
2102    -0.15386991360858354,
2103    -0.10280693796673702,
2104    -0.0514718425553177,
2105    0.0,
2106    0.0514718425553177,
2107    0.10280693796673702,
2108    0.15386991360858354,
2109    0.20452511668230988,
2110    0.25463692616788985,
2111    0.30407320227362505,
2112    0.3527047255308781,
2113    0.4004012548303944,
2114    0.44703376953808915,
2115    0.49248046786177857,
2116    0.5366241481420199,
2117    0.5793452358263617,
2118    0.6205261829892429,
2119    0.6600610641266269,
2120    0.6978504947933158,
2121    0.7337900624532268,
2122    0.7677774321048262,
2123    0.799727835821839,
2124    0.8295657623827684,
2125    0.8572052335460612,
2126    0.8825605357920527,
2127    0.9055733076999078,
2128    0.9262000474292743,
2129    0.94437444474856,
2130    0.9600218649683075,
2131    0.9731163225011262,
2132    0.9836681232797472,
2133    0.9916309968704046,
2134    0.9968934840746495,
2135    0.9994844100504906,
2136];
2137
2138const KRONROD_WEIGHTS_61: [f64; 61] = [
2139    0.001389013698677008,
2140    0.003890461127099884,
2141    0.006630703915931292,
2142    0.009273279659517762,
2143    0.011823015253496341,
2144    0.0143697295070458,
2145    0.016920889189053274,
2146    0.019414141193942382,
2147    0.02182803582160919,
2148    0.0241911620780806,
2149    0.0265099548823331,
2150    0.02875404876504129,
2151    0.030907257562387762,
2152    0.03298144705748373,
2153    0.034979338028060025,
2154    0.03688236465182123,
2155    0.038678945624727595,
2156    0.04037453895153596,
2157    0.041969810215164244,
2158    0.04345253970135607,
2159    0.04481480013316266,
2160    0.04605923827100699,
2161    0.04718554656929915,
2162    0.048185861757087126,
2163    0.04905543455502978,
2164    0.0497956834270742,
2165    0.05040592140278235,
2166    0.0508817958987496,
2167    0.051221547849258774,
2168    0.05142612853745902,
2169    0.05149472942945157,
2170    0.05142612853745902,
2171    0.051221547849258774,
2172    0.0508817958987496,
2173    0.05040592140278235,
2174    0.0497956834270742,
2175    0.04905543455502978,
2176    0.048185861757087126,
2177    0.04718554656929915,
2178    0.04605923827100699,
2179    0.04481480013316266,
2180    0.04345253970135607,
2181    0.041969810215164244,
2182    0.04037453895153596,
2183    0.038678945624727595,
2184    0.03688236465182123,
2185    0.034979338028060025,
2186    0.03298144705748373,
2187    0.030907257562387762,
2188    0.02875404876504129,
2189    0.0265099548823331,
2190    0.0241911620780806,
2191    0.02182803582160919,
2192    0.019414141193942382,
2193    0.016920889189053274,
2194    0.0143697295070458,
2195    0.011823015253496341,
2196    0.009273279659517762,
2197    0.006630703915931292,
2198    0.003890461127099884,
2199    0.001389013698677008,
2200];
peroxide/numerical/integral.rs

peroxide/numerical/
integral.rs
