From a8f91a7527f9df770377ff4a11275a0eb61c34ab Mon Sep 17 00:00:00 2001 From: Marc Mezzarobba Date: Wed, 6 May 2026 12:22:49 +0200 Subject: [PATCH] =?UTF-8?q?[WIP]=20Solutions=20of=20linear=20ODEs=20with?= =?UTF-8?q?=20coeffs=20in=20=E2=84=82[x]?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile.in | 2 +- dev/check_examples.sh | 3 + doc/source/acb.rst | 8 + doc/source/acb_ode.rst | 654 ++++++++++++++++++ doc/source/acb_poly.rst | 22 + doc/source/index.rst | 1 + examples/ode_basic.c | 69 ++ .../ode_fundamental_matrix_manual_exponents.c | 144 ++++ examples/ode_low_level.c | 197 ++++++ src/acb.h | 2 + src/acb/vec_get_mag.c | 34 + src/acb/vec_get_mag_lower.c | 34 + src/acb_ode.h | 489 +++++++++++++ src/acb_ode/apply_diffop_basecase.c | 167 +++++ src/acb_ode/apply_diffop_polmul.c | 106 +++ src/acb_ode/bound_init_clear.c | 32 + src/acb_ode/bound_precompute.c | 194 ++++++ src/acb_ode/bound_precompute_group.c | 41 ++ src/acb_ode/bound_precompute_integrand.c | 44 ++ src/acb_ode/bound_rat_ordinary_vec.c | 169 +++++ src/acb_ode/bound_rat_ref_vec.c | 66 ++ src/acb_ode/bound_rat_vec.c | 54 ++ src/acb_ode/choose_prec.c | 85 +++ src/acb_ode/cvest.c | 73 ++ src/acb_ode/exponents.c | 252 +++++++ src/acb_ode/fundamental_matrix.c | 329 +++++++++ src/acb_ode/ind_lbound.c | 160 +++++ src/acb_ode/inlines.c | 3 + src/acb_ode/poly_negdivrevhigh.c | 64 ++ src/acb_ode/poly_taylor_shift_aps_trunc.c | 36 + src/acb_ode/sol.c | 117 ++++ src/acb_ode/sol_estimate.c | 56 ++ src/acb_ode/sol_jet.c | 96 +++ src/acb_ode/solution_growth.c | 89 +++ src/acb_ode/stairs.c | 66 ++ src/acb_ode/sum_context.c | 170 +++++ src/acb_ode/sum_divconquer.c | 46 ++ src/acb_ode/sum_done.c | 295 ++++++++ src/acb_ode/sum_fix.c | 41 ++ src/acb_ode/sum_forward_1.c | 206 ++++++ src/acb_ode/sum_forward_divconquer.c | 169 +++++ src/acb_ode/tail_bound_jet.c | 272 ++++++++ src/acb_ode/test/main.c | 12 + src/acb_ode/test/t-apply_diffop.c | 98 +++ src/acb_ode/test/t-fundamental_matrix.c | 258 +++++++ src/acb_ode/test/t-sum_divconquer.c | 124 ++++ src/acb_poly.h | 13 + src/acb_poly/gr_stuff.c | 70 ++ src/acb_poly/vec.c | 89 +++ 49 files changed, 5820 insertions(+), 1 deletion(-) create mode 100644 doc/source/acb_ode.rst create mode 100644 examples/ode_basic.c create mode 100644 examples/ode_fundamental_matrix_manual_exponents.c create mode 100644 examples/ode_low_level.c create mode 100644 src/acb/vec_get_mag.c create mode 100644 src/acb/vec_get_mag_lower.c create mode 100644 src/acb_ode.h create mode 100644 src/acb_ode/apply_diffop_basecase.c create mode 100644 src/acb_ode/apply_diffop_polmul.c create mode 100644 src/acb_ode/bound_init_clear.c create mode 100644 src/acb_ode/bound_precompute.c create mode 100644 src/acb_ode/bound_precompute_group.c create mode 100644 src/acb_ode/bound_precompute_integrand.c create mode 100644 src/acb_ode/bound_rat_ordinary_vec.c create mode 100644 src/acb_ode/bound_rat_ref_vec.c create mode 100644 src/acb_ode/bound_rat_vec.c create mode 100644 src/acb_ode/choose_prec.c create mode 100644 src/acb_ode/cvest.c create mode 100644 src/acb_ode/exponents.c create mode 100644 src/acb_ode/fundamental_matrix.c create mode 100644 src/acb_ode/ind_lbound.c create mode 100644 src/acb_ode/inlines.c create mode 100644 src/acb_ode/poly_negdivrevhigh.c create mode 100644 src/acb_ode/poly_taylor_shift_aps_trunc.c create mode 100644 src/acb_ode/sol.c create mode 100644 src/acb_ode/sol_estimate.c create mode 100644 src/acb_ode/sol_jet.c create mode 100644 src/acb_ode/solution_growth.c create mode 100644 src/acb_ode/stairs.c create mode 100644 src/acb_ode/sum_context.c create mode 100644 src/acb_ode/sum_divconquer.c create mode 100644 src/acb_ode/sum_done.c create mode 100644 src/acb_ode/sum_fix.c create mode 100644 src/acb_ode/sum_forward_1.c create mode 100644 src/acb_ode/sum_forward_divconquer.c create mode 100644 src/acb_ode/tail_bound_jet.c create mode 100644 src/acb_ode/test/main.c create mode 100644 src/acb_ode/test/t-apply_diffop.c create mode 100644 src/acb_ode/test/t-fundamental_matrix.c create mode 100644 src/acb_ode/test/t-sum_divconquer.c create mode 100644 src/acb_poly/gr_stuff.c create mode 100644 src/acb_poly/vec.c diff --git a/Makefile.in b/Makefile.in index d0d990c632..c5eff3e9b1 100644 --- a/Makefile.in +++ b/Makefile.in @@ -217,7 +217,7 @@ HEADER_DIRS := \ double_interval dlog \ arb_fmpz_poly arb_fpwrap \ acb_dft acb_elliptic acb_modular \ - acb_dirichlet acb_theta \ + acb_dirichlet acb_theta acb_ode \ dirichlet bernoulli hypgeom \ bool_mat partitions \ \ diff --git a/dev/check_examples.sh b/dev/check_examples.sh index bb548facd6..e48a1e2bc0 100755 --- a/dev/check_examples.sh +++ b/dev/check_examples.sh @@ -169,6 +169,9 @@ then elif test "$1" = "hilbert_matrix_ca"; then echo "hilbert_matrix_ca....SKIPPED" +elif test "$1" = "ode"; +then + echo "ode....SKIPPED" elif test "$1" = "huge_expr"; then echo "huge_expr....SKIPPED" diff --git a/doc/source/acb.rst b/doc/source/acb.rst index a4d0f54136..8aec7e2217 100644 --- a/doc/source/acb.rst +++ b/doc/source/acb.rst @@ -1206,6 +1206,14 @@ Vector functions Swaps the entries of *vec1* and *vec2*. +.. function:: void _acb_vec_get_mag(mag_t bound, acb_srcptr vec, slong len) + + Sets *bound* to an upper bound for the entries in *vec*. + +.. function:: void _acb_vec_get_mag_lower(mag_t bound, acb_srcptr vec, slong len) + + Sets *bound* to an lower bound for the absolute values of the entries in *vec*. + .. function:: void _acb_vec_get_real(arb_ptr re, acb_srcptr vec, slong len) .. function:: void _acb_vec_get_imag(arb_ptr im, acb_srcptr vec, slong len) diff --git a/doc/source/acb_ode.rst b/doc/source/acb_ode.rst new file mode 100644 index 0000000000..84089cd8d2 --- /dev/null +++ b/doc/source/acb_ode.rst @@ -0,0 +1,654 @@ +.. _acb-ode: + +**acb_ode.h** -- linear differential equations with polynomial coefficients +================================================================================== + +.. note:: + + This module is experimental; the API is subject to change. + +XXX discuss the general form of logarithmic series solutions etc. + +Fundamental matrices +------------------------------------------------------------------------------- + +See also the :ref:`usage example ` below. + +.. type:: acb_ode_basis_t + + Represents one of the following standard bases of solutions: + + .. macro:: ACB_ODE_BASIS_CASCADE + + XXX + + .. macro:: ACB_ODE_BASIS_ECHELON + + XXX + + + In both cases, the basis is the concatenation of blocks + corresponding to groups of solutions. + Within each block, the solutions are sorted by increasing shift and then + by decreasing degree in `\log(x)` of the coefficient of `x` of minimal + valuation. + + XXX Specify the order of the blocks (as in *expos*?...), + provide more flexibility (individual solutions ordered by asymptotic + behavior, breaking the block structure). + +.. function:: int acb_ode_fundamental_matrix(acb_mat_t mat, const gr_ore_poly_t dop, gr_ore_poly_ctx_t Dop, const acb_ode_exponents_t expos, acb_srcptr lcrt, acb_srcptr pt, acb_ode_basis_t basis, slong prec) + + Sets *mat* to the value at *pt* of the fundamental matrix in the + basis *basis* of the differential operator *dop*. + + This function is currently limited to operators *dop* of algebra type + :macro:`GR_ORE_POLY_EULER_DERIVATIVE`. + + The origin may be an ordinary point or a regular singular point of *dop*. + The point *pt* must be closer to the origin than any nonzero singular + points of *dop*. + + The *expos* and *lcrt* parameters respectively denote the exponent structure + of *dop* at the origin and the roots of its leading coefficient. + They may be set to ``NULL`` when the corresponding quantities can be + automatically determined from *dop* using :func:`acb_ode_exponents`, + resp. :func:`gr_poly_roots_other`. + (This is typically the case over an exact field of constant supporting a + complete equality test.) + When specified, *expos* must be set to the exponent structure at the origin + in the format returned by :func:`acb_ode_exponents` + (see :ref:`example `). + When specified, *lcrt* must be a vector of enclosures of the roots of the + leading coefficient of *dop*, with multiple roots repeated. + + The return value is *gr*-style status code. A return value other that + ``GR_SUCCESS`` indicates a failure in the algebraic part of the computation. + It may happen, however, that the return value is ``GR_SUCCESS`` but the + output matrices are indeterminate. + +.. function:: void _acb_ode_fundamental_matrix_vec(acb_mat_struct * mat, const acb_poly_struct * dop, slong dop_len, const acb_ode_exponents_t expos, acb_srcptr lcrt, acb_srcptr pts, slong npts, acb_ode_basis_t basis, acb_ode_sum_worker_t sum_worker, slong prec) + int acb_ode_fundamental_matrix_vec(acb_mat_struct * mat, const gr_ore_poly_t dop, gr_ore_poly_ctx_t Dop, const acb_ode_exponents_t expos, acb_srcptr lcrt, acb_srcptr pts, slong npts, acb_ode_basis_t basis, slong prec) + + Like :func:`acb_ode_fundamental_matrix` but evaluates the fundamental + matrix at multiple points. + + The underscore method takes the differential operator as an array of + ``acb_poly`` coefficients instead of an Ore polynomial and assumes that it + has order at least one. It also requires the *expos* and *lcrt* parameters + to be specified, Additionally, it accepts an additional parameter + *sum_worker* that specifies the summation algorithm to be used. + +Solution groups and exponents +------------------------------------------------------------------------------- + +See below for an :ref:`example ` of how to use these +functions in combination with :func:`acb_ode_fundamental_matrix`. + +XXX define + +.. type:: acb_ode_shift_struct + acb_ode_shift_t + acb_ode_group_struct + acb_ode_group_t + acb_ode_exponents_struct + acb_ode_exponents_t + + Types used to represent the structure of a local basis of solutions, that + is, the roots of the indicial polynomial and the integer differences between + these roots. + + An *acb_ode_exponents_struct* contains an array *groups* of length *ngroups* + of *acb_ode_group_struct*. + Each *acb_ode_group_struct* contains an exponent *leader* of type + :type:`acb_t` and an array of length *nshifts* of *acb_ode_shift_struct*. + An *acb_ode_shift_struct* is a pair of an integer shift *n* and a + multiplicity *mult*. + The shifts are sorted in increasing order. + +.. function:: void acb_ode_group_init(acb_ode_group_t group, slong len) + void acb_ode_group_clear(acb_ode_group_t group) + void acb_ode_exponents_init(acb_ode_exponents_t expos) + void acb_ode_exponents_clear(acb_ode_exponents_t expos) + +.. function:: void acb_ode_group_set(acb_ode_group_t dest, const acb_ode_group_t src) + +.. function:: void acb_ode_exponents_ordinary(acb_ode_exponents_t expos, slong order) + + Sets *expos* to the exponent structure of an operator of order *order* at + an ordinary point. + +.. function:: int acb_ode_exponents(acb_ode_exponents_t expos, const gr_ore_poly_t dop, gr_ctx_t Dop, gr_ctx_t CC) + + Sets *expos* to the exponent structure of an operator of *dop* at the + origin, computed as elements of *CC*. The only supported *CC* are currently + complex ball fields. + +.. function:: slong acb_ode_group_length(const acb_ode_group_t group) + + Returns the number of solutions of *group*, i.e., the sum of + multiplicities of its shifts. + +.. function:: slong acb_ode_group_multiplicity (const acb_ode_group_t group, slong n) + + Returns the multiplicity of *n* as a shift in group *group*. + +.. function:: slong acb_ode_group_nlogs (const acb_ode_group_t group, slong n) + + Returns a bound on the length of the term of index *n* of a solution of the + group *group* when that term is viewed as a polynomial in log(*x*). + +Partial sums of logarithmic series solutions +------------------------------------------------------------------------------- + +This section describes low-level functions for computing series expansions +and/or values of partial sums or derivatives of partial sums of solutions. +These functions can handle several solutions from the same group (or possibly, +in some cases, from conjugate groups) in parallel, and can evaluate them at +zero, one, or more points. + +The general usage pattern is as follows: + +1. initialize a context object of type *acb_ode_sum_t*; +2. configure it by setting the differential operator + (*acb_ode_sum_set_diffop*), solution group (*acb_ode_sum_set_group* or + similar) and initial values (*acb_ode_sum_set_ini_\**) of interest, in + this order; +3. optionally set some evaluation points and other options; +4. call one or more summation functions such as *acb_ode_sum_divconquer* to + compute the partial sum; +5. extract the value and/or coefficients of the solution (XXX this step + currently requires using *sol_* functions). + +.. type:: acb_ode_sum_struct + acb_ode_sum_t + + Context object. + + The following options can be set by writing to the corresponding fields + after initializing the object: + + .. member:: ulong flags + + .. macro:: ACB_ODE_WANT_SERIES + + Keep series coefficients. + + .. macro:: ACB_ODE_APPROX + + If set, skip the computation of rigorous error bounds. + +.. function:: void acb_ode_sum_init(acb_ode_sum_t sum, slong dop_len, slong npts, slong nsols, slong nder) + + Initializes a context object suitable for computing the derivatives of order + 0 to *nder - 1* of *nsols* solutions of an equation *dop_len - 1* at + *npts* points. + When only the coefficients of the series are of interest, *npts* can be set + to zero. + + XXX Clarify the case *nder = 0*. It makes sense to have *npts = 0* with + *nder = 1* if one wants to compute the coefficients of a rigorous + approximation valid in a certain disk. For computing a fixed number of + coefficients without any kind of tail bound, though, we may want to allow + *nder = 0*. + +.. function:: void acb_ode_sum_clear(acb_ode_sum_t sum) + + Clears the context object. + +.. function:: void acb_ode_sum_set_diffop(acb_ode_sum_t sum, const acb_poly_t dop, slong dop_len, const mag_t cvrad) + + Sets the differential operator whose solutions are to be computed. + The operator must be nonzero. + The *cvrad* parameter is a lower bound on the radius of convergence of the + solutions and will limit the magnitude of the evaluation points. + +.. function:: void acb_ode_sum_set_group(acb_ode_sum_t sum, const acb_ode_group_t grp) + void acb_ode_sum_set_ordinary(acb_ode_sum_t sum) + + Sets the solution group. + +.. function:: void acb_ode_sum_set_ini_echelon(acb_ode_sum_t sum) + void acb_ode_sum_set_ini_highest(acb_ode_sum_t sum) + + Sets initial values defining the solutions to be computed. + The solution group must have been set before. + +.. function:: void acb_ode_sum_set_points(acb_ode_sum_t sum, acb_srcptr pts, slong npts) + + Sets the evaluation points. + +.. function:: slong acb_ode_sum_max_nlogs(acb_ode_sum_struct * sum) + +.. function:: slong acb_ode_sum_max_nlogs_xn(acb_ode_sum_struct * sum, slong off) + +.. function:: void acb_ode_sum_precompute(acb_ode_sum_t sum) + + Performs various precomputations. + +.. type:: acb_ode_sum_worker_t + + Pointer to a worker function + ``void worker(acb_ode_sum_t sum, slong nterms, slong prec)``. + + XXX maybe attach some data + +.. function:: void acb_ode_sum_divconquer(acb_ode_sum_t sum, slong nterms, slong prec) + + Sums the series expansions of the solutions contained in *sum*, updating the + *sum* object with their coefficients and/or partial sums. + The *nterms* parameter can be used to limit the number of terms (-1 means + unlimited). + This function computes the sum from scratch even if *sum* already contains a + partial sum from a previous run. + + XXX Clarify what happens with no evaluation points and how to evaluate the + *coefficients* to a given precision. + + Updates *sum* with the coefficients and/or partial sums of the solutions. + +The *_acb_ode_sum_forward_\** functions all assume that +:func:`acb_ode_sum_precompute` has been applies and that the *sum* object's +component *acb_ode_sol_t* all have space for the new coefficients. +Additionally, they may leave some of these components in a state that requires +:func:`_acb_ode_sum_fix` to be called before using other functions on them. + +.. function:: void _acb_ode_sum_forward_1(acb_ode_sum_struct * sum) + + Adds one term to *sum*, whose components are assumed to have space for the + update. + +.. function:: void _acb_ode_sum_forward_divconquer(XXX) + + Extends the sum until the partial sum of order *high* is reached, some error + tolerance based on *prec* is met, or the working precision and/or tail bound + quality are determined to be insufficient (whichever comes first). + + The algorithm computes strides of *stride* terms (checking for convergence + and making some internal adjustments after each stride) by chaining blocks + of *block_len* terms that are computed using a divide-and-conquer method. + It is assumed that *block_len* is at least the degree of the coefficients of + the differential operator written in terms the Euler derivation. + Currently it is also assumed that *stride* is a multiple of *block_len* and + that the starting index (value of ``sum->n``) of a top-level recursive call + is a multiple of *stride*. + +.. function:: void _acb_ode_sum_fix(acb_ode_sum_struct * sum) + + Restores some invariants of *sum* broken by the *sum_foward_\** functions. + After a sequence of calls to those functions, *_acb_ode_sum_fix* must be + called before accessing the computed sum. + +.. function:: int acb_ode_sum_done(acb_ode_sum_struct * sum, slong stride, slong prec) + + Returns nonzero when either the sum has converged (within a tolerance based + on *prec*) or the current working precision is insufficient to significantly + improve the current estimate. + Updates some heuristic convergence estimates stored in *sum*, assuming that + *stride* terms have been added since the last update. + When returning nonzero, also updates the rigorous tail bounds (unless the + *ACB_ODE_APPROX* flag is set in *sum*). + May be called on a *sum* object modified by *sum_foward_\** without first + calling *sum_fix*. + +Partial sums of individual solutions +------------------------------------------------------------------------------- + +.. type:: acb_ode_sol_struct + acb_ode_sol_t + + Represents a logarithmic series solution in the process of being computed. + Usually makes sense only as a component of an :type:`acb_ode_sum_t` object. + + XXX Initial values, `x`-slice of coefficients, partial sums at zero or more + points. + +.. function:: void acb_ode_sol_init(acb_ode_sol_t sol, slong nshifts, slong nlogs, slong npts, slong nder) + + Initializes an :type:`acb_ode_sol_t` object suitable for a solution + belonging to a solution group with *nshifts* shifts, involving `\log(x)` to + powers less than *nlogs*, and to be evaluated at *npts* points. + +.. function:: void acb_ode_sol_clear(acb_ode_sol_t sol) + + Clears *sol*. + +.. function:: void acb_ode_sol_fit_length(acb_ode_sol_t sol, slong len) + + Ensures that the *series* and *sums* fields of *sol* have space for *len* + coefficients. + +.. function:: void acb_ode_sol_zero(acb_ode_sol_t sol) + + Zeroes the series coefficients and partial sums (without touching the + initial conditions). + +.. function:: acb_ptr acb_ode_sol_sum_ptr(const acb_ode_sol_t sol, slong j, slong k, slong i) + + Returns a pointer to the coefficient of index *i* in the Taylor expansion + at the evaluation point of index *j* of the partial sum of the coefficient + of `\log(x)^k/k!`. + +.. function:: void _acb_ode_sol_jet(XXX) + void acb_ode_sol_jet(acb_poly_struct * val, acb_srcptr expo, const acb_ode_sol_t sol, slong j, acb_srcptr pt, slong nder, slong nfrobshifts, slong prec) + + Computes the jets of order *nder* at the evaluation point of index *pt* of + the Frobenius shifts of order zero to *nfrobshifts* of the solution *sol*. + The integer *j* is the index of the partial sum at *pt* in *sol*. + The Frobenius shift of order *p* of `y(x)` is defined as the series obtained + by replacing `\log(x)^k/k!` for `k ≥ p` by `\log(x)^{k-p}/(k-p)!` in the + series expansion and dropping the terms with `k < p`. + + The underscore version assumes *nder* > 0. + +.. function:: slong acb_ode_sol_estimate_terms(mag_t est, const acb_ode_sol_t sol, slong off, slong len, const mag_t radpow) + + Sets *est* to heuristic estimate of the magnitude of the coefficients at + positions *off* to *off + len* in the coefficient buffer of *sol*, + multiplied by *radpow*. + Returns an estimate of the relative accuracy in bits of these coefficients. + +.. function:: slong acb_ode_sol_estimate_sums(mag_t mag, mag_t rad, const acb_ode_sol_t sol) + + Sets *mag* to a heuristic estimate of the magnitude of the partial sums + encoded by *sol* and *rad* to an estimate of their radii. + Returns an estimate of the minimum relative accuracy in bits of the sums. + +Functions for working with differential operators +------------------------------------------------------------------------------- + +.. function:: void _acb_ode_apply_diffop_basecase_precomp(acb_poly_struct * g, slong goff, acb_srcptr weights, slong weights_nlogs, const acb_poly_struct * f, slong foff, slong flen, slong nlogs, slong start, slong len, slong prec) + void acb_ode_apply_diffop_basecase(acb_poly_struct * g, const acb_poly_struct * dop, slong dop_len, acb_srcptr expo, slong offset, const acb_poly_struct * f, slong nlogs, slong start, slong len, slong prec) + void _acb_ode_apply_diffop_polmul(acb_poly_struct * g, slong goff, const acb_poly_struct * dop, slong dop_len, acb_srcptr expo, slong offset, const acb_poly_struct * f, slong foff, slong flen, slong nlogs, slong start, slong len, slong prec) + void acb_ode_apply_diffop_polmul(acb_poly_struct * g, const acb_poly_struct * dop, slong dop_len, acb_srcptr expo, slong offset, const acb_poly_struct * f, slong nlogs, slong start, slong len, slong prec) + + Computes a slice of the coefficients of the image of a series by a + differential operator. + + XXX :: + + We are considering an operator + \sum_{i,j} a[i,j] x^j (x d/dx)^i + applied to a series + f = x^{expo + offset} \sum_{p,k} f[k,p] x^p log(x)^k/k!, + resulting in + g = x^{expo + offset} \sum{m1,k1} g[k1,m1] x^m1 log(x)^k1/k1!. + We want to compute the terms of index m1 = start + p1 of g for 0 <= p1 < len. + + The non-underscore versions represent *f* and *g* as arrays of polynomials + corresponding to the coefficients of `\log(x)^k/k!`. + The underscore versions read the coefficients of `f[k]` (resp. write those + of `g[k]`) starting at offset *offset* in the coefficient array of the + element of index `k` of the array *f* (resp. *g*) + + The *polmul* version uses fast polynomial multiplication. The *basecase* version XXX. The *precomp* version is similar to the *basecase* version but takes as input, instead of the differential operator, an array of weights that can be precomputed using :func:`_acb_ode_apply_diffop_basecase_weights`. + +.. function:: void _acb_ode_apply_diffop_basecase_weights(acb_ptr weights, const acb_poly_struct * dop, slong dop_len, acb_srcptr expo, slong offset, slong flen, slong nlogs, slong start, slong len, slong prec) + + Sets *weights* to an array of *len* × *nlogs* × *flen* weights that can be + used to apply the differential operator *dop* to multiple series using + :func:`_acb_ode_apply_diffop_basecase_precomp`. + +Other summation utilities +------------------------------------------------------------------------------- + +.. function:: void _acb_ode_poly_negdivrevhigh(acb_ptr res, acb_srcptr a, acb_srcptr cst, acb_srcptr b, slong len, slong prec) + + Solves the equation *a(1/X)·res(X) + cst·b(X) = (terms of negative degree)*, + where *a*, *res*, *b* are polynomials of length *len*. The *cst* pointer may + be null; this is understood as *cst = 1*. + +.. function:: void acb_ode_poly_taylor_shift_aps_trunc(acb_poly_t g, const acb_poly_t f, acb_srcptr a, slong n, slong len, slong prec); + + Shifts the variable by *a + n*, truncates to length *len*. + +Error bounds +------------------------------------------------------------------------------- + +The functions in this section operate simultaneously on one or more logarithmic +series solutions of the differential equation `L(y) = 0` belonging to the same +*group*, that is, of the form + +.. math:: + + y(x) = x^λ \sum_{k=0}^{K-1} y_k(x) \frac{\log(x)^k}{k!} + = x^λ \sum_{k=0}^{K-1} \sum_{j=0}^{∞} c_{k,j} x^j \frac{\log(x)^k}{k!}. + +with the same `λ`. +Given one such solution, the main goal is to bound the tails + +.. math:: + + y_{k,n:}(a) = \sum_{j=n}^{∞} c_{k,j} a^j, \qquad 0 ≤ k < K, + +of the power series `y_k` truncated to some order `n` and evaluated at `a ∈ ℂ`. + +The tail bound is computed from the *normalized residual* `q(x)` defined as +follows. Let `Q_0` be the monic indicial polynomial of `L` at the origin. +The equation + +.. math:: + + (Q_0(x·\mathrm d/\mathrm dx))(q) = L(y_{n:}) + +has a unique solution of the form + +.. math:: + + q(x) = x^{λ+n} \sum_{k} q_{k}(x) \frac{\log(x)^k}{k!} + = x^{λ+n} \sum_{k,j} q_{k,j} x^j \frac{\log(x)^k}{k!} + +where the sums are finite and `q_{k,j} = 0` whenever `λ + n + j` is a root +of multiplicity strictly more than `k` of `Q_0`. + +.. type:: acb_ode_bound_struct + acb_ode_bound_t + + Stores precomputed bounds depending neither on the particular solution in + the group nor on the truncation order. + +.. function:: void acb_ode_bound_init(acb_ode_ind_lbound_t ind_lbound) + void acb_ode_bound_clear(acb_ode_ind_lbound_t ind_lbound) + +.. function:: void acb_ode_bound_precompute(acb_ode_bound_t bound, const acb_poly_struct *dop, slong dop_len, acb_srcptr lcrt, slong pol_part_len, slong prec) + + Fills *bound* with data used when computing bounds on the tails of series + solutions of a given equation. + + The input consists of the operator *dop* of length *dop_len*, + a vector *lcrt* containing enclosures of the roots of the leading + coefficient (with multiple roots repeated), + a minimum number *pol_part_len* of initial terms of the local series + expansion of *dop* to be considered individually, + and a working precision *prec*. + + As a general rule, increasing *pol_part_len* leads to sharper bounds but + makes both the precomputation and subsequent evaluation of the bound slower. + Values of *prec* beyond 30–50 bits are rarely useful except for equations + with extremely large coefficients and/or clusters of singular points. + + When this function is called several times with the same destination object, + *dop* must correspond to the same operator (possibly represented at + different precision) for all calls. + In addition, calling this function again with a different input invalidates + any group-dependent data stored in *bound* (meaning that + :func:`acb_ode_bound_precompute_group` must be called again). + +.. function:: void acb_ode_bound_precompute_group(acb_ode_bound_t bound, const acb_poly_struct * dop, slong dop_len, const acb_ode_exponents_t expos, slong grp, slong prec) + + Takes a *bound* object partially filled with group-independent data and + updates it for computing bounds on solutions from a given group. + This may be done several times with different groups over the lifetime + of *bound*. + + The input consists of the operator *dop* of length *dop_len*, + its local exponent structure *expos*, + the index *grp* of the group in the exponent structure, + and a numeric working precision *prec*. + +.. function:: void acb_ode_tail_bound_jet(arb_poly_t res, const acb_ode_group_t group, const acb_ode_bound_t bound, slong n, slong nlogs, const arb_poly_t nres_maj, const arb_t rad, slong ord, slong prec) + void acb_ode_tail_bound_jet_precomp(arb_poly_t res, const acb_ode_bound_t bound, slong n, const arb_poly_t itg_pol, const arb_poly_t itg_num, const arb_poly_t nres_maj, const arb_t rad, slong ord, slong prec) + + Computes upper bounds on the first *ord* coefficients of the Taylor + expansion with respect to `z` of `y_{k,n:}(a + z)` where `|a|` is bounded + by *rad*. + The returned bounds are valid uniformly for all *k*, for all evaluation + points `a` in the disk of radius *rad* around zero, and for all + solutions `y` that belong to the given *group* and satisfy the following + conditions: + + * The polynomial *nres_maj* is a common majorant of the + coefficients `q_k(x)` of the normalized residual, i.e., for all `j` + and `k`, the absolute value `|q_{k,j}|` is bounded by the coefficient + of `x^j` in *nres_maj*. + + * The coefficients `c_{k,j}` with `j ≥ n` such that `λ + j` is a root of the + indicial polynomial `Q_0` of multiplicity *k* or more are zero. + + * The degree in `\log(x)` of the truncation at order `x^n` of `y(x)` is + strictly smaller than *nlogs*. In other words, `c_{k,j} = 0` when *j < n* + and *k ≥ nlogs*. + + The caller must supply a bound object *bound* initialized with precomputed + data associated with the differential equation and solution group of + interest. + + The *precomp* variant is suitable for repeated calls with the same *group* + and *nlogs* and values of *n* bounded from below. + It takes two polynomials *itg_pol* and *itg_num* obtained by calling + :func:`acb_ode_bound_precompute_integrand` with *n0 ≤ n*. + Typically, values of *n0* closer to *n* yield sharper bounds, but the effect + becomes less pronounced as *n* grows. + +.. function:: void acb_ode_bound_precompute_integrand(arb_poly_t itg_pol, arb_poly_t itg_num, const acb_ode_group_t group, const acb_ode_bound_t bound, slong n0, slong ord, slong prec) + + Computes polynomials *itg_pol* and *itg_num* that can be passed to + :func:`acb_ode_tail_bound_jet_precomp` to obtain bounds on the tails of + solutions from the group *group* at truncation orders ≥ *n0*. + + The remaining input arguments have the same meaning as the corresponding + arguments of :func:`acb_ode_tail_bound_jet`. + +Bounds on rational sequences +------------------------------------------------------------------------------- + +The functions in this section compute upper bounds of the expression + +.. math:: + + F_T(n) = n \cdot \sum_{t=0}^{T-1} \left| [X^t] \frac{p(n+X)}{X^{-m(n)} q(n+X)} \right| + +where `q` is a monic polynomial, `m(n)` denotes the multiplicity of `n` as a +root of `q`, and `p` is a polynomial of strictly smaller degree than `q`. +All functions work with a vector of numerators `p` and a single denominator `q`. + +The function and parameter names reflect that, in the intended application, +`q` is the monic indicial polynomial of the differential equation being solved. + +Some functions accept a working precision *prec* for intermediate computations +even though their outputs are represented in fixed precision. + +.. type:: acb_ode_ind_lbound_struct + acb_ode_ind_lbound_t + + A precomputed lower bound on `q(n)/n^r` where `q` is a monic polynomial + of degree `r`. + The data structure contains an integer *length* and a vector *r* of + length *length* of structures each containing + a root *root* of `q`, + its multiplicity *mult*, + an integer *n_min* such that *|1-root/n|* is nondecreasing for *n ≥ nmin*, + and a lower bound *global_lbound* of *|1-root/n|^mult* for *n ≥ nmin*. + +.. function:: void acb_ode_ind_lbound_init(acb_ode_ind_lbound_t ind_lbound) + void acb_ode_ind_lbound_clear(acb_ode_ind_lbound_t ind_lbound) + +.. function:: void acb_ode_ind_lbound_precompute(acb_ode_ind_lbound_t ind_lbound, const acb_ode_exponents_t expos, slong grp, slong prec) + + Fills *ind_lbound* with data corresponding to `q(X) = q_0(λ + X)` where + `q_0` is the monic polynomial whose roots are given by *expos* and `λ` is + the leader of the group of index *grp*. + +In the following descriptions `\tau(n) = \sum_{k=0}^{n} m(k)`. + +.. type:: acb_ode_stairs_struct + acb_ode_stairs_t + + A vector *h* of length *length* of *mag_struct*. + +.. function:: void acb_ode_stairs_init(acb_ode_stairs_t stairs) + void acb_ode_stairs_clear(acb_ode_stairs_t stairs) + +.. function:: void acb_ode_stairs_precompute(acb_ode_stairs_t stairs, const acb_poly_struct * num, slong len, const acb_poly_t ind, const acb_ode_group_t group, const acb_ode_ind_lbound_t ind_lbound, slong prec) + + For each `p` in the vector *num* of length *len* and for each shift `s` in + *group*, computes an upper bound of `F_{\tau(n)}(n)` valid for *n ≥ s*. + The parameter *ind* corresponds to `q`. + The parameter *group* must describe the integer roots of `q` and their + multiplicities (its *leader* field is ignored). + Additionally, the caller must supply an *ind_lbound* structure precomputed + as described above. + +.. function:: void acb_ode_bound_rat_vec(mag_ptr res, const acb_poly_struct * num, slong len, const acb_poly_t ind, const acb_ode_group_t group, const acb_ode_ind_lbound_t ind_lbound, const acb_ode_stairs_t stairs, slong n0, slong ord, slong prec) + + For each `p` in the vector *num* of length *len*, computes a bound on + `F_{\rho(n)}(n)` valid for all `n \geq n_0`, where `\rho(n) = \tau(n)` when + there is no integer root of `q` between `n_0` (exclusive) and `n` + (inclusive), and `\rho(n)` is equal to *ord* otherwise. + The parameter *ind* specifies the common denominator `q`. + + The caller must supply the integer roots of `q` in the parameter *group* + (the *leader* field is ignored), as well as two data structures *ind_lbound* + and *stairs* precomputed as described above. + +.. function:: void acb_ode_bound_rat_ordinary_vec(mag_ptr res, const acb_poly_struct * num, slong len, const acb_poly_t ind, const acb_ode_ind_lbound_t ind_lbound, slong n0, slong ord, slong prec) + + Similar but assumes `m(n_0) = 0` and computes a bound on `F_{T}(n)` valid + for all `n \geq n_0` with `m(n) = 0`. + The value of `T` is given by the parameter *ord*. + +.. function:: void acb_ode_bound_rat_ref_vec(mag_ptr res, const acb_poly_struct * num, slong len, const acb_poly_t ind, slong n, slong mult, slong ord, slong prec) + + Similar but assumes that `m(n)` is equal to *mult* and computes a bound on + `F_{T}(n)` at the given *n*. + +Miscellaneous bounds and estimates +------------------------------------------------------------------------------- + +.. function:: void _acb_ode_solution_growth(mag_t order, mag_t base, const acb_poly_struct * dop, slong dop_len) + + Estimates parameters *order* and *base* such that the coefficient sequences + of series solutions of *dop* are + `\mathrm O(\mathit{base}^n/n!^(1/\mathit{order}))`. + Assumes that the leading coefficient of *dop* is a constant. + The return values are currently not rigorous bounds. + +.. function:: slong acb_ode_choose_prec(slong * rec_prec, const acb_poly_struct * dop, slong dop_len, mag_srcptr rad, mag_srcptr cvrad, slong tgt_prec) + +.. type:: acb_ode_cvest_struct + acb_ode_cvest_t + +.. function:: void acb_ode_cvest_init(acb_ode_cvest_t cvest) + void acb_ode_cvest_clear(acb_ode_cvest_t cvest) + +.. function:: void acb_ode_cvest_update(acb_ode_cvest_t cvest, const acb_ode_cvest_t old, mag_t est, slong acc, slong stride, slong prec, slong work_prec) + +Examples +------------------------------------------------------------------------------- + +XXX simplify examples + +.. _ode-example-basic: + +Basic usage for computing a fundamental matrix: + +.. literalinclude:: ../../examples/ode_basic.c + :language: c + +.. _ode-example-exponents: + +In the following example, the coefficients of the operator are approximate and +we need to specify the exponents: + +.. literalinclude:: ../../examples/ode_fundamental_matrix_manual_exponents.c + :language: c diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst index 17d0e6c852..40d432ccb2 100644 --- a/doc/source/acb_poly.rst +++ b/doc/source/acb_poly.rst @@ -1205,3 +1205,25 @@ Root-finding If this function returns zero, then the signs of the imaginary parts are not known for certain, based on the accuracy of the inputs and the working precision *prec*. + + +Vector functions +------------------------------------------------------------------------------- + +.. function:: acb_poly_struct * _acb_poly_vec_init(slong n) + +.. function:: void _acb_poly_vec_clear(acb_poly_struct * vec, slong n) + +.. function:: void _acb_poly_vec_set(acb_poly_struct * dest, const acb_poly_struct * src, slong n) + +.. function:: void _acb_poly_vec_set_block(acb_poly_struct * dest, const acb_poly_struct * src, slong n, slong base, slong len) + +.. function:: slong _acb_poly_vec_length(const acb_poly_struct * vec, slong n) + +.. function:: void _acb_poly_vec_fit_length(acb_poly_struct * vec, slong n, slong len) + +.. function:: void _acb_poly_vec_set_length(acb_poly_struct * vec, slong n, slong len) + +.. function:: void _acb_poly_vec_normalise(acb_poly_struct * vec, slong n) + +.. function:: int _acb_poly_vec_overlaps(acb_poly_struct * vec1, acb_poly_struct * vec2, slong n) diff --git a/doc/source/index.rst b/doc/source/index.rst index e433f11521..e82a4c6873 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -209,6 +209,7 @@ Real and complex numbers acb_modular.rst acb_theta.rst acb_dirichlet.rst + acb_ode.rst bernoulli.rst hypgeom.rst partitions.rst diff --git a/examples/ode_basic.c b/examples/ode_basic.c new file mode 100644 index 0000000000..0a76c8035e --- /dev/null +++ b/examples/ode_basic.c @@ -0,0 +1,69 @@ +#include + +#include "acb_types.h" +#include "acb.h" +#include "acb_mat.h" +#include "acb_ode.h" +#include "gr.h" +#include "gr_ore_poly.h" + +int main(void) +{ + gr_ctx_t ZZ, Pol, Dop; + gr_ptr dop; + + int status = GR_SUCCESS; + + slong prec = 50; + + // Differential operator + + // XXX fmpz_poly? fmpq_poly?? + // XXX std derivative + gr_ctx_init_fmpz(ZZ); + gr_ctx_init_gr_poly(Pol, ZZ); + gr_ctx_init_gr_ore_poly(Dop, Pol, 0, ORE_ALGEBRA_EULER_DERIVATIVE); + status |= gr_ctx_set_gen_name(Pol, "z"); + status |= gr_ctx_set_gen_name(Dop, "Tz"); + GR_TMP_INIT(dop, Dop); + status |= gr_ore_poly_set_str(dop, "(z^2 + 1)*Tz^2 + (z^2 - 1)*Tz", Dop); + // status |= gr_ore_poly_set_str(dop, "4*Tz^2 - 4*Tz - z^2 + 8*z - 11", Dop); + // status |= gr_ore_poly_set_str(dop, "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2", Dop); + // status |= gr_ore_poly_set_str(dop, "4*Tz^2 - 4*Tz - z^2 + 8*z - 11", Dop); + GR_MUST_SUCCEED(status); + + // flint_printf("dop = %{gr}\n", dop, Dop); + + // Evaluation point + + acb_t pt; + acb_init(pt); + // acb_set_d(pt, .5); + acb_const_pi(pt, prec + 4); + arb_one(acb_imagref(pt)); + acb_inv(pt, pt, prec + 4); + + // Output matrix + + acb_mat_t mat; + slong dop_order = gr_ore_poly_length(dop, Dop) - 1; + acb_mat_init(mat, dop_order, dop_order); + + // Numerical solution + + status |= acb_ode_fundamental_matrix(mat, dop, Dop, NULL, NULL, pt, 0, prec); + GR_MUST_SUCCEED(status); + + flint_printf("\n------\n"); + acb_mat_printd(mat, prec+4); + flint_printf("\n%f\n", prec/log2(10.)); + + acb_clear(pt); + acb_mat_clear(mat); + GR_TMP_CLEAR(dop, Dop); + + gr_ctx_clear(Dop); + gr_ctx_clear(Pol); + + flint_cleanup_master(); +} diff --git a/examples/ode_fundamental_matrix_manual_exponents.c b/examples/ode_fundamental_matrix_manual_exponents.c new file mode 100644 index 0000000000..dbc38da4a3 --- /dev/null +++ b/examples/ode_fundamental_matrix_manual_exponents.c @@ -0,0 +1,144 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_mat.h" +#include "acb_ode.h" +#include "gr.h" +#include "gr_ore_poly.h" + +// XXX also manual sing? + +void +fundamental_matrix(const char * dop_str, + const acb_ode_exponents_struct * expos, + double pt_d) +{ + gr_ctx_t CC, Pol, Dop; + gr_ptr dop; + + slong prec = 30; + + int status = GR_SUCCESS; + + gr_ctx_init_complex_acb(CC, prec); + gr_ctx_init_gr_poly(Pol, CC); + gr_ctx_init_gr_ore_poly(Dop, Pol, 0, ORE_ALGEBRA_EULER_DERIVATIVE); + + GR_TMP_INIT(dop, Dop); + + status |= gr_ctx_set_gen_name(Pol, "z"); + status |= gr_ctx_set_gen_name(Dop, "Tz"); + status |= gr_ore_poly_set_str(dop, dop_str, Dop); + + status |= gr_println(dop, Dop); + flint_printf("\n"); + + GR_MUST_SUCCEED(status); + + slong dop_order = gr_ore_poly_length(dop, Dop) - 1; + + acb_mat_t mat; + acb_mat_init(mat, dop_order, dop_order); + + acb_t pt; + acb_init(pt); + acb_set_d(pt, pt_d); + + GR_MUST_SUCCEED(acb_ode_fundamental_matrix(mat, dop, Dop, expos, NULL, pt, 0, prec)); + + flint_printf("%{acb_mat}\n\n", mat); + flint_printf("--------\n\n"); + + acb_mat_clear(mat); + GR_TMP_CLEAR(dop, Dop); + gr_ctx_clear(Dop); + gr_ctx_clear(Pol); + gr_ctx_clear(CC); + acb_clear(pt); +} + + +void +apery(void) +{ + acb_ode_shift_struct shift[1] = {{ .n = 0, .mult = 3 }}; + acb_ode_group_struct grp[1] = {{ .nshifts = 1, .shifts = shift }}; + acb_init(grp->leader); + acb_zero(grp->leader); + acb_ode_exponents_struct expos[1] = {{ .ngroups = 1, .groups = grp }}; + + fundamental_matrix( + "(z^2 - 34*z + 1)*Tz^3 + (3*z^2 - 51*z)*Tz^2 + (3*z^2 - 27*z)*Tz + z^2 - 5*z", + expos, + 0.015625); + + acb_clear(grp->leader); +} + + +void +multiple_shifts(void) +{ + const char * dop = "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2"; + + acb_ode_shift_struct shift[3] = { + { .n = 0, .mult = 2 }, + { .n = 1, .mult = 3 }, + { .n = 3, .mult = 1 }, + }; + acb_ode_group_struct grp[1] = {{ .nshifts = 3, .shifts = shift }}; + acb_init(grp->leader); + acb_zero(grp->leader); + acb_ode_exponents_struct expos[1] = {{ .ngroups = 1, .groups = grp }}; + + fundamental_matrix(dop, expos, 2.); + + acb_clear(grp->leader); +} + + +void +whittaker(void) +{ + slong prec = 30; + + acb_t kappa, mu, half; + acb_init(kappa); + acb_init(mu); + acb_init(half); + + const char * dop = "4*Tz^2 - 4*Tz - z^2 + 8*z - 11"; + acb_set_si(kappa, 2); + acb_set_si(mu, 3); + acb_sqrt(mu, mu, prec); + acb_set_d(half, .5); + + acb_ode_shift_struct shift[1] = { { .n = 0, .mult = 1 } }; + acb_ode_group_struct grp[2] = { + { .nshifts = 1, .shifts = shift }, + { .nshifts = 1, .shifts = shift }, + }; + acb_init(grp[0].leader); + acb_sub(grp[0].leader, half, mu, prec); + acb_init(grp[1].leader); + acb_add(grp[1].leader, half, mu, prec); + acb_ode_exponents_struct expos[1] = {{ .ngroups = 2, .groups = grp }}; + + fundamental_matrix(dop, expos, 2.); + + acb_clear(grp[0].leader); + acb_clear(grp[1].leader); + acb_clear(kappa); + acb_clear(mu); + acb_clear(half); +} + + +int +main(void) +{ + apery(); + multiple_shifts(); + whittaker(); + + flint_cleanup_master(); +} diff --git a/examples/ode_low_level.c b/examples/ode_low_level.c new file mode 100644 index 0000000000..f476958fed --- /dev/null +++ b/examples/ode_low_level.c @@ -0,0 +1,197 @@ +#include "acb_types.h" +#include "acb_poly.h" +#include "acb_mat.h" +#include "acb_ode.h" + + +// XXX get rid of this file? + +void +_acb_ode_sum_swap_mat_ordinary( + acb_mat_t mat, + const acb_ode_sum_struct * sum, slong s) +{ + for (slong j = 0; j < sum->nsols; j++) + { + FLINT_ASSERT(sum->sol[j].nlogs <= 1); + if (sum->sol[j].nlogs < 1) + continue; + for (slong i = 0; i < sum->nder; i++) + acb_swap(acb_mat_entry(mat, i, j), + acb_ode_sol_sum_ptr(sum->sol + j, s, 0, i)); + } +} + + +void +ordinary(void) +{ + acb_ode_sum_t sum; + + /* (x^2 + 1)*Dx^3 + 7*x + * = (x^2 + 1)*Tx^3 + (-3*x^2 - 3)*Tx^2 + (2*x^2 + 2)*Tx + 7*x^4 */ + acb_ode_sum_init(sum, 4, 2, 3, 3); + acb_poly_set_coeff_si(sum->dop + 3, 2, 1); + acb_poly_set_coeff_si(sum->dop + 3, 0, 1); + acb_poly_set_coeff_si(sum->dop + 2, 2, -3); + acb_poly_set_coeff_si(sum->dop + 2, 0, -3); + acb_poly_set_coeff_si(sum->dop + 1, 2, 2); + acb_poly_set_coeff_si(sum->dop + 1, 0, 2); + acb_poly_set_coeff_si(sum->dop + 0, 4, 7); + + /* sum->flags |= ACB_ODE_WANT_SERIES; */ + sum->flags |= ACB_ODE_APPROX; + + acb_ode_sum_set_ordinary(sum); + acb_ode_sum_set_ini_echelon(sum); + + acb_set_d_d(sum->pts, 0.25, 0.25); + acb_set_si(sum->pts + 1, 0); + mag_set_d(&((sum->pts + 1)->real.rad), 1e-8); + + acb_ode_sum_divconquer(sum, 20, 64); + + acb_mat_t mat; + acb_mat_init(mat, sum->nder, sum->nsols); + + for (slong i = 0; i < sum->npts; i++) + { + _acb_ode_sum_swap_mat_ordinary(mat, sum, i); + acb_mat_printd(mat, 8); + /* flint_printf("%{acb_mat}\n\n", mat); */ + } + + acb_ode_sum_clear(sum); + acb_mat_clear(mat); +} + + +void +series(void) +{ + acb_ode_sum_t sum; + + acb_ode_sum_init(sum, 2, 0, 1, 1); + acb_poly_set_coeff_si(sum->dop + 1, 0, 1); + acb_poly_set_coeff_si(sum->dop + 0, 1, -1); + + acb_ode_sum_set_ordinary(sum); + acb_ode_sum_set_ini_echelon(sum); + sum->flags |= ACB_ODE_WANT_SERIES; + + slong len = 5; + acb_ode_sum_divconquer(sum, len, 64); + /* Clear the high part reserved for the residual. (In this special case, the + * high part is zero because block_length = 1.) */ + acb_poly_truncate(sum->sol[0].series, len); + + flint_printf("%{acb_poly}\n", sum->sol[0].series); + + acb_ode_sum_clear(sum); +} + + +void +bessel_j0(void) +{ + acb_ode_sum_t sum; + + slong dop_order = 2; + acb_ode_sum_init(sum, dop_order + 1, 1, dop_order, dop_order); + acb_poly_set_coeff_si(sum->dop + 2, 0, 1); + acb_poly_set_coeff_si(sum->dop + 0, 2, 1); + + sum->group->shifts[0].n = 0; + sum->group->shifts[0].mult = 2; + + acb_ode_sum_set_ini_echelon(sum); + + acb_set_d(sum->pts, 0.25); + + acb_ode_sum_divconquer(sum, 10, 64); + + for (slong m = 0; m < sum->nsols; m++) + { + acb_ode_sol_struct * sol = sum->sol + m; + /* series in x */ + flint_printf("f%wd =", m); + for (slong k = 0; k < sol->nlogs; k++) + { + flint_printf(" + (%{acb*})*log(%{acb} + x)^%wd/%wd!", + acb_ode_sol_sum_ptr(sol, 0, k, 0), sum->nder, + sum->pts, k, k); + } + flint_printf("\n"); + } + + acb_ode_sum_clear(sum); +} + + +void +whittaker_m(void) /* non-integer exponent, no logs */ +{ + acb_t kappa, mu2, half; + acb_poly_t val; + + acb_init(kappa); + acb_init(mu2); + acb_init(half); + acb_poly_init(val); + + acb_ode_sum_t sum; + + acb_set_si(kappa, 2); + acb_set_si(mu2, 3); + acb_set_d(half, .5); + + slong prec = 64; + slong dop_order = 2; + slong len = dop_order; + + acb_ode_sum_init(sum, dop_order + 1, 1, 1, len); + + acb_poly_set_coeff_si(sum->dop + 2, 0, 4); + acb_poly_set_coeff_si(sum->dop + 1, 0, -4); + acb_poly_set_coeff_si(sum->dop + 0, 2, -1); + acb_mul_si((sum->dop + 0)->coeffs + 1, kappa, 4, prec); + acb_poly_set_coeff_si(sum->dop + 0, 0, 1); + acb_addmul_si((sum->dop + 0)->coeffs + 0, mu2, -4, prec); + + acb_sqrt(sum->group->leader, mu2, prec); + acb_add(sum->group->leader, sum->group->leader, half, prec); /* other leader = 1/2 - mu */ + + sum->group->shifts[0].n = 0; + sum->group->shifts[0].mult = 1; + + acb_ode_sum_set_ini_echelon(sum); + + acb_set_d(sum->pts, 1.4242); + + acb_ode_sum_divconquer(sum, -1, prec); + + flint_printf("(%{acb} + x)^(%{acb}) * (%{acb*})\n", sum->pts, sum->group->leader, + sum->sol->sums, sum->nder); + + acb_ode_sol_jet(val, sum->group->leader, sum->sol, 0, sum->pts, sum->nder, 1, prec); + + flint_printf("M(%{acb} + x) = %{acb_poly} + O(x^%wd)\n", sum->pts, val, len); + + acb_ode_sum_clear(sum); + acb_clear(kappa); + acb_clear(mu2); + acb_clear(half); + acb_poly_clear(val); +} + + +int +main(void) +{ + ordinary(); + series(); + bessel_j0(); + whittaker_m(); + + flint_cleanup_master(); +} diff --git a/src/acb.h b/src/acb.h index 2b6069a6a8..55502d8837 100644 --- a/src/acb.h +++ b/src/acb.h @@ -903,6 +903,8 @@ int _acb_vec_is_finite(acb_srcptr vec, slong len); int _acb_vec_equal(acb_srcptr vec1, acb_srcptr vec2, slong len); int _acb_vec_overlaps(acb_srcptr vec1, acb_srcptr vec2, slong len); int _acb_vec_contains(acb_srcptr vec1, acb_srcptr vec2, slong len); +void _acb_vec_get_mag(mag_t bound, acb_srcptr vec, slong len); +void _acb_vec_get_mag_lower(mag_t bound, acb_srcptr vec, slong len); void _acb_vec_get_real(arb_ptr re, acb_srcptr vec, slong len); void _acb_vec_get_imag(arb_ptr im, acb_srcptr vec, slong len); void _acb_vec_set_real_imag(acb_ptr vec, arb_srcptr re, arb_srcptr im, slong len); diff --git a/src/acb/vec_get_mag.c b/src/acb/vec_get_mag.c new file mode 100644 index 0000000000..4e6f8e1510 --- /dev/null +++ b/src/acb/vec_get_mag.c @@ -0,0 +1,34 @@ +/* + Copyright (C) 2014 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "acb.h" + +void +_acb_vec_get_mag(mag_t bound, acb_srcptr vec, slong len) +{ + if (len < 1) + { + mag_zero(bound); + } + else + { + mag_t t; + slong i; + acb_get_mag(bound, vec); + mag_init(t); + for (i = 1; i < len; i++) + { + acb_get_mag(t, vec + i); + mag_max(bound, bound, t); + } + mag_clear(t); + } +} diff --git a/src/acb/vec_get_mag_lower.c b/src/acb/vec_get_mag_lower.c new file mode 100644 index 0000000000..7ccc9468d8 --- /dev/null +++ b/src/acb/vec_get_mag_lower.c @@ -0,0 +1,34 @@ +/* + Copyright (C) 2014 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "acb.h" + +void +_acb_vec_get_mag_lower(mag_t bound, acb_srcptr vec, slong len) +{ + if (len < 1) + { + mag_inf(bound); + } + else + { + mag_t t; + slong i; + acb_get_mag(bound, vec); + mag_init(t); + for (i = 1; i < len; i++) + { + acb_get_mag_lower(t, vec + i); + mag_min(bound, bound, t); + } + mag_clear(t); + } +} diff --git a/src/acb_ode.h b/src/acb_ode.h new file mode 100644 index 0000000000..b8b9513923 --- /dev/null +++ b/src/acb_ode.h @@ -0,0 +1,489 @@ +#ifndef ACB_ODE_H +#define ACB_ODE_H + +#ifdef ACB_ODE_INLINES_C +#define ACB_ODE_INLINE +#else +#define ACB_ODE_INLINE static inline +#endif + +// /* temporary (for sage) */ +// #include "gmp.h" +// typedef mp_limb_signed_t slong; +// /* */ + +#include "acb_types.h" +#include "gr_types.h" +#include "gr_ore_poly.h" + +// #ifdef __cplusplus +// extern "C" { +// #endif + +/****************************************************************************/ + +// probably want: +// - a data structure for storing “local” information (shifted singularities, +// local exponent structure...), +// - a function (acb_ode_focus?) for “moving it” to another point +typedef struct +{ + // XXX sing_gr? + acb_ptr sing_acb; + slong sing_len; + + // XXX maybe something to track the precision at which transition matrices + // are known?? +} +acb_ode_context_struct; + +typedef acb_ode_context_struct acb_ode_context_t[1]; + +/****************************************************************************/ + +typedef struct +{ + slong n; + slong mult; +} +acb_ode_shift_struct; + +typedef acb_ode_shift_struct acb_ode_shift_t[1]; + + +typedef struct +{ + acb_t leader; + acb_ode_shift_struct * shifts; + slong nshifts; +} +acb_ode_group_struct; + +typedef acb_ode_group_struct acb_ode_group_t[1]; + + +typedef struct +{ + /* Memory managed by the exponent object, partial sharing across groups */ + acb_ode_group_struct * groups; + slong ngroups; +} +acb_ode_exponents_struct; + +typedef acb_ode_exponents_struct acb_ode_exponents_t[1]; + + +typedef enum +{ + ACB_ODE_BASIS_ECHELON = 0, + ACB_ODE_BASIS_CASCADE, /* XXX: same as Frobenius? */ + + ACB_ODE_NUM_BASES +} +acb_ode_basis_t; + +void acb_ode_group_init(acb_ode_group_t group, slong len); +void acb_ode_group_clear(acb_ode_group_t group); + +void acb_ode_group_set(acb_ode_group_t dest, const acb_ode_group_t src); +slong acb_ode_group_length(const acb_ode_group_t grp); +slong acb_ode_group_multiplicity (const acb_ode_group_t group, slong n); +slong acb_ode_group_nlogs (const acb_ode_group_t group, slong n); + +void acb_ode_exponents_init(acb_ode_exponents_t expos); +void acb_ode_exponents_clear(acb_ode_exponents_t expos); + +void acb_ode_exponents_ordinary(acb_ode_exponents_t expos, slong dop_order); +int acb_ode_exponents(acb_ode_exponents_t expos, const gr_ore_poly_t dop, gr_ctx_t Dop, gr_ctx_t CC); + +/**************************************************************************** + * Error bounds + ****************************************************************************/ + +typedef struct +{ + slong length; + struct + { + acb_t root; + mag_t global_lbound; + slong mult; + slong n_min; + } + * r; +} +acb_ode_ind_lbound_struct; + +typedef acb_ode_ind_lbound_struct acb_ode_ind_lbound_t[1]; + +typedef struct +{ + slong length; + mag_ptr h; +} +acb_ode_stairs_struct; + +typedef acb_ode_stairs_struct acb_ode_stairs_t[1]; + +typedef struct +{ + slong prec; + + arb_t cst; /* 1/abs(lc_x(lc_θ(dop)) */ + arb_poly_t den_lbound; + arb_ptr den_rt; /* roots(lc_θ(dop)), with multiple roots repeated */ + slong den_rt_len; + + acb_poly_struct * all_nums; + slong all_nums_len; + slong pol_part_len; + + acb_poly_t ind; /* MONIC indicial polynomial shifted so that 0 corresponds + to the current group leader */ + acb_ode_ind_lbound_t ind_lbound; + acb_ode_stairs_t stairs; +} +acb_ode_bound_struct; + +typedef acb_ode_bound_struct acb_ode_bound_t[1]; + +void acb_ode_bound_init(acb_ode_bound_t bound); +void acb_ode_bound_clear(acb_ode_bound_t bound); +void acb_ode_bound_precompute(acb_ode_bound_t bound, + const acb_poly_struct *dop, slong dop_len, acb_srcptr lcrt, + slong pol_part_len, slong prec); +void acb_ode_bound_precompute_group(acb_ode_bound_t bound, + const acb_poly_struct * dop, slong dop_len, + const acb_ode_exponents_t expos, slong grp, slong prec); + +void acb_ode_bound_precompute_integrand(arb_poly_t itg_pol, arb_poly_t itg_num, + const acb_ode_group_t group, const acb_ode_bound_t bound, + slong n0, slong ord, slong prec); +void acb_ode_tail_bound_jet_precomp(arb_poly_t res, + const acb_ode_bound_t bound, slong n, + const arb_poly_t itg_pol, const arb_poly_t itg_num, + const arb_poly_t nres_maj, + const arb_t rad, slong ord, slong prec); +void acb_ode_tail_bound_jet(arb_poly_t res, + const acb_ode_group_t group, const acb_ode_bound_t bound, + slong n, slong nlogs, const arb_poly_t nres_maj, + const arb_t rad, slong ord, slong prec); + +void acb_ode_ind_lbound_init(acb_ode_ind_lbound_t ind_lbound); +void acb_ode_ind_lbound_clear(acb_ode_ind_lbound_t ind_lbound); +void acb_ode_ind_lbound_precompute(acb_ode_ind_lbound_t ind_lbound, + const acb_ode_exponents_t expos, slong grp, slong prec); + +void acb_ode_stairs_init(acb_ode_stairs_t stairs); +void acb_ode_stairs_clear(acb_ode_stairs_t stairs); +void acb_ode_stairs_precompute(acb_ode_stairs_t stairs, + const acb_poly_struct * num, slong len, + const acb_poly_t ind, const acb_ode_group_t group, + const acb_ode_ind_lbound_t ind_lbound, slong prec); + +void acb_ode_bound_rat_vec(mag_ptr res, + const acb_poly_struct * num, slong len, + const acb_poly_t ind, const acb_ode_group_t group, + const acb_ode_ind_lbound_t ind_lbound, const acb_ode_stairs_t stairs, + slong n0, slong ord, slong prec); +void acb_ode_bound_rat_ordinary_vec(mag_ptr res, + const acb_poly_struct * num, slong len, + const acb_poly_t ind, const acb_ode_ind_lbound_t ind_lbound, + slong n0, slong ord, slong prec); +void acb_ode_bound_rat_ref_vec(mag_ptr res, + const acb_poly_struct *num, slong len, + const acb_poly_t ind, slong n, slong mult, slong ord, slong prec); + +/****************************************************************************/ + +typedef struct +{ + double neglogterm; + double cvg_rate; + slong accuracy; + double loss_rate; + double terms_wanted; + double terms_full_prec; + slong prec_wanted; +} +acb_ode_cvest_struct; + +typedef acb_ode_cvest_struct acb_ode_cvest_t[1]; + +void acb_ode_cvest_init(acb_ode_cvest_t cvest); +void acb_ode_cvest_clear(acb_ode_cvest_t cvest); +void acb_ode_cvest_update(acb_ode_cvest_t cvest, const acb_ode_cvest_t old, mag_t est, slong accuracy, slong stride, slong prec, slong work_prec); + +/****************************************************************************/ + +/* todo: better vectors of polynomials / polynomial with vector coefficients + (single buffer, row stride, essentially self-resizing matrices) */ + +typedef struct +{ + /* vector of polynomials in x, entries = coefficients of log(x)^k/k! of + * a slice of solution and/or rhs */ + acb_poly_struct * series; + + /* coefficients of Taylor expansions (jets) of the partial sums + * dim: sums[npts][alloc_logs][nder] */ + acb_ptr sums; + acb_mat_t extini; + + /* bounds on the tails of successive derivatives, valid uniformly for all + * coefficients of log(x)^k/k! and all evaluation points */ + // XXX move into cvest? + mag_ptr tb; + + /* convergence estimates */ + acb_ode_cvest_t cvest; + + slong alloc_logs; + slong nlogs; + slong npts; + slong nder; + + int done; +} +acb_ode_sol_struct; + +typedef acb_ode_sol_struct acb_ode_sol_t[1]; + +ACB_ODE_INLINE acb_ptr +acb_ode_sol_sum_ptr(const acb_ode_sol_t sol, + slong j, slong k, slong i) +{ + return sol->sums + (j * sol->alloc_logs + k) * sol->nder + i; +} + +void acb_ode_sol_init(acb_ode_sol_t sol, slong nshifts, + slong nlogs, slong npts, slong nder); + +void acb_ode_sol_clear(acb_ode_sol_t sol); + +void acb_ode_sol_zero(acb_ode_sol_t sol); + +void acb_ode_sol_fit_length(acb_ode_sol_t sol, slong len); + + +void +acb_ode_sol_set_ini(acb_ode_sol_t sol, acb_srcptr ini, + const acb_ode_shift_t shifts); // unused +void acb_ode_sol_unit_ini(acb_ode_sol_t sol, slong i0, + const acb_ode_shift_t shifts); + +void +_acb_ode_sol_add_term( + acb_ode_sol_t sol, slong base, slong n, + acb_poly_struct * rhs, slong rhs_nlogs, + /* XXX define a separate struct with the next three? */ + slong mult, slong ini_i, const acb_poly_t ind_n, + acb_srcptr pows, const char * shifted, slong npts, + const fmpz * binom_n, + slong nder, slong prec, slong sums_prec); + +void _acb_ode_sol_jet(acb_poly_struct * val, acb_srcptr expo, + acb_srcptr sums, mag_srcptr tb, acb_srcptr pt, + slong nlogs, slong nder, slong nfrobshifts, + slong prec); +void acb_ode_sol_jet(acb_poly_struct * val, acb_srcptr expo, const acb_ode_sol_t + sol, slong j, acb_srcptr pt, slong nder, slong nfrobshifts, + slong prec); + +slong +acb_ode_sol_estimate_terms(mag_t est, + const acb_ode_sol_t sol, slong off, slong len, + const mag_t radpow); +slong acb_ode_sol_estimate_sums(mag_t est, mag_t rad, const acb_ode_sol_t sol); + + +/****************************************************************************/ + +#define ACB_ODE_WANT_SERIES 1 +#define ACB_ODE_APPROX 2 + +/* to be adapted for compatibility with other summation algorithms */ + +typedef struct +{ + /* operator */ + acb_poly_struct * dop; /* XXX FJ: move outside the struct? */ + slong dop_len; + slong dop_clen; + mag_t cvrad; + + /* solution group */ // XXX full exponent structure??? + acb_ode_group_t group; + acb_poly_t ind; + + /* solutions */ + /* XXX possible alternative design: move out of sum struct, use sum as a + context when manipulating sols, replace some sum_* functions by sol_*_vec + functions, some of which would also update the context object */ + acb_ode_sol_struct * sol; + slong nsols; + slong n0; /* index of the first stored coefficient */ + slong n; /* current truncation order */ + + /* evaluation points */ + acb_ptr pts; /* evaluation points x[j] for 0 < j < npts */ + acb_ptr pows; /* x[j]^{n-i} for i < nder (i-major, cyclic on i); + only the block i=0 is really used for nonzero points */ + char * shifted_sums; /* nonzero if the i-th derivatives of the sums at x[j] + stored in sol have extra x[j]^i factors */ + slong npts; + slong nder; + + /* working precisions */ + slong wp; + slong sum_wp; + + /* error bounds + XXX at least the pointer probably should be a separate argument */ + acb_ode_bound_struct * bound; + mag_t mag; + mag_t magpow; + arb_poly_t itg_pol, itg_num; // itg_n ? + + /* options + XXX at least some of these would make better sense outside the sum object + since one may want to pass them to higher-level functions than sum_* */ + ulong flags; + // slong bound_pol_part_len? + void * data; /* user data for sum workers */ + + /* miscellaneous state */ + fmpz * binom_n; /* binom(n, j) for j < nder */ + int have_precomputed; +} +acb_ode_sum_struct; + +typedef acb_ode_sum_struct acb_ode_sum_t[1]; + +void acb_ode_sum_init(acb_ode_sum_t sum, slong dop_len, slong npts, slong nsols, slong nder); +void acb_ode_sum_clear(acb_ode_sum_t sum); + +void acb_ode_sum_set_diffop(acb_ode_sum_t sum, const acb_poly_t dop, slong dop_len, const mag_t cvrad); + +void acb_ode_sum_set_group(acb_ode_sum_t sum, const acb_ode_group_t grp); +void acb_ode_sum_set_ordinary(acb_ode_sum_t sum); + +void acb_ode_sum_set_ini_echelon(acb_ode_sum_t sum); +void acb_ode_sum_set_ini_highest(acb_ode_sum_t sum); + +void acb_ode_sum_set_points(acb_ode_sum_t sum, acb_srcptr pts, slong npts); + +slong acb_ode_sum_max_nlogs(acb_ode_sum_struct * sum); +slong acb_ode_sum_max_nlogs_xn(acb_ode_sum_struct * sum, slong off); + +void acb_ode_sum_precompute(acb_ode_sum_t sum); + +void _acb_ode_sum_forward_1(acb_ode_sum_struct * sum); +int _acb_ode_sum_forward_divconquer(acb_ode_sum_struct * sum, slong high, slong block_len, slong stride, slong prec); +int acb_ode_sum_done(acb_ode_sum_struct * sum, slong stride, slong prec); +void _acb_ode_sum_fix(acb_ode_sum_struct * sum); + +void acb_ode_sum_divconquer(acb_ode_sum_t sum, slong nterms, slong prec); + +void _acb_ode_poly_negdivrevhigh(acb_ptr res, acb_srcptr a, acb_srcptr cst, + acb_srcptr b, slong len, slong prec); +void acb_ode_poly_taylor_shift_aps_trunc(acb_poly_t g, const acb_poly_t f, acb_srcptr a, slong n, slong len, slong prec); + + +/****************************************************************************/ + +/* XXX struct with a data field? currently the data is part of sum, but sum may + not have the right lifespan */ +typedef void (* acb_ode_sum_worker_t)(acb_ode_sum_t sum, slong nterms, slong prec); + +void _acb_ode_fundamental_matrix_vec( + acb_mat_struct * mat, + const acb_poly_struct * dop, slong dop_len, + const acb_ode_exponents_t expos, + acb_srcptr lcrt, + acb_srcptr pts, slong npts, + acb_ode_basis_t basis, + acb_ode_sum_worker_t sum_worker, + slong prec); + +WARN_UNUSED_RESULT int acb_ode_fundamental_matrix_vec( + acb_mat_struct * mat, + const gr_ore_poly_t dop, gr_ore_poly_ctx_t dop_ctx, + const acb_ode_exponents_t expos, + acb_srcptr lcrt, + acb_srcptr pts, slong npts, + acb_ode_basis_t basis, + slong prec); + +WARN_UNUSED_RESULT int acb_ode_fundamental_matrix( + acb_mat_t mat, + const gr_ore_poly_t dop, gr_ore_poly_ctx_t dop_ctx, + const acb_ode_exponents_t expos, + acb_srcptr lcrt, + acb_srcptr pt, + acb_ode_basis_t basis, + slong prec); + +/****************************************************************************/ + +void _acb_ode_apply_diffop_basecase_weights( + acb_ptr weights, + const acb_poly_struct * dop, slong dop_len, + acb_srcptr expo, slong offset, + slong flen, slong nlogs, slong start, slong len, + slong prec); + +void _acb_ode_apply_diffop_basecase_precomp( + acb_poly_struct * g, slong goff, + acb_srcptr weights, slong weights_nlogs, + const acb_poly_struct * f, slong foff, slong flen, + slong nlogs, + slong start, slong len, + slong prec); + +void acb_ode_apply_diffop_basecase( + acb_poly_struct * g, + const acb_poly_struct * dop, slong dop_len, // gr_ore_poly? + acb_srcptr expo, slong offset, + const acb_poly_struct * f, + slong nlogs, + slong start, slong len, + slong prec); + +void _acb_ode_apply_diffop_polmul( + acb_poly_struct * g, slong goff, + const acb_poly_struct * dop, slong dop_len, + acb_srcptr expo, slong offset, + const acb_poly_struct * f, slong foff, slong flen, + slong nlogs, + slong start, slong len, + slong prec); + +void acb_ode_apply_diffop_polmul( + acb_poly_struct * g, + const acb_poly_struct * dop, slong dop_len, // gr_ore_poly? + acb_srcptr expo, slong offset, + const acb_poly_struct * f, + slong nlogs, + slong start, slong len, + slong prec); + +// todo: acb_ode_apply_diffop + +/****************************************************************************/ + +void _acb_ode_solution_growth(mag_t order, mag_t base, const acb_poly_struct * dop, slong dop_len); +slong acb_ode_choose_prec(slong * rec_prec, const acb_poly_struct * dop, slong dop_len, mag_srcptr rad, mag_srcptr cvrad, slong tgt_prec); + +/* utilities probably to be moved to gr_ore_poly */ + +int gr_ore_poly_ctx_over_gr_poly_base_ptrs(gr_ctx_struct ** Scalars, gr_ctx_struct ** Pol, const gr_ctx_t Ore); +int _gr_ore_poly_indicial_polynomial_euler_derivative(gr_ptr ind, gr_srcptr op, slong len, gr_ctx_t Ore); +int _gr_ore_poly_indicial_polynomial(gr_ptr ind, gr_srcptr op, slong len, gr_ctx_t Ore); +int gr_ore_poly_indicial_polynomial(gr_poly_t ind, const gr_ore_poly_t op, gr_ctx_t Ore); + +// #ifdef __cplusplus +// } +// #endif + +#endif diff --git a/src/acb_ode/apply_diffop_basecase.c b/src/acb_ode/apply_diffop_basecase.c new file mode 100644 index 0000000000..76ca0a2f22 --- /dev/null +++ b/src/acb_ode/apply_diffop_basecase.c @@ -0,0 +1,167 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_poly.h" +#include "acb_ode.h" +#include "fmpz_mat.h" +#include "gr.h" +#include "gr_mat.h" + + +/* Notation: + * + * We are considering an operator + * + * \sum_{i,j} a[i,j] x^j (x d/dx)^i + * + * applied to a series + * + * f = x^{expo + offset} \sum_{p,k} f[k,p] x^p log(x)^k/k!, + * + * resulting in + * + * g = x^{expo + offset} \sum{m1,k1} g[k1,m1] x^m1 log(x)^k1/k1!. + * + * We want to compute the terms of index m1 = start + p1 of g for 0 <= p1 < len. + */ + + +/* TODO specialized version for dop with fmpz coefficients (and expo == 0) */ + + +static void +acb_addmul_binom(acb_ptr c, acb_srcptr b, fmpz_mat_t binom, slong i, slong t, + slong prec) +{ + if (t == 0) + acb_add(c, c, b, prec); + else if (t == 1) + acb_addmul_si(c, b, i, prec); + else + acb_addmul_fmpz(c, b, fmpz_mat_entry(binom, i, t), prec); +} + + +void +_acb_ode_apply_diffop_basecase_weights( + acb_ptr weights, /* len * nlogs * flen */ + const acb_poly_struct * dop, slong dop_len, + acb_srcptr expo, slong offset, + slong flen, slong nlogs, slong start, slong len, + slong prec) +{ + + fmpz_mat_t binom; + if (nlogs >= 2) /* XXX: worth caching? */ + { + gr_ctx_t fmpz_ctx; + gr_ctx_init_fmpz(fmpz_ctx); + fmpz_mat_init(binom, dop_len, dop_len); + GR_MUST_SUCCEED(gr_mat_pascal((gr_mat_struct *) binom, -1, fmpz_ctx)); + gr_ctx_clear(fmpz_ctx); + } + + acb_t expo1; + acb_init(expo1); + + if (expo != NULL && acb_is_zero(expo)) + expo = NULL; + + for (slong p1 = 0; p1 < len; p1++) /* n1 = offset + start + p1 */ + { + for (slong t = 0; t < nlogs; t++) /* t = k - k1 */ + { + for (slong p = 0; p < flen; p++) + { + slong j = start + p1 - p; + slong n = offset + p; + + if (j < 0) + break; + + if (expo != NULL) + acb_add_si(expo1, expo, n, prec); + + acb_ptr c = weights + p1 * nlogs * flen + t * flen + p; + acb_zero(c); + for (slong i = dop_len - 1; i >= t; i--) /* Horner */ + { + if (expo == NULL) + acb_mul_si(c, c, n, prec); + else + acb_mul(c, c, expo1, prec); + if (j >= (dop + i)->length) + continue; + acb_ptr aij = (dop + i)->coeffs + j; + acb_addmul_binom(c, aij, binom, i, t, prec); + } + } + } + } + + if (nlogs >= 2) + fmpz_mat_clear(binom); + + acb_clear(expo1); +} + + +void +_acb_ode_apply_diffop_basecase_precomp( + acb_poly_struct * g, slong goff, + acb_srcptr weights, slong weights_nlogs, + const acb_poly_struct * f, slong foff, slong flen, + slong nlogs, + slong start, slong len, + slong prec) +{ + for (slong p1 = 0; p1 < len; p1++) /* n1 = offset + start + p1 */ + { + for (slong k1 = 0; k1 < nlogs; k1++) + { + acb_ptr dest = (g + k1)->coeffs + goff + p1; + + /* XXX integrate the loop on k in acb_dot? */ + for (slong k = k1; k < nlogs; k++) + { + acb_srcptr src = (f + k)->coeffs + foff; + acb_srcptr cofac = weights + p1 * weights_nlogs * flen + (k - k1) * flen; + slong fklen = FLINT_MIN(flen, (f + k)->length - foff); + /* loop on p */ + acb_dot(dest, dest, 0, cofac, 1, src, 1, fklen, prec); + } + } + } +} + + +void +acb_ode_apply_diffop_basecase( + acb_poly_struct * g, + const acb_poly_struct * dop, slong dop_len, + acb_srcptr expo, slong offset, + const acb_poly_struct * f, + slong nlogs, + slong start, slong len, + slong prec) +{ + slong flen = 0; + for (int k = 0; k < nlogs; k++) + flen = FLINT_MAX(flen, (f + k)->length); + /* XXX could we take flen = FLINT_MIN(flen, len) or similar? */ + + acb_ptr weights = _acb_vec_init(len * nlogs * flen); + + _acb_ode_apply_diffop_basecase_weights( + weights, dop, dop_len, expo, offset, + flen, nlogs, start, len, prec); + + _acb_poly_vec_fit_length(g, nlogs, len); + _acb_ode_apply_diffop_basecase_precomp( + g, 0, + weights, nlogs, f, 0, flen, + nlogs, start, len, prec); + _acb_poly_vec_set_length(g, nlogs, len); + _acb_poly_vec_normalise(g, nlogs); + + _acb_vec_clear(weights, len * nlogs * flen); +} diff --git a/src/acb_ode/apply_diffop_polmul.c b/src/acb_ode/apply_diffop_polmul.c new file mode 100644 index 0000000000..6fc043d5e3 --- /dev/null +++ b/src/acb_ode/apply_diffop_polmul.c @@ -0,0 +1,106 @@ +#include "acb_types.h" +#include "acb_poly.h" +#include "acb_ode.h" + + +/* XXX take f + k instead of (f, k), handle the test on k in some other way? */ +static void +xD_logcoeff_inplace(acb_poly_struct * f, acb_srcptr expo, slong offset, slong k, + slong nlogs, slong prec) +{ + acb_poly_t tmp; + acb_poly_init(tmp); + + /* XXX handle rational expo more efficiently (-> gr?) */ + if (expo != NULL) + acb_poly_scalar_mul(tmp, f + k, expo, prec); + for (slong j = 0; j < (f + k)->length; j++) + acb_mul_ui((f + k)->coeffs + j, (f + k)->coeffs + j, offset + j, prec); + acb_poly_add(f + k, f + k, tmp, prec); + if (k + 1 < nlogs) + acb_poly_add(f + k, f + k, f + k + 1, prec); + + acb_poly_clear(tmp); +} + + +/* Computes the coefficients of x^start to x^{start+len-1} (inclusive) of the + * log-polynomial g such that dop(x^{expo+offset}·f) = x^{expo+offset}·g. + * + * expo may be NULL (treated as zero) + * + * The underscore version operates on the chunk of coefficients of length + * flen starting at offset foff in each of the components of the input vector, + * and _adds_ each component of g at offset goff in the corresponding component + * of the output vector. + */ + +/* XXX Take acb_struct**s instead of acb_poly_struct*s? More generally, how can + * we make the interface more natural / consistent with FLINT? */ + +void +_acb_ode_apply_diffop_polmul( + acb_poly_struct * g, slong goff, + const acb_poly_struct * dop, slong dop_len, + acb_srcptr expo, slong offset, + const acb_poly_struct * f, slong foff, slong flen, + slong nlogs, + slong start, slong len, + slong prec) +{ + acb_poly_t tmp; + acb_poly_init(tmp); + acb_poly_fit_length(tmp, start + len); + + /* XXX leave it to the caller to do this, and destroy f? */ + acb_poly_struct * curder = _acb_poly_vec_init(nlogs); + _acb_poly_vec_set_block(curder, f, nlogs, foff, flen); + + for (slong i = 0; i < dop_len; i++) + { + for (slong k = 0; k < nlogs; k++) + { + /* should be a mulmid in the typical case (and should accept a + * pretransformed first operand) */ + acb_poly_mullow(tmp, dop + i, curder + k, start + len, prec); + acb_poly_shift_right(tmp, tmp, start); + + FLINT_ASSERT (tmp->length <= len); + _acb_poly_add((g + k)->coeffs + goff, + (g + k)->coeffs + goff, len, + tmp->coeffs, tmp->length, + prec); + + /* curder[k] ← (x·d/dx(prev curder))[k] */ + xD_logcoeff_inplace(curder, expo, offset, k, nlogs, prec); + } + } + + _acb_poly_vec_clear(curder, nlogs); + acb_poly_clear(tmp); +} + + +/* XXX support aliasing */ +void +acb_ode_apply_diffop_polmul( + acb_poly_struct * g, + const acb_poly_struct * dop, slong dop_len, + acb_srcptr expo, slong offset, + const acb_poly_struct * f, + slong nlogs, + slong start, slong len, + slong prec) +{ + _acb_poly_vec_fit_length(g, nlogs, len); + _acb_ode_apply_diffop_polmul( + g, 0, + dop, dop_len, + expo, offset, + f, 0, start + len, + nlogs, + start, len, + prec); + _acb_poly_vec_set_length(g, nlogs, len); + _acb_poly_vec_normalise(g, nlogs); +} diff --git a/src/acb_ode/bound_init_clear.c b/src/acb_ode/bound_init_clear.c new file mode 100644 index 0000000000..34aff3741f --- /dev/null +++ b/src/acb_ode/bound_init_clear.c @@ -0,0 +1,32 @@ +#include "acb_ode.h" +#include "acb_poly.h" +#include "arb_poly.h" + + +void +acb_ode_bound_init(acb_ode_bound_t bound) +{ + bound->prec = -1; + arb_init(bound->cst); + arb_poly_init(bound->den_lbound); + bound->den_rt = NULL; + bound->den_rt_len = 0; + bound->all_nums = NULL; + bound->all_nums_len = 0; + bound->pol_part_len = 0; + acb_poly_init(bound->ind); + acb_ode_ind_lbound_init(bound->ind_lbound); + acb_ode_stairs_init(bound->stairs); +} + +void +acb_ode_bound_clear(acb_ode_bound_t bound) +{ + arb_clear(bound->cst); + arb_poly_clear(bound->den_lbound); + _arb_vec_clear(bound->den_rt, bound->den_rt_len); + _acb_poly_vec_clear(bound->all_nums, bound->all_nums_len); + acb_poly_clear(bound->ind); + acb_ode_ind_lbound_clear(bound->ind_lbound); + acb_ode_stairs_clear(bound->stairs); +} diff --git a/src/acb_ode/bound_precompute.c b/src/acb_ode/bound_precompute.c new file mode 100644 index 0000000000..c2016f17bc --- /dev/null +++ b/src/acb_ode/bound_precompute.c @@ -0,0 +1,194 @@ +#include "acb.h" +#include "acb_ode.h" +#include "acb_poly.h" +#include "arb_poly.h" +#include "fmpz_mat.h" +#include "gr_mat.h" + + +/* Write dop ∈ K[x][θ] as rcoeff[0] + θ·rcoeff[1] + ··· */ + +static void +dop_rcoeffs(acb_poly_struct * rcoeff, const acb_poly_struct * dop, slong len, + slong prec) +{ + fmpz_t c, pow; + fmpz_mat_t binom; + + fmpz_init(c); + fmpz_init(pow); + fmpz_mat_init(binom, len, len); + + gr_ctx_t fmpz_ctx; + gr_ctx_init_fmpz(fmpz_ctx); + GR_MUST_SUCCEED(gr_mat_pascal((gr_mat_struct *) binom, -1, fmpz_ctx)); + gr_ctx_clear(fmpz_ctx); + + slong xlen = 0; + for (slong i = 0; i < len; i++) + xlen = FLINT_MAX(xlen, (dop + i)->length); + + for (slong k = 0; k < len; k++) + { + acb_poly_struct * pol = rcoeff + k; + acb_poly_zero(pol); + acb_poly_fit_length(pol, xlen); + + for (slong j = 0; j < xlen; j++) + { + fmpz_one(pow); + for (slong i = 0; i < len - k; i++) + { + if ((dop + k + i)->length > j) + { + fmpz_mul(c, pow, fmpz_mat_entry(binom, k + i, i)); + acb_addmul_fmpz(pol->coeffs + j, (dop + k + i)->coeffs + j, + c, prec); + } + if (j == 0) + break; + fmpz_mul_si(pow, pow, -j); + } + } + + _acb_poly_set_length(pol, xlen); + _acb_poly_normalise(pol); + } + + fmpz_mat_clear(binom); + fmpz_clear(c); + fmpz_clear(pow); +} + +/* Write (rcoeffs[0] + θ·rcoeffs[1] + ··· + θ^{dop_len-1}·lc)·lc⁻¹ as + initial[0] + θ·initial[1] + ··· + + (rem_num[0] + θ·rem_num[1] + ···)·x^trunc·lc⁻¹ + where the polynomials initial[i] have degree < trunc. + We assume trunc >= 1, so that rem_num has length < dop_len. */ + +static void +split_dop(acb_poly_struct * initial, acb_poly_struct * rem_num, + const acb_poly_struct * rcoeffs, slong dop_len, + slong trunc, slong prec) +{ + acb_poly_t inv, rev; + + FLINT_ASSERT(trunc >= 1); + + acb_poly_init(inv); + acb_poly_init(rev); + + /* todo: when deg(lc) < trunc, we should only invert to order ~deg(lc)... */ + const acb_poly_struct * lc = rcoeffs + dop_len - 1; + acb_poly_inv_series(inv, lc, trunc, prec); + + for (slong i = 0; i < dop_len - 1; i++) + { + acb_poly_mullow(initial + i, rcoeffs + i, inv, trunc, prec); + acb_poly_mul(rem_num + i, initial + i, lc, prec); + acb_poly_sub(rem_num + i, rcoeffs + i, rem_num + i, prec); + /* drop known zeroes (todo: avoid computing them) */ + acb_poly_shift_right(rem_num + i, rem_num + i, trunc); + } + acb_poly_one(initial + dop_len - 1); + + acb_poly_clear(inv); + acb_poly_clear(rev); +} + +static void +swap_vars(acb_poly_struct * b, acb_poly_struct * a, slong len) +{ + for (slong i = 0; i < len; i++) + { + for (slong j = 0; j < (a + i)->length; j++) + { + acb_poly_fit_length(b + j, i + 1); + acb_swap((b + j)->coeffs + i, (a + i)->coeffs + j); + _acb_poly_set_length(b + j, i + 1); + _acb_poly_normalise(b + j); + } + } +} + +void +acb_ode_bound_precompute(acb_ode_bound_t bound, + const acb_poly_struct * dop, slong dop_len, + acb_srcptr lcrt, + slong pol_part_len, slong prec) +{ + const acb_poly_struct * den = dop + dop_len - 1; + + /* Denominator bound: 1/lc_Θ(dop) << cst/prod(rt-x) = cst/den_lbound */ + + if (prec > bound->prec) + { + bound->prec = prec; + + acb_abs(bound->cst, den->coeffs + den->length - 1, prec); + arb_inv(bound->cst, bound->cst, prec); + + slong rtlen = den->length - 1; + if (bound->den_rt != NULL) + _arb_vec_clear(bound->den_rt, bound->den_rt_len); + bound->den_rt = _arb_vec_init(rtlen); + bound->den_rt_len = rtlen; + + for (slong j = 0; j < rtlen; j++) + acb_get_abs_lbound_arf(arb_midref(bound->den_rt + j), + lcrt + j, prec); + + arb_poly_product_roots(bound->den_lbound, + bound->den_rt, den->length - 1, prec); + } + + /* Write dop with θ on the left, right-divide by lc, isolate pol_part_len + initial terms in the series expansion wrt x of the resulting fractions. + + We have order(initial) = order(dop), order(rem_num) < order(dop) (because + the leading coefficient of dop/lc is 1 and pol_part_len >= 1). + + todo: Avoid recomputations when refining a bound computed for a different + pol_part_len. */ + + if (pol_part_len <= 0) + pol_part_len = 1; + + if (prec > bound->prec || pol_part_len > bound->pol_part_len) + { + acb_poly_struct * rcoeffs = _acb_poly_vec_init(3*dop_len - 1); + acb_poly_struct * initial = rcoeffs + dop_len; + acb_poly_struct * rem_num = rcoeffs + 2*dop_len; + + dop_rcoeffs(rcoeffs, dop, dop_len, prec); + + /* Our pol_part_len = (that of ore_algebra) + 1, see comment below */ + split_dop(initial, rem_num, rcoeffs, dop_len, pol_part_len, prec); + + /* Rewrite the truncated part and the numerator of the remainder as + elements of K[θ][x], which may also be viewed as elements of + K[n][S⁻¹]. + + Note: here initial_n[0] = indicial polynomial, whereas in ore_algebra + this first element is dropped (so that only those that contribute to + the integrand in Algo 6.11, step 3 of [M19] remain) and all indices + are shifted by one. */ + + slong rem_num_clen = 0; + for (slong i = 0; i < dop_len - 1; i++) + rem_num_clen = FLINT_MAX(rem_num_clen, (rem_num + i)->length); + + if (bound->all_nums != NULL) + _acb_poly_vec_clear(bound->all_nums, bound->all_nums_len); + bound->pol_part_len = pol_part_len; + bound->all_nums_len = pol_part_len + rem_num_clen; + bound->all_nums = _acb_poly_vec_init(bound->all_nums_len); + + acb_poly_struct * initial_n = bound->all_nums; + acb_poly_struct * rem_num_n = bound->all_nums + pol_part_len; + swap_vars(initial_n, initial, dop_len); + swap_vars(rem_num_n, rem_num, dop_len - 1); + + _acb_poly_vec_clear(rcoeffs, 3*dop_len - 1); + } +} diff --git a/src/acb_ode/bound_precompute_group.c b/src/acb_ode/bound_precompute_group.c new file mode 100644 index 0000000000..a2029e6770 --- /dev/null +++ b/src/acb_ode/bound_precompute_group.c @@ -0,0 +1,41 @@ +#include "acb.h" +#include "acb_ode.h" +#include "acb_poly.h" + +/* todo: incremental version for increased pol_part_len (and maybe prec) but + * same group */ + +void +acb_ode_bound_precompute_group(acb_ode_bound_t bound, + const acb_poly_struct * dop, slong dop_len, + const acb_ode_exponents_t expos, + slong grp, // group index in expos + slong prec) +{ + /* flint_printf("== acb_ode_bound_precompute_group %{acb} ==\n", expos->groups[grp].leader); */ + + /* Monic indicial polynomial, centered at the current group leader */ + acb_poly_struct * ind = bound->all_nums; + acb_poly_taylor_shift(bound->ind, ind, expos->groups[grp].leader, prec); + FLINT_ASSERT(expos->groups[grp].shifts[0].n == 0); + _acb_vec_zero(bound->ind->coeffs, expos->groups[grp].shifts[0].mult); + FLINT_ASSERT(bound->ind->length == dop_len); + _acb_vec_scalar_div(bound->ind->coeffs, bound->ind->coeffs, dop_len - 1, + bound->ind->coeffs + dop_len - 1, prec); + acb_one(bound->ind->coeffs + dop_len - 1); + + /* flint_printf("ind=%{acb_poly}\n", bound->ind); */ + + /* Rat(n) denominator data */ + + acb_ode_ind_lbound_precompute(bound->ind_lbound, expos, grp, prec); + + /* Ignore the first entry of all_nums, anticipating the factor w⁻¹ under the + * integral in Algo 6.11, step 3 of [M19] */ + acb_ode_stairs_precompute(bound->stairs, + bound->all_nums + 1, bound->all_nums_len - 1, bound->ind, + expos->groups + grp, bound->ind_lbound, prec); + + /* flint_printf("stairs = %{mag*}\n", bound->stairs->h, bound->stairs->length); */ +} + diff --git a/src/acb_ode/bound_precompute_integrand.c b/src/acb_ode/bound_precompute_integrand.c new file mode 100644 index 0000000000..54fb2c0f27 --- /dev/null +++ b/src/acb_ode/bound_precompute_integrand.c @@ -0,0 +1,44 @@ +#include "arf.h" +#include "acb_ode.h" +#include "arb_poly.h" + + +/* [M19], Algorithm 6.1, step 4. */ + +void +acb_ode_bound_precompute_integrand(arb_poly_t itg_pol, arb_poly_t itg_num, + const acb_ode_group_t group, + const acb_ode_bound_t bound, + slong n0, slong nlogs, slong prec) +{ + /* flint_printf("== acb_ode_bound_precompute_integrand n0=%ld ==\n", n0); */ + + mag_ptr maj_all_nums = _mag_vec_init(bound->all_nums_len); + + /* In view of [M19], Eq. (5.9), one can take here for nlogs the actual + * length wrt log(x) of the solutions of interest truncated after their term + * in x^n, instead of τ(n0) as Algo 6.1, step 3 would prescribe. */ + acb_ode_bound_rat_vec(maj_all_nums, bound->all_nums + 1, bound->all_nums_len - 1, + bound->ind, group, bound->ind_lbound, bound->stairs, + n0, nlogs, prec); + + /* flint_printf("n0=%ld maj_all_nums = %{mag*}\n", n0, maj_all_nums, bound->all_nums_len - 1); */ + + arb_poly_fit_length(itg_pol, bound->pol_part_len - 1); + for (slong i = 0; i < bound->pol_part_len - 1; i++) + arf_set_mag(arb_midref(itg_pol->coeffs + i), maj_all_nums + i); + _arb_poly_set_length(itg_pol, bound->pol_part_len - 1); + _arb_poly_normalise(itg_pol); + + slong num_len = bound->all_nums_len - bound->pol_part_len; + arb_poly_fit_length(itg_num, num_len); + _arb_poly_set_length(itg_num, num_len); + for (slong i = 0; i < num_len; i++) + arf_set_mag(arb_midref(itg_num->coeffs + i), + maj_all_nums + bound->pol_part_len - 1 + i); + _arb_poly_normalise(itg_num); + + /* flint_printf("n0=%ld itg_pol=%{arb_poly} itg_num=%{arb_poly}\n", n0, itg_pol, itg_num); */ + + _mag_vec_clear(maj_all_nums, bound->all_nums_len); +} diff --git a/src/acb_ode/bound_rat_ordinary_vec.c b/src/acb_ode/bound_rat_ordinary_vec.c new file mode 100644 index 0000000000..c4e82edf02 --- /dev/null +++ b/src/acb_ode/bound_rat_ordinary_vec.c @@ -0,0 +1,169 @@ +#include "arf.h" +#include "acb.h" +#include "acb_ode.h" +#include "acb_poly.h" +#include "arb_poly.h" + + +/* A lower bound on prod[ind(α) = 0](|1-α/k|), where ind is the shifted indicial + * polynomial, valid for all k ≥ n with n, k ∈ ℕ ∖ exponents. + * Cf. [M19], Lemma 7.2 */ + +static void +ind_lbound_eval(mag_t res, const acb_ode_ind_lbound_t ind_lbound, slong n, + slong prec) +{ + acb_t a; + mag_t m; + + if (n == 0) + { + mag_zero(res); + return; + } + + acb_init(a); + mag_init(m); + + mag_one(res); + for (slong i = 0; i < ind_lbound->length; i++) + { + if (n < ind_lbound->r[i].n_min) + { + mag_mul_lower(res, res, ind_lbound->r[i].global_lbound); + } + else + { + acb_div_si(a, ind_lbound->r[i].root, n, prec); + acb_sub_si(a, a, 1, prec); + acb_get_mag_lower(m, a); + mag_pow_ui_lower(m, m, ind_lbound->r[i].mult); + + mag_mul_lower(res, res, m); + } + + /* flint_printf("root=%{acb} n=%ld n_min=%ld lb=%{mag} res=%{mag}\n", + ind_lbound->r[i].root, n, ind_lbound->r[i].n_min, + ind_lbound->r[i].global_lbound, res); */ + } + + mag_clear(m); + acb_clear(a); +} + +/* [M19], Algorithm 7.1 */ + +void +acb_ode_bound_rat_ordinary_vec(mag_ptr res, + const acb_poly_struct * num, slong len, + const acb_poly_t ind, + const acb_ode_ind_lbound_t ind_lbound, + slong n0, slong ord, slong prec) +{ + FLINT_ASSERT(n0 >= 0); + + acb_t invcst, range_acb; + mag_t a, invcstmag; + arb_poly_t range_jet0, range_tail; + acb_poly_t range_jet0_acb, range_tail_acb; + acb_poly_t rat_jet, ind_jet, common_jet; + + acb_init(invcst); + acb_init(range_acb); + mag_init(a); + mag_init(invcstmag); + arb_poly_init(range_jet0); + arb_poly_init(range_tail); + acb_poly_init(range_jet0_acb); + acb_poly_init(range_tail_acb); + acb_poly_init(rat_jet); + acb_poly_init(ind_jet); + acb_poly_init(common_jet); + + arb_ptr range = acb_realref(range_acb); + + /* [n0,∞]⁻¹ */ + mag_set_ui(arb_radref(range), 2*n0); + mag_inv(arb_radref(range), arb_radref(range)); + arf_set_mag(arb_midref(range), arb_radref(range)); + + /* 1/(1+[n0,∞]⁻¹ε) */ + arb_poly_one(range_jet0); + arb_poly_set_coeff_arb(range_jet0, 1, range); + arb_poly_inv_series(range_jet0, range_jet0, ord, prec); + acb_poly_set_arb_poly(range_jet0_acb, range_jet0); + + /* 1/([n0,∞]+ε) = range + range_tail + O(x^ord) */ + arb_poly_fit_length(range_tail, range_jet0->length); + _arb_vec_scalar_mul(range_tail->coeffs + 1, range_jet0->coeffs + 1, + range_jet0->length - 1, range, prec); + _arb_poly_set_length(range_tail, range_jet0->length); + _arb_poly_normalise(range_tail); + acb_poly_set_arb_poly(range_tail_acb, range_tail); + + /* flint_printf("n0=%ld ord=%ld range_acb=%{acb} range_jet0_acb=%{acb_poly} range_tail_acb=%{acb_poly}\n0", n0, ord, range_acb, range_jet0_acb, range_tail_acb); */ + + /* 1/([n0,∞]+ε)^r · ind([n0,∞]+ε) */ + // could precompute the reciprocal polynomial (or use a shallow copy) + slong ind_jet_len0 = FLINT_MAX(ind->length, ord); + acb_poly_fit_length(ind_jet, ind_jet_len0); + _acb_poly_reverse(ind_jet->coeffs, ind->coeffs, ind->length, ind->length); + _acb_poly_set_length(ind_jet, ind->length); + _acb_poly_normalise(ind_jet); + /* flint_printf("rcpq_den=%{acb_poly}\n0", ind_jet); */ + /* Use the factored form of ind? But ind_lbound currently stores only a + * subset of the roots. */ + _acb_poly_taylor_shift(ind_jet->coeffs, range_acb, ind_jet_len0, prec); + acb_poly_compose_series(ind_jet, ind_jet, range_tail_acb, ord, prec); + + /* flint_printf("ind_jet=%{acb_poly}\n0", ind_jet); */ + + ind_lbound_eval(invcstmag, ind_lbound, n0, prec); + mag_inv(invcstmag, invcstmag); + acb_add_error_mag(invcst, invcstmag); /* XXX overkill? */ + _acb_vec_scalar_mul(ind_jet->coeffs + 1, ind_jet->coeffs + 1, + ind_jet->length - 1, invcst, prec); + acb_one(ind_jet->coeffs); + + /* flint_printf("invcstmag=%{mag} invcst=%{acb} improved ind_jet=%{acb_poly}\n0", invcstmag, invcst, ind_jet); */ + + /* 1/[n0,∞]·([n0,∞]+ε)^{r-1}/ind([n0,∞]+ε) */ + acb_poly_inv_series(common_jet, ind_jet, ord, prec); + acb_poly_mullow(common_jet, common_jet, range_jet0_acb, ord, prec); + + for (slong i = 0; i < len; i++) + { + acb_poly_fit_length(rat_jet, FLINT_MAX(ind->length, ord)); + _acb_poly_reverse(rat_jet->coeffs, (num + i)->coeffs, + FLINT_MIN(ind->length - 1, (num + i)->length), + ind->length - 1); + _acb_poly_set_length(rat_jet, ind->length); + _acb_poly_normalise(rat_jet); + /* flint_printf("i=%ld rcpqnum=%{acb_poly} range_tail_acb=%{acb_poly}\n0", i, rat_jet, range_tail_acb); */ + _acb_poly_taylor_shift(rat_jet->coeffs, range_acb, + FLINT_MAX(rat_jet->length, ord), prec); + acb_poly_compose_series(rat_jet, rat_jet, range_tail_acb, ord, prec); + /* flint_printf("num=%{acb_poly}\n0", rat_jet); */ + acb_poly_mullow(rat_jet, rat_jet, common_jet, ord, prec); + mag_zero(res + i); + for (slong j = 0; j < rat_jet->length; j++) + { + acb_get_mag(a, rat_jet->coeffs + j); + mag_add(res + i, res + i, a); + } + mag_mul(res + i, res + i, invcstmag); + } + + acb_clear(invcst); + acb_clear(range_acb); + mag_clear(a); + mag_clear(invcstmag); + arb_poly_clear(range_jet0); + arb_poly_clear(range_tail); + acb_poly_clear(range_jet0_acb); + acb_poly_clear(range_tail_acb); + acb_poly_clear(rat_jet); + acb_poly_clear(ind_jet); + acb_poly_clear(common_jet); +} + diff --git a/src/acb_ode/bound_rat_ref_vec.c b/src/acb_ode/bound_rat_ref_vec.c new file mode 100644 index 0000000000..e446ab879d --- /dev/null +++ b/src/acb_ode/bound_rat_ref_vec.c @@ -0,0 +1,66 @@ +#include "acb.h" +#include "acb_ode.h" +#include "acb_poly.h" + + +static void +acb_poly_taylor_shift_si_trunc(acb_poly_t res, const acb_poly_t pol, + const slong n, slong ord, slong prec) +{ + /* todo: faster code, especially for small ord */ + acb_t a; + acb_init(a); + + acb_set_si(a, n); + acb_poly_taylor_shift(res, pol, a, prec); + acb_poly_truncate(res, ord); + + acb_clear(a); +} + + +void +acb_ode_bound_rat_ref_vec(mag_ptr res, + const acb_poly_struct * num, slong len, const acb_poly_t ind, + slong n, slong mult, slong ord, slong prec) +{ + acb_poly_t inv, ser; + mag_t m; + + if (n == 0) + { + for (slong i = 0; i < len; i++) + mag_zero(res + i); + return; + } + + acb_poly_init(inv); + acb_poly_init(ser); + mag_init(m); + + acb_poly_taylor_shift_si_trunc(inv, ind, n, ord + mult, prec); + acb_poly_shift_right(inv, inv, mult); + acb_poly_inv_series(inv, inv, ord, prec); + + /* flint_printf("ref n=%ld ord=%ld mult=%ld inv=%{acb_poly}\n", n, ord, mult, inv); */ + + for (slong i = 0; i < len; i++) + { + acb_poly_taylor_shift_si_trunc(ser, num + i, n, ord, prec); + acb_poly_mullow(ser, ser, inv, ord, prec); + + /* flint_printf("n=%ld i=%ld ser=%{acb_poly}\n", n, i, ser); */ + + mag_zero(res + i); + for (slong j = 0; j < ser->length; j++) + { + acb_get_mag(m, ser->coeffs + j); + mag_add(res + i, res + i, m); + } + mag_mul_ui(res + i, res + i, n); + } + + acb_poly_clear(inv); + acb_poly_clear(ser); + mag_clear(m); +} diff --git a/src/acb_ode/bound_rat_vec.c b/src/acb_ode/bound_rat_vec.c new file mode 100644 index 0000000000..a6d2862bbd --- /dev/null +++ b/src/acb_ode/bound_rat_vec.c @@ -0,0 +1,54 @@ +#include "acb_ode.h" +#include "mag.h" + + +static void +acb_ode_bound_rat_exceptional_vec(mag_ptr res, slong len, + const acb_ode_group_t group, + const acb_ode_stairs_t stairs, slong n) +{ + for (slong i = 0; i < len; i++) + { + slong s = 0; + while (s < group->nshifts && group->shifts[s].n < n) + s++; + if (s < group->nshifts) + mag_set(res + i, stairs->h + i * group->nshifts + s); + else + mag_zero(res + i); + } +} + +/* [M19], Algorithm 7.4, steps 2 to 4, with the slight refinement that + * τ(n₀) = ord can differ from sum(m(k), k=0..n₀). */ + +void +acb_ode_bound_rat_vec(mag_ptr res, + const acb_poly_struct * num, slong len, + const acb_poly_t ind, const acb_ode_group_t group, + const acb_ode_ind_lbound_t ind_lbound, + const acb_ode_stairs_t stairs, + slong n0, slong ord, slong prec) +{ + /* Bound at indices >= (smallest exceptional index > n0, if any). This term + * does not depend on ord (possible values are precomputed in stairs using + * worst-case bounds for the degrees in log(x) of the solutions). */ + acb_ode_bound_rat_exceptional_vec(res, len, group, stairs, n0); + + for (slong s = 0; s < group->nshifts; s++) + if (group->shifts[s].n == n0) + return; + + mag_ptr ordinary = _mag_vec_init(len); + + /* Bound at indices up < (smallest exceptional index > n0, if any). + * Past the last exceptional index, this is the only contribution. */ + acb_ode_bound_rat_ordinary_vec(ordinary, num, len, ind, ind_lbound, n0, ord, + prec); + + for (slong i = 0; i < len; i++) + mag_max(res + i, res + i, ordinary + i); + + _mag_vec_clear(ordinary, len); +} + diff --git a/src/acb_ode/choose_prec.c b/src/acb_ode/choose_prec.c new file mode 100644 index 0000000000..4c1818e4a6 --- /dev/null +++ b/src/acb_ode/choose_prec.c @@ -0,0 +1,85 @@ +#include + +#include "acb_types.h" +#include "acb_poly.h" +#include "acb_ode.h" +#include "ulong_extras.h" + + +/* XXX + * - include overestimation of tail bound? how? + * - cap to input accuracy? */ +slong +acb_ode_choose_prec(slong * rec_prec, const acb_poly_struct * dop, slong dop_len, + mag_srcptr rad, mag_srcptr cvrad, slong tgt_prec) +{ + slong sum_prec, prec0, prec; + double lgmag, nterms; + mag_t ratio; + + mag_init(ratio); + + prec0 = tgt_prec + 2*dop_len; + + /* rough estimate of bits lost to overestimation when unrolling a recurrence + * in interval arithmetic (todo: linear error analysis) */ + slong dop_clen = _acb_poly_vec_length(dop, dop_len); + prec = 8 + prec0 * (1 + n_clog(FLINT_MAX(1, dop_clen - 3), 2)); + + lgmag = 0.; + if (mag_is_zero(rad)) + { + nterms = 1; //??? + } + else if (mag_is_finite(cvrad)) + { + mag_div_lower(ratio, cvrad, rad); + nterms = prec / mag_get_d_log2_approx(ratio); + flint_printf("rad=%{mag} cvrad=%{mag} prec0=%ld nterms=%ld\n", rad, cvrad, prec0, nterms); + } + else + { + /* estimate hump height based on worst-case asymptotic behavior */ + mag_t order, base; + + mag_init(order); + mag_init(base); + + _acb_ode_solution_growth(order, base, dop, dop_len); + flint_printf("order=%{mag} base=%{mag}\n", order, base); + if (mag_is_zero(order)) + nterms = 1; + else + { + /* mag so we don't have to worry about overflows */ + mag_mul(base, base, rad); + + double base_d = mag_get_d(base); + double order_d = mag_get_d(order); + + double hump = exp(log(base_d) * order_d + 1.); + /* cap the cancellation we'll attempt to absorb */ + hump = FLINT_MIN(hump, 100. + prec * log2(prec)); + double den = log2(prec) / order_d; + if (base_d < 1.) + den = FLINT_MAX(den, - log2(base_d)); + nterms = hump + prec / den; + lgmag = FLINT_MAX(0., order_d * hump * log2(hump)); + + flint_printf("prec0=%ld lgmag=%f nterms=%f base_d=%f hump=%f\n", prec0, lgmag, nterms, base_d, hump); + } + } + + sum_prec = 8 + 1.125 * (prec0 + 2 * lgmag + log2(nterms)); + prec = FLINT_MAX(prec, sum_prec); + + flint_printf("initial prec=%ld sum_prec=%ld\n", prec, sum_prec); + + mag_clear(ratio); + + if (rec_prec != NULL) + * rec_prec = prec; + return sum_prec; +} + + diff --git a/src/acb_ode/cvest.c b/src/acb_ode/cvest.c new file mode 100644 index 0000000000..b8dbc648d5 --- /dev/null +++ b/src/acb_ode/cvest.c @@ -0,0 +1,73 @@ +#include + +#include "acb_ode.h" +#include "arf.h" + +void +acb_ode_cvest_init(acb_ode_cvest_t cvest) +{ + cvest->neglogterm = 0.; + cvest->cvg_rate = 0.; + cvest->terms_wanted = INFINITY; + cvest->accuracy = WORD_MAX; // (?) + cvest->loss_rate = 1.; // (?) + cvest->terms_full_prec = INFINITY; + cvest->prec_wanted = ARF_PREC_EXACT; +} + +void +acb_ode_cvest_clear(acb_ode_cvest_t cvest) +{ + return; +} + +void +acb_ode_cvest_update(acb_ode_cvest_t cvest, /* aliasing allowed */ + const acb_ode_cvest_t old, + mag_t est, slong accuracy, + slong stride, + slong prec, slong work_prec) +{ + acb_ode_cvest_struct _old = * old; + + /* estimated absolute accuracy of the sum in bits based on the magnitude of + the last few terms */ + cvest->neglogterm = - mag_get_d_log2_approx(est); + + /* running estimate of convergence rate in bits/term, + exponential discounting */ + double cvg_rate_stride = (cvest->neglogterm - _old.neglogterm)/stride; + cvest->cvg_rate = _old.cvg_rate/2 + cvg_rate_stride/2; + + /* number of terms still needed for full accuracy at the current rate of + convergence */ + cvest->terms_wanted = cvest->cvg_rate > 0. + ? (prec - cvest->neglogterm)/cvest->cvg_rate + : INFINITY; + + /* relative accuracy of the last few terms */ + cvest->accuracy = accuracy; + + /* running estimate of interval growth (lost bits/term), + exponential discounting + (work_prec may have decreased since _old.accuracy was computed) */ + double loss_stride = FLINT_MIN(_old.accuracy, work_prec) - cvest->accuracy; + double loss_rate_stride = FLINT_MAX(0., loss_stride)/stride; + cvest->loss_rate = _old.loss_rate/2 + loss_rate_stride/2; + + /* estimated number of terms that can be computed before interval radii + exceed the target error */ + double abs_acc = cvest->neglogterm + cvest->accuracy; + // todo: double abs_tol = prec - log(sum_est); ? + cvest->terms_full_prec = (abs_acc - prec)/cvest->loss_rate; + + /* suggested working precision */ + if (cvest->terms_full_prec <= cvest->terms_wanted + || cvest->loss_rate >= _old.loss_rate * 1.0125) + cvest->prec_wanted = ARF_PREC_EXACT; + else if (cvest->terms_wanted > 0) + cvest->prec_wanted = prec - 0.875 * cvest->neglogterm + + cvest->terms_wanted * cvest->loss_rate; + else + cvest->prec_wanted = prec + 4; +} diff --git a/src/acb_ode/exponents.c b/src/acb_ode/exponents.c new file mode 100644 index 0000000000..907c64cfce --- /dev/null +++ b/src/acb_ode/exponents.c @@ -0,0 +1,252 @@ +#include +#include +#include "acb.h" +#include "acb_ode.h" +#include "fmpz.h" +#include "gr_poly.h" +#include "gr_vec.h" + +void +acb_ode_group_init(acb_ode_group_t group, slong len) +{ + acb_init(group->leader); + group->shifts = flint_malloc(sizeof(acb_ode_shift_struct) * len); + group->nshifts = 0; +} + +void +acb_ode_group_clear(acb_ode_group_t group) +{ + acb_clear(group->leader); + flint_free(group->shifts); +} + +void +acb_ode_group_set(acb_ode_group_t dest, const acb_ode_group_t src) +{ + acb_set(dest->leader, src->leader); + memcpy(dest->shifts, src->shifts, + src->nshifts * sizeof(acb_ode_shift_struct)); + dest->nshifts = src->nshifts; +} + +slong +acb_ode_group_length(const acb_ode_group_t grp) +{ + slong len = 0; + + for (slong s = 0; s < grp->nshifts; s++) + len += grp->shifts[s].mult; + + return len; +} + +slong +acb_ode_group_multiplicity (const acb_ode_group_t group, slong n) +{ + for (slong s = 0; s < group->nshifts; s++) + if (group->shifts[s].n == n) + return group->shifts[s].mult; + return 0; +} + +slong +acb_ode_group_nlogs (const acb_ode_group_t group, slong n) +{ + slong nlogs = 0; + for (slong s = 0; s < group->nshifts; s++) { + if (group->shifts[s].n > n) + break; + nlogs += group->shifts[s].mult; + } + return nlogs; +} + +void +acb_ode_exponents_init(acb_ode_exponents_t expos) +{ + expos->ngroups = 0; + expos->groups = NULL; +} + +void +acb_ode_exponents_clear(acb_ode_exponents_t expos) +{ + for (slong g = 0; g < expos->ngroups; g++) + /* group objects referenced by exponents objects share data, so we + cannot call group_clear */ + acb_clear(expos->groups[g].leader); + flint_free(expos->groups); +} + + +void +acb_ode_exponents_ordinary(acb_ode_exponents_t expos, slong order) +{ + slong off = sizeof(acb_ode_group_struct); + slong bufsz = off + sizeof(acb_ode_shift_struct) * order; + unsigned char * buf = flint_realloc(expos->groups, bufsz); + expos->groups = (acb_ode_group_struct *) buf; + expos->groups->shifts = (acb_ode_shift_struct *) (buf + off); + + acb_zero(expos->groups->leader); + for (slong n = 0; n < order; n++) + { + expos->groups->shifts[n].n = n; + expos->groups->shifts[n].mult = 1; + } + + expos->ngroups = 1; +} + + +/* todo: also provide a function for computing the exponents starting from the + indicial polynomial */ +int +acb_ode_exponents(acb_ode_exponents_t expos, const gr_ore_poly_t dop, + gr_ctx_t Dop, gr_ctx_t CC) +{ + int status = GR_SUCCESS; + + gr_ctx_struct * Scalars, * Pol; + gr_ctx_t ZZ, ZZvec; + gr_poly_t ind, polc; + gr_vec_t fac, mult, slfac, slshifts, slmult, roots, rtmult; + + if (gr_ore_poly_ctx_over_gr_poly_base_ptrs(&Scalars, &Pol, Dop) != GR_SUCCESS) + flint_abort(); + /* return GR_UNABLE; */ + + /* For now exponents are of type acb_t, but this might change */ + if (CC->which_ring != GR_CTX_CC_ACB) + flint_abort(); + /* return GR_UNABLE; */ + + gr_ctx_init_fmpz(ZZ); + gr_ctx_init_vector_gr_vec(ZZvec, ZZ); + + gr_poly_init(ind, Scalars); + gr_poly_init(polc, Scalars); /* unused */ + gr_vec_init(fac, 0, Pol); + gr_vec_init(mult, 0, ZZ); + gr_vec_init(slfac, 0, Pol); + gr_vec_init(slshifts, 0, ZZvec); + gr_vec_init(slmult, 0, ZZvec); + gr_vec_init(roots, 0, CC); + gr_vec_init(rtmult, 0, ZZ); /* unused */ + + /* Compute the local exponents (indicial roots) and organize them into shift + equivalence classes */ + + status |= gr_ore_poly_indicial_polynomial(ind, dop, Dop); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + // GR_MUST_SUCCEED(status); // tmp + + /* XXX In some cases this should use the irreducible factorization for + representing the exponents exactly; in others, this should compute the + roots of the shiftless factors without further decomposing them. + + On pourrait (au moins dans la plupart des cas) continuer même si + gr_factor échoue, mais dans l'état actuel des choses on n'a pas de + décomposition shiftless qui fonctionne sans factorisation. */ + + status |= gr_factor(polc, fac, mult, ind, 0, Pol); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + // GR_MUST_SUCCEED(status); // tmp + status |= gr_poly_shiftless_decomposition_from_factors( + slfac, slshifts, slmult, fac, mult, Scalars); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + // GR_MUST_SUCCEED(status); // tmp + + /* + gr_ctx_t Polvec, ZZvecvec; + gr_ctx_init_vector_gr_vec(Polvec, Pol); + gr_ctx_init_vector_gr_vec(ZZvecvec, ZZvec); + flint_printf("sldec: %{gr} %{gr} %{gr}\n", slfac, Polvec, slshifts, + ZZvecvec, slmult, ZZvecvec); + gr_ctx_clear(ZZvecvec); + gr_ctx_clear(Polvec); + */ + + /* Allocate space for sum(deg(slfac)) exponents and sum(len(shiftset)) + shifts; exponents belonging to the same irred factor will share a shift + set */ + + expos->ngroups = 0; + for (slong i = 0; i < slfac->length; i++) + { + gr_poly_struct * indfactor = gr_vec_entry_ptr(slfac, i, Pol); + expos->ngroups += indfactor->length - 1; + } + + size_t shifts_offset = expos->ngroups * sizeof(acb_ode_group_struct); + size_t defect = shifts_offset % alignof(acb_ode_shift_struct); + if (defect != 0) + shifts_offset += alignof(acb_ode_shift_struct) - defect; + + size_t bufsz = shifts_offset; + for (slong i = 0; i < slshifts->length; i++) + { + gr_vec_struct * shi = GR_ENTRY(slshifts->entries, i, sizeof(gr_vec_struct)); + bufsz += shi->length * sizeof(acb_ode_shift_struct); + } + + unsigned char * buf = flint_realloc(expos->groups, bufsz); + + expos->groups = (acb_ode_group_struct *) buf; + acb_ode_shift_struct * shifts = (acb_ode_shift_struct *) (buf + shifts_offset); + + /* Populate the data structure */ + + for (slong i = 0, g = 0; i < slfac->length; i++) + { + gr_poly_struct * indfactor = gr_vec_entry_ptr(slfac, i, Pol); + gr_vec_struct * shi = gr_vec_entry_ptr(slshifts, i, ZZvec); + gr_vec_struct * mi = gr_vec_entry_ptr(slmult, i, ZZvec); + + status |= gr_poly_roots_other(roots, rtmult, indfactor, Scalars, 0, CC); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + FLINT_ASSERT(roots->length + 1 == indfactor->length || status != GR_SUCCESS); + + slong len = shi->length; + slong maxshift = fmpz_get_si((fmpz *) shi->entries + len - 1); + for (slong s = 0; s < len; s++) + { + /* The shifts in the shiftless decomposition apply to the variable + of the indicial polynomial; we want shifts on the roots */ + shifts[s].n = maxshift - fmpz_get_si((fmpz *) shi->entries + len - 1 - s); + shifts[s].mult = fmpz_get_si((fmpz *) mi->entries + len - 1 - s); + } + + for (slong j = 0; j < roots->length; j++, g++) + { + gr_init(expos->groups[g].leader, CC); + status |= gr_sub_si(expos->groups[g].leader, + GR_ENTRY(roots->entries, j, CC->sizeof_elem), + maxshift, CC); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + /* XXX exponents share shift pointers but not shift counts... */ + expos->groups[g].nshifts = shi->length; + expos->groups[g].shifts = shifts; + } + + shifts += shi->length; + } + + FLINT_ASSERT((unsigned char *) shifts == buf + bufsz); + + gr_vec_clear(rtmult, ZZ); + gr_vec_clear(roots, CC); + gr_vec_clear(slmult, ZZvec); + gr_vec_clear(slshifts, ZZvec); + gr_vec_clear(slfac, Pol); + gr_vec_clear(mult, ZZ); + gr_vec_clear(fac, Pol); + gr_poly_clear(polc, Scalars); + gr_poly_clear(ind, Scalars); + + gr_ctx_clear(ZZvec); + gr_ctx_clear(ZZ); + + return status; +} diff --git a/src/acb_ode/fundamental_matrix.c b/src/acb_ode/fundamental_matrix.c new file mode 100644 index 0000000000..669f8f7d7b --- /dev/null +++ b/src/acb_ode/fundamental_matrix.c @@ -0,0 +1,329 @@ +#include "acb_types.h" +#include "acb_poly.h" +#include "acb_mat.h" +#include "acb_ode.h" + +#include "gr.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" +#include "gr_vec.h" + + +static slong +col(const acb_ode_group_struct * grp, slong s, slong k) +{ + slong j = 0; + for (slong s1 = 0; s1 <= s; s1++) + j += grp->shifts[s1].mult; + return j - 1 - k; +} + + +static void +fill_column(acb_mat_t mat, + const acb_ode_group_struct * grp, slong s, slong k, + const acb_poly_struct * val) +{ + slong j = col(grp, s, k); + slong len = FLINT_MIN(acb_mat_nrows(mat), val->length); + + slong i = 0; + for (; i < len; i++) + { + acb_ptr a = acb_mat_entry(mat, i, j); + acb_swap(a, val->coeffs + i); + } + for (; i < acb_mat_nrows(mat); i++) + acb_zero(acb_mat_entry(mat, i, j)); +} + + +static void +fix_column_echelon(acb_mat_struct * mat, + const acb_ode_group_struct * grp, slong s, slong k, + const acb_mat_t extini, slong prec) +{ + slong mult = grp->shifts[s].mult; + slong delta = mult - 1 - k; + slong j = col(grp, s, k); + + for (slong s1 = s + 1; s1 < grp->nshifts; s1++) + { + slong j1 = j1 = col(grp, s1, 0); + slong mult1 = grp->shifts[s1].mult; + for (slong k1 = FLINT_MAX(0, mult1 - delta); k1 < mult1; k1++) + { + acb_ptr cc = acb_mat_entry(extini, s1, k1 + delta); + /* flint_printf("fix: s=%wd k=%wd s1=%wd k1=%wd j=%wd j1=%d cc=%{acb}\n", s, k, s1, k1, j, j1, cc); */ + for (slong i = 0; i < acb_mat_nrows(mat); i++) + { + acb_ptr a = acb_mat_entry(mat, i, j); + acb_ptr a1 = acb_mat_entry(mat, i, j1); + /* flint_printf("i=%wd %{acb} - %{acb}", i, a, a1); */ + acb_submul(a, cc, a1, prec); + /* flint_printf(" = %{acb}\n", a); */ + } + j1++; + } + } + +} + + +static void +fill_group(acb_mat_t mat, const acb_ode_sum_t sum, + const acb_ode_group_struct * grp, slong p, + acb_ode_basis_t basis, slong prec) +{ + for (slong s = grp->nshifts - 1; s >= 0; s--) + { + acb_ode_sol_struct * sol = sum->sol + s; + acb_poly_struct * val = _acb_poly_vec_init(sol->nlogs); + + slong mult = grp->shifts[s].mult; + acb_ode_sol_jet(val, sum->group->leader, sol, p, sum->pts + p, + acb_mat_nrows(mat), mult, prec); + // flint_printf("s=%wd mult=%wd val=%{acb_poly}\n\n", s, mult, val); + + for (slong k = 0; k < mult; k++) + { + slong delta = mult - 1 - k; + fill_column(mat, grp, s, k, val + delta); + + // flint_printf("s=%wd k=%wd after fill \n%{acb_mat}\n\n", s, k, mat); + + switch (basis) + { + case ACB_ODE_BASIS_CASCADE: + break; + case ACB_ODE_BASIS_ECHELON: + fix_column_echelon(mat, grp, s, k, sum->sol[s].extini, prec); + break; + default: + FLINT_ASSERT(0); + } + + // flint_printf("s=%wd k=%wd after fix \n%{acb_mat}\n\n", s, k, mat); + } + + _acb_poly_vec_clear(val, mult); + } +} + + +void +_acb_ode_fundamental_matrix_vec( + acb_mat_struct * mat, + const acb_poly_struct * dop, slong dop_len, + const acb_ode_exponents_t expos, + acb_srcptr lcrt, + acb_srcptr pts, slong npts, + acb_ode_basis_t basis, + acb_ode_sum_worker_t sum_worker, + slong prec) +{ + mag_t cvrad, mag; + acb_ode_bound_t bound; + + mag_init(cvrad); + mag_init(mag); + acb_ode_bound_init(bound); + + slong dop_clen = _acb_poly_vec_length(dop, dop_len); + + slong nder = 0; + for (slong p = 0; p < npts; p++) + nder = FLINT_MAX(nder, acb_mat_nrows(mat + p)); + + if (nder == 0) + return; + + // XXX somehow make pol_part_len and bounds_prec customizable + slong pol_part_len = dop_clen / 2 + 3; + acb_ode_bound_precompute(bound, dop, dop_len, lcrt, pol_part_len, MAG_BITS); + + _acb_vec_get_mag_lower(cvrad, lcrt, dop[dop_len - 1].length - 1); + + /* XXX maybe not the best place to test this; + redundant with sum_set_points */ + _acb_vec_get_mag(mag, pts, npts); + flint_printf("mag=%{mag} cvrad=%{mag}\n", mag, cvrad); + if (mag_cmp(mag, cvrad) >= 0) + { + for (slong p = 0; p < npts; p++) + acb_mat_indeterminate(mat + p); + goto cleanup; + } + + for (slong g = 0, j = 0; g < expos->ngroups; g++) + { + acb_ode_sum_t sum; + + acb_ode_group_struct * group = expos->groups + g; + slong glen = acb_ode_group_length(expos->groups + g); + + acb_ode_bound_precompute_group(bound, dop, dop_len, expos, g, MAG_BITS); + + acb_ode_sum_init(sum, dop_len, npts, group->nshifts, nder); + acb_ode_sum_set_diffop(sum, dop, dop_len, cvrad); + sum->bound = bound; /* XXX handle in set_diffop or move out of sum */ + acb_ode_sum_set_group(sum, group); + /* todo: This only avoids part of the redundant computations between + solutions of the same group in the presence of logs. */ + acb_ode_sum_set_ini_highest(sum); + acb_ode_sum_set_points(sum, pts, npts); + + sum_worker(sum, -1, prec); + + for (slong p = 0; p < npts; p++) + { + acb_mat_t win; + acb_mat_window_init(win, mat + p, 0, j, acb_mat_nrows(mat + p), j + glen); + + fill_group(win, sum, expos->groups + g, p, basis, prec); + + acb_mat_window_clear(win); + } + + acb_ode_sum_clear(sum); + + j += glen; + } + +cleanup: + acb_ode_bound_clear(bound); + mag_clear(cvrad); + mag_clear(mag); +} + + +int +acb_ode_fundamental_matrix_vec( + acb_mat_struct * mat, + const gr_ore_poly_t dop, gr_ore_poly_ctx_t Dop, + const acb_ode_exponents_t expos, + acb_srcptr lcrt, + acb_srcptr pts, slong npts, + acb_ode_basis_t basis, + slong prec) +{ + gr_ctx_struct * Scalars, * Pol; + gr_ctx_t CC, CCx, CCxT, bCC; + acb_ode_exponents_t _expos; + gr_ore_poly_struct * _dop; + gr_vec_t sing; + acb_struct * _lcrt = NULL; + + if (gr_ore_poly_ctx_over_gr_poly_base_ptrs(&Scalars, &Pol, Dop) != GR_SUCCESS) + // flint_abort(); + return GR_UNABLE; + + /* todo: support d/dx (handle lcrt) */ + if (GR_ORE_POLY_CTX(Dop)->which_algebra != ORE_ALGEBRA_EULER_DERIVATIVE) + // flint_abort(); + return GR_UNABLE; + + if (dop->length == 0) + return GR_DOMAIN; + if (dop->length == 1) + return GR_SUCCESS; + if (gr_ore_poly_is_zero(dop, Dop) != T_FALSE) + return GR_UNABLE; + + gr_ctx_init_complex_acb(bCC, MAG_BITS); // XXX prec choice? + gr_vec_init(sing, 0, bCC); + + gr_poly_struct * lc = (gr_poly_struct *) (dop->coeffs) + dop->length - 1; + + int status = GR_SUCCESS; + + gr_ctx_init_complex_acb(CC, prec); + gr_ctx_init_gr_poly(CCx, CC); + gr_ore_poly_ctx_init(CCxT, CCx, 0, ORE_ALGEBRA_EULER_DERIVATIVE); + + GR_TMP_INIT(_dop, CCxT); + + char * gen = NULL; + status |= gr_ctx_gen_name(&gen, 0, Pol); + status |= gr_ctx_set_gen_name(CCx, gen); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + flint_free(gen); + + status |= gr_ctx_set_gen_name(CCxT, GR_ORE_POLY_CTX(Dop)->var); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + + /* XXX reorganize to do this at the *working* precision */ + /* todo: use gr_set_other once its gr_ore_poly version is powerful enough */ + status |= gr_poly_set_gr_poly_other((gr_ptr) _dop, (gr_ptr) dop, Pol, CCx); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + + if (expos == NULL) + { + acb_ode_exponents_init(_expos); + status |= acb_ode_exponents(_expos, dop, Dop, CC); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + expos = _expos; + } + + if (lcrt == NULL) + { + gr_ctx_t ZZ; + gr_vec_t sing_mult; + gr_ctx_init_fmpz(ZZ); + gr_vec_init(sing_mult, 0, ZZ); + + /* todo: should use (multiple) root enclosures supporting polynomials + with interval coefficients */ + status |= gr_poly_roots_other(sing, sing_mult, lc, + Scalars, 0, bCC); + if(status) flint_printf("%s:%d status=%d\n", __FILE__, __LINE__, status); + + _lcrt = flint_malloc(sizeof(acb_struct) * (lc->length - 1)); + for (slong j0 = 0; j0 < sing->length; j0++) + { + slong m = fmpz_get_si((fmpz *) (sing_mult->entries) + j0); + for (slong j1 = 0; j1 < m; j1++) + _lcrt[j0 + j1] = ((acb_ptr) sing->entries)[j0]; + } + lcrt = _lcrt; + + gr_vec_clear(sing_mult, ZZ); + gr_ctx_clear(ZZ); + } + + flint_printf("lcrt=%{acb*}\n", lcrt, sing->length); + + if (status == GR_SUCCESS) + _acb_ode_fundamental_matrix_vec(mat, _dop->coeffs, _dop->length, + expos, lcrt, + pts, npts, + basis, acb_ode_sum_divconquer, + prec); + + flint_free(_lcrt); + gr_vec_clear(sing, bCC); + if (expos == _expos) + acb_ode_exponents_clear(_expos); + GR_TMP_CLEAR(_dop, CCxT); + gr_ctx_clear(bCC); + gr_ctx_clear(CCxT); + gr_ctx_clear(CCx); + gr_ctx_clear(CC); + + return status; +} + + +int +acb_ode_fundamental_matrix( + acb_mat_t mat, + const gr_ore_poly_t dop, gr_ore_poly_ctx_t Dop, + const acb_ode_exponents_t expos, + acb_srcptr lcrt, + acb_srcptr pt, + acb_ode_basis_t basis, + slong prec) +{ + return acb_ode_fundamental_matrix_vec(mat, dop, Dop, expos, lcrt, pt, 1, + basis, prec); +} diff --git a/src/acb_ode/ind_lbound.c b/src/acb_ode/ind_lbound.c new file mode 100644 index 0000000000..c68b39f220 --- /dev/null +++ b/src/acb_ode/ind_lbound.c @@ -0,0 +1,160 @@ +#include "arf.h" +#include "acb.h" +#include "acb_ode.h" + +void +acb_ode_ind_lbound_init(acb_ode_ind_lbound_t ind_lbound) +{ + ind_lbound->length = 0; + ind_lbound->r = NULL; +} + +void +acb_ode_ind_lbound_clear(acb_ode_ind_lbound_t ind_lbound) +{ + for (slong i = 0; i < ind_lbound->length; i++) + { + acb_clear(ind_lbound->r[i].root); + mag_clear(ind_lbound->r[i].global_lbound); + } + flint_free(ind_lbound->r); +} + +static void +find_min(slong * n_min, mag_ptr lb, const acb_ode_group_t group, + const acb_t rt, slong mult, slong prec) +{ + acb_t tmpacb; + arb_t crit_n; + arf_t x0, x1; + mag_t tmpmag; + + acb_init(tmpacb); + arb_init(crit_n); + arf_init(x0); + arf_init(x1); + mag_init(tmpmag); + + arb_sqr(crit_n, acb_realref(rt), prec); + arb_addmul(crit_n, acb_imagref(rt), acb_imagref(rt), prec); + arb_div(crit_n, crit_n, acb_realref(rt), prec); + + if (mag_cmp_2exp_si(arb_radref(crit_n), 3) > 0) + { + * n_min = WORD_MAX; + /* Lower bound valid for all real n > 0 */ + arb_get_mag_lower(lb, acb_imagref(rt)); + acb_get_mag(tmpmag, rt); + mag_div(lb, lb, tmpmag); + } + else + { + arb_get_interval_arf(x0, x1, crit_n, prec); + slong n0 = arf_get_si(x0, ARF_RND_DOWN); + slong n1 = arf_get_si(x1, ARF_RND_UP); + + /* Extend the interval so that the endpoints are ordinary indices */ + for (slong s = 0; s < group->nshifts;) + { + slong exn = group->shifts[s].n; + if (exn == n0) + n0--, s = FLINT_MAX(0, s-1); + else if (exn == n1) + n1++, s++; + else + s++; + } + + n0 = FLINT_MAX(n0, 1); + * n_min = n1; + + mag_one(lb); + for (slong n = n0, s = 0; n <= n1; n++) + { + /* Skip exceptional indices (which would lead to lb = 0 and will be + * handled separately) */ + while (s + 1 < group->nshifts && group->shifts[s].n < n) + s++; + if (n == group->shifts[s].n) + continue; + + /* Evaluate |1-rt/n| */ + acb_div_si(tmpacb, rt, n, prec); + acb_sub_si(tmpacb, tmpacb, 1, prec); + acb_get_mag_lower(tmpmag, tmpacb); + mag_min(lb, lb, tmpmag); + } + mag_pow_ui_lower(lb, lb, mult); + } + + acb_clear(tmpacb); + arb_clear(crit_n); + arf_clear(x0); + arf_clear(x1); + mag_clear(tmpmag); +} + +void +acb_ode_ind_lbound_precompute(acb_ode_ind_lbound_t ind_lbound, + const acb_ode_exponents_t expos, + slong grp, // group index in expos + slong prec) +{ + acb_t leader, rt; + /* flint_printf("== acb_ode_ind_lbound_precompute ==\n"); */ + + acb_init(leader); + acb_init(rt); + + slong length = 0; + for (slong g = 0; g < expos->ngroups; g++) + length += expos->groups[g].nshifts; + + if (ind_lbound->r != NULL) + acb_ode_ind_lbound_clear(ind_lbound); + ind_lbound->r = flint_malloc(sizeof(*(ind_lbound->r)) * length); + + slong i = 0; + for (slong g = 0; g < expos->ngroups; g++) + { + if (g == grp) + acb_zero(leader); + else + acb_sub(leader, expos->groups[g].leader, expos->groups[grp].leader, prec); + + for (slong s = 0; s < expos->groups[g].nshifts; s++) + { + /* A root of the shifted indicial polynomial */ + acb_add_si(rt, leader, expos->groups[g].shifts[s].n, prec); + + /* flint_printf("g=%ld s=%ld rt=%{acb}\n", g, s, rt); */ + + /* When Re(α) ≤ 0, the sequence |1-α/n| decreases to 1. */ + if (arb_is_nonpositive(acb_realref(rt))) + continue; + + /* Otherwise, it first decreases to its minimum (which may be 0 + * if α is an integer), then increases to 1. We compute the + * minimum and a value of n after which the sequence is + * nondecreasing. */ + acb_init(ind_lbound->r[i].root); + mag_init(ind_lbound->r[i].global_lbound); + + acb_set(ind_lbound->r[i].root, rt); + ind_lbound->r[i].mult = expos->groups[g].shifts[s].mult; + + find_min(&ind_lbound->r[i].n_min, ind_lbound->r[i].global_lbound, + expos->groups + g, rt, ind_lbound->r[i].mult, prec); + + /* flint_printf("mult=%ld n_min=%ld lbound=%{mag}\n", + ind_lbound->r[i].mult, ind_lbound->r[i].n_min, + ind_lbound->r[i].global_lbound + i); */ + + i++; + } + } + ind_lbound->length = i; + + acb_clear(rt); + acb_clear(leader); +} diff --git a/src/acb_ode/inlines.c b/src/acb_ode/inlines.c new file mode 100644 index 0000000000..52d1062eaa --- /dev/null +++ b/src/acb_ode/inlines.c @@ -0,0 +1,3 @@ +#define ACB_ODE_INLINES_C +#include "acb_ode.h" +#undef ACB_ODE_INLINES_C diff --git a/src/acb_ode/poly_negdivrevhigh.c b/src/acb_ode/poly_negdivrevhigh.c new file mode 100644 index 0000000000..9627acd70d --- /dev/null +++ b/src/acb_ode/poly_negdivrevhigh.c @@ -0,0 +1,64 @@ +#include "acb.h" +#include "acb_ode.h" + +typedef enum +{ + LC_UNKNOWN, + LC_ACB, + LC_FMPZ +} +lc_status_t; + +// XXX maybe accept separate lengths +void +_acb_ode_poly_negdivrevhigh(acb_ptr res, acb_srcptr a, acb_srcptr cst, acb_srcptr b, + slong len, slong prec) +{ + fmpz_t lc_fmpz; + acb_t invlc; + lc_status_t lc_status = LC_UNKNOWN; + + fmpz_init(lc_fmpz); + acb_init(invlc); + + for (slong k = len - 1; k >= 0; k--) + { + /* cst*b[k] + sum(a[1+j] * res[k+1+j], j≥0) */ + if (cst == NULL) + acb_dot(res + k, b + k, 0, a + 1, 1, res + k + 1, 1, len - k - 1, prec); + else + { + acb_mul(res + k, cst, b + k, prec); + acb_dot(res + k, res + k, 0, a + 1, 1, res + k + 1, 1, len - k - 1, prec); + } + + if (acb_is_zero(res + k)) + continue; + + /* only compute invlc if needed */ + // XXX share between calls? + if (lc_status == LC_UNKNOWN) + { + if (acb_is_exact(a) && acb_get_unique_fmpz(lc_fmpz, a)) + { + lc_status = LC_FMPZ; + acb_indeterminate(invlc); + } + else + { + lc_status = LC_ACB; + acb_inv(invlc, a, prec); + } + } + + if (lc_status == LC_FMPZ) + acb_div_fmpz(res + k, res + k, lc_fmpz, prec); + else + acb_mul(res + k, res + k, invlc, prec); + + acb_neg(res + k, res + k); + } + + acb_clear(invlc); + fmpz_clear(lc_fmpz); +} diff --git a/src/acb_ode/poly_taylor_shift_aps_trunc.c b/src/acb_ode/poly_taylor_shift_aps_trunc.c new file mode 100644 index 0000000000..fe00654adf --- /dev/null +++ b/src/acb_ode/poly_taylor_shift_aps_trunc.c @@ -0,0 +1,36 @@ +#include "acb_poly.h" + +void +acb_ode_poly_taylor_shift_aps_trunc(acb_poly_t g, + const acb_poly_t f, acb_srcptr a, + slong n, slong len, slong prec) +{ + acb_poly_fit_length(g, len); + + if (len == 1 && acb_is_zero(a)) /* covers ordinary points */ + { + acb_ptr c = g->coeffs; + acb_zero(c); + + for (slong i = acb_poly_degree(f); i >= 0; i--) + { + acb_mul_si(c, c, n, prec); + acb_add(c, c, f->coeffs + i, prec); + } + + _acb_poly_set_length(g, 1); + _acb_poly_normalise(g); + } + else + { + acb_t b; + acb_init(b); + + acb_add_si(b, a, n, prec); + acb_poly_taylor_shift(g, f, b, prec); /* wasteful */ + acb_poly_truncate(g, len); + + acb_clear(b); + } +} + diff --git a/src/acb_ode/sol.c b/src/acb_ode/sol.c new file mode 100644 index 0000000000..95d707a852 --- /dev/null +++ b/src/acb_ode/sol.c @@ -0,0 +1,117 @@ +#include "acb_mat.h" +#include "acb_poly.h" +#include "acb_ode.h" + + +void +acb_ode_sol_init(acb_ode_sol_t sol, slong nshifts, + slong nlogs, slong npts, slong nder) +{ + sol->series = _acb_poly_vec_init(nlogs); + sol->sums = _acb_vec_init(npts * nlogs * nder); + acb_mat_init(sol->extini, nshifts, nlogs); + sol->tb = _mag_vec_init(nder); + + acb_ode_cvest_init(sol->cvest); + + sol->alloc_logs = nlogs; + sol->nlogs = 0; + sol->npts = npts; + sol->nder = nder; +} + + +void +acb_ode_sol_set_ini(acb_ode_sol_t sol, acb_srcptr ini, + const acb_ode_shift_t shifts) // unused +{ + slong nshifts = acb_mat_nrows(sol->extini); + slong i = 0; + for (slong j = 0; j < nshifts; j++) + { + slong mult = shifts ? shifts[j].mult : 1; + for (slong k = 0; k < mult; k++) + acb_set(acb_mat_entry(sol->extini, j, k), ini + i++); + } +} + + +void +acb_ode_sol_unit_ini(acb_ode_sol_t sol, slong i0, + const acb_ode_shift_t shifts) +{ + slong nshifts = acb_mat_nrows(sol->extini); + for (slong i = 0, j = 0; j < nshifts; j++) + { + slong mult = shifts ? shifts[j].mult : 1; + for (slong k = 0; k < mult; k++, i++) + acb_set_si(acb_mat_entry(sol->extini, j, k), i == i0 ? 1 : 0); + } +} + + +void +acb_ode_sol_clear(acb_ode_sol_t sol) +{ + acb_mat_clear(sol->extini); + _acb_vec_clear(sol->sums, sol->npts * sol->alloc_logs * sol->nder); + _acb_poly_vec_clear(sol->series, sol->alloc_logs); + _mag_vec_clear(sol->tb, sol->nder); + acb_ode_cvest_clear(sol->cvest); +} + + +void +acb_ode_sol_fit_length(acb_ode_sol_t sol, slong len) +{ + for (slong k = 0; k < sol->alloc_logs; k++) + { + acb_poly_struct * f = sol->series + k; + acb_poly_fit_length(f, len); + _acb_poly_set_length(f, f->alloc); /* for printing */ + } +} + + +void +acb_ode_sol_zero(acb_ode_sol_t sol) +{ + for (slong k = 0; k < sol->alloc_logs; k++) + acb_poly_zero(sol->series + k); + + _acb_vec_zero(sol->sums, sol->npts * sol->alloc_logs * sol->nder); + + sol->done = 0; + // XXX reset tail bounds, cvest?? +} + + + +slong +acb_ode_sum_max_nlogs(acb_ode_sum_struct * sum) +{ + slong nlogs = 0; + for (slong m = 0; m < sum->nsols; m++) + nlogs = FLINT_MAX(nlogs, sum->sol[m].nlogs); + return nlogs; +} + + +/* max degree in log at offset off in any of the solutions */ +slong +acb_ode_sum_max_nlogs_xn(acb_ode_sum_struct * sum, slong off) +{ + for (slong nlogs = acb_ode_sum_max_nlogs(sum); nlogs > 0; nlogs--) + { + for (slong m = 0; m < sum->nsols; m++) + { + if (nlogs <= sum->sol[m].nlogs) + { + acb_poly_struct * p = sum->sol[m].series + nlogs - 1; + if (off < p->length && !acb_is_zero(p->coeffs + off)) + return nlogs; + } + } + } + return 0; +} diff --git a/src/acb_ode/sol_estimate.c b/src/acb_ode/sol_estimate.c new file mode 100644 index 0000000000..f2d3ba3bde --- /dev/null +++ b/src/acb_ode/sol_estimate.c @@ -0,0 +1,56 @@ +#include "acb.h" +#include "acb_ode.h" +#include "arf.h" +#include "mag.h" + +slong +acb_ode_sol_estimate_terms(mag_t est, + const acb_ode_sol_t sol, slong off, slong len, + const mag_t radpow) +{ + mag_t m; + + mag_init(m); + + slong accuracy = ARF_PREC_EXACT; + mag_zero(est); + + for (slong k = 0; k < sol->nlogs; k++) + { + for (slong i = 0; i < len; i++) + { + acb_ptr c = (sol->series + k)->coeffs + off + i; + + acb_get_mag(m, c); + mag_add(est, est, m); + + accuracy = FLINT_MIN(accuracy, acb_rel_accuracy_bits(c)); + } + } + + mag_mul(est, est, radpow); + + mag_clear(m); + return accuracy; +} + +slong +acb_ode_sol_estimate_sums(mag_t mag, mag_t rad, const acb_ode_sol_t sol) +{ + slong len = sol->npts * sol->nlogs * sol->nder; + + _acb_vec_get_mag(mag, sol->sums, len); /* XXX good enough? */ + + mag_zero(rad); + for (slong i = 0; i < len; i++) + { + mag_max(rad, rad, &sol->sums[i].real.rad); + mag_max(rad, rad, &sol->sums[i].imag.rad); + } + + slong accuracy = ARF_PREC_EXACT; + for (slong i = 0; i < len; i++) + accuracy = FLINT_MIN(accuracy, acb_rel_accuracy_bits(sol->sums + i)); + + return accuracy; +} diff --git a/src/acb_ode/sol_jet.c b/src/acb_ode/sol_jet.c new file mode 100644 index 0000000000..c42256275b --- /dev/null +++ b/src/acb_ode/sol_jet.c @@ -0,0 +1,96 @@ +#include "acb_types.h" +#include "acb_poly.h" +#include "acb_ode.h" + +/* XXX should val be a plain acb array? acb functions that compute jets usually + * return them as acb_ptrs, not acb_polys, but here we are computing several + * jets at once, so the user cannot have these functions write into an acb_poly + * for further processing */ + +void +_acb_ode_sol_jet(acb_poly_struct * val, acb_srcptr expo, + acb_srcptr sums, mag_srcptr tb, acb_srcptr pt, + slong nlogs, slong nder, slong nfrobshifts, + slong prec) +{ + acb_poly_t expt, exptpow, log, logpow, term; + + if (nlogs <= 0) + { + _acb_poly_vec_zero(val, nfrobshifts); + return; + } + + acb_poly_init(expt); + acb_poly_init(exptpow); + acb_poly_init(log); + acb_poly_init(logpow); + acb_poly_init(term); + + acb_poly_fit_length(term, nder); + + acb_poly_set_coeff_si(expt, 1, 1); + acb_poly_set_coeff_acb(expt, 0, pt); + + /* XXX does this work when pt contains 0 but expo is a nonneg int? */ + acb_poly_pow_acb_series(exptpow, expt, expo, nder, prec); + acb_poly_log_series(log, expt, nder, prec); + + acb_poly_one(logpow); + + slong len = FLINT_MIN(nlogs, nfrobshifts); + for (slong d = 0; d < len; d++) + { + acb_poly_fit_length(val + d, nder); + _acb_vec_set((val + d)->coeffs, sums + d * nder, nder); + for (slong i = 0; i < nder; i++) + acb_add_error_mag((val + d)->coeffs + i, tb + i); + _acb_poly_set_length(val + d, nder); + _acb_poly_normalise(val + d); + } + _acb_poly_vec_zero(val + len, nfrobshifts - len); + + for (slong p = 1; p < nlogs; p++) + { + acb_poly_mullow(logpow, logpow, log, nder, prec); + _acb_vec_scalar_div_ui(logpow->coeffs, logpow->coeffs, logpow->length, + p, prec); + + for (slong d = 0; d < FLINT_MIN(nlogs - p, nfrobshifts); d++) + { + _acb_poly_mullow(term->coeffs, + sums + (d + p) * nder, nder, + logpow->coeffs, nder, + nder, prec); + _acb_poly_set_length(term, nder); + _acb_poly_normalise(term); + acb_poly_add(val + d, val + d, term, prec); + } + } + + for (slong d = 0; d < nfrobshifts; d++) + acb_poly_mullow(val + d, val + d, exptpow, nder, prec); + + acb_poly_clear(expt); + acb_poly_clear(exptpow); + acb_poly_clear(log); + acb_poly_clear(logpow); + acb_poly_clear(term); +} + + +void +acb_ode_sol_jet(acb_poly_struct * val, + acb_srcptr expo, const acb_ode_sol_t sol, slong j, + acb_srcptr pt, slong nder, slong nfrobshifts, slong prec) +{ + FLINT_ASSERT(nder <= sol->nder); + + if (nder == 0) + return; + + _acb_ode_sol_jet(val, expo, + acb_ode_sol_sum_ptr(sol, j, 0, 0), sol->tb, pt, + sol->nlogs, nder, nfrobshifts, + prec); +} diff --git a/src/acb_ode/solution_growth.c b/src/acb_ode/solution_growth.c new file mode 100644 index 0000000000..3225fbfe6f --- /dev/null +++ b/src/acb_ode/solution_growth.c @@ -0,0 +1,89 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_poly.h" + +/* max growth of + * - coefficient sequence ≈ base^n / n!^(1/order) + * - function at infinity (if entire) ≈ exp((1/order)*(base*x)^order) */ + +void +_acb_ode_solution_growth(mag_t order, mag_t base, + const acb_poly_struct * dop, slong dop_len) +{ + FLINT_ASSERT(acb_poly_length(dop + dop_len - 1) == 1); + + /* Newton polygon: If the general term of f(x) is ≈ base^n/n!^(1/order)·x^n, + applying x^j·(x·d/dx)^i to f(x) multiplies the general term by + ≈ base^(-j)·n^(i+j/order). So the possible values for order are those + such that i+j/order is reached twice or more, and we want the largest of + these. */ + + slong i0 = -1; + slong j0 = -1; + for (slong i = 0; i < dop_len; i++) + for (slong j = 0; j < dop[i].length; j++) + if (!acb_is_zero(dop[i].coeffs + j)) + if (i > i0 || i == i0 && j > j0) + i0 = i, j0 = j; + + slong p = -1, q = 0; + for (slong i = 0; i < dop_len; i++) + for (slong j = 0; j < dop[i].length; j++) + { + if (j > j0 && i < i0 && !acb_is_zero(dop[i].coeffs + j)) + { + /* (j, i) belongs to the Newton polygon */ + // flint_printf("(%ld, %ld) ", j, i); + slong p1 = i - i0, q1 = j - j0; + if (p1 * q > p * q1) + p = p1, q = q1; + } + } + // flint_printf("\n"); + + FLINT_ASSERT(p <= 0); + FLINT_ASSERT(q >= 0); + + mag_set_ui(order, q); + if (q == 0) + { + mag_zero(base); + return; + } + mag_div_ui(order, order, -p); + + /* estimate the roots of the corresponding characteristic equation */ + + acb_poly_t eqn; + mag_t mag; + acb_poly_init(eqn); + mag_init(mag); + + for (slong i = 0; i <= i0; i++) + for (slong j = j0; j < dop[i].length; j++) + if ((i - i0) * q == p * (j - j0)) + acb_poly_set_coeff_acb(eqn, i0 - i, dop[i].coeffs + j); + + /* change to a cheap but rigorous lower bound? (not needed at the moment but + * may be useful for other applications) */ + acb_ptr rts = _acb_vec_init(eqn->length - 1); + _acb_poly_find_roots_double(rts, eqn->coeffs, NULL, eqn->length, 4, + MAG_BITS); + + // flint_printf("eqn=%{acb_poly} rts=%{acb*}\n", eqn, rts, eqn->length - 1); + + mag_inf(base); + for (slong i = 0; i < eqn->length - 1; i++) + { + acb_get_mag_lower(mag, rts + i); + mag_min(base, base, mag); + } + + mag_inv(base, base); + mag_root(base, base, q); + mag_pow_ui(base, base, -p); + + _acb_vec_clear(rts, eqn->length - 1); + mag_clear(mag); + acb_poly_clear(eqn); +} diff --git a/src/acb_ode/stairs.c b/src/acb_ode/stairs.c new file mode 100644 index 0000000000..9e87171e02 --- /dev/null +++ b/src/acb_ode/stairs.c @@ -0,0 +1,66 @@ +#include "acb_ode.h" +#include "mag.h" + + +void +acb_ode_stairs_init(acb_ode_stairs_t stairs) +{ + stairs->length = 0; + stairs->h = NULL; +} + +void +acb_ode_stairs_clear(acb_ode_stairs_t stairs) +{ + _mag_vec_clear(stairs->h, stairs->length); +} + +/* [M19], Algorithm 7.4, step 1 */ + +void +acb_ode_stairs_precompute(acb_ode_stairs_t stairs, + const acb_poly_struct * num, slong len, + const acb_poly_t ind, const acb_ode_group_t group, + const acb_ode_ind_lbound_t ind_lbound, + slong prec) +{ + /* flint_printf("== _acb_ode_stairs_precompute %{acb} ==\n", group->leader); */ + + if (stairs->h != NULL) + acb_ode_stairs_clear(stairs); + stairs->length = len /* +1 ? */ * group->nshifts; + stairs->h = _mag_vec_init(stairs->length); + + mag_ptr bound = _mag_vec_init(2*len); + mag_ptr bound1 = bound + len; + + /* This could in principle be refined when it is known that the degree in + * log(x) of the solutions will not reach the bound. */ + slong ord = 0; + for (slong s = 0; s < group->nshifts; s++) + ord += group->shifts[s].mult; + + for (slong s = group->nshifts - 1; s >= 0; s--) { + slong n = group->shifts[s].n; + slong mult = acb_ode_group_multiplicity(group, n); + /* flint_printf("s=%ld n=%ld mult=%ld ord=%ld\n", s, n, mult, ord); */ + acb_ode_bound_rat_ref_vec(bound, num, len, ind, n, mult, ord, prec); + /* flint_printf(" n=%ld ref = %{mag*}\n", n, bound + i, len); */ + if (s == group->nshifts - 1 || group->shifts[s + 1].n != n + 1) { + acb_ode_bound_rat_ordinary_vec(bound1, num, len, ind, ind_lbound, + n + 1, ord, prec); + /* flint_printf(" n=%ld bound1 = %{mag*}\n", n, bound1 + i, len); */ + for (slong i = 0; i < len; i++) + mag_max(bound + i, bound + i, bound1 + i); + } + for (slong i = 0; i < len; i++) { + slong j = i * group->nshifts + s; + mag_swap(stairs->h + j, bound + i); + if (s + 1 < group->nshifts) + mag_max(stairs->h + j, stairs->h + j, stairs->h + j + 1); + } + ord -= group->shifts[s].mult; + } + + _mag_vec_clear(bound, 2*len); +} diff --git a/src/acb_ode/sum_context.c b/src/acb_ode/sum_context.c new file mode 100644 index 0000000000..89ceb88eb9 --- /dev/null +++ b/src/acb_ode/sum_context.c @@ -0,0 +1,170 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_mat.h" +#include "acb_poly.h" +#include "acb_ode.h" +#include "arb_poly.h" +#include "arf.h" +#include "fmpz_vec.h" + + +void +acb_ode_sum_init(acb_ode_sum_t sum, slong dop_len, slong npts, + slong nsols, slong nder) +{ + slong dop_order = dop_len - 1; + + sum->dop = _acb_poly_vec_init(dop_len); + sum->dop_len = dop_len; + + acb_ode_group_init(sum->group, dop_order); + for (slong s = 0; s < dop_order; s++) + sum->group->shifts[s].n = -1; + acb_poly_init(sum->ind); + + sum->sol = flint_malloc(nsols * sizeof(acb_ode_sol_struct)); + /* using dop_order as a bound + * - for max possible log prec (will need updating to support inhomogeneous + * equations), + * - for number of initial value positions */ + for (slong m = 0; m < nsols; m++) + acb_ode_sol_init(sum->sol + m, dop_order, dop_order, npts, nder); + sum->nsols = nsols; + + sum->pts = _acb_vec_init(npts); + sum->pows = _acb_vec_init(npts * nder); + sum->shifted_sums = flint_malloc(npts); + sum->npts = npts; + sum->nder = nder; + + mag_init(sum->mag); + mag_init(sum->magpow); + mag_init(sum->cvrad); + arb_poly_init(sum->itg_pol); + arb_poly_init(sum->itg_num); + + sum->binom_n = _fmpz_vec_init(nder); + + sum->flags = 0; + sum->have_precomputed = 0; +} + +void +acb_ode_sum_clear(acb_ode_sum_t sum) +{ + _acb_poly_vec_clear(sum->dop, sum->dop_len); + + acb_ode_group_clear(sum->group); + acb_poly_clear(sum->ind); + + for (slong m = 0; m < sum->nsols; m++) + acb_ode_sol_clear(sum->sol + m); + flint_free(sum->sol); + + _acb_vec_clear(sum->pts, sum->npts); + _acb_vec_clear(sum->pows, sum->npts * sum->nder); + flint_free(sum->shifted_sums); + + mag_clear(sum->mag); + mag_clear(sum->magpow); + mag_clear(sum->cvrad); + arb_poly_clear(sum->itg_pol); + arb_poly_clear(sum->itg_num); + + _fmpz_vec_clear(sum->binom_n, sum->nder); +} + +/* Operator */ + +void +acb_ode_sum_set_diffop(acb_ode_sum_t sum, const acb_poly_t dop, slong dop_len, + const mag_t cvrad) +{ + FLINT_ASSERT(dop_len > 0); + _acb_poly_vec_set(sum->dop, dop, dop_len); + mag_set(sum->cvrad, cvrad); +} + +/* Solution group */ + +void +acb_ode_sum_set_ordinary(acb_ode_sum_t sum) +{ + acb_zero(sum->group->leader); + for (slong n = 0; n < sum->dop_len - 1; n++) + { + sum->group->shifts[n].n = n; + sum->group->shifts[n].mult = 1; + } +} + +void +acb_ode_sum_set_group(acb_ode_sum_t sum, const acb_ode_group_t group) +{ + acb_ode_group_set(sum->group, group); +} + +/* Initial values */ + +void +acb_ode_sum_set_ini_echelon(acb_ode_sum_t sum) +{ + for (int m = 0; m < sum->nsols; m++) + acb_ode_sol_unit_ini(sum->sol + m, m, sum->group->shifts); +} + + +void +acb_ode_sum_set_ini_highest(acb_ode_sum_t sum) +{ + for (int m = 0; m < sum->nsols; m++) + { + acb_mat_zero(sum->sol[m].extini); + acb_one(acb_mat_entry(sum->sol[m].extini, m, + sum->group->shifts[m].mult - 1)); + } +} + +/* Evaluation points */ + +void +acb_ode_sum_set_points(acb_ode_sum_t sum, acb_srcptr pts, slong npts) +{ + _acb_vec_set(sum->pts, pts, npts); +} + +void +acb_ode_sum_precompute(acb_ode_sum_t sum) +{ + if (sum->have_precomputed) + return; + + slong dop_len = sum->dop_len; + + arf_t abspt; + arf_init(abspt); + + sum->dop_clen = _acb_poly_vec_length(sum->dop, dop_len); + + /* XXX _gr_ore_poly_indicial_polynomial? */ + acb_poly_fit_length(sum->ind, dop_len); + _acb_poly_set_length(sum->ind, dop_len); + for (slong i = 0; i < dop_len; i++) + acb_poly_get_coeff_acb((sum->ind)->coeffs + i, sum->dop + i, 0); + _acb_poly_normalise(sum->ind); + + for (slong j = 0; j < sum->npts; j++) /* XXX move? */ + { + /* We can trade some full-width muls for muls by small integers by + multiplying everyone by ξ^n, but this optimization is only + valid when ξ does not contain 0. */ + sum->shifted_sums[j] = !acb_contains_zero(sum->pts + j); + } + + _acb_vec_get_mag(sum->mag, sum->pts, sum->npts); + FLINT_ASSERT(mag_cmp(sum->mag, sum->cvrad) <= 0); + + sum->have_precomputed = 1; + + arf_clear(abspt); +} diff --git a/src/acb_ode/sum_divconquer.c b/src/acb_ode/sum_divconquer.c new file mode 100644 index 0000000000..ee70a6a7df --- /dev/null +++ b/src/acb_ode/sum_divconquer.c @@ -0,0 +1,46 @@ +#include "acb.h" +#include "acb_ode.h" + +static void +acb_ode_sum_reset(acb_ode_sum_t sum) +{ + sum->n = 0; + sum->n0 = 0; + + for (slong m = 0; m < sum->nsols; m++) + { + acb_ode_sol_zero(sum->sol + m); + + sum->sol[m].cvest->accuracy = sum->wp; + sum->sol[m].cvest->loss_rate = sum->dop_len * sum->dop_clen; + } + + for (slong i = 0; i < sum->npts; i++) + acb_one(sum->pows + i); + + mag_one(sum->magpow); +} + +/* XXX make it possible to customize the working precision--how exactly? */ +void +acb_ode_sum_divconquer(acb_ode_sum_t sum, slong nterms, slong prec) +{ + if (nterms < 0) + nterms = WORD_MAX; + + acb_ode_sum_precompute(sum); + acb_ode_sum_reset(sum); + + /* todo: take into account nterms and the APPROX flag */ + sum->sum_wp = acb_ode_choose_prec(&sum->wp, sum->dop, sum->dop_len, + sum->mag, sum->cvrad, prec); + + slong block_len = FLINT_MAX(1, sum->dop_clen - 1); + slong stride = block_len; /* todo: less when block_len is large */ + + for (slong m = 0; m < sum->nsols; m++) + acb_ode_sol_fit_length(sum->sol + m, 2 * block_len); + + _acb_ode_sum_forward_divconquer(sum, nterms, block_len, stride, prec); + _acb_ode_sum_fix(sum); +} diff --git a/src/acb_ode/sum_done.c b/src/acb_ode/sum_done.c new file mode 100644 index 0000000000..12fc06ae3b --- /dev/null +++ b/src/acb_ode/sum_done.c @@ -0,0 +1,295 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_poly.h" +#include "acb_ode.h" +#include "arb_poly.h" + +void +acb_ode_cvest_println(const acb_ode_cvest_t cvest) +{ + flint_printf("neglogterm=%f cvg_rate=%f terms_wanted=%f accuracy=%ld " + "loss_rate=%f terms_full_prec=%f prec_wanted=%ld\n", + cvest->neglogterm, cvest->cvg_rate, cvest->terms_wanted, + cvest->accuracy, cvest->loss_rate, cvest->terms_full_prec, + cvest->prec_wanted); +} + + +/* TODO support the case when only part of the residual is known */ +static void +bound_normalized_residuals(arb_poly_struct * nres_maj, + acb_ode_sum_struct * sum, + slong prec) +{ + acb_t cst; + arf_t tmparf; + acb_poly_t ind_nd; + TMP_INIT; + TMP_START; + + acb_init(cst); + arf_init(tmparf); + + acb_poly_struct * dop_lc = sum->dop + sum->dop_len - 1; + + acb_poly_init(ind_nd); + slong nlogs0 = acb_ode_sum_max_nlogs(sum); + acb_ptr rhs_term = TMP_ALLOC(sizeof(acb_struct) * nlogs0); + acb_ptr nres_term = _acb_vec_init(nlogs0); + + for (slong m = 0; m < sum->nsols; m++) + { + arb_poly_fit_length(nres_maj + m, sum->dop_clen); + _arb_poly_set_length(nres_maj + m, sum->dop_clen); + } + + for (slong d = 0, s = 0; d < sum->dop_clen - 1; d++) + { + slong n = sum->n; + + /* Multiplicity of λ+n+d as an indicial root */ + slong mult = 0; + for (; s < sum->dop_len - 1; s++) { + if (sum->group->shifts[s].n == n) + { + mult = sum->group->shifts[s].mult; + break; + } + } + + /* ind_nd/cst = jet of monic indicial polynomial at λ+n+d */ + // XXX change acb_ode_sum_max_nlogs to take k (=n+d here) or k - n instead of k - n0 + slong nlogs_d = acb_ode_sum_max_nlogs_xn(sum, n - sum->n0 + d); + + acb_ode_poly_taylor_shift_aps_trunc(ind_nd, sum->ind, sum->group->leader, + n, nlogs_d + mult, prec); + acb_set_round(cst, dop_lc->coeffs, prec); + acb_neg(cst, cst); + + acb_poly_fit_length(ind_nd, nlogs_d + mult); + for (slong m = 0; m < sum->nsols; m++) + { + // XXX skip solutions that have converged? + + // XXX tighten? + slong nlogs_dm = FLINT_MIN(sum->sol[m].nlogs, nlogs_d); + + /* nres_term = coeff of x^{λ+n+d} in normalized residual of sol m */ + for (slong k = 0; k < nlogs_dm; k++) + rhs_term[k] = (sum->sol[m].series + k)->coeffs[n + d - sum->n0]; + _acb_ode_poly_negdivrevhigh(nres_term, ind_nd->coeffs + mult, + cst, rhs_term, nlogs_dm, prec); + + /* nres_maj[m] = common majorant of the coefficients of log(x)^k/k! + * in the normalized residual of solution m */ + for (slong k = 0; k < nlogs_dm; k++) + { + arf_ptr c = arb_midref((nres_maj + m)->coeffs + d); + acb_get_abs_ubound_arf(tmparf, nres_term + k, prec); + arf_max(c, c, tmparf); + } + } + } + + for (slong m = 0; m < sum->nsols; m++) + _arb_poly_normalise(nres_maj + m); + + _acb_vec_clear(nres_term, nlogs0); + acb_poly_clear(ind_nd); + acb_clear(cst); + arf_clear(tmparf); + TMP_END; +} + + +static void +update_tail_bounds(acb_ode_sum_struct * sum, slong prec) +{ + arb_t rad; + arb_poly_t jet_tb; + TMP_INIT; + TMP_START; + + arb_init(rad); + arb_poly_init(jet_tb); + arb_poly_struct * nres_maj = TMP_ALLOC(sizeof(arb_poly_t) * sum->nsols); + for (slong m = 0; m < sum->nsols; m++) + arb_poly_init(nres_maj + m); + + /* ore_algebra uses a bound on nlogs: + slong nlogs = acb_ode_group_nlogs(sum->group, sum->n); + */ + + // TODO Only call precompute_integrand if the current bounds are not tight + // enough or nlogs has changed. (We can also use the larger nlogs from the + // start, but this leads to worse bounds before nlogs increases.) + acb_ode_bound_precompute_integrand(sum->itg_pol, sum->itg_num, sum->group, + sum->bound, sum->n, acb_ode_sum_max_nlogs(sum), MAG_BITS); + + bound_normalized_residuals(nres_maj, sum, prec); + + /* ore_algebra uses a single nres_maj for all solutions: + + for (slong m = 1; m < sum->nsols; m++) + { + arb_poly_fit_length(nres_maj, nres_maj[m].length); + for (slong i = 0; i < nres_maj[m].length; i++) + arb_max(nres_maj->coeffs + i, nres_maj->coeffs + i, + nres_maj[m].coeffs + i, prec); + nres_maj->length = FLINT_MAX(nres_maj->length, nres_maj[m].length); + } + _arb_poly_normalise(nres_maj); + */ + + for (slong m = 0; m < sum->nsols; m++) + { + slong i; + arf_set_mag(arb_midref(rad), sum->mag); + acb_ode_tail_bound_jet_precomp(jet_tb, + sum->bound, sum->n, sum->itg_pol, sum->itg_num, + nres_maj + m, rad, sum->nder, MAG_BITS); + for (i = 0; i < jet_tb->length; i++) + arb_get_mag((sum->sol + m)->tb + i, jet_tb->coeffs + i); + for (; i < sum->nder; i++) + mag_zero((sum->sol + m)->tb + i); + } + + arb_poly_clear(jet_tb); + arb_clear(rad); + for (slong m = 0; m < sum->nsols; m++) + arb_poly_clear(nres_maj + m); + TMP_END; +} + +/* stride = terms since last check */ +int +acb_ode_sum_done(acb_ode_sum_struct * sum, slong stride, slong prec) +{ + mag_t eps, est, sum_mag, sum_rad; + mag_init(eps); + mag_init(est); + mag_init(sum_mag); + mag_init(sum_rad); + + int stop = 1; + int have_tail_bounds = 0; + + slong deg = sum->dop_clen - 1; + FLINT_ASSERT(sum->n >= sum->n0 + deg); /* todo: raise this limitation */ + + /* detect solutions that have converged or will not improve */ + + for (slong m = 0; m < sum->nsols; m++) + { + acb_ode_sol_struct * sol = sum->sol + m; + + /* todo: non-rigorous check with less than deg terms? + (maybe using known terms from the previous residuals in addition to + the available coefficients) */ + slong off = sum->n - deg - sum->n0; + slong accuracy = acb_ode_sol_estimate_terms(est, sol, off, deg, + sum->magpow); + slong sum_accuracy = acb_ode_sol_estimate_sums(sum_mag, sum_rad, sol); + acb_ode_cvest_update(sol->cvest, sol->cvest, + est, accuracy, stride, prec, sum->wp); + + /* Attempt to detect cases where it is better to bail out even though, + technically, we may still be able to improve the results (typically + entire functions with a large hump to climb). XXX Not sure if this + belongs here; maybe find a better criterion. */ + if (mag_cmp_2exp_si(sum_rad, 4 * prec) >= 0 + && (sum_accuracy < 0 || sum->n > prec * FLINT_BIT_COUNT(prec))) + { + mag_inf(eps); + } + /* sum_accuracy may be <= 0 and still have a chance to improve, + typically when the terms are still increasing. + + If accuracy <= 0, the approximation of the sum is no longer going to + get better, but the error bound still might (e.g., if the evaluation + point contains zero), roughly until the radius of the partial sum + exceeds the tail bound. */ + else if (accuracy >= 0 || sum_accuracy >= 0) + { + /* stopping condition: ~prec correct bits OR working precision + insufficient to do better */ + mag_mul_2exp_si(eps, sum_mag, -prec); + mag_max(eps, eps, sum_rad); + mag_mul_2exp_si(eps, eps, -8); + } + else + { + /* todo: The point about the bound having a chance to improve still + holds. A better heuristic (maybe based on the rate at which the + tail bounds are improving?) should improve evaluations on wide + intervals (typically containing 0). */ + mag_inf(eps); + } + + // flint_printf("n=%ld ", sum->n); acb_ode_cvest_println(sol->cvest); + // if (sum->n % (10 * stride) == 0) + // flint_printf("n=%ld prec=%ld m=%ld est=%{mag} accuracy=%ld sum_mag=%{mag} sum_rad=%{mag} sum_accuracy=%ld eps=%{mag}\n", sum->n, prec, m, est, accuracy, sum_mag, sum_rad, sum_accuracy, eps); + + /* quick non-rigorous check */ + if (mag_cmp(est, eps) > 0 && accuracy >= 0) + { + stop = 0; + // flint_printf("cont\n"); + continue; + } + + if (sum->flags & ACB_ODE_APPROX) + { + sol->done = 1; + // flint_printf("done (approx)\n"); + continue; + } + + /* rigorous check */ + if (!have_tail_bounds) + { + // XXX wasteful when only some sols have converged + update_tail_bounds(sum, MAG_BITS); + have_tail_bounds = 1; + } + + slong max_shift = sum->group->shifts[sum->group->nshifts - 1].n; + + // XXX also check derivatives? + if (sum->nder > 0 && mag_cmp(sol->tb + 0, eps) <= 0) + { + // flint_printf("done\n"); + sol->done = 1; // XXX currently unused + } + else if (sum_accuracy == -ARF_PREC_EXACT) + { + sol->done = 1; // XXX currently unused + } + /* Crude heuristic to detect situations where the tail bound might as + well be infinite. XXX Should somehow take into account the rate at + which the bounds improve, and maybe also indicial roots from other + groups. */ + else if (sum->n > max_shift + sum->dop_clen + && mag_get_d_log2_approx(sol->tb + 0) + > mag_get_d_log2_approx(eps) + prec) + { + sol->done = 1; // XXX currently unused + } + else + { + // flint_printf("cont\n"); + stop = 0; + } + + if (stop || sum->n % (10 * stride) == 0) + // flint_printf("n=%ld m=%ld tb=%{mag*}\n", sum->n, m, sum->sol[m].tb, sum->nder); + flint_printf("n=%ld m=%ld tb=%{mag}\n", sum->n, m, sum->sol[m].tb); + } + + mag_clear(eps); + mag_clear(est); + mag_clear(sum_mag); + mag_clear(sum_rad); + + return stop; +} diff --git a/src/acb_ode/sum_fix.c b/src/acb_ode/sum_fix.c new file mode 100644 index 0000000000..040c69ffef --- /dev/null +++ b/src/acb_ode/sum_fix.c @@ -0,0 +1,41 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_ode.h" + + +void +_acb_ode_sum_fix(acb_ode_sum_struct * sum) +{ + acb_t inv, invpow; + acb_init(inv); + acb_init(invpow); + + for (slong j = 0; j < sum->npts; j++) + { + if (!sum->shifted_sums[j]) + continue; + + acb_inv(inv, sum->pts + j, sum->sum_wp); + + acb_one(invpow); + for (slong i = 1; i < sum->nder; i++) + { + acb_mul(invpow, invpow, inv, sum->sum_wp); + for (slong m = 0; m < sum->nsols; m++) + { + for (slong k = 0; k < sum->sol[m].nlogs; k++) + { + acb_ptr c = acb_ode_sol_sum_ptr(sum->sol + m, j, k, i); + acb_mul(c, c, invpow, sum->sum_wp); + } + } + } + + sum->shifted_sums[j] = 0; + } + + /* XXX in WANT_SERIES mode, set series lengths as required */ + + acb_clear(invpow); + acb_clear(inv); +} diff --git a/src/acb_ode/sum_forward_1.c b/src/acb_ode/sum_forward_1.c new file mode 100644 index 0000000000..a726d4436f --- /dev/null +++ b/src/acb_ode/sum_forward_1.c @@ -0,0 +1,206 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_mat.h" +#include "acb_poly.h" +#include "acb_ode.h" +#include "fmpz.h" +#include "fmpz_vec.h" + +static void +next_binom_n(fmpz * binom_n, slong n, slong len) +{ + if (len == 0) + return; + + if (n == 0) + { + _fmpz_vec_zero(binom_n, len); + fmpz_one(binom_n); + } + else + { + for (slong i = len - 1; i > 0; i--) + fmpz_add(binom_n + i, binom_n + i, binom_n + i - 1); + } +} + +static void +next_coeff(acb_ode_sol_t sol, slong base, slong n, + acb_poly_struct * rhs, slong rhs_nlogs, + slong mult, slong ini_i, const acb_poly_t ind_n, + int approx, slong prec) +{ + slong k; + TMP_INIT; + TMP_START; + + while (rhs_nlogs > 0 + && acb_is_zero((rhs + rhs_nlogs - 1)->coeffs + n - base)) + rhs_nlogs--; + + /* Compute the high-log part of the new term from the coefficient of x^{λ+n} + * on the rhs and the indicial polynomial */ + + acb_ptr new_term = TMP_ALLOC(sizeof(acb_t) * rhs_nlogs); + acb_ptr rhs_term = TMP_ALLOC(sizeof(acb_t) * rhs_nlogs); + + for (slong k = 0; k < rhs_nlogs; k++) { + acb_init(new_term + k); + rhs_term[k] = (rhs + k)->coeffs[n - base]; + } + + _acb_ode_poly_negdivrevhigh(new_term, ind_n->coeffs + mult, NULL, rhs_term, + rhs_nlogs, prec); + + /* todo: also do this in rigorous mode but keep the radii for linear error + analysis */ + if (approx) + { + for (slong k = 0; k < rhs_nlogs; k++) + { + mag_zero(arb_radref(acb_realref(new_term + k))); + mag_zero(arb_radref(acb_imagref(new_term + k))); + } + } + + /* Write the new term to sol, store new extini */ + + for (k = 0; k < mult; k++) /* initial conditions */ + { + acb_set((sol->series + k)->coeffs + n - base, + acb_mat_entry(sol->extini, ini_i, k)); + } + for (; k < rhs_nlogs + mult; k++) + { + acb_swap((sol->series + k)->coeffs + n - base, new_term + k - mult); + + if (mult > 0) + acb_set(acb_mat_entry(sol->extini, ini_i, k), + (sol->series + k)->coeffs + n - base); + } + for (; k < sol->nlogs; k++) + { + acb_zero((sol->series + k)->coeffs + n - base); + } + + /* Update log-degree and used initial values */ + + sol->nlogs = FLINT_MAX(sol->nlogs, rhs_nlogs + mult); + + for (slong k = 0; k < rhs_nlogs; k++) + acb_clear(new_term + k); + TMP_END; +} + +static void +next_sums(acb_ode_sol_t sol, slong base, slong n, + acb_srcptr pows, const char * shifted, slong npts, + const fmpz * binom_n, + slong nder, slong prec) +{ + acb_t tmp; + acb_init(tmp); + + for (slong j = 0; j < npts; j++) + { + for (slong k = 0; k < sol->nlogs; k++) + { + acb_ptr c = (sol->series + k)->coeffs + n - base; + acb_ptr s = acb_ode_sol_sum_ptr(sol, j, k, 0); + + if (shifted[j]) + { + slong pows_offset = (n % nder) * npts; + acb_srcptr ptpow = pows + pows_offset + j; + acb_mul(tmp, c, ptpow, prec); + for (slong i = 0; i < nder; i++) + acb_addmul_fmpz(s + i, tmp, binom_n + i, prec); + } + else + { + for (slong i = 0; i < nder; i++) + { + slong pows_offset = ((n + nder - i) % nder) * npts; + acb_srcptr ptpow = pows + pows_offset + j; + acb_mul(tmp, c, ptpow, prec); + acb_addmul_fmpz(s + i, tmp, binom_n + i, prec); + } + } + } + } + + acb_clear(tmp); +} + +void +_acb_ode_sum_forward_1(acb_ode_sum_struct * sum) +{ + slong ini_i, mult, len; + acb_poly_t ind_n; + + acb_poly_init(ind_n); + + slong n = sum->n; + + /* find the position and multiplicity of the current shift as an indicial + root */ + + ini_i = -1; + mult = 0; + for (slong i = 0; i < sum->dop_len - 1; i++) + { + /* flint_printf("i=%wd shift=%wd mult=%wd\n", i, sum->group->shifts[i].n, sum->group->shifts[i].mult); */ + if (sum->group->shifts[i].n == n) + { + ini_i = i; + mult = sum->group->shifts[i].mult; + break; + } + } + + /* find the maximum log(x)-length of the corresponding coefficient of x on + the rhs */ + + len = acb_ode_sum_max_nlogs_xn(sum, n - sum->n0); + + /* evaluate the indicial polynomial */ + + acb_ode_poly_taylor_shift_aps_trunc(ind_n, sum->ind, sum->group->leader, n, + len + mult, sum->wp); + + /* flint_printf("fwd1: n=%wd mult=%wd len=%wd ind_n=%{acb_poly}\n", n, mult, len, ind_n); */ + + /* advance solutions */ + + next_binom_n(sum->binom_n, n, sum->nder); + + for (slong m = 0; m < sum->nsols; m++) + { + acb_ode_sol_struct * sol = sum->sol + m; + /* XXX when rewriting with a better data structure, pass n0 - n instead + of (n0, n) and an rhs pointer/offset separate from the sum + pointer/offset */ + next_coeff(sol, sum->n0, n, sol->series, sol->nlogs, mult, ini_i, + ind_n, sum->flags & ACB_ODE_APPROX, sum->wp); + next_sums(sol, sum->n0, n, sum->pows, sum->shifted_sums, sum->npts, + sum->binom_n, sum->nder, sum->sum_wp); + } + + /* compute the next power of each evaluation point */ + /* XXX wasteful near the end */ + + for (slong j = 0; j < sum->npts; j++) + { + slong prec = FLINT_MIN(sum->wp, sum->sum_wp); + acb_mul(sum->pows + ((n + 1) % sum->nder) * sum->npts + j, + sum->pows + (n % sum->nder) * sum->npts + j, + sum->pts + j, + prec); + } + + mag_mul(sum->magpow, sum->magpow, sum->mag); + + sum->n++; + + acb_poly_clear(ind_n); +} diff --git a/src/acb_ode/sum_forward_divconquer.c b/src/acb_ode/sum_forward_divconquer.c new file mode 100644 index 0000000000..d60bea31ce --- /dev/null +++ b/src/acb_ode/sum_forward_divconquer.c @@ -0,0 +1,169 @@ +#include "acb_types.h" +#include "acb.h" +#include "acb_poly.h" +#include "acb_ode.h" + + +/* todo: document how solutions are represented */ + +/* XXX better decouple rhs from sol */ + + +/* reads series[low-n0:mid-n0], writes series[mid-n0:high-n0] */ +static void +apply_diffop(acb_ode_sum_struct * sum, slong low, slong mid, slong high) +{ + /* flint_printf("apply_diffop: n0=%wd low=%wd mid=%wd high=%wd\n"); */ + + if (high - mid >= 7 * sum->dop_len) /* threshold from ore_algebra */ + { + for (slong m = 0; m < sum->nsols; m++) + _acb_ode_apply_diffop_polmul( + sum->sol[m].series, mid - sum->n0, + sum->dop, sum->dop_len, + sum->group->leader, low, + sum->sol[m].series, low - sum->n0, mid - low, + sum->sol[m].nlogs, + mid - low, high - mid, + sum->wp); + } + else + { + slong nlogs = acb_ode_sum_max_nlogs(sum); + + slong weights_len = (high - mid) * nlogs * (mid - low); + acb_ptr weights = _acb_vec_init(weights_len); + + _acb_ode_apply_diffop_basecase_weights( + weights, sum->dop, sum->dop_len, sum->group->leader, + low, mid - low, nlogs, mid - low, high - mid, + sum->wp); + + for (slong m = 0; m < sum->nsols; m++) + _acb_ode_apply_diffop_basecase_precomp( + sum->sol[m].series, mid - sum->n0, + weights, nlogs, sum->sol[m].series, low - sum->n0, mid - low, + sum->sol[m].nlogs, mid - low, high - mid, sum->wp); + + _acb_vec_clear(weights, weights_len); + } + + /* todo: also do this in rigorous mode but keep the radii for linear error + analysis */ + if (sum->flags & ACB_ODE_APPROX) + { + for (slong m = 0; m < sum->nsols; m++) + { + for (slong n = mid; n < high; n++) + { + acb_ptr c = sum->sol[m].series->coeffs + (n - sum->n0); + mag_zero(arb_radref(acb_realref(c))); + mag_zero(arb_radref(acb_imagref(c))); + } + } + } +} + + +static void +discard_block(acb_ode_sum_struct * sum, slong block_len) +{ + for (slong m = 0; m < sum->nsols; m++) + { + for (slong k = 0; k < sum->sol[m].nlogs; k++) + { + acb_poly_struct * f = sum->sol[m].series + k; + /* the sage version uses slightly different sizes near max_terms */ + _acb_poly_shift_right(f->coeffs, f->coeffs, 2 * block_len, block_len); + _acb_vec_zero(f->coeffs + block_len, block_len); /* XXX useful? */ + } + } + sum->n0 += block_len; +} + + +// assumes all sol have enough space for 2*block_len coefficients +int +_acb_ode_sum_forward_divconquer(acb_ode_sum_struct * sum, slong high, + slong block_len, slong stride, slong prec) +{ + slong low = sum->n; + FLINT_ASSERT(low <= high); + + if (high == low) + return 0; + else if (high == low + 1) { + _acb_ode_sum_forward_1(sum); + return 0; + } + + slong mid = high <= low + block_len ? (high + low)/2 : low + block_len; + + int done = _acb_ode_sum_forward_divconquer(sum, mid, block_len, stride, prec); + + /* + acb_poly_t tmp; + acb_poly_init(tmp); + flint_printf("divconquer: block_len=%wd n0=%wd low=%wd mid=%wd high=%wd\n", block_len, n0, low, mid, high); + for (slong m = 0; m < sum->nsols; m++) + { + acb_poly_set(tmp, sum->sol[m].series); + acb_poly_truncate(tmp, mid - n0); + flint_printf("series[%wd]=%{acb_poly}\n", m, tmp); + } + acb_poly_clear(tmp); + */ + + if (done) + return 1; + + slong resid_len = FLINT_MIN(high - mid, block_len); + apply_diffop(sum, low, mid, mid + resid_len); + + /* + for (slong m = 0; m < sum->nsols; m++) + flint_printf("series+res[%wd]=%{acb_poly}\n", m, sum->sol[m].series); + for (slong j = 0; j < sum->npts; j++) + for (slong m = 0; m < sum->nsols; m++) + flint_printf("sums[%wd][%wd]=%{acb_poly}\n", m, j, sum->sol[m].sums + j); + */ + + /* Doing the convergence test here allows the call to apply_diffop to be + shared, but otherwise it would make sense as well to do it in or just + after sum_forward_1. */ + + if (mid % stride == 0) + { + /* TODO we currently rely on the choice of block_len and stride to + ensure that sum_done has access to (1) enough terms for a robust + quick convergence check and (2) the full residual */ + FLINT_ASSERT(sum->n >= sum->n0 + sum->dop_clen - 1); /* (1) */ + FLINT_ASSERT(resid_len >= sum->dop_clen - 1); /* (2) */ + if (acb_ode_sum_done(sum, stride, prec)) + return 1; + + /* gradually decrease the working precision when it appears wasteful */ + slong prec_wanted = 0; + for (slong m = 0; m < sum->nsols; m++) + prec_wanted = FLINT_MAX(prec_wanted, sum->sol[m].cvest->prec_wanted); + if (prec_wanted < sum->wp) + sum->wp -= (sum->wp - prec_wanted)/2; + } + + if (mid - low >= block_len) + { + /* extend or shift buffer + XXX TBI once we have a better representation of vectors of + polynomials */ + if(sum->flags & ACB_ODE_WANT_SERIES) + /* support returning series coefficients even when the number of + terms is not known in advance */ + for (slong m = 0; m < sum->nsols; m++) + acb_ode_sol_fit_length(sum->sol + m, high + block_len); + else + discard_block(sum, block_len); + } + + /* XXX terminal recursion --> loop for -O0? */ + return _acb_ode_sum_forward_divconquer(sum, high, block_len, stride, prec); +} diff --git a/src/acb_ode/tail_bound_jet.c b/src/acb_ode/tail_bound_jet.c new file mode 100644 index 0000000000..5ef95631e7 --- /dev/null +++ b/src/acb_ode/tail_bound_jet.c @@ -0,0 +1,272 @@ +#include "acb_ode.h" +#include "arb_poly.h" + + +static void +arb_poly_taylor_shift_trunc(arb_poly_t res, const arb_poly_t pol, + const arb_t a, slong ord, slong prec) +{ + /* todo: ord is typically small, evaluate derivatives... */ + arb_poly_taylor_shift(res, pol, a, prec); + arb_poly_truncate(res, ord); +} + +static void +_arb_poly_product_roots_trunc(arb_ptr res, arb_ptr rts, slong rts_len, + slong trunc, slong prec) +{ + FLINT_ASSERT(trunc >= 1); + + if (rts_len == 1) + arb_neg(res, rts); + else if (rts_len > 1) + { + slong mid = (rts_len + 1) / 2; + slong len1 = FLINT_MIN(mid + 1, trunc); + slong len2 = FLINT_MIN(rts_len - mid + 1, trunc); + + arb_ptr scratch = _arb_vec_init(len1 + len2); + + _arb_poly_product_roots_trunc(scratch, rts, mid, trunc, prec); + _arb_poly_product_roots_trunc(scratch + len1, rts + mid, rts_len - mid, + trunc, prec); + _arb_poly_mullow(res, scratch, len1, scratch + len1, len2, + FLINT_MIN(rts_len, trunc), prec); + + _arb_vec_clear(scratch, len1 + len2); + } + + if (trunc > rts_len) + arb_one(res + rts_len); +} + +static void +hexp_series(arb_poly_t res, const arb_poly_t itg_pol, slong shift, + const arb_poly_t itg_num, const arb_poly_t itg_den, + slong ord, slong prec) +{ + if (ord > shift + 1) { + /* XXX Not yet tested(?) + * 1/itg_den could be cached. + * And couldn't we just choose pol_part_len so that this case never + * occurs? */ + arb_poly_div_series(res, itg_num, itg_den, ord - 1 - shift, prec); + arb_poly_shift_left(res, res, shift); + } + arb_poly_add_series(res, res, itg_pol, ord - 1, prec); + + arb_poly_neg(res, res); + arb_poly_integral(res, res, prec); + arb_poly_exp_series(res, res, ord, prec); +} + +/* [M19], Algorithm 6.11, step 5, + * specialized to w⁻¹·â(w) = itg_pol + w^itg_rat_shift·itg_num/itg_den, + * pre = g with an implicit x^n factor. + * + * (To reiterate: itg_rat_shift is about the expression of \hat{f} and has + * nothing to do with the w^{n-1} factor under the outer integral.) */ + +static void +maj_pre_num(arb_poly_t pre, const arb_poly_t nres_maj, slong n, + const arb_poly_t itg_pol, slong itg_rat_shift, + const arb_poly_t itg_num, const arb_poly_t itg_den, + slong prec) +{ + arb_poly_t inv_exp; + arb_poly_init(inv_exp); + + arb_poly_fit_length(pre, nres_maj->length); + _arb_poly_set_length(pre, nres_maj->length); + for (slong j = 0; j < pre->length; j++) + arb_mul_si(pre->coeffs + j, nres_maj->coeffs + j, n + j, prec); + + hexp_series(inv_exp, itg_pol, itg_rat_shift, itg_num, + itg_den, nres_maj->length, prec); + + /* polynomial factor in rhs of majorizing equation */ + arb_poly_mullow(pre, pre, inv_exp, nres_maj->length, prec); + /* flint_printf("rhs=%{arb_poly}\n", pre); */ + + for (slong j = 0; j < pre->length; j++) { + arb_ptr c = pre->coeffs + j; + if (arb_contains_positive(c)) + arb_div_si(c, c, n + j, prec); + else + arb_zero(c); + } + _arb_poly_normalise(pre); + + arb_poly_clear(inv_exp); +} + +/* Bound the series expansion at rad of + * + * x^n·pre_num/den × exp ∫(itg_pol + t^itg_rat_shift·itg_num/den). + * + * [M19], Algorithm 8.1, specialized to this case and with a slightly different + * output format. + */ + +static void +maj_jet(arb_poly_t res, slong n, const arb_poly_t pre_num, + const arb_poly_t itg_pol, slong itg_rat_shift, const arb_poly_t itg_num, + const arb_t cst, arb_ptr den_rt, slong den_rt_len, + const arb_t rad, slong ord, slong prec) +{ + /* flint_printf("== maj_jet n=%ld ord=%ld rad=%{arb} ==\n", n, ord, rad); */ + + arb_poly_t pre_ser, invden_ser, shx_ser, int_pol_ser, int_rat_ser, int_ser; + arb_ptr rt_rad = _arb_vec_init(den_rt_len); + + arb_poly_init(pre_ser); + arb_poly_init(invden_ser); + arb_poly_init(shx_ser); + arb_poly_init(int_pol_ser); + arb_poly_init(int_rat_ser); + arb_poly_init(int_ser); + + /* Denominator: (1/den)(rad+ε) */ + + for (slong i = 0; i < den_rt_len; i++) + arb_sub(rt_rad + i, den_rt + i, rad, prec); + arb_poly_fit_length(invden_ser, ord); + _arb_poly_product_roots_trunc(invden_ser->coeffs, + rt_rad, den_rt_len, ord, prec); + _arb_poly_set_length(invden_ser, ord); + _arb_poly_normalise(invden_ser); + + if (arb_contains_zero(invden_ser->coeffs)) + { + arb_poly_fit_length(res, ord); + _arb_vec_indeterminate(res->coeffs, ord); + _arb_poly_set_length(res, ord); + goto cleanup; + } + arb_poly_inv_series(invden_ser, invden_ser, ord, prec); + // XXX not 100% sure this cst should be there + arb_poly_scalar_mul(invden_ser, invden_ser, cst, prec); + + /* flint_printf("invden_ser=%{arb_poly}\n", invden_ser); */ + + /* Rational part */ + + arb_poly_set_coeff_si(shx_ser, 1, 1); + arb_poly_set_coeff_arb(shx_ser, 0, rad); + arb_poly_pow_ui_trunc_binexp(shx_ser, shx_ser, n, ord, prec); + + arb_poly_taylor_shift_trunc(pre_ser, pre_num, rad, ord, prec); + arb_poly_mullow(pre_ser, pre_ser, shx_ser, ord, prec); + /* pre_ser corresponds to rat_ser in ore_algebra */ + /* subtract den' from the numerator of itg_rat below instead? */ + arb_poly_mullow(pre_ser, pre_ser, invden_ser, ord, prec); + + /* flint_printf("pre_ser=%{arb_poly}\n", pre_ser); */ + + /* Exponential part. Up to the sign in the exponential, this is the + * same expression as in hexp_series, but now we want a bound on the series + * expansion at rad instead of a series expansion at zero. For this we bound + * ∫(f+p/q) by ∫f + (∫p)/q, then compose with rad+ε. */ + + arb_poly_integral(int_pol_ser, itg_pol, prec); + arb_poly_taylor_shift_trunc(int_pol_ser, int_pol_ser, rad, ord, prec); + + arb_poly_fit_length(int_rat_ser, itg_num->length); + _arb_poly_set_length(int_rat_ser, itg_num->length); + itg_rat_shift++; + for (slong i = 0; i < itg_num->length; i++) + arb_div_si(int_rat_ser->coeffs + i, itg_num->coeffs + i, + itg_rat_shift + i, prec); + arb_poly_taylor_shift_trunc(int_rat_ser, int_rat_ser, rad, ord, prec); + + arb_poly_set_arb(shx_ser, rad); + arb_poly_set_coeff_si(shx_ser, 1, 1); + arb_poly_pow_ui_trunc_binexp(shx_ser, shx_ser, itg_rat_shift, ord, prec); + + arb_poly_mullow(int_rat_ser, int_rat_ser, shx_ser, ord, prec); + arb_poly_mullow(int_rat_ser, int_rat_ser, invden_ser, ord, prec); + + /* flint_printf("int_rat_ser=%{arb_poly}\n", int_ser); */ + + arb_poly_add(int_ser, int_pol_ser, int_rat_ser, prec); + + /* flint_printf("int_ser=%{arb_poly}\n", int_ser); */ + + arb_poly_exp_series(res, int_ser, ord, prec); + + /* flint_printf("exp_ser=%{arb_poly}\n", res); */ + + /* Final series */ + + arb_poly_mullow(res, res, pre_ser, ord, prec); + +cleanup: + arb_poly_clear(pre_ser); + arb_poly_clear(invden_ser); + arb_poly_clear(shx_ser); + arb_poly_clear(int_pol_ser); + arb_poly_clear(int_rat_ser); + arb_poly_clear(int_ser); + _arb_vec_clear(rt_rad, den_rt_len); +} + +void +acb_ode_tail_bound_jet_precomp(arb_poly_t res, + const acb_ode_bound_t bound, slong n, + const arb_poly_t itg_pol, const arb_poly_t itg_num, /* n0 <= n */ + const arb_poly_t nres_maj, /* implicit x^n factor */ + const arb_t rad, slong ord, slong prec) +{ + // flint_printf("== acb_ode_tail_bound_jet_precomp n=%ld ==\n", n); + + arb_poly_t pre; + + // flint_printf("nres_maj=%{arb_poly}\n", nres_maj); + + if (arb_poly_is_zero(nres_maj)) + { + arb_poly_zero(res); + return; + } + + arb_poly_init(pre); + + /* The solution of the homogeneous part of the majorizing equation obtained + * after specialization at n is exp(∫w⁻¹·â(w)) where + * w⁻¹·â(w) = itg_pol + w^{pol_part_len-1} · itg_num/itg_den. + * In other words, the itg_pol term already incorporates the w⁻¹ factor, but + * itg_num does not. */ + slong rat_shift = bound->pol_part_len - 1; + + /* Variation of constants */ + maj_pre_num(pre, nres_maj, n, itg_pol, rat_shift, itg_num, + bound->den_lbound, prec); + + /* Evaluate the resulting majorant */ + maj_jet(res, n, pre, itg_pol, rat_shift, itg_num, bound->cst, + bound->den_rt, bound->den_rt_len, rad, ord, prec); + + // flint_printf("n=%ld pre=%{arb_poly}\n", n, pre); + // flint_printf("n=%ld maj_jet=%{arb_poly}\n", n, res); + + arb_poly_clear(pre); +} + +void acb_ode_tail_bound_jet(arb_poly_t res, + const acb_ode_group_t group, const acb_ode_bound_t bound, slong n, + slong nlogs, const arb_poly_t nres_maj, const arb_t rad, + slong ord, slong prec) +{ + arb_poly_t itg_pol, itg_num; + + arb_poly_init(itg_pol); + arb_poly_init(itg_num); + + acb_ode_bound_precompute_integrand(itg_pol, itg_num, group, bound, n, nlogs, + prec); + acb_ode_tail_bound_jet_precomp(res, bound, n, itg_pol, itg_num, nres_maj, + rad, ord, prec); + + arb_poly_clear(itg_pol); + arb_poly_clear(itg_num); +} diff --git a/src/acb_ode/test/main.c b/src/acb_ode/test/main.c new file mode 100644 index 0000000000..afc74ae27a --- /dev/null +++ b/src/acb_ode/test/main.c @@ -0,0 +1,12 @@ +#include "t-apply_diffop.c" +#include "t-fundamental_matrix.c" +#include "t-sum_divconquer.c" + +test_struct tests[] = +{ + // TEST_FUNCTION(acb_ode_apply_diffop), + TEST_FUNCTION(acb_ode_fundamental_matrix), + TEST_FUNCTION(acb_ode_sum_divconquer) +}; + +TEST_MAIN(tests) diff --git a/src/acb_ode/test/t-apply_diffop.c b/src/acb_ode/test/t-apply_diffop.c new file mode 100644 index 0000000000..de003c98e8 --- /dev/null +++ b/src/acb_ode/test/t-apply_diffop.c @@ -0,0 +1,98 @@ +#include "test_helpers.h" +#include "acb_types.h" +#include "acb_poly.h" +#include "fmpq_types.h" +#include "fmpq_poly.h" +#include "acb_ode.h" + +TEST_FUNCTION_START(acb_ode_apply_diffop, state) +{ + fmpq_poly_t tmp; + fmpq_poly_init(tmp); + + for (slong iter = 0; iter < 10000 * flint_test_multiplier(); iter++) + { + slong dop_len = n_randint(state, 3); + acb_poly_struct * dop = _acb_poly_vec_init(dop_len); + for (slong i = 0; i < dop_len; i++) + { + acb_poly_randtest(dop + i, state, 1 + n_randint(state, 10), 1 + + n_randint(state, 200), 5); + /* fmpq_poly_randtest(tmp, state, 1 + n_randint(state, 5), 1 + n_randint(state, 5)); */ + /* flint_printf("dop[%wd] = %{fmpq_poly}\n", i, tmp); */ + /* acb_poly_set_fmpq_poly(dop + i, tmp, 100); */ + } + + acb_t expo; + acb_init(expo); + acb_randtest(expo, state, 1 + n_randint(state, 200), 5); + slong offset = n_randint(state, 1000); + + slong logs = n_randint(state, 4); + acb_poly_struct * f = _acb_poly_vec_init(logs); + for (slong k = 0; k < logs; k++) + { + acb_poly_randtest(f + k, state, 1 + n_randint(state, 20), 1 + + n_randint(state, 200), 5); + /* fmpq_poly_randtest(tmp, state, 1 + n_randint(state, 10), 1 + n_randint(state, 5)); */ + /* flint_printf("f[%wd] = %{fmpq_poly}\n", k, tmp); */ + /* acb_poly_set_fmpq_poly(f + k, tmp, 100); */ + } + + slong len1 = 1 + n_randint(state, 5); + slong len2 = n_randint(state, 5); + slong len = len1 + len2; + slong start = n_randint(state, len); + + slong prec = 2 + n_randint(state, 100); + + acb_poly_struct * g1 = _acb_poly_vec_init(logs); + acb_poly_struct * g2 = _acb_poly_vec_init(logs); + acb_poly_struct * g3 = _acb_poly_vec_init(logs); + + acb_ode_apply_diffop_polmul(g1, dop, dop_len, expo, offset, f, logs, start, len, prec); + + _acb_poly_vec_fit_length(g2, logs, len1 + len2); + /* TODO play with foff and offset too */ + _acb_ode_apply_diffop_polmul(g2, 0, dop, dop_len, expo, offset, f, 0, start + len1, logs, start, len1, prec); + _acb_ode_apply_diffop_polmul(g2, len1, dop, dop_len, expo, offset, f, 0, start + len, logs, start + len1, len2, prec); + _acb_poly_vec_set_length(g2, logs, len); + _acb_poly_vec_normalise(g2, logs); + + acb_ode_apply_diffop_basecase(g3, dop, dop_len, expo, offset, f, logs, start, len, prec); + + if (!_acb_poly_vec_overlaps(g1, g2, logs) + || !_acb_poly_vec_overlaps(g1, g3, logs) + || !_acb_poly_vec_overlaps(g2, g3, logs)) + { + flint_printf("FAIL\n\n"); + for (slong i = 0; i < dop_len; i++) + flint_printf("dop[%wd] = %{acb_poly}\n", i, dop + i); + for (slong k = 0; k < logs; k++) + flint_printf("f[%wd] = %{acb_poly}\n", k, f + k); + flint_printf("expo = %{acb}\n", expo); + flint_printf("offset = %wd\n", offset); + flint_printf("start = %wd\n", start); + flint_printf("len1 = %wd\n", len1); + flint_printf("len2 = %wd\n", len2); + for (slong k = 0; k < logs; k++) + { + flint_printf("g1[%wd] = %{acb_poly}\n", k, g1 + k); + flint_printf("g2[%wd] = %{acb_poly}\n", k, g2 + k); + flint_printf("g3[%wd] = %{acb_poly}\n", k, g3 + k); + } + flint_abort(); + } + + /* todo: left or right-multiply by something; aliasing tests */ + + _acb_poly_vec_clear(dop, dop_len); + acb_clear(expo); + _acb_poly_vec_clear(f, logs); + _acb_poly_vec_clear(g1, logs); + _acb_poly_vec_clear(g2, logs); + _acb_poly_vec_clear(g3, logs); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/acb_ode/test/t-fundamental_matrix.c b/src/acb_ode/test/t-fundamental_matrix.c new file mode 100644 index 0000000000..14d11a488a --- /dev/null +++ b/src/acb_ode/test/t-fundamental_matrix.c @@ -0,0 +1,258 @@ +#include "test_helpers.h" +#include "acb_types.h" +#include "acb_mat.h" +#include "acb_ode.h" +#include "acb.h" +#include "arf.h" +#include "gr_poly.h" + +FLINT_DLL extern gr_static_method_table _ca_methods, _qqbar_methods; + +TEST_FUNCTION_START(acb_ode_fundamental_matrix, state) +{ + gr_ctx_t ZZ, QQ, ZZi, QQbar, CaCC; + gr_ctx_init_fmpz(ZZ); + gr_ctx_init_fmpq(QQ); + gr_ctx_init_fmpzi(ZZi); + gr_ctx_init_complex_qqbar(QQbar); + gr_ctx_init_complex_ca(CaCC); + + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + int status = GR_SUCCESS; + + gr_ctx_t CC, Pol; + gr_ore_poly_ctx_t Dop; + gr_ptr pert; + gr_ore_poly_t dop, dop1; + acb_t r; + arf_t frac; + + acb_init(r); + arf_init(frac); + + slong prec = 2 + n_randint(state, 100); + slong prec1 = 2 + n_randint(state, prec); + gr_ctx_init_complex_acb(CC, prec); + + int expect_success = 0; + if (!n_randint(state, 8)) + gr_ore_poly_ctx_init_randtest2(Pol, Dop, state); + else + { + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random_poly(Pol, state); + break; + case 1: + gr_ctx_init_gr_poly(Pol, CC); + break; + case 2: + gr_ctx_init_gr_poly(Pol, QQ); + break; + case 3: + gr_ctx_init_gr_poly(Pol, ZZi); + break; + case 4: + gr_ctx_init_gr_poly(Pol, QQbar); + break; + case 5: + gr_ctx_init_gr_poly(Pol, CaCC); + break; + default: + gr_ctx_init_gr_poly(Pol, ZZ); + expect_success = 1; + break; + } + ore_algebra_t alg = n_randint(state, 4) + ? ORE_ALGEBRA_EULER_DERIVATIVE + : ORE_ALGEBRA_DERIVATIVE; + if (alg != ORE_ALGEBRA_EULER_DERIVATIVE) // XXX + expect_success = 0; + gr_ore_poly_ctx_init(Dop, Pol, 0, alg); + } + + GR_TMP_INIT(pert, Pol); + gr_ore_poly_init(dop, Dop); + gr_ore_poly_init(dop1, Dop); + + /* todo: Test expos != NULL and lcrt != NULL (by generating the leading + coefficient and/or indicial polynomial from its roots, cf. + t-shiftless_decomposition). Update expect_success as needed. + + todo: With reasonable probability, generate evaluation points in the + disk of convergence. */ + + /* also generate operators with coefficients of higher degree? */ + slong dop_len_bound = 20; + if (Pol->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(Pol)->methods == _ca_methods + || POLYNOMIAL_ELEM_CTX(Pol)->methods == _qqbar_methods)) + dop_len_bound = 4; + status |= gr_ore_poly_randtest(dop, state, + n_randint(state, dop_len_bound), Dop); + status |= gr_randtest_not_zero(pert, state, Pol); + /* dop1 may have extra singular points; ensure however that there is no + extra singular point at zero in interesting cases */ + if (Pol->which_ring == GR_CTX_GR_POLY) + status |= gr_randtest_not_zero( + gr_poly_coeff_ptr(pert, 0, POLYNOMIAL_ELEM_CTX(Pol)), + state, + POLYNOMIAL_ELEM_CTX(Pol)); + status |= gr_ore_poly_other_mul(dop1, pert, Pol, dop, Dop); + slong order = dop->length - 1; + + // XXX large npts hits hard entire function case too often + // slong npts = n_randint(state, n_randint(state, 8) ? 3 : 100); + slong npts = n_randint(state, 3); + acb_ptr pts = _acb_vec_init(npts * 2); + acb_ptr pts1 = pts + npts; + for (slong p = 0; p < npts; p++) + { + acb_randtest_precise(pts + p, state, prec + 30, 3); + if (n_randint(state, 2)) + acb_set(pts1 + p, pts + p); + else + { + acb_randtest_special(r, state, prec, 3); + acb_add(pts1 + p, pts + p, r, prec); + acb_sub(pts1 + p, pts1 + p, r, prec); + } + } + + acb_ode_basis_t basis = n_randint(state, ACB_ODE_NUM_BASES); + + acb_mat_struct * mat = flint_malloc(sizeof(acb_mat_struct) * npts * 3); + acb_mat_struct * mat1 = mat + npts; + acb_mat_struct * mat2 = mat1 + npts; + for (slong p = 0; p < npts; p++) + { + slong nrows = n_randint(state, order + 2); + acb_mat_init(mat + p, nrows, order); + acb_mat_init(mat1 + p, nrows, order); + acb_mat_init(mat2 + p, nrows, order); + } + + if (status != GR_SUCCESS) + goto cleanup; + + flint_printf("******** testing: ********\n", Dop); + flint_printf("Dop = %{gr_ctx}\n", Dop); + flint_printf("dop = %{gr}\n", dop, Dop); + flint_printf("dop1 = %{gr}\n", dop1, Dop); + flint_printf("pts = %{acb*}\n", pts, npts); + flint_printf("pts1 = %{acb*}\n", pts1, npts); + flint_printf("prec = %wd, prec1 = %wd\n", prec, prec1); + + status |= acb_ode_fundamental_matrix_vec(mat, dop, Dop, NULL, NULL, + pts, npts, basis, prec + 30); + flint_printf("status0 = %wd\n", status); + + if (expect_success && status == GR_UNABLE) + { + flint_printf("FAIL: unable\n"); + flint_abort(); + } + + if (status != GR_SUCCESS) + { + goto cleanup; + } + + int test0_is_accurate = 1; + for (slong p = 0; p < npts; p++) + for (slong i = 0; i < acb_mat_nrows(mat + p); i++) + for (slong j = 0; j < order; j++) + { + slong acc = acb_rel_accuracy_bits( + acb_mat_entry_ptr(mat + p, i, j)); + if (acc < prec) + { + flint_printf("poor accuracy @pt%wd (%wd,%wd) = %wd\n", + p, i, j, acc); + test0_is_accurate = 0; + } + } + + for (slong p = 0; p < npts; p++) // tmp + { + flint_printf("res@%{acb} =\n", pts + p); + acb_mat_printd(mat + p, prec+4); + } + + status |= acb_ode_fundamental_matrix_vec(mat1, dop, Dop, NULL, NULL, + pts1, npts, basis, prec1); + flint_printf("status1 = %wd\n", status); + + if (status == GR_SUCCESS) // tmp + for (slong p = 0; p < npts; p++) + { + flint_printf("res1@%{acb} =\n", pts1 + p); + acb_mat_printd(mat1 + p, prec+4); + } + + if (status == GR_SUCCESS) + { + for (slong p = 0; p < npts; p++) + if (test0_is_accurate) + { + if (!acb_mat_contains(mat1 + p, mat + p)) + { + flint_printf("FAIL: containment (@pt%ld)\n\n", p); + flint_abort(); + } + } + else if (!acb_mat_overlaps(mat1 + p, mat + p)) + { + flint_printf("FAIL: overlap 1 (@pt%ld)\n\n", p); + flint_abort(); + } + } + + if (Pol->which_ring == GR_CTX_GR_POLY + && POLYNOMIAL_ELEM_CTX(Pol)->methods == _ca_methods) + goto cleanup; + + status |= acb_ode_fundamental_matrix_vec(mat2, dop1, Dop, NULL, NULL, + pts, npts, basis, prec); + flint_printf("status2 = %wd\n", status); + + if (status == GR_SUCCESS) + { + for (slong p = 0; p < npts; p++) + if (!acb_mat_overlaps(mat2 + p, mat + p)) + { + flint_printf("FAIL: overlap 2\n\n"); + flint_abort(); + } + } + +cleanup: + for (slong p = 0; p < npts; p++) + { + acb_mat_clear(mat2 + p); + acb_mat_clear(mat1 + p); + acb_mat_clear(mat + p); + } + flint_free(mat); + _acb_vec_clear(pts, npts * 2); + gr_ore_poly_clear(dop1, Dop); + gr_ore_poly_clear(dop, Dop); + GR_TMP_CLEAR(pert, Pol); + gr_ctx_clear(Dop); + gr_ctx_clear(Pol); + gr_ctx_clear(CC); + arf_clear(frac); + acb_clear(r); + } + + gr_ctx_clear(CaCC); + gr_ctx_clear(QQbar); + gr_ctx_clear(ZZi); + gr_ctx_clear(QQ); + gr_ctx_clear(ZZ); + + TEST_FUNCTION_END(state); +} + diff --git a/src/acb_ode/test/t-sum_divconquer.c b/src/acb_ode/test/t-sum_divconquer.c new file mode 100644 index 0000000000..c365ee61c7 --- /dev/null +++ b/src/acb_ode/test/t-sum_divconquer.c @@ -0,0 +1,124 @@ +#include "test_helpers.h" +#include "acb_types.h" +#include "acb_poly.h" +#include "acb_ode.h" + +// XXX very very incomplete and too involved + +TEST_FUNCTION_START(acb_ode_sum_divconquer, state) +{ + + for (slong iter = 0; iter < 10 * flint_test_multiplier(); iter++) + { + /* classical functions -- well for now just atan */ + + acb_ode_sum_t sum; + + slong dop_order = 2; + acb_ode_sum_init(sum, dop_order + 1, + n_randint(state, 5), + dop_order, + 1 + n_randint(state, 3)); + + /* sum_set_diffop */ + acb_poly_set_coeff_si(sum->dop + 2, 2, 1); + acb_poly_set_coeff_si(sum->dop + 2, 0, 1); + acb_poly_set_coeff_si(sum->dop + 1, 2, 1); + acb_poly_set_coeff_si(sum->dop + 1, 0, -1); + mag_one(sum->cvrad); + + acb_ode_sum_set_ordinary(sum); + acb_ode_sum_set_ini_echelon(sum); + + slong prec = n_randint(state, 1000); + /* XXX at the moment nterms must be a multiple of the degrees due to + incomplete termination logic in sum_divconquer in other cases */ + slong nterms = n_randint(state, 500) * 2; + + for (slong i = 0; i < sum->npts; i++) + { + acb_ptr a = sum->pts + i; + acb_randtest_precise(a, state, prec, 0); + acb_div_si(a, a, 4, prec); + } + + /* note that sum->bound is null */ + sum->flags = ACB_ODE_APPROX + | (n_randint(state, 2) * ACB_ODE_WANT_SERIES); + + acb_ode_sum_divconquer(sum, nterms, prec); + + acb_poly_struct * ref = _acb_poly_vec_init(dop_order); + acb_poly_one(ref); /* atan: first sol=1 */ + + if (prec <= 12) + goto cleanup; + + if (nterms >= prec) /* todo: test rigorous version */ + { + for (slong i = 0; i < sum->npts; i++) + { + acb_poly_zero(ref + 1); /* 2nd sol = atan */ + acb_poly_set_coeff_acb(ref + 1, 0, sum->pts + i); + acb_poly_set_coeff_si(ref + 1, 1, 1); + /* APPROX => compute the ref to lower precision */ + acb_poly_atan_series(ref + 1, ref + 1, sum->nder, prec - 8); + + for (slong m = 0; m < dop_order; m++) + { + if (!_acb_poly_overlaps( + acb_ode_sol_sum_ptr(sum->sol + m, i, 0, 0), + sum->sol->nder, + (ref + m)->coeffs, + (ref + m)->length)) + { + flint_printf("FAIL (nonrigorous partial sum)\n\n"); + flint_printf("prec = %wd, nterms = %wd\n", prec, nterms); + flint_printf("i = %wd, m = %wd, a = %{acb}\n", i, m, sum->pts + i); + for (slong j = 0; j < 2; j++) + flint_printf("sum[%wd] = %{acb*}\n", j, acb_ode_sol_sum_ptr(sum->sol + m, i, 0, 0), sum->sol->nder); + + flint_printf("ref = %{acb_poly}\n", ref + m); + flint_abort(); + } + } + } + } + + if (sum->flags & ACB_ODE_WANT_SERIES) + { + if (sum->npts == 0 && sum->n != nterms) + flint_printf("FAIL (wrong number of terms)\n"); + + acb_poly_zero(ref + 1); + acb_poly_set_coeff_si(ref + 1, 1, 1); + /* APPROX => compute the ref to lower precision */ + acb_poly_atan_series(ref + 1, ref + 1, nterms, prec - 8); + + for (slong m = 0; m < dop_order; m++) + { + /* in general, only the first sum->n terms are correct */ + for (slong n = 0; n < sum->n; n++) + { + if (n >= (ref + m)->length) + break; // for now + if (!acb_overlaps((sum->sol[m].series)->coeffs + n, + (ref + m)->coeffs + n)) + { + flint_printf("FAIL (series)\n\n"); + flint_printf("m = %wd, nterms = %wd, sum->n = %wd, n=%wd\n", m, nterms, sum->n, n); + flint_printf("series = %{acb_poly}\n", sum->sol[m].series); + flint_printf("ref = %{acb_poly}\n", ref + m); + flint_abort(); + } + } + } + } + +cleanup: + _acb_poly_vec_clear(ref, dop_order); + acb_ode_sum_clear(sum); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/acb_poly.h b/src/acb_poly.h index 2682619ec2..1c860def55 100644 --- a/src/acb_poly.h +++ b/src/acb_poly.h @@ -716,6 +716,19 @@ void cd_poly_horner(double* results_r, double* results_i, double cd_poly_find_roots(double * z, const double * p, const double * z0, slong n, slong num_iter, double reltol); double _acb_poly_find_roots_double(acb_ptr roots, acb_srcptr poly, acb_srcptr initial, slong len, slong maxiter, slong prec); +/* Vector functions */ + +acb_poly_struct * _acb_poly_vec_init(slong n); +void _acb_poly_vec_clear(acb_poly_struct *vec, slong n); +void _acb_poly_vec_zero(acb_poly_struct *dest, slong n); +void _acb_poly_vec_set(acb_poly_struct *dest, const acb_poly_struct *src, slong n); +void _acb_poly_vec_set_block(acb_poly_struct *dest, const acb_poly_struct *src, + slong n, slong base, slong len); +void _acb_poly_vec_fit_length(acb_poly_struct *vec, slong n, slong len); +void _acb_poly_vec_set_length(acb_poly_struct *vec, slong n, slong len); +void _acb_poly_vec_normalise(acb_poly_struct *vec, slong n); +int _acb_poly_vec_overlaps(acb_poly_struct *vec1, acb_poly_struct *vec2, slong n); +slong _acb_poly_vec_length(const acb_poly_struct * vec, slong n); #ifdef __cplusplus } diff --git a/src/acb_poly/gr_stuff.c b/src/acb_poly/gr_stuff.c new file mode 100644 index 0000000000..0ca79fab2a --- /dev/null +++ b/src/acb_poly/gr_stuff.c @@ -0,0 +1,70 @@ +#include "gr_ore_poly.h" +#include "gr_poly.h" + +int +gr_ore_poly_ctx_over_gr_poly_base_ptrs(gr_ctx_struct ** Scalars, + gr_ctx_struct ** Pol, const gr_ctx_t Ore) +{ + gr_ctx_struct * _Pol = GR_ORE_POLY_ELEM_CTX(Ore); + + ulong which_base = _Pol->which_ring; + if (which_base != GR_CTX_GR_POLY) + return GR_UNABLE; + + if (Pol != NULL) + * Pol = _Pol; + if (Scalars != NULL) + * Scalars = POLYNOMIAL_ELEM_CTX(_Pol); + + return GR_SUCCESS; +} + +int +_gr_ore_poly_indicial_polynomial_euler_derivative(gr_ptr ind, gr_srcptr op, + slong len, gr_ctx_t Ore) +{ + gr_ctx_struct * Scalars, * Pol; + + int status = GR_SUCCESS; + + GR_MUST_SUCCEED (gr_ore_poly_ctx_over_gr_poly_base_ptrs(&Scalars, &Pol, Ore)); + + for (slong i = 0; i < len; i++) + status = gr_poly_get_coeff_scalar( + GR_ENTRY(ind, i, Scalars->sizeof_elem), + GR_ENTRY(op, i, Pol->sizeof_elem), + 0, Scalars); + + return status; +} + +int +_gr_ore_poly_indicial_polynomial(gr_ptr ind, gr_srcptr op, slong len, gr_ctx_t Ore) +{ + switch (GR_ORE_POLY_CTX(Ore)->which_algebra) + { + case ORE_ALGEBRA_EULER_DERIVATIVE: + return _gr_ore_poly_indicial_polynomial_euler_derivative(ind, op, len, Ore); + default: + return GR_UNABLE; + } +} + +int +gr_ore_poly_indicial_polynomial(gr_poly_t ind, const gr_ore_poly_t op, gr_ctx_t Ore) +{ + gr_ctx_struct * Scalars; + + if (gr_ore_poly_ctx_over_gr_poly_base_ptrs(&Scalars, NULL, Ore) != GR_SUCCESS) + return GR_UNABLE; + + // XXX move to level depending on operator type? + gr_poly_fit_length(ind, op->length, Scalars); + + int status = _gr_ore_poly_indicial_polynomial(ind->coeffs, op->coeffs, op->length, Ore); + + _gr_poly_set_length(ind, op->length, Scalars); + _gr_poly_normalise(ind, Scalars); + + return status; +} diff --git a/src/acb_poly/vec.c b/src/acb_poly/vec.c new file mode 100644 index 0000000000..2af10b1e9b --- /dev/null +++ b/src/acb_poly/vec.c @@ -0,0 +1,89 @@ +#include "acb_poly.h" + +acb_poly_struct * +_acb_poly_vec_init(slong n) +{ + acb_poly_struct *vec = (acb_poly_struct *) flint_malloc(sizeof(acb_poly_struct) * n); + + for (slong i = 0; i < n; i++) + acb_poly_init(vec + i); + + return vec; +} + +void +_acb_poly_vec_clear(acb_poly_struct *vec, slong n) +{ + for (slong i = 0; i < n; i++) + acb_poly_clear(vec + i); + flint_free(vec); +} + +void +_acb_poly_vec_zero(acb_poly_struct *dest, slong n) +{ + for (slong i = 0; i < n; i++) + acb_poly_zero(dest + i); +} + +void +_acb_poly_vec_set(acb_poly_struct *dest, const acb_poly_struct *src, slong n) +{ + for (slong i = 0; i < n; i++) + acb_poly_set(dest + i, src + i); +} + +void +_acb_poly_vec_set_block(acb_poly_struct *dest, const acb_poly_struct *src, + slong n, slong start, slong len) +{ + for (slong k = 0; k < n; k++) + { + if (start >= (src + k)->length) + continue; + slong len1 = FLINT_MIN(len, (src + k)->length - start); + acb_poly_fit_length(dest + k, len1); + _acb_vec_set((dest + k)->coeffs, (src + k)->coeffs + start, len1); + _acb_poly_set_length(dest + k, len1); + _acb_poly_normalise(dest + k); + } +} + +void +_acb_poly_vec_fit_length(acb_poly_struct *vec, slong n, slong len) +{ + for (slong i = 0; i < n; i++) + acb_poly_fit_length(vec + i, len); +} + +void +_acb_poly_vec_set_length(acb_poly_struct *vec, slong n, slong len) +{ + for (slong i = 0; i < n; i++) + _acb_poly_set_length(vec + i, len); +} + +void +_acb_poly_vec_normalise(acb_poly_struct *vec, slong n) +{ + for (slong i = 0; i < n; i++) + _acb_poly_normalise(vec + i); +} + +int +_acb_poly_vec_overlaps(acb_poly_struct *vec1, acb_poly_struct *vec2, slong n) +{ + for (slong i = 0; i < n; i++) + if (!acb_poly_overlaps(vec1 + i, vec2 + i)) + return 0; + return 1; +} + +slong +_acb_poly_vec_length(const acb_poly_struct * vec, slong n) +{ + slong len = 0; + for (slong i = 0; i < n; i++) + len = FLINT_MAX(len, vec[i].length); + return len; +}