Wishlist for D-finite/holonomic functions

[rburing]

Based on discussions with @mezzarobba and @fredrik-johansson:

• A module for generic Ore polynomials: gr_ore_poly

  • differential operators in d/dz and in z*d/dz, with gr_polys as coefficients
  • difference operators (forward and backward shift)
  • q-difference operators
  • Conversion functions, including also Ore polynomial to (companion) matrix of polynomials.

• Implementation of numerical algorithms: acb_holonomic or acb_ode

  • Functions can take Ore polynomials or matrices of polynomials, whichever is most natural for the algorithm.
  • Can have specializations for high precision and low precision.
  • For many applications you know the path you want to take for analytic continuation, so it can be an input.
  • Subdivision of a path should ideally be done automatically; often you notice the behavior of a solution as you expand, so you need to subdivide more.
  • For singular case, some functions need additional input to specify the structure of the local solutions.
  • Evaluation of holonomic function at a point.
  • Implementation of the bit-burst algorithm, e.g. for evaluation of Bessel functions.

• Generic implementations of D-finite function expansion/evaluation can go in gr_poly

  • Ideally expansion and evaluation would share as much code as possible.
    • Should support evaluating at a point directly (or e.g. computing a single coefficient), without keeping the whole polynomial in memory.
    • Dynamically decide when to stop expanding a power series, to get the needed precision to evaluate at a point.
  • Newton iteration
  • Divide and conquer
  • Naive algorithm using the recurrence directly
  • Faster algorithms for higher precision: fast computation of a single coefficient, or a single partial sum.

[rburing]

To give an update on this, at https://gitlab.inria.fr/ricardo-thomas.buring/d-finite-fun I am working on generic implementations of Newton iteration and divide-and-conquer for first-order systems (to obtain a basis of solutions), that I plan to contribute to FLINT. When the first-order system comes from a higher-order scalar equation, some matrices have a special form that allows some optimizations. I am currently working on making it possible to continue expanding a series further without re-doing a lot of the work that was already done (which will be exposed in some function starting with _), so this can be used for dynamically deciding when to stop expanding (useful for evaluation / analytic continuation), and/or for switching strategies (e.g. if you want to add just one term or a few terms). Since this (as well as other possible optimizations) might affect the API I have not made the PR to FLINT yet.

cc @lairez

[rburing]

Here is a proposal for a part of the "high-level" API, biased by what I'm working on. I would welcome some feedback on this, what you would change/rename/remove/add, if anything (when you have time).

In gr_ore_poly.h (depending on #2299):

int gr_ore_poly_fundamental_matrix_nonsingular_newton(gr_mat_t res, const gr_ore_poly_t poly, gr_srcptr point, slong len, int flags, gr_ctx_t sol_poly_ctx, gr_ore_poly_ctx_t ctx);
int gr_ore_poly_fundamental_matrix_nonsingular_divide_and_conquer(gr_mat_t res, const gr_ore_poly_t poly, gr_srcptr point, slong len, int flags, gr_ctx_t sol_poly_ctx, gr_ore_poly_ctx_t ctx);

These would call slightly lower level methods with the initial conditions Y0 set to the identity matrix, and also pass the flags to the lower level methods.

In gr_mat.h (which has the precedent gr_mat_gr_poly_evaluate), functions taking matrices of univariate polynomials:

int gr_mat_gr_poly_solve_first_order_homogeneous_ode_nonsingular_newton(gr_mat_t res, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_poly_ctx, gr_ctx_t sol_poly_ctx);
int gr_mat_gr_poly_solve_first_order_homogeneous_ode_nonsingular_divide_and_conquer(gr_mat_t res, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_poly_ctx, gr_ctx_t sol_poly_ctx);

These would produce matrices of truncated series. One of the flags would be to indicate that A is a companion matrix.

In gr_poly.h, analogous functions taking univariate polynomials with matrix coefficients:

int gr_poly_gr_mat_solve_first_order_homogeneous_ode_nonsingular_newton(gr_poly_t res, const gr_poly_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_mat_ctx, gr_ctx_t sol_mat_ctx);
int gr_poly_gr_mat_solve_first_order_homogeneous_ode_nonsingular_divide_and_conquer(gr_poly_t res, const gr_poly_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_mat_ctx, gr_ctx_t sol_mat_ctx);

In addition to the above, there would be functions for evaluation at (a list of) points / transition matrices, and also a high-level API for analytic continuation along a path (that can be shared by different methods).

[mezzarobba]

In gr_ore_poly.h (depending on #2299):

int gr_ore_poly_fundamental_matrix_nonsingular_newton(gr_mat_t res, const gr_ore_poly_t poly, gr_srcptr point, slong len, int flags, gr_ctx_t sol_poly_ctx, gr_ore_poly_ctx_t ctx);
int gr_ore_poly_fundamental_matrix_nonsingular_divide_and_conquer(gr_mat_t res, const gr_ore_poly_t poly, gr_srcptr point, slong len, int flags, gr_ctx_t sol_poly_ctx, gr_ore_poly_ctx_t ctx);

From the arguments I assume this is for computing partial sums evaluated at a point of fundamental solutions. Is that correct? Maybe this should somehow appear in the name.

Or is it for truncated solutions? But then, what is the role of point? If it is an expansion point, I think it is probably better to have a separate function for shifting the operator and then always expand at 0.

Also, I am not sure I understand what context corresponds to what object(s). Is the idea that sol_poly_ctx is the parent of res and point? But then what does its name refer to?

Maybe drop nonsingular_ from the names and return GR_UNABLE in singular cases? (The same signature can work at singular points, at least in some cases.)

In gr_mat.h (which has the precedent gr_mat_gr_poly_evaluate), functions taking matrices of univariate polynomials:

int gr_mat_gr_poly_solve_first_order_homogeneous_ode_nonsingular_newton(gr_mat_t res, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_poly_ctx, gr_ctx_t sol_poly_ctx);
int gr_mat_gr_poly_solve_first_order_homogeneous_ode_nonsingular_divide_and_conquer(gr_mat_t res, const gr_mat_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_poly_ctx, gr_ctx_t sol_poly_ctx);

I would call these gr_mat_gr_poly_solve_lode_newton and gr_mat_gr_poly_solve_lode_divconquer.

Same questions/remarks as above otherwise.

These would produce matrices of truncated series. One of the flags would be to indicate that A is a companion matrix.

Why not detect it at run time?

In gr_poly.h, analogous functions taking univariate polynomials with matrix coefficients:

int gr_poly_gr_mat_solve_first_order_homogeneous_ode_nonsingular_newton(gr_poly_t res, const gr_poly_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_mat_ctx, gr_ctx_t sol_mat_ctx);
int gr_poly_gr_mat_solve_first_order_homogeneous_ode_nonsingular_divide_and_conquer(gr_poly_t res, const gr_poly_t A_numerator, const gr_poly_t A_denominator, const gr_mat_t Y0, gr_srcptr z0, slong len, int flags, gr_ctx_t A_mat_ctx, gr_ctx_t sol_mat_ctx);

If these are meant to behave the same as the gr_mat_gr_poly functions above, I am not sure it is worth implementing both.

In addition to the above, there would be functions for evaluation at (a list of) points / transition matrices, and also a high-level API for analytic continuation along a path (that can be shared by different methods).

Maybe this part should go in a separate ode (or gr_ode) support module? (Edit: or simply in acb_ode, since the use cases for analytic continuation paths outside the context of acb are a bit limited.)