BESSELI

What it does

BESSELI(x, n)

BESSELI(x, n) returns I_n(x), the modified Bessel function of the first kind, of integer order n at the real point x. Microsoft describes it as the Bessel function evaluated for a purely imaginary argument, which is the standard relation

I_n(x)  =  i^(−n) · J_n(i·x)

I_n is the solution that stays finite at the origin of the modified Bessel equation

x² w″ + x w′ − (x² + n²) w  =  0

It has an everywhere-convergent power series and, for integer order, a short integral form:

I_n(x)  =  Σ_{k≥0}  (x/2)^(2k+n) / ( k! · (k+n)! )
I_n(x)  =  (1/π) ∫₀^π  e^(x cos θ) cos(nθ) dθ

Three facts about its shape matter below.

• Growth. For large x, I_n(x) ≈ e^x / √(2πx), and at that scale every order looks alike: the order only appears as a factor close to e^(−n²/(2x)). I_0 is still representable in binary64 up to about x = 713.9.
• The origin. I_0(0) = 1 and I_n(0) = 0 for n ≥ 1, with I_n(x) ≈ (x/2)^n / n! for small x.
• Parity. I_n(−x) = (−1)^n · I_n(x): even orders are even functions and odd orders odd.

For integer order the function is real on the whole line, positive for x > 0, and has no poles and no branch points. Everything difficult about BESSELI in Excel comes from how it is approximated.

Arguments

Truncation and the negative test. Microsoft states both that a non-integer n is truncated and that n < 0 is #NUM!, and does not say which comes first. The recorded answers do: the sign test sees the value before truncation. All 64 recorded calls with n strictly between −1 and 0 returned #NUM! (Excel 16.0 build 20430, Windows 64-bit, CPU not recorded), although truncating any of them gives 0. Positive fractions behave as documented: BESSELI(-1, 0.9) and BESSELI(-1, 0) read the same.

BESSELI is not a function of real order. The notes below compute their "correct value" at the integer part of n for that reason, and BESSELI(x, 3.99) is order 3.

The slots are numeric. Beyond Microsoft's statement that a nonnumeric argument gives #VALUE!, no reading of text, a logical, an empty cell or an array in either slot is on record for this function; see What has not been checked.

Try it

Quirks

Result and edge cases

Returns Number. Unless stated otherwise the behaviour below is read from recorded Excel answers (builds 20430 and 20527, Windows 64-bit).

• x = 0 gives 1 at order 0 and 0 at every higher order read (1, 2 and 10), as the mathematics requires.
• Negative x is admitted. Even orders return the value at |x|, and odd orders return it negated, which was read at ±1 and at ±709.7.
• The documented example BESSELI(1.5, 1) is printed by Microsoft as 0.981666428. Nine digits is about what the function offers: the answers are accurate to roughly seven or eight significant figures at ordinary arguments, not to fifteen. BESSELI(-1, 0) is the plainest example in the notes below.
• |x| just above 709.7 is #NUM! at every order read, although the true value is far inside the binary64 range there. The notes below record the boundary to the adjacent double. The failure is abrupt: BESSELI(709.7, 0) is a number near 2.48E306, and the next double up is an error.
• Orders 2 and above can be badly wrong at large x. At x = 709.7, order 3 is 53% below the true value. This is a real wrong answer with no error raised, recorded in the notes below and discussed under Numerical notes.
• Very high orders can collapse to 0. At x = ±709.7, orders 1568 to 1572 all read 0 (build 20527), although the true values, from about 1.5E-299 down to 3.1E-302, are ordinary binary64 numbers. Order 10000 also read 0, which there is the right rounded answer: the true value is near 2.9E-10154.

Errors

As documented on Microsoft's page:

As recorded, in addition to the documented #NUM! for n ≤ −1:

How to implement

Numerical notes

These are low-order polynomials, by design. Orders 0 and 1 are the polynomial approximations printed in Abramowitz & Stegun as 9.8.1 to 9.8.4: one polynomial in (x/3.75)² for |x| < 3.75 and one in 3.75/|x|, multiplied by e^|x|/√|x|, from 3.75 up. The coefficients carry seven or eight decimal places, and the printed error bounds reach 1.9E-7 for I_0 and 2.2E-7 for I_1 on each form's own scale. In the large-x form that scale is a factor near 0.4, so the relative error of I_0 can reach about 2.5 times the bound; the largest recorded at order 0 is 4.7E-7, at x = 100. The notes below measure the result: correctly rounded on 2 of 1,163 recorded rows. Treat the answer as a seven-digit figure.

There is a seam at 3.75, and Excel shows it. Two different polynomials meet there and they do not agree exactly. Build 20527 reads I_0 at adjacent doubles: the double below 3.75 is on the small polynomial, 3.75 itself is on the large one, and the answer steps up by 2.1E-8 relative between them. A numerical derivative taken across 3.75 is not to be trusted. The printed coefficients for I_1 imply a smaller step, about 1.1E-8, which has not been read.

Orders 2 and above come from a backward recurrence. Stepping upward in n from I_0 and I_1 is unstable for this function, so the standard method, due to Miller, starts high with an arbitrary value, runs the recurrence downward where it is stable, and rescales the whole sequence at the end by the known I_0(x). The notes below give Excel's exact staging. Two consequences are worth knowing:

1. The error of I_0 is inherited. The normalisation uses the polynomial value of I_0, so a healthy high-order answer carries order 0's relative error. The recorded rows show it plainly: at x = 12.51 orders 0, 2 and 7 are all low by 3.34E-8, and at x = 709.7 orders 0 and 100 are both low by 1.45E-7. The step at the seam reaches orders 2 and 17 by the same route.
2. The starting point ignores x. The recurrence starts at j = 2(n + ⌊√(40n)⌋), a function of the order alone. That start is high enough while x is small compared with it. When x is much larger, the sequence has too few steps to settle, and the answer is simply wrong.

How wrong, and where. The recorded values do not change smoothly with the order. At x near 709, order 2 is high by 0.20%, order 3 is low by 53%, order 10 is high by 0.08%, and order 100 is as good as order 0. The 1,236-row recording behind the notes stops at x = 50, apart from six rows at x = 100. Its worst row at x ≤ 50 is about 9E-6 relative, at order 2, and the 5.0E-4 quoted in the notes below is BESSELI(100, 2). So order 2 degrades steadily with x over the range read. Order 3 is about 1.6E-6 off near x = 50 and 53% low at 709, and the path between them has not been read.

The large-x failure is the one defect on this page that matters to an ordinary user: it raises no error, and the wrong answer looks entirely plausible.

The collapse to 0 has a likely cause in the same recurrence. It computes the ratio I_n/I_0 and only then multiplies by I_0. At x = 709.7 and high order that ratio is below the binary64 range although the product is not, so the ratio underflows to 0 and the multiplication cannot recover it. A plain binary64 transcription of the published body does this somewhere between orders 1100 and 1200 at this x. That is a computation on the published body, and the cause is an inference: Excel itself has been read only at order 100 and then from order 1568.

The cut at 709.7. e^709.7 is about 1.65E308, just inside the binary64 range, and the exponential itself overflows near 709.78. The polynomial's e^|x| factor is therefore the natural reason for a guard at about this size, but the guard then gives up about four units of argument over which I_0 is still finite. The reason is an inference; only the boundary is measured.

Rust code

60 besseli "BESSELI" (2, 2) [V "x", V "n"] |s, a, sv| { bessel(s, a, k::bessel::besseli_kernel, false) }

/// `BESSELI(x, n)`.
pub fn besseli_kernel(x: f64, n: f64) -> Result<f64, WorksheetErrorCode> {
    // Use one admission policy across orders, extending the adjacent-input
    // measurements at orders 0, 1, 2, 3, 10 and 100. Higher orders still need
    // Excel measurements; a representable series value alone cannot settle
    // their admission. EXP itself admits larger x.
    if !x.is_finite() || x.abs() > 709.7 { return Err(WorksheetErrorCode::Num); }
    let n = order(n)?;
    let x = D::new(x);
    finish(Ok(match n {
        0 => D::new(super::bessel_rational::i0(x.0, &I0_SMALL, &I0_LARGE)),
        1 => i1(x),
        _ if x.0 == 0.0 => D(0.0),
        _ => {
            let ax = x.abs();
            let normalizer = super::bessel_rational::i0(ax.0, &I0_SMALL, &I0_LARGE);
            odd_sign(x, n, D::new(super::bessel_rational::in_(n, ax.0, normalizer)))
        },
    }))
}

pub(super) fn i0(x: f64, small: &[f64], large: &[f64]) -> f64 {
    let magnitude = x.abs();
    if magnitude < 3.75 {
        let reduced = excel_x87_div(x, 3.75);
        let square = excel_x87_mul(reduced, reduced);
        polynomial_with(square, small, excel_x87_mul, excel_x87_add)
    } else {
        let reduced = excel_x87_div(3.75, magnitude);
        let scale = excel_x87_div(excel_exp(magnitude), excel_sqrt_staged(magnitude));
        excel_x87_mul(scale, polynomial_with(reduced, large, excel_x87_mul, excel_x87_add))
    }
}

/// The published higher-order body; its caller owns order admission and parity.
pub(super) fn in_(order: u32, magnitude: f64, normalizer: f64) -> f64 {
    if magnitude == 0.0 {
        return 0.0;
    }
    let twice_over_x = excel_x87_div(2.0, magnitude);
    let mut next = 0.0;
    let mut current = 1.0;
    let mut answer = 0.0;
    let start = 2 * (order + (40.0 * f64::from(order)).sqrt() as u32);
    for j in (1..=start).rev() {
        let coefficient = excel_x87_mul(f64::from(j), twice_over_x);
        let previous = excel_x87_add(next, excel_x87_mul(coefficient, current));
        next = current;
        current = previous;
        if current.abs() > 1e10 {
            answer = excel_x87_mul(answer, 1e-10);
            current = excel_x87_mul(current, 1e-10);
            next = excel_x87_mul(next, 1e-10);
        }
        if j == order {
            answer = next;
        }
    }
    excel_x87_mul(answer, excel_x87_div(normalizer, current))
}

No implementation has been published with a Handbook test suite for this function.

Version differences

The recorded environments below distinguish a reading from a same-input comparison. Counts describe captured inputs in each environment; overlapping captures and comparisons are not totals of unique worksheet calls.

The record contains no same-input comparison between these environments. Reading different inputs on different builds does not establish either a version change or cross-build agreement.

Recorded environments and comparison details.

Alternatives

Relationships

• BESSELK — the other modified Bessel function, decaying where this one grows. The two are tied by an exact cross-product, I_n(x)·K_(n+1)(x) + I_(n+1)(x)·K_n(x) = 1/x, which needs no oracle and is a useful sanity test on any pair of answers. BESSELK shares this page's accuracy note but none of the large-x trouble below, because it steps upward.
• EXP — every answer at |x| ≥ 3.75 carries a factor e^|x|, which is the likely reason for the cut near 709.7.
• The recurrence. I_(n−1)(x) − I_(n+1)(x) = (2n/x)·I_n(x) links neighbouring orders. It gives a workaround for the large-x defect at low order: BESSELI(x,0) - 2/x*BESSELI(x,1) is I_2(x) from the two orders that do not use the recurrence. That identity is mathematics; it has not been read in Excel as a repair.
• Where readers meet it. I_0(κ) is the normalising constant of the von Mises distribution, and I_0 and I_1 appear in the Rician distribution and in radially symmetric diffusion problems.

Named variants

No named variants are recorded yet.