gkquad 0.0.4

Numerical integration library for Rust
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

use super::super::common::{Integrand, Range};
use super::super::qk::QKResult;
use super::super::util::{rescale_error, Array};

#[cfg(not(feature = "std"))]
use crate::float::Float;

/// perform Gauss-Kronrod integration with custom points
///
/// # Parameters
///
/// * xgk: abscissae of the Gauss-Kronrod rule
/// * wg: weight of the Gauss rule
/// * wgk: weight of the Kronrod rule
/// * wck: Kronrod weight for center point
/// * buf: buffer to hold return value of integrand, may be not initialized.
///
/// # Precondition
///
/// The arguments must satisfy the following conditions, otherwise this
/// function causes an undefined behaviour.
///
/// * xgk.len() > 0
/// * xgk.len() / 2 == wg.len()
/// * buf.len() >= xgk.len() * 2 - 1
///
/// # TODO
///
/// - SIMD implementation
pub unsafe fn qk<F, K, G>(
    f: &mut F,
    range: &Range,
    xgk: &K,
    wg: &G,
    wgk: &K,
    wck: f64,
    buf: &mut [f64],
) -> QKResult
where
    F: Integrand,
    K: Array<Item = f64>,
    G: Array<Item = f64>,
{
    let xgk = xgk.as_slice();
    let wg = wg.as_slice();
    let wgk = wgk.as_slice();

    debug_assert!(!xgk.is_empty());
    debug_assert!(xgk.len() == wg.len() * 2);
    debug_assert!(buf.len() > xgk.len() * 2);

    let n = K::CAPACITY;
    let center = 0.5 * (range.begin + range.end);
    let half_length = 0.5 * (range.end - range.begin);
    let abs_half_length = half_length.abs();

    *buf.get_unchecked_mut(n << 1) = center;

    for j in 0..n {
        let abscissa = half_length * xgk.get_unchecked(j);
        *buf.get_unchecked_mut(j) = center - abscissa;
        *buf.get_unchecked_mut(j + n) = center + abscissa;
    }

    f.apply_to_slice(buf);

    let f_center = buf.get_unchecked(n << 1);
    let mut result_gauss = 0.;
    let mut result_kronrod = f_center * wck;
    let mut result_abs = result_kronrod.abs();

    for j in 0..n {
        let fval1 = *buf.get_unchecked(j);
        let fval2 = *buf.get_unchecked(j + n);
        let fsum = fval1 + fval2;
        result_kronrod += wgk.get_unchecked(j) * fsum;
        result_abs += wgk.get_unchecked(j) * (fval1.abs() + fval2.abs());

        if j % 2 == 0 {
            result_gauss += *wg.get_unchecked(j / 2) * fsum;
        }
    }

    let mean = result_kronrod * 0.5;
    let mut result_asc = wck * (f_center - mean).abs();

    for j in 0..n {
        result_asc += wgk.get_unchecked(j)
            * ((buf.get_unchecked(j) - mean).abs() + (buf.get_unchecked(j + n) - mean).abs());
    }

    let err = (result_kronrod - result_gauss) * half_length;
    result_kronrod *= half_length;
    result_abs *= abs_half_length;
    result_asc *= abs_half_length;

    QKResult {
        estimate: result_kronrod,
        delta: rescale_error(err, result_abs, result_asc),
        absvalue: result_abs,
        asc: result_asc,
    }
}