Implement acb_ode_solution_evaluate

doc/source/acb_ode_solution.rst

@@ -50,3 +50,8 @@ There should rarely be a need to call these directly.
 .. function:: void _acb_ode_solution_extend (acb_ode_solution_t sol, slong nu, acb_poly_t g_nu, slong prec)
 
 	Set the *nu*-th coefficients of the internal power series to the *rho*-derivatives of *g_nu*.
+
+.. function:: void acb_ode_solution_evaluate (acb_t out, acb_ode_solution_t sol, acb_t a, slong prec)
+
+	Evaluate the solution stored in *sol* at the point a.
+	Because any solution may contain logarithms, a must not be zero.

src/acb_ode_solution.c

@@ -21,6 +21,42 @@ void acb_ode_solution_clear (acb_ode_solution_t sol)
 	flint_free(sol->gens);
 }
 
+void acb_ode_solution_evaluate (acb_t out, acb_ode_solution_t sol, acb_t a, slong prec)
+{
+	acb_t l, p, res;
+	slong binom = 1;
+
+	if (acb_is_zero(a))
+	{
+		acb_indeterminate(out);
+		return;
+	}
+
+	acb_init(l);
+	acb_init(p);
+	acb_init(res);
+
+	acb_log(l, a, prec);
+	acb_poly_evaluate(res, sol->gens, a, prec);
+	for (slong i = 1; i < sol->M; i++)
+	{
+		acb_mul(res, res, l, prec);
+
+		acb_poly_evaluate(p, sol->gens + i, a, prec);
+		binom = (binom * (sol->M - i + 1)) / i;
+		acb_mul_si(p, p, binom, prec);
+
+		acb_add(res, res, p, prec);
+	}
+
+	acb_pow(p, a, sol->rho, prec);
+	acb_mul(out, res, p, prec);
+
+	acb_clear(res);
+	acb_clear(l);
+	acb_clear(p);
+}
+
 void _acb_ode_solution_update (acb_ode_solution_t sol, acb_poly_t f, slong prec)
 {
 	acb_struct *F;

tests/CMakeLists.txt

@@ -9,6 +9,7 @@ set(Tests
 	indicial_polynomial
 	solution_extend
 	solution_update
+	solution_eval
 	singleton_frobenius
 	frobenius
 )

tests/solution_eval.c

@@ -0,0 +1,71 @@
+#include "acb_ode.h"
+
+int main ()
+{
+	int return_value = EXIT_SUCCESS;
+	slong prec;
+	fmpz_t binom;
+	acb_t rho, exp, num;
+	acb_t temp1, temp2;
+	acb_ode_solution_t sol;
+	flint_rand_t state;
+
+	/* Initialization */
+	flint_randinit(state);
+
+	fmpz_init(binom);
+	acb_init(rho);
+	acb_init(exp);
+	acb_init(num);
+	acb_init(temp1);
+	acb_init(temp2);
+
+	for (slong iter = 0; iter < 100; iter++)
+	{
+		prec = 30 + n_randint(state, 128);
+
+		acb_randtest(rho, state, prec, 10);
+		acb_ode_solution_init(sol, rho, 1 + n_randint(state, 10), 0);
+
+		for (slong i = 0; i < sol->M; i++)
+			acb_poly_randtest(sol->gens, state, 5, prec, 10);
+
+		acb_randtest(rho, state, prec, 10);
+		acb_ode_solution_evaluate(exp, sol, rho, prec);
+
+		acb_zero(num);
+		for (slong i = 0; i < sol->M; i++)
+		{
+			acb_poly_evaluate(temp1, sol->gens + i, rho, prec);
+
+			fmpz_bin_uiui(binom, sol->M - 1, i);
+			acb_mul_fmpz(temp2, temp1, binom, prec);
+
+			acb_log(temp1, rho, prec);
+			acb_pow_si(temp1, temp1, sol->M - i - 1, prec);
+			acb_mul(temp2, temp2, temp1, prec);
+
+			acb_add(num, num, temp2, prec);
+		}
+		acb_pow(temp1, rho, sol->rho, prec);
+		acb_mul(num, num, temp1, prec);
+
+		acb_ode_solution_clear(sol);
+
+		if (!acb_overlaps(exp, num))
+		{
+			return_value = EXIT_FAILURE;
+			break;
+		}
+	}
+
+	acb_clear(rho);
+	acb_clear(num);
+	acb_clear(exp);
+	acb_clear(temp1);
+	acb_clear(temp2);
+	fmpz_clear(binom);
+	flint_randclear(state);
+	flint_cleanup();
+	return return_value;
+}