pub(crate) fn erf(x: f64) -> f64 {
    let cof = vec![
        -1.3026537197817094,
        6.419_697_923_564_902e-1,
        1.9476473204185836e-2,
        -9.561_514_786_808_63e-3,
        -9.46595344482036e-4,
        3.66839497852761e-4,
        4.2523324806907e-5,
        -2.0278578112534e-5,
        -1.624290004647e-6,
        1.303655835580e-6,
        1.5626441722e-8,
        -8.5238095915e-8,
        6.529054439e-9,
        5.059343495e-9,
        -9.91364156e-10,
        -2.27365122e-10,
        9.6467911e-11,
        2.394038e-12,
        -6.886027e-12,
        8.94487e-13,
        3.13092e-13,
        -1.12708e-13,
        3.81e-16,
        7.106e-15,
        -1.523e-15,
        -9.4e-17,
        1.21e-16,
        -2.8e-17,
    ];

    let mut d = 0.0;
    let mut dd = 0.0;

    let x_abs = x.abs();

    let t = 2.0 / (2.0 + x_abs);
    let ty = 4.0 * t - 2.0;

    for j in (1..=cof.len() - 1).rev() {
        let tmp = d;
        d = ty * d - dd + cof[j];
        dd = tmp;
    }

    let res = t * f64::exp(-x_abs * x_abs + 0.5 * (cof[0] + ty * d) - dd);
    if x < 0.0 {
        res - 1.0
    } else {
        1.0 - res
    }
}