fundamental matrices

Author: mezzarobba

doc/source/gr.rst

@@ -175,8 +175,8 @@ Flags can be OR'ed and checked only at the top level of a computation
 to avoid complex control flow::
 
     status = GR_SUCCESS;
-    gr |= gr_add(res, a, b, ctx);
-    gr |= gr_pow_ui(res, res, 2, ctx);
+    status |= gr_add(res, a, b, ctx);
+    status |= gr_pow_ui(res, res, 2, ctx);
     ...
 
 If we do not care about recovering from *undefined*/*unknown* results,

examples/holonomic.c

@@ -1,9 +1,9 @@
 #include "acb_types.h"
 #include "acb_poly.h"
 #include "acb_mat.h"
-#include "fmpq_types.h"
-#include "fmpq_poly.h"
 #include "acb_holonomic.h"
+#include "gr.h"
+#include "gr_poly.h"
 
 
 void
@@ -43,7 +43,7 @@ ordinary(void)
     acb_poly_set_coeff_si(ctx->dop + 1, 0, 2);
     acb_poly_set_coeff_si(ctx->dop + 0, 4, 7);
 
-    ctx->flags |= ACB_HOLONOMIC_WANT_SERIES;
+    /* ctx->flags |= ACB_HOLONOMIC_WANT_SERIES; */
 
     acb_holonomic_sum_ordinary(ctx);
     acb_holonomic_sum_canonical_basis(ctx);
@@ -142,9 +142,14 @@ 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_holonomic_sum_context_t ctx;
 
     acb_set_si(kappa, 2);
     acb_set_si(mu2, 3);
@@ -154,8 +159,6 @@ whittaker_m(void)  /* non-integer exponent, no logs */
     slong dop_order = 2;
     slong len = dop_order;
 
-    acb_holonomic_sum_context_t ctx;
-
     acb_holonomic_sum_context_init(ctx, dop_order + 1, 1, 1, len);
 
     acb_poly_set_coeff_si(ctx->dop + 2, 0, 4);
@@ -166,7 +169,7 @@ whittaker_m(void)  /* non-integer exponent, no logs */
     acb_addmul_si((ctx->dop + 0)->coeffs + 0, mu2, -4, prec);
 
     acb_sqrt(ctx->expo, mu2, prec);
-    acb_add(ctx->expo, ctx->expo, half, prec);  /* sub for other expo */
+    acb_add(ctx->expo, ctx->expo, half, prec);  /* other expo = mu2 - 1/2 */
 
     ctx->sing_shifts[0].n = 0;
     ctx->sing_shifts[0].mult = 1;
@@ -180,38 +183,120 @@ whittaker_m(void)  /* non-integer exponent, no logs */
 
     acb_holonomic_sum_divconquer(ctx, prec);
 
-    flint_printf("(%{acb} + x)^(%{acb}) * (%{acb_poly})\n",
-                 ctx->pts, ctx->expo,
-                 acb_holonomic_sol_sum_ptr(ctx->sol, 0, 0));
+    acb_poly_struct * f = acb_holonomic_sol_sum_ptr(ctx->sol, 0, 0);
 
-    acb_poly_t tmp;
-    acb_poly_init(tmp);
-    acb_poly_set_coeff_si(tmp, 1, 1);
-    acb_poly_set_coeff_acb(tmp, 0, ctx->pts);
-    acb_poly_pow_acb_series(tmp, tmp, ctx->expo, len, prec);
-    acb_poly_mullow(tmp, tmp, acb_holonomic_sol_sum_ptr(ctx->sol, 0, 0), len, prec);
+    flint_printf("(%{acb} + x)^(%{acb}) * (%{acb_poly})\n",
+                 ctx->pts, ctx->expo, f);
 
-    flint_printf("M(%{acb} + x) = %{acb_poly} + O(x^%wd)\n", ctx->pts, tmp, len);
+    _acb_holonomic_sol_value(val, ctx->expo, f, ctx->sol[0].nlogs, ctx->pts,
+                             ctx->nder, 1, ctx->prec);
 
-    acb_poly_clear(tmp);
+    flint_printf("M(%{acb} + x) = %{acb_poly} + O(x^%wd)\n", ctx->pts, val, len);
 
     acb_holonomic_sum_context_clear(ctx);
     acb_clear(kappa);
     acb_clear(mu2);
     acb_clear(half);
+    acb_poly_clear(val);
 }
 
 
-/* TODO once we have more support for extracting the results: Apéry */
+void
+fundamental_matrix(const char * dop_str,
+                   const acb_holonomic_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_poly(Dop, Pol);  /* should be Ore poly */
+
+    GR_TMP_INIT(dop, Dop);
+
+    status |= gr_ctx_set_gen_name(Pol, "z");
+    status |= gr_ctx_set_gen_name(Dop, "Tz");
+    status |= gr_set_str(dop, dop_str, Dop);
+    status |= gr_println(dop, Dop);
+    GR_MUST_SUCCEED(status);
+
+    slong dop_order = gr_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);
+
+    acb_holonomic_fundamental_matrix(mat, dop, Dop, expos, pt, 1, 0, 8, prec);
+
+    flint_printf("%{acb_mat}\n", mat);
+
+    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_holonomic_shift_struct shift[1] = {{ .n = 0, .mult = 3 }};
+    acb_holonomic_group_struct grp[1] = {{ .nshifts = 1, .shifts = shift }};
+    acb_init(grp->expo);
+    acb_zero(grp->expo);
+    acb_holonomic_exponents_struct expos[1] = {{ .len = 1, .grps = 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->expo);
+}
+
+
+void
+multiple_shifts(void)
+{
+    /* const char * dop = "Tz^4 - 4*Tz^3 + 3*Tz^2 - z"; */
+    const char * dop = "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2";
+
+    acb_holonomic_shift_struct shift[3] = {
+        { .n = 0, .mult = 2 },
+        { .n = 1, .mult = 3 },
+        { .n = 3, .mult = 1 },
+    };
+    acb_holonomic_group_struct grp[1] = {{ .nshifts = 3, .shifts = shift }};
+    acb_init(grp->expo);
+    acb_zero(grp->expo);
+    acb_holonomic_exponents_struct expos[1] = {{ .len = 1, .grps = grp }};
+
+    fundamental_matrix(dop, expos, 2.);
+
+    acb_clear(grp->expo);
+}
 
 
 int
 main(void)
 {
-    ordinary();
-    series();
-    bessel_j0();
-    whittaker_m();
+    /* ordinary(); */
+    /* series(); */
+    /* bessel_j0(); */
+    /* whittaker_m(); */
+
+    /* apery(); */
+    multiple_shifts();
 
     flint_cleanup_master();
 }

src/acb_holonomic.h

@@ -7,6 +7,7 @@ typedef mp_limb_signed_t slong;
 /* */
 
 #include "acb_types.h"
+#include "gr_types.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -21,6 +22,36 @@ typedef struct
 }
 acb_holonomic_shift_struct;
 
+
+typedef struct
+{
+    acb_t expo;
+    slong nshifts;
+    acb_holonomic_shift_struct * shifts;
+}
+acb_holonomic_group_struct;
+
+
+typedef struct
+{
+    slong len;
+    acb_holonomic_group_struct * grps;
+}
+acb_holonomic_exponents_struct;
+
+
+typedef enum
+{
+    ACB_HOLONOMIC_BASIS_ECHELON = 0,
+    ACB_HOLONOMIC_BASIS_CASCADE  /* XXX: same as Frobenius? */
+}
+acb_holonomic_basis_t;
+
+
+slong acb_holonomic_group_length(const acb_holonomic_group_struct * grp);
+
+/****************************************************************************/
+
 typedef struct
 {
     acb_mat_t extini;
@@ -70,6 +101,11 @@ _acb_holonomic_sol_add_term(
         const fmpz * binom_n,
         slong nder, slong prec, slong sums_prec);
 
+void _acb_holonomic_sol_value(acb_poly_struct * val, acb_srcptr expo,
+                              const acb_poly_struct * sums, slong nlogs,
+                              acb_srcptr pt, slong nder,
+                              slong nshifts, slong prec);
+
 /****************************************************************************/
 
 /* XXX see if this can be reused by other summation algorithms (probably not as
@@ -80,7 +116,7 @@ _acb_holonomic_sol_add_term(
 typedef struct
 {
     /* operator */
-    acb_poly_struct * dop;
+    acb_poly_struct * dop;  /* XXX FJ: move outside the struct? */
     slong dop_len;
     slong dop_clen;
 
@@ -127,15 +163,44 @@ void acb_holonomic_sum_context_init(
 void acb_holonomic_sum_context_clear(acb_holonomic_sum_context_struct * ctx);
 
 void acb_holonomic_sum_ordinary(acb_holonomic_sum_context_struct * ctx);
-void acb_holonomic_sum_mum(acb_holonomic_sum_context_struct * ctx);
 void acb_holonomic_sum_canonical_basis(acb_holonomic_sum_context_struct * ctx);
+void acb_holonomic_sum_highest(acb_holonomic_sum_context_struct * ctx);
+void acb_holonomic_sum_group(acb_holonomic_sum_context_struct * ctx, const acb_holonomic_group_struct * grp);
+
 
 void _acb_holonomic_sum_precompute(acb_holonomic_sum_context_struct * ctx);
 void _acb_holonomic_sum_reset(acb_holonomic_sum_context_struct * ctx);
 
 void acb_holonomic_sum_divconquer(
         acb_holonomic_sum_context_struct * ctx, slong nterms);
 
+void acb_holonomic_sum_value(acb_poly_struct * val, slong nfrobshifts,
+                             const acb_holonomic_sum_context_struct * ctx,
+                             slong i, slong j, slong prec);
+
+/****************************************************************************/
+
+/* XXX void * data or fix args */
+typedef void (* acb_holonomic_sum_worker_t)(
+        acb_holonomic_sum_context_struct * ctx, slong nterms);
+
+void _acb_holonomic_fundamental_matrix(
+        acb_mat_struct * mat,
+        const acb_poly_struct * dop, slong dop_len,
+        const acb_holonomic_group_struct * groups, slong ngroups,
+        acb_srcptr pts, slong npts,
+        acb_holonomic_basis_t basis,
+        acb_holonomic_sum_worker_t sum_worker,
+        slong nterms, slong prec);
+
+int acb_holonomic_fundamental_matrix(
+        acb_mat_struct * mat,
+        gr_srcptr dop, gr_ctx_t dop_ctx,
+        const acb_holonomic_exponents_struct * expos,
+        acb_srcptr pts, slong npts,
+        acb_holonomic_basis_t basis,
+        slong nterms, slong prec);
+
 /****************************************************************************/
 
 void _acb_holonomic_apply_diffop_basecase_weights(
@@ -147,7 +212,7 @@ void _acb_holonomic_apply_diffop_basecase_weights(
 
 void _acb_holonomic_apply_diffop_basecase_precomp(
         acb_poly_struct * g, slong goff,
-        acb_srcptr weights,
+        acb_srcptr weights, slong weights_nlogs,
         const acb_poly_struct * f, slong foff, slong flen,
         slong nlogs,
         slong start, slong len,

src/acb_holonomic/apply_diffop_basecase.c

@@ -105,7 +105,7 @@ _acb_holonomic_apply_diffop_basecase_weights(
 void
 _acb_holonomic_apply_diffop_basecase_precomp(
         acb_poly_struct * g, slong goff,
-        acb_srcptr weights,
+        acb_srcptr weights, slong weights_nlogs,
         const acb_poly_struct * f, slong foff, slong flen,
         slong nlogs,
         slong start, slong len,
@@ -121,7 +121,7 @@ _acb_holonomic_apply_diffop_basecase_precomp(
             for (slong k = k1; k < nlogs; k++)
             {
                 acb_srcptr src = (f + k)->coeffs + foff;
-                acb_srcptr cofac = weights + p1 * nlogs * flen + (k - k1) * flen;
+                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);
@@ -155,7 +155,7 @@ acb_holonomic_apply_diffop_basecase(
     _acb_poly_vec_fit_length(g, nlogs, len);
     _acb_holonomic_apply_diffop_basecase_precomp(
             g, 0,
-            weights, f, 0, flen,
+            weights, nlogs, f, 0, flen,
             nlogs, start, len, prec);
     _acb_poly_vec_set_length(g, nlogs, len);
     _acb_poly_vec_normalise(g, nlogs);

src/acb_holonomic/exponents.c

@@ -0,0 +1,32 @@
+#include "acb_holonomic.h"
+
+
+slong
+acb_holonomic_group_length(const acb_holonomic_group_struct * grp)
+{
+    slong len = 0;
+
+    for (slong s = 0; s < grp->nshifts; s++)
+        len += grp->shifts[s].mult;
+
+    return len;
+}
+
+
+void
+acb_holonomic_exponents_init(acb_holonomic_exponents_struct * expos)
+{
+    expos->len = 0;
+    expos->grps = NULL;
+}
+
+
+void
+acb_holonomic_exponents_free(acb_holonomic_exponents_struct * expos)
+{
+    flint_free(expos->grps);
+}
+
+
+
+