GR series methods for Bessel functions

[Joel-Dahne]

I have been looking at adding adding support for Bessel functions to the gr_series module. It seems like the interface has stabilized enough for this to be reasonable for me to try.

This can implemented by computing the first two terms in the expansion explicitly and then using the recurrence that can be derived from the differential equation. I have already implemented a version of this in ArbExtras.jl (that one computes the first four terms explicitly).

The ODE has a singular point at zero, which means that a different recurrence is needed there. One could just implement this as a special case. However, that will not directly work for balls overlapping zero but not identically zero. I don't know if there is any good way around this? When I have encountered this case in some of my work I have just explicitly computed the formula for the required derivatives in terms of Bessel functions (which do support evaluation around zero). I don't think handling this case is necessary for an initial implementation, but it is a use case I have encountered in practice.

Similar recurrences could of course be implemented for any special function satisfying a linear ODE. It seems like this could be handled automatically by some of the functionality discussed in #1881? However, I don't know what the timeline for anything like that looks like. Would a specialized implementation still make sense?

[mezzarobba]

FWIW, I have some reasonably general code for computing series solutions at regular singular points in #2186. However:

• At the moment, it only works over acb. There is no real obstacle to making it more generic, but that's too much for my poor little brain to do while also working on other issues.
• There is no high-level function that returns the series coefficients in the form you would expect (only low-level functions that store them in some internal format).

The second point is on my todo-list for this week—which does not mean I will actually manage to work on it though!

[fredrik-johansson]

Good idea. I guess these are exactly the recurrences listed here: https://fungrim.org/topic/Recurrence_relations_for_Bessel_functions/

For the location of the actual algorithms I'm inclined to have _jet variants in the gr_special module, so there could be gr_bessel_j_jet(gr_ptr res, gr_srcptr nu, gr_srcptr z, slong len, gr_ctx_t ctx), etc.

Then there could be quite generic wrappers of such algorithms for gr_poly, gr_series and gr_series_mod. Compare with the existing gr_mat_func_jordan which takes such a jet function as input. Indeed -- once you have the Bessel J jet, you could immediately try it on matrices too via gr_mat_func_param_jordan (you just need a small wrapper to call with nu in the right parameter slot -- and similarly for other multivariate functions).

I don't know how to treat intervals around zero easily. I think this can't be done purely "algebraically" but needs some analysis. Over $\mathbb{R}$ and $\mathbb{C}$ the most straightforward solution is probably to implement special-purpose algorithms using arb and acb, convert input to these types, and convert back. More interesting would be to have some kind of construction for evaluating analytic functions over generic $\mathbb{R}$-algebras using something like a truncated power series + norm bound of geometric tail (here you'd use the truncated power series at 0). But I'm not really sure what that would look like.