From cd4d473139f4aedc90cb8ae78341b091dc3b7ddc Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Fri, 20 Dec 2024 23:39:32 +0100 Subject: [PATCH] fixing mysterious power series bug --- src/ODE.jl | 2 +- src/parametrizations.jl | 6 +----- src/power_series_utils.jl | 13 ++++++------- src/wronskian.jl | 2 +- 4 files changed, 9 insertions(+), 14 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index cecf0d08..f8d29318 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -241,7 +241,7 @@ end Input: ode is an ODE over QQ, p is a prime number Output: the reduction mod p, throws an exception if p divides one of the denominators """ -function reduce_ode_mod_p(ode::ODE, p::Int) +function reduce_ode_mod_p(ode::ODE{<:MPolyRingElem{Nemo.QQFieldElem}}, p::Int) new_ring, new_vars = Nemo.polynomial_ring(Nemo.Native.GF(p), map(var_to_str, gens(ode.poly_ring))) new_type = typeof(new_vars[1]) diff --git a/src/parametrizations.jl b/src/parametrizations.jl index 6b873ac5..951b0b62 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -218,11 +218,7 @@ $sat_string $tagged_mqs Monom ordering: $(ord)""" - tagged_mqs_gb = groebner( - tagged_mqs, - ordering = ord, - homogenize = :no, - ) + tagged_mqs_gb = groebner(tagged_mqs, ordering = ord, homogenize = :no) # Relations between tags in K[T] relations_between_tags = filter( poly -> isempty(intersect(vars(poly), vcat(sat_var, orig_vars))), diff --git a/src/power_series_utils.jl b/src/power_series_utils.jl index 3054e503..848b20f8 100644 --- a/src/power_series_utils.jl +++ b/src/power_series_utils.jl @@ -52,7 +52,7 @@ function ps_matrix_inv( power_series_ring = base_ring(parent(M)) result = map(a -> power_series_ring(a), Base.inv(const_term)) cur_prec = 1 - while cur_prec <= prec + while cur_prec < prec result = _matrix_inv_newton_iteration(M, result) cur_prec *= 2 end @@ -161,7 +161,7 @@ function ps_matrix_homlinear_de( (one(parent(A)) + gen(ps_ring) * truncate_matrix(A, cur_prec)) * map(e -> truncate(ps_ring(e), cur_prec), Y0) Z = map(e -> truncate(ps_ring(e), cur_prec), Base.inv(Y0)) - while cur_prec <= prec + while cur_prec < prec matrix_set_precision!(Y, 2 * cur_prec) matrix_set_precision!(Z, cur_prec) Y, Z = _matrix_homlinear_de_newton_iteration(A, Y, Z, cur_prec) @@ -239,8 +239,7 @@ function ps_ode_solution( Svconst = AbstractAlgebra.matrix_space(base_ring(ring), n, 1) eqs = Sv(equations) - x_vars = filter(v -> ("$(v)_dot" in map(string, gens(ring))), gens(ring)) - x_vars = [x for x in x_vars] + x_vars = collect(filter(v -> ("$(v)_dot" in map(var_to_str, gens(ring))), gens(ring))) x_dot_vars = [str_to_var(var_to_str(x) * "_dot", ring) for x in x_vars] Jac_dots = S([derivative(p, xd) for p in equations, xd in x_dot_vars]) @@ -260,14 +259,14 @@ function ps_ode_solution( end cur_prec = 1 - while cur_prec <= prec + while cur_prec < prec @debug "\t Computing power series solution, currently at precision $cur_prec" - new_prec = min(prec + 1, 2 * cur_prec) + new_prec = min(prec, 2 * cur_prec) for i in 1:length(x_vars) set_precision!(solution[x_vars[i]], new_prec) set_precision!(solution[x_dot_vars[i]], new_prec) end - eval_point = [solution[v] for v in gens(ring)] + eval_point = [copy(solution[v]) for v in gens(ring)] map(ps -> set_precision!(ps, 2 * cur_prec), eval_point) eqs_eval = map(p -> evaluate(p, eval_point), eqs) J_eval = map(p -> evaluate(p, eval_point), Jac_xs) diff --git a/src/wronskian.jl b/src/wronskian.jl index 96450961..acf941ec 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -237,7 +237,7 @@ Computes the Wronskians of io_equations @debug "Constructing Wronskians" result = [] - for (i, tlist) in enumerate(termlists) + for tlist in termlists n = length(tlist) evaled = massive_eval(tlist, ps_ext) S = Nemo.matrix_space(F, n, n)