ported/bessel

fn abs(self) -> Self {
    D(self.0.abs())
}

fn sqrt(self) -> Self {
    D::new(self.0.sqrt())
}

fn exp<const S: bool>(x: D<S>) -> D<S> {
    D::new(excel_exp(x.0))
}

/// `c[0] + y (c[1] + y (c[2] + ...))`.
fn horner<const S: bool>(c: &[f64], y: D<S>) -> D<S> {
    let mut r = D(c[c.len() - 1]);
    for &k in c[..c.len() - 1].iter().rev() {
        r = D(k) + y * r;
    }
    r
}

fn i1<const S: bool>(x: D<S>) -> D<S> {
    let ax = x.abs();
    if ax.0 < 3.75 {
        let t = x / D(3.75);
        x * horner(&I1_SMALL, t * t)
    } else {
        let r = exp(ax) / ax.sqrt() * horner(&I1_LARGE, D(3.75) / ax);
        if x.0 < 0.0 { -r } else { r }
    }
}

fn admit_phase(x: D, phase: f64) -> Result<(), WorksheetErrorCode> {
    // R16/R19 Value2 observations, including adjacent input pairs at eleven
    // function/order/sign combinations. Compare the stored published phase;
    // no answer, fitted numerical coefficient or argument-specific selector
    // enters the calculation. Higher forward orders inherit the J0/Y0 bound.
    if x.0.abs() - phase > ((1_u32 << 27) - 1) as f64 {
        Err(WorksheetErrorCode::Num)
    } else {
        Ok(())
    }
}

/// The order as Excel reads it: truncated; #NUM! when negative or above [`MAX_ORDER`].
fn order(n: f64) -> Result<u32, WorksheetErrorCode> {
    if !(n >= 0.0) {
        return Err(WorksheetErrorCode::Num);
    }
    let n = n.trunc();
    if n > MAX_ORDER {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(n as u32)
}

/// The forward recurrence from orders 0 and 1 up to `n`: `next = j (2/x) c - previous` for J and
/// Y, `next = previous + j (2/x) c` for K.
fn forward<const S: bool>(x: D<S>, n: u32, c0: D<S>, c1: D<S>, plus: bool) -> D<S> {
    let tox = D(2.0) / x;
    let (mut before, mut c) = (c0, c1);
    for j in 1..n {
        let step = D(j as f64) * tox * c;
        let next = if plus { before + step } else { step - before };
        before = c;
        c = next;
    }
    c
}

fn odd_sign(x: D, n: u32, v: D) -> D {
    if x.0 < 0.0 && n % 2 == 1 { -v } else { v }
}

fn finish(v: R) -> Result<f64, WorksheetErrorCode> {
    let v = v?.0;
    if v.is_finite() { Ok(v) } else { Err(WorksheetErrorCode::Num) }
}