diff --git a/examples/Sage-ACSV-Tests.ipynb b/examples/Sage-ACSV-Tests.ipynb index 2a0814c..9ee549e 100644 --- a/examples/Sage-ACSV-Tests.ipynb +++ b/examples/Sage-ACSV-Tests.ipynb @@ -49,7 +49,7 @@ "source": [ "# Binomial coefficients\n", "F1 = 1/(1-x-y)\n", - "diagonal_asy(F1, as_symbolic=True)" + "diagonal_asy(F1, as_symbolic=True, use_msolve=False)" ] }, { @@ -60,7 +60,7 @@ { "data": { "text/plain": [ - "4*4^n/(pi^1.0*n^1.0)" + "4^n/(sqrt(pi)*sqrt(n))" ] }, "execution_count": 4, @@ -68,6 +68,28 @@ "output_type": "execute_result" } ], + "source": [ + "# Binomial coefficients\n", + "F1 = 1/(1-x-y)\n", + "diagonal_asy(F1, as_symbolic=True, use_msolve=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4*4^n/(pi*n)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "# Lattice paths enumeration sequence (on main diagonal)\n", "F2 = (1+x)*(1+y)/(1-w*x*y*(x+y+1/x+1/y))\n", @@ -76,16 +98,16 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "1.225275868941647?*33.97056274847714?^n/(pi^1.5*n^1.5)" + "1.225275868941647?*33.97056274847714?^n/(pi^(3/2)*n^(3/2))" ] }, - "execution_count": 5, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -98,19 +120,41 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.225275868941647?*33.97056274847714?^n/(pi^(3/2)*n^(3/2))" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Apéry sequence (on main diagonal)\n", + "F3 = 1/(1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1))\n", + "diagonal_asy(F3, as_symbolic=True, use_msolve=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ - "\\(\\displaystyle \\frac{{\\left(12 \\, \\sqrt{2} + 17\\right)}^{n} \\sqrt{\\frac{17}{2} \\, \\sqrt{2} + 12}}{4 \\, \\pi^{1.5} n^{1.5}}\\)" + "\\(\\displaystyle \\frac{{\\left(12 \\, \\sqrt{2} + 17\\right)}^{n} \\sqrt{\\frac{17}{2} \\, \\sqrt{2} + 12}}{4 \\, \\pi^{\\frac{3}{2}} n^{\\frac{3}{2}}}\\)" ], "text/latex": [ - "$\\displaystyle \\frac{{\\left(12 \\, \\sqrt{2} + 17\\right)}^{n} \\sqrt{\\frac{17}{2} \\, \\sqrt{2} + 12}}{4 \\, \\pi^{1.5} n^{1.5}}$" + "$\\displaystyle \\frac{{\\left(12 \\, \\sqrt{2} + 17\\right)}^{n} \\sqrt{\\frac{17}{2} \\, \\sqrt{2} + 12}}{4 \\, \\pi^{\\frac{3}{2}} n^{\\frac{3}{2}}}$" ], "text/plain": [ - "1/4*(12*sqrt(2) + 17)^n*sqrt(17/2*sqrt(2) + 12)/(pi^1.5*n^1.5)" + "1/4*(12*sqrt(2) + 17)^n*sqrt(17/2*sqrt(2) + 12)/(pi^(3/2)*n^(3/2))" ] }, "metadata": {}, @@ -125,16 +169,17 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "0.09677555757474702?*5.884442204019508?^n/(sqrt(pi)*sqrt(n))" + "(0.09677555757474702?*5.88444220401951?^n/(sqrt(pi)*sqrt(n)),\n", + " [[0.548232473567013?, 0.3099773361396642?]])" ] }, - "execution_count": 7, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } @@ -142,12 +187,35 @@ "source": [ "# Example with two critical points with positive coords, neither of which is obviously non-minimal\n", "F4 = -1/(1-(1-x-y)*(20-x-40*y))\n", - "diagonal_asy(F4, as_symbolic=True)" + "diagonal_asy(F4, as_symbolic=True, return_points=True)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.09677555757474702?*5.884442204019508?^n/(sqrt(pi)*sqrt(n)),\n", + " [[0.548232473567013?, 0.3099773361396642?]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Example with two critical points with positive coords, neither of which is obviously non-minimal\n", + "F4 = -1/(1-(1-x-y)*(20-x-40*y))\n", + "diagonal_asy(F4, as_symbolic=True, use_msolve=True, return_points=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, "metadata": {}, "outputs": [ { @@ -169,7 +237,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -191,7 +259,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -200,7 +268,7 @@ "1.015051765128218?*5.828427124746190?^n/(sqrt(pi)*sqrt(n))" ] }, - "execution_count": 10, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -213,7 +281,29 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.015051765128218?*5.828427124746190?^n/(sqrt(pi)*sqrt(n))" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Delannoy generating function\n", + "F7 = 1/(1 - x - y - x*y)\n", + "diagonal_asy(F7, as_symbolic=True, use_msolve=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, "metadata": {}, "outputs": [ { @@ -222,7 +312,7 @@ "1/7*(0.6234898018587335? + 0.7818314824680299?*I)^n + 1/7*(0.6234898018587335? - 0.7818314824680299?*I)^n + 1/7*(-0.2225209339563144? + 0.9749279121818236?*I)^n + 1/7*(-0.2225209339563144? - 0.9749279121818236?*I)^n + 1/7*(-0.9009688679024191? + 0.4338837391175581?*I)^n + 1/7*(-0.9009688679024191? - 0.4338837391175581?*I)^n + 1/7" ] }, - "execution_count": 11, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -235,7 +325,29 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1/7*(0.6234898018587335? + 0.7818314824680299?*I)^n + 1/7*(0.6234898018587335? - 0.7818314824680299?*I)^n + 1/7*(-0.2225209339563144? + 0.9749279121818236?*I)^n + 1/7*(-0.2225209339563144? - 0.9749279121818236?*I)^n + 1/7*(-0.9009688679024191? + 0.4338837391175581?*I)^n + 1/7*(-0.9009688679024191? - 0.4338837391175581?*I)^n + 1/7" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# One variable example with many minimal critical points on same torus\n", + "F8 = 1/(1 - x^7)\n", + "diagonal_asy(F8, as_symbolic=True, use_msolve=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, "metadata": {}, "outputs": [ { @@ -244,7 +356,7 @@ "([(1.285654384750451?, 1, 1, 0.03396226416457560?)], [[0.7778140158516262?]])" ] }, - "execution_count": 12, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -257,7 +369,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 18, "metadata": {}, "outputs": [ { @@ -266,7 +378,7 @@ "0.9430514023983397?*4.518911369262258?^n/(sqrt(pi)*sqrt(n))" ] }, - "execution_count": 13, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -278,12 +390,47 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 19, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "0.9430514023983397?*4.518911369262258?^n/(sqrt(pi)*sqrt(n))" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "F10 = (x^2*y^2 - x*y + 1)/(1-x-y-x*y+x*y^2+x^2*y-x^2*y^3-x^3*y^2)\n", + "diagonal_asy(F10, as_symbolic=True, use_msolve=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.068693934883829?*3.910193204291776?^n/(sqrt(pi)*sqrt(n))" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# Does not finish computing (hangs on computing GB for critical point system)\n", - "#F = (1-x^3*y^6+x^3*y^4+x^2*y^4+x^2*y^3)/(1-x-y+x^2*y^3-x^3*y^3-x^4*y^4-x^3*y^6+x^4*y^6)" + "#Does not finish computing (hangs on computing GB for critical point system)\n", + "from sage_acsv import diagonal_asy\n", + "var('x,y,w,z,n')\n", + "F = (1-x^3*y^6+x^3*y^4+x^2*y^4+x^2*y^3)/(1-x-y+x^2*y^3-x^3*y^3-x^4*y^4-x^3*y^6+x^4*y^6)\n", + "diagonal_asy(F, as_symbolic=True, use_msolve=True)" ] }, { @@ -307,7 +454,7 @@ "lastKernelId": null }, "kernelspec": { - "display_name": "SageMath 9.7", + "display_name": "SageMath 10.2", "language": "sage", "name": "sagemath" }, @@ -321,7 +468,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.1" } }, "nbformat": 4, diff --git a/sage_acsv/asymptotics.py b/sage_acsv/asymptotics.py index 8743e60..ea20543 100644 --- a/sage_acsv/asymptotics.py +++ b/sage_acsv/asymptotics.py @@ -2,17 +2,17 @@ of multivariate rational functions. """ -from sage.all import AA, PolynomialRing, QQ, QQbar, SR, gcd, prod, pi +from sage.all import AA, PolynomialRing, QQ, QQbar, SR, RIF, RealIntervalField, gcd, prod, pi, xgcd -from sage_acsv.kronecker import _kronecker_representation -from sage_acsv.helpers import ACSVException, RationalFunctionReduce, DetHessianWithLog, OutputFormat +from sage_acsv.kronecker import _kronecker_representation, _msolve_kronecker_representation +from sage_acsv.helpers import ACSVException, RationalFunctionReduce, DetHessianWithLog, OutputFormat, CountRootsInInterval from sage_acsv.debug import Timer, acsv_logger MAX_MIN_CRIT_RETRIES = 3 -def diagonal_asy(F, r=None, linear_form=None, return_points=False, output_format=None, as_symbolic=False): +def diagonal_asy(F, r=None, linear_form=None, return_points=False, output_format=None, as_symbolic=False, use_msolve=False): r"""Asymptotics in a given direction r of the multivariate rational function F. INPUT: @@ -137,14 +137,14 @@ def diagonal_asy(F, r=None, linear_form=None, return_points=False, output_format # Initialize variables vs = list(H.variables()) - RR, (t, lambda_, u_) = PolynomialRing(QQ, 't, lambda_, u_').objgens() - expanded_R, _ = PolynomialRing(QQ, len(vs)+3, vs + [t, lambda_, u_]).objgens() + RR, (lambda_, u_) = PolynomialRing(QQ, 'lambda_, u_').objgens() + expanded_R, _ = PolynomialRing(QQ, len(vs)+2, vs + [lambda_, u_]).objgens() vs = [expanded_R(v) for v in vs] - t, lambda_, u_ = expanded_R(t), expanded_R(lambda_), expanded_R(u_) - vsT = vs + [t, lambda_] + lambda_, u_ = expanded_R(lambda_), expanded_R(u_) + vsT = vs + [lambda_] - all_variables = (vs, lambda_, t, u_) + all_variables = (vs, lambda_, u_) d = len(vs) rd = r[-1] @@ -161,7 +161,8 @@ def diagonal_asy(F, r=None, linear_form=None, return_points=False, output_format min_crit_pts = MinimalCriticalCombinatorial( G, H, all_variables, r=r, - linear_form=linear_form + linear_form=linear_form, + use_msolve=use_msolve ) break except Exception as e: @@ -178,7 +179,7 @@ def diagonal_asy(F, r=None, linear_form=None, return_points=False, output_format timer = Timer() # Find det(zH_z Hess) where Hess is the Hessian of z_1...z_n * log(g(z_1, ..., z_n)) - Det = DetHessianWithLog(H, vsT[0:-2], r) + Det = DetHessianWithLog(H, vsT[0:-1], r) # Find exponential growth T = prod([vs[i]**r[i] for i in range(d)]) @@ -245,7 +246,7 @@ def diagonal_asy(F, r=None, linear_form=None, return_points=False, output_format return result -def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None): +def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None, use_msolve=False): r"""Compute minimal critical points of a combinatorial multivariate rational function F=G/H admitting a finite number of critical points. @@ -255,7 +256,7 @@ def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None): * ``G, H`` -- Coprime polynomials with `F = G/H` * ``variables`` -- Tuple of variables of ``G`` and ``H``, followed - by ``lambda_, t, u_`` + by ``lambda_, u_`` * ``r`` -- (Optional) Length `d` vector of positive integers * ``linear_form`` -- (Optional) A linear combination of the input variables that separates the critical point solutions @@ -275,11 +276,11 @@ def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None): Examples:: sage: from sage_acsv import MinimalCriticalCombinatorial - sage: R. = QQ[] + sage: R. = QQ[] sage: pts = MinimalCriticalCombinatorial( ....: 1, ....: 1 - w*(y + x + x^2*y + x*y^2), - ....: ([w, x, y], lambda_, t, u_) + ....: ([w, x, y], lambda_, u_) ....: ) sage: sorted(pts) [[-1/4, -1, -1], [1/4, 1, 1]] @@ -288,8 +289,8 @@ def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None): timer = Timer() # Fetch the variables we need - vs, lambda_, t, u_ = variables - vsT = vs + [t, lambda_] + vs, lambda_, u_ = variables + vsT = vs + [lambda_] # If direction r is not given, default to the diagonal if r is None: @@ -300,44 +301,32 @@ def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None): system = [ vsH[i]*H.derivative(vsH[i]) - r[i]*lambda_ for i in range(len(vsH)) - ] + [H, H.subs({z: z*t for z in vsH})] + ] + [H] # Compute the Kronecker representation of our system timer.checkpoint() - P, Qs = _kronecker_representation(system, u_, vsT, lambda_, linear_form) + + P, Qs = _msolve_kronecker_representation(system, u_, vsT) if use_msolve else \ + _kronecker_representation(system, u_, vsT, lambda_, linear_form) timer.checkpoint("Kronecker") - Qt = Qs[-2] # Qs ordering is H.variables() + [t, lambda_] Pd = P.derivative() + pos_minimals = [] - # Solutions to Pt are solutions to the system where t is not 1 - one_minus_t = gcd(Pd - Qt, P) - Pt, _ = P.quo_rem(one_minus_t) - rts_t_zo = list( - filter( - lambda k: (Qt/Pd).subs(u_=k) > 0 and (Qt/Pd).subs(u_=k) < 1, - Pt.roots(AA, multiplicities=False) - ) - ) - non_min = [[(q/Pd).subs(u_=u) for q in Qs[0:-2]] for u in rts_t_zo] + Rt, t_ = PolynomialRing(AA, ['t_']).objgen() - # Filter the real roots for minimal points with positive coords - pos_minimals = [] - for u in one_minus_t.roots(AA, multiplicities=False): - is_min = True - v = [(q/Pd).subs(u_=u) for q in Qs[0:-2]] - if any([k <= 0 for k in v]): + for u in P.roots(AA, multiplicities=False): + pt = [(q/Pd).subs(u_=u) for q in Qs] + if any([k <= 0 for k in pt[:-1]]): continue - for pt in non_min: - if all([a == b for (a, b) in zip(v, pt)]): - is_min = False - break - if is_min: - pos_minimals.append(u) + ptt = [x * t_ for x in pt] + Htt = Rt(H.subs({v:xt for (v, xt) in zip(vs, ptt)}).radical()/(t_-1)) + if CountRootsInInterval(Htt, t_, 0, 1) == 0: + pos_minimals.append(pt) # Remove non-smooth points and points with zero coordinates (where lambda=0) for i in range(len(pos_minimals)): - x = (Qs[-1]/Pd).subs(u_=pos_minimals[i]) + x = pos_minimals[i][-1] if x == 0: acsv_logger.warning( f"Removing critical point {pos_minimals[i]} because it either " @@ -349,22 +338,23 @@ def MinimalCriticalCombinatorial(G, H, variables, r=None, linear_form=None): if len(pos_minimals) == 0: raise ACSVException("No smooth minimal critical points found.") elif len(pos_minimals) > 1: + print(pos_minimals) raise ACSVException( "More than one minimal point with positive real coordinates found." ) # Find all minimal critical points - minCP = [(q/Pd).subs(u_=pos_minimals[0]) for q in Qs[0:-2]] + minCP = pos_minimals[0][:-1] minimals = [] - for u in one_minus_t.roots(QQbar, multiplicities=False): - v = [(q/Pd).subs(u_=u) for q in Qs[0:-2]] + for u in P.roots(QQbar, multiplicities=False): + v = [(q/Pd).subs(u_=u) for q in Qs[0:-1]] if all([a.abs() == b.abs() for (a, b) in zip(minCP, v)]): minimals.append(u) # Get minimal point coords, and make exact if possible - minimal_coords = [[(q/Pd).subs(u_=u) for q in Qs[0:-2]] for u in minimals] + minimal_coords = [[(q/Pd).subs(u_=u) for q in Qs[0:-1]] for u in minimals] [[a.exactify() for a in b] for b in minimal_coords] timer.checkpoint("Minimal Points") - return [[(q/Pd).subs(u_=u) for q in Qs[0:-2]] for u in minimals] + return [[(q/Pd).subs(u_=u) for q in Qs[0:-1]] for u in minimals] diff --git a/sage_acsv/helpers.py b/sage_acsv/helpers.py index 9df9fb9..b72cfd0 100644 --- a/sage_acsv/helpers.py +++ b/sage_acsv/helpers.py @@ -1,6 +1,8 @@ from enum import Enum -from sage.all import QQ, ceil, gcd, matrix, randint + +from sage.all import QQ, AA, RIF, CIF, PolynomialRing, RealIntervalField, ComplexIntervalField +from sage.all import ceil, gcd, matrix, randint, sqrt, log, factorial class OutputFormat(Enum): @@ -112,6 +114,26 @@ def DetHessianWithLog(H, vs, r): # Return determinant return matrix(Hess).determinant() +def CountRootsInInterval(F, t_, start, end): + R = PolynomialRing(AA, t_) + F_prev = R(F.radical()) + F_cur = F_prev.derivative(t_) + + ss = [F_prev.subs({t_:start})] + es = [F_prev.subs({t_:end})] + while F_cur != 0: + ss.append(F_cur.subs({t_:start})) + es.append(F_cur.subs({t_:end})) + + _, rem = F_prev.quo_rem(F_cur) + F_prev, F_cur = F_cur, -rem + + ss = list(filter(lambda x: x != 0, ss)) + es = list(filter(lambda x: x != 0, es)) + alt_s = len([_ for i in range(1, len(ss)) if ss[i] * ss[i-1] < 0 ]) + alt_e = len([_ for i in range(1, len(es)) if es[i] * es[i-1] < 0 ]) + + return alt_s - alt_e class ACSVException(Exception): def __init__(self, message, retry=False): diff --git a/sage_acsv/kronecker.py b/sage_acsv/kronecker.py index f23e8ef..da738e2 100644 --- a/sage_acsv/kronecker.py +++ b/sage_acsv/kronecker.py @@ -3,6 +3,7 @@ from sage_acsv.helpers import ACSVException, GenerateLinearForm from sage_acsv.debug import acsv_logger +from sage_acsv.msolve import get_parametrization def _kronecker_representation(system, u_, vs, lambda_=None, linear_form=None): @@ -130,8 +131,38 @@ def _kronecker_representation(system, u_, vs, lambda_=None, linear_form=None): return P, Qs +def _msolve_kronecker_representation(system, u_, vs): + result = get_parametrization(vs, system) + _, nvars, _, msvars, _, param = result[1] + + # msolve may reorder the variables, so order them back + Qparams = param[1][2] + vsExt = [str(v) for v in vs] + # Check if no new variable was created by msolve + # If so, the linear_form used is just zd + # i.e. u_ = zd and Qd = zd * P'(zd) + if nvars == len(vs): + Pdz = [0] + [-c for c in param[1][1][1]] + Qparams.append([[1, Pdz], 1]) + pidx = [msvars.index(v) for v in vsExt] + Qparams = [Qparams[i] for i in pidx] -def kronecker(system, vs, linear_form=None): + R = PolynomialRing(QQ, u_) + u_ = R(u_) + + P_coeffs = param[1][0][1] + P = sum([c * u_**i for (i, c) in enumerate(P_coeffs)]) + + Qs = [] + for Q_param in Qparams: + Q_coeffs = Q_param[0][1] + c_div = Q_param[1] + Q = -sum([c * u_**i for (i, c) in enumerate(Q_coeffs)])/c_div + Qs.append(Q) + + return P, Qs + +def kronecker(system, vs, linear_form=None, use_msolve=False): r"""Computes the Kronecker Representation of a system of polynomials INPUT: @@ -159,4 +190,6 @@ def kronecker(system, vs, linear_form=None): system = [R(f) for f in system] vs = [R(v) for v in vs] u_ = R(u_) + if use_msolve: + return _msolve_kronecker_representation(system, u_, vs) return _kronecker_representation(system, u_, vs, linear_form=linear_form) diff --git a/sage_acsv/msolve.py b/sage_acsv/msolve.py new file mode 100644 index 0000000..7999ca7 --- /dev/null +++ b/sage_acsv/msolve.py @@ -0,0 +1,34 @@ +import os +import tempfile +import subprocess + +from sage.misc.sage_eval import sage_eval +from sage.features.msolve import msolve +from sage_acsv.helpers import ACSVException + +def get_parametrization(vs, system): + filename = msolve().absolute_filename() + msolve_in = tempfile.NamedTemporaryFile(mode='w', + encoding='ascii', delete=False) + command = [filename, "-f", msolve_in.name, "-P", "2"] + + system = list(str(e) for e in system) + try: + print(",".join([str(v) for v in vs]), file=msolve_in) + print(0, file=msolve_in) + print(*(pol.replace(" ", "") for pol in system), + sep=',\n', file=msolve_in) + msolve_in.close() + msolve_out = subprocess.run(command, capture_output=True, text=True) + finally: + os.unlink(msolve_in.name) + + msolve_out.check_returncode() + + result = msolve_out.stdout + result = sage_eval(result[:-2]) + + if result[0] != 0: + raise ACSVException("Issue with msolve parametrization - system does not have finitely many solutions") + + return result \ No newline at end of file