use nalgebra::{ComplexField, RealField};
use num_traits::{FromPrimitive, One, Zero};
use std::f64;
mod gaussian;
pub use gaussian::*;
mod tables;
use tables::WEIGHTS_DE;
fn integrate_core<N: ComplexField + FromPrimitive + Copy, F: FnMut(N::RealField) -> N>(
mut f: F,
tol: N::RealField,
) -> Result<N, String>
where
<N as ComplexField>::RealField: FromPrimitive + Copy,
{
let mut error_estimate = N::RealField::one() + tol;
let mut num_function_evaluations = 1;
let mut current_delta = N::RealField::zero();
let half = N::from_f64(0.5).unwrap();
let one_point_nine = N::from_f64(1.9).unwrap().real();
let two_point_one = N::from_f64(2.1).unwrap().real();
let pi = N::from_f64(f64::consts::PI).unwrap();
let mut integral = pi * f(N::RealField::zero());
for &weight in &WEIGHTS_DE {
let new_contribution = weight
.iter()
.map(|&(w, x)| {
let x = N::from_f64(x).unwrap().real();
N::from_f64(w).unwrap() * (f(x) + f(-x))
})
.fold(N::zero(), |sum, x| sum + x);
num_function_evaluations += 2 * weight.len();
let previous_delta_ln = current_delta.ln();
current_delta = (half * integral - new_contribution).abs();
integral = half * integral + new_contribution;
if num_function_evaluations <= 13 {
continue; }
if current_delta == N::RealField::zero() {
error_estimate = N::RealField::zero();
break;
}
let r = current_delta.ln() / previous_delta_ln;
if r > one_point_nine && r < two_point_one {
error_estimate = current_delta * current_delta;
} else {
error_estimate = current_delta;
}
if error_estimate < tol {
break;
}
}
if error_estimate < tol {
Ok(integral)
} else {
Err("integrate: maximum iterations exceeded".to_owned())
}
}
pub fn integrate<N: ComplexField + FromPrimitive + Copy, F: FnMut(N::RealField) -> N>(
left: N::RealField,
right: N::RealField,
mut f: F,
tol: N::RealField,
) -> Result<N, String>
where
<N as ComplexField>::RealField: FromPrimitive + Copy,
{
if left >= right {
return Err("integrate: left must be less than right".to_owned());
}
if !tol.is_sign_positive() {
return Err("integrate: tolerance must be positive".to_owned());
}
let half = N::from_f64(0.5).unwrap().real();
let scale = (right - left) * half;
let shift = (right + left) * half;
let scale_cmplx = N::from_real(scale);
let fun = |x: N::RealField| -> N {
let out = f(scale * x + shift);
if out.is_finite() {
out
} else {
N::zero()
}
};
Ok(integrate_core(fun, tol)? * scale_cmplx)
}
pub fn integrate_simpson<N: ComplexField + FromPrimitive + Copy, F: FnMut(N::RealField) -> N>(
left: N::RealField,
right: N::RealField,
mut f: F,
tol: N::RealField,
n_max: usize,
) -> Result<N, String>
where
<N as ComplexField>::RealField: FromPrimitive + Copy,
{
if left >= right {
return Err("integrate: left must be less than right".to_owned());
}
if !tol.is_sign_positive() {
return Err("integrate: tolerance must be positive".to_owned());
}
let sixth = N::from_f64(1.0 / 6.0).unwrap();
let third = N::from_f64(1.0 / 3.0).unwrap();
let half_real = N::from_f64(0.5).unwrap().real();
let one_and_a_half_real = N::from_f64(1.5).unwrap().real();
let four = N::from_i32(4).unwrap();
let mut area = N::zero();
let mut i = 1;
let mut tol_i = vec![N::from_i32(10).unwrap().real() * tol];
let mut left_i = vec![left];
let mut step_i = vec![(right - left) * half_real];
let mut f_ai = vec![f(left)];
let mut f_ci = vec![f(left + step_i[0])];
let mut f_bi = vec![f(right)];
let mut sum_i = vec![N::from_real(step_i[0]) * (f_ai[0] + four * f_ci[0] + f_bi[0]) * third];
let mut l_i = vec![1];
while i > 0 {
let f_d = f(left_i[i - 1] + half_real * step_i[i - 1]);
let f_e = f(left_i[i - 1] + one_and_a_half_real * step_i[i - 1]);
let s1 = N::from_real(step_i[i - 1]) * (f_ai[i - 1] + four * f_d + f_ci[i - 1]) * sixth;
let s2 = N::from_real(step_i[i - 1]) * (f_ci[i - 1] + four * f_e + f_bi[i - 1]) * sixth;
let v_1 = left_i[i - 1];
let v_2 = f_ai[i - 1];
let v_3 = f_ci[i - 1];
let v_4 = f_bi[i - 1];
let v_5 = step_i[i - 1];
let v_6 = tol_i[i - 1];
let v_7 = sum_i[i - 1];
let v_8 = l_i[i - 1];
i -= 1;
if (s1 + s2 - v_7).abs() < v_6 {
area += s1 + s2;
} else {
if v_8 >= n_max {
return Err("integrate: maximum iterations exceeded".to_owned());
}
i += 1;
if i > left_i.len() {
left_i.push(v_1 + v_5);
f_ai.push(v_3);
f_ci.push(f_e);
f_bi.push(v_4);
step_i.push(half_real * v_5);
tol_i.push(half_real * v_6);
sum_i.push(s2);
l_i.push(v_8 + 1);
} else {
left_i[i - 1] = v_1 + v_5;
f_ai[i - 1] = v_3;
f_ci[i - 1] = f_e;
f_bi[i - 1] = v_4;
step_i[i - 1] = half_real * v_5;
tol_i[i - 1] = half_real * v_6;
sum_i[i - 1] = s2;
l_i[i - 1] = v_8 + 1;
}
i += 1;
if i > left_i.len() {
left_i.push(v_1);
f_ai.push(v_2);
f_ci.push(f_d);
f_bi.push(v_3);
step_i.push(step_i[i - 2]);
tol_i.push(tol_i[i - 2]);
sum_i.push(s1);
l_i.push(l_i[i - 2]);
} else {
left_i[i - 1] = v_1;
f_ai[i - 1] = v_2;
f_ci[i - 1] = f_d;
f_bi[i - 1] = v_3;
step_i[i - 1] = step_i[i - 2];
tol_i[i - 1] = tol_i[i - 2];
sum_i[i - 1] = sum_i[i - 2];
l_i[i - 1] = l_i[i - 2];
}
}
}
Ok(area)
}
pub fn integrate_fixed<N: ComplexField + FromPrimitive + Copy, F: FnMut(N::RealField) -> N>(
left: N::RealField,
right: N::RealField,
mut f: F,
n: usize,
) -> Result<N, String>
where
<N as ComplexField>::RealField: FromPrimitive + Copy,
{
if left >= right {
return Err("integrate_fixed: left must be less than right".to_owned());
}
let half = N::from_f64(0.5).unwrap();
let half_real = N::from_f64(0.5).unwrap().real();
let four = N::from_i32(4).unwrap();
let mut h = right - left;
let mut prev_rows = vec![N::zero(); n];
prev_rows[0] = N::from_real(h) * half * (f(left) + f(right));
let mut next = vec![N::zero(); n];
for i in 2..=n {
let mut acc = N::zero();
for k in 1..=(1 << (i - 2)) {
acc += f(left + N::from_f64(k as f64 - 0.5).unwrap().real() * h);
}
acc *= N::from_real(h);
acc += prev_rows[0];
acc *= half;
next[0] = acc;
for j in 2..=i {
next[j - 1] = next[j - 2]
+ (next[j - 2] - prev_rows[j - 2]) / (four.powi(j as i32 - 1) - N::one());
}
h *= half_real;
prev_rows[..i].clone_from_slice(&next[..i]);
}
Ok(prev_rows[n - 1])
}