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