Skip to content

Instantly share code, notes, and snippets.

@behackl

behackl/2021-06 minimal critical points.ipynb

Last active June 9, 2021 00:04
Show Gist options
  • Save behackl/4b9168b014c5b17d299958a103d44906 to your computer and use it in GitHub Desktop.
Save behackl/4b9168b014c5b17d299958a103d44906 to your computer and use it in GitHub Desktop.
Download ZIP
SageMath code (translated from Maple) for finding minimal critical points for multivariate combinatorial generating functions.
Display the source blob
Display the rendered blob
Raw
2021-06 minimal critical points.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"This worksheet contains a translation of the Maple code from [a recent paper by Steve Melczer and Bruno Salvy](https://www.sciencedirect.com/science/article/abs/pii/S0747717120300079) for the computation of minimal critical points of multivariate combinatorial generating functions.\n",
"\n",
"This worksheet has been authored by Benjamin Hackl, Jesse Selover, and Elaine Wong in the context of the MRC on *Combinatorial Applications of Computational Geometry and Algebraic Topology* which took place from May 30, 2021 to June 5, 2021. Work on including this algorithm (in a more SageMath compatible way) into SageMath itself is underway, see [trac ticket #31908](https://trac.sagemath.org/ticket/31908)."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"########################################\n",
"#\n",
"# Kronecker (A translation of Maple Code)\n",
"#\n",
"# Given: A (regular reduced) sequence of polynomials F.\n",
"# \n",
"# Goal: Achieve a symbolic Kronecker representation of the solutions to the system {f = 0 : f in F}.\n",
"#\n",
"# INPUT:\n",
"# - system: {f = 0 : f in F}\n",
"# - u: a new variable not contained in the system\n",
"# - lambda_: parameter introduced for critical point computation\n",
"# - vs: variables of the system\n",
"# - Optional: A linear form in the given variables and u\n",
"#\n",
"# OUTPUT:\n",
"# - P: a polynomial in u\n",
"# - Qs: a list of polynomials in u\n",
"#\n",
"########################################\n",
"\n",
"def Kronecker(system, u, lambda_, vs, linf):\n",
" # set up everything to compute the groebner basis of the ideals generated\n",
" # by everything contained in system, the specified linear form, and\n",
" # additionally 1 - lambda_*x0, where x0 is a new variable (Rabinowitsch's trick)\n",
" expanded_R = lambda_.parent()\n",
" rabinowitsch_R.<x0> = PolynomialRing(expanded_R, 1)\n",
" rabinowitsch_system = [rabinowitsch_R(f) for f in system] + [rabinowitsch_R(linf), 1-lambda_*x0]\n",
" flat = rabinowitsch_R.flattening_morphism()\n",
" Flat_Ring = flat.codomain()\n",
" n = len(Flat_Ring.gens())\n",
" swapXandU = matrix(Permutation((n-1,n)))\n",
" fvs = [Flat_Ring(v) for v in vs]\n",
" Order = TermOrder(swapXandU)\n",
" Ordered_Ring = Flat_Ring.change_ring(order=Order)\n",
" ovs = [Ordered_Ring(v) for v in fvs]\n",
" unflat = flat.inverse()\n",
" osystem = [Ordered_Ring(flat(f)) for f in rabinowitsch_system]\n",
" I = ideal(osystem)\n",
" big_gb = I.groebner_basis()\n",
" gb = [p for p in I.groebner_basis() if unflat(p).is_constant()] # the ones without x0\n",
" Ps = [p for p in gb if not any([z in ovs for z in p.variables()])]\n",
"\n",
" if len(Ps) != 1:\n",
" raise ValueError(\"The P polynomial could not be found.\")\n",
"\n",
" P = Ps[0]\n",
" squarefree, _ = P.quo_rem(gcd(P, P.derivative(Ordered_Ring(flat(u)))))\n",
" P = squarefree/squarefree.content()\n",
" Pd = P.derivative(Ordered_Ring(flat(u)))\n",
"\n",
" # Find Q_i for each variable\n",
" Qs = []\n",
" for z in ovs:\n",
" eqns = [f for f in gb if z in f.variables()]\n",
" if len(eqns) != 1:\n",
" raise ValueError(\"Linear form was not separating, try again.\")\n",
"\n",
" eq = eqns[0].polynomial(z)\n",
"\n",
" if eq.degree() != 1:\n",
" raise ValueError(\"Linear form was not separating, try again.\")\n",
"\n",
" _, rem = (Pd*eq.roots()[0][0]).quo_rem(P)\n",
" Qs.append(rem)\n",
"\n",
" return P, Qs\n",
"\n",
"\n",
"def NewtonRefine(F, p, digits):\n",
" r\"\"\"\n",
" Finds and returns a root p0 of F to `digits` binary digits of precision.\n",
" Input:\n",
" F: a polynomial\n",
" p: a starting approximation to the root p0. Make sure it is not an interval.\n",
" digits: number of binary digits of precision needed\n",
" \"\"\"\n",
" CF = RealIntervalField(2*digits) # Computation field. TODO pick field better\n",
" x = CF(p)\n",
" prec = CF(2^(-digits))\n",
" MAX_STEPS = 2*digits\n",
" delta = 1\n",
" for _ in range(MAX_STEPS):\n",
" y = x\n",
" delta = -F(y)/F.derivative()(y)\n",
" x = y + delta\n",
" if abs(delta) < prec:\n",
" # Assuming Newton's method is converging, the value of the root p0 is within delta of x\n",
" return CF(x-delta, x+delta)\n",
" print(\"Did not converge\")\n",
" return CF(x-delta, x+delta)\n",
"\n",
"########################################\n",
"#\n",
"# Numerical Kronecker (A translation of Maple Code)\n",
"#\n",
"# Given: The symbolic Kronecker representation encoded by a polynomial in the variable u and a list of polynomials in u.\n",
"# \n",
"# Goal: Achieve the isolating intervals (resp. disks) for the real (resp. complex) roots of the polynomial up to a certain precision.\n",
"#\n",
"# INPUT:\n",
"# - P: polynomial in the variable u\n",
"# - Qs: list of polynomials in the variable u\n",
"#\n",
"# OUTPUT:\n",
"# - precise_roots: list of approximated roots\n",
"# - N: ???\n",
"# - kappa: precision in number of digits\n",
"#\n",
"########################################\n",
"\n",
"def NumericalKronecker(P, Qs):\n",
" # Does not support complex roots yet\n",
"\n",
" O = P.parent()\n",
" x0 = O.gens()[-1]\n",
"\n",
" u = P.variables()[0] # ordered ring\n",
" Pd = P.derivative(u)\n",
" deg = P.polynomial(u).degree()\n",
" norm = sqrt(sum([c^2 for c in P.dict().values()]))\n",
" precision = deg^((deg+2)/2)*norm^(deg-1)\n",
"\n",
" roots = [root[0] for root in P.polynomial(u).roots(RR)] # no multiplicities\n",
"\n",
" precisionByQ = []\n",
" for Q in Qs:\n",
" resultant = P.resultant(Pd*x0-Q, u) # TODO: primitive part?\n",
" r_squarefree, _ = resultant.quo_rem(gcd(resultant, resultant.derivative(x0)))\n",
" res = r_squarefree/r_squarefree.content()\n",
" res_norm = sqrt(sum([c^2 for c in res.dict().values()]))\n",
" precisionByQ.append(deg^((deg+2)/2)*res_norm^(deg-1) + 5)\n",
"\n",
" # Forget base ring, move to univariate polynomial ring in u over a field\n",
" # TODO do this before returning from Kronecker\n",
" R = PolynomialRing(QQ, u)\n",
" P = R(P.polynomial(u).map_coefficients(lambda x: x.constant_coefficient()))\n",
" Pd = R(Pd.polynomial(u).map_coefficients(lambda x: x.constant_coefficient()))\n",
" Qs = [R(Q.polynomial(u).map_coefficients(lambda x: x.constant_coefficient())) for Q in Qs]\n",
"\n",
" # Number of binary digits needed in approximation of Q_j/P'\n",
" digits = ceil(log(max(precisionByQ+[precision]), 2))\n",
"\n",
" # TODO revisit\n",
" approx_Pds = [Pd(p) for p in roots]\n",
" approx_Qs = [[Q(p) for p in roots] for Q in Qs]\n",
"\n",
" Pdm = floor(min(approx_Pds))\n",
" Qm = ceil(max([max(values_of_one_Q) for values_of_one_Q in approx_Qs]))\n",
"\n",
" digits = digits + ceil(log(1/Pdm + 2*Qm/(Pdm^2), 2))\n",
"\n",
" Umax = max([abs(r) for r in roots]+[1])\n",
" Pd_max = sum([abs(coeff)*(Umax+1)^k for (k,coeff) in Pd.dict().items()])\n",
" Q_max = max([sum([abs(coeff)*(Umax+1)^k for (k,coeff) in Q.dict().items()]) for Q in Qs])\n",
" kappa = ceil(log(max(1, Pd_max, Q_max), 2)) + digits\n",
"\n",
" # TODO: investigate precision assumptions\n",
" precise_roots = [NewtonRefine(P, r, kappa) for r in roots]\n",
" N = []\n",
" for root in precise_roots:\n",
" # Calculate the values of the variables at the solution corresponding to this root\n",
" pd_val = Pd(root)\n",
" N.append([Q(root)/pd_val for Q in Qs])\n",
"\n",
" # For each variable, sanitize N so that values within the precision bound are returned equal (including 0s)\n",
" for var_i in range(len(Qs)):\n",
" for root_j in range(len(precise_roots)):\n",
" # TODO actually check that the diameter of these intervals is less than precision\n",
" # Zeros\n",
" if N[root_j][var_i].contains_zero():\n",
" N[root_j][var_i] = 0\n",
" # Equal values\n",
" for root_k in range(root_j+1, len(precise_roots)):\n",
" if N[root_j][var_i].overlaps(N[root_k][var_i]):\n",
" N[root_k][var_i] = N[root_j][var_i]\n",
"\n",
" return precise_roots, N, kappa\n",
"\n",
"\n",
"########################################\n",
"#\n",
"# MinimalCriticalCombinatorial (A translation of Maple Code)\n",
"#\n",
"# Given: A combinatorial multivariate rational function F=G/H admitting a finite number of \n",
"# critical points with non-degenerate minimal critical points.\n",
"# \n",
"# Goal: Find the minimial critical points of F.\n",
"#\n",
"# INPUT:\n",
"# - Two polynomials G and H.\n",
"# - Optional: A linear form in the given variables and u.\n",
"#\n",
"# OUTPUT:\n",
"# - P: A polynomial in u representing the minimal polynomial of the algebraic\n",
"# numbers in dominant asymptotics (OR: set of roots of P corresponding to the\n",
"# minimal critical points)\n",
"# - Qs: Polynomials Q[i] for the variables in the Kronecker representation defined by P.\n",
"#\n",
"# Intermediate Functions:\n",
"# - Kronecker\n",
"# - NumericalKronecker\n",
"#\n",
"########################################\n",
"\n",
"# TODO: can't specify useful linf_constant if we create t and lambda_\n",
"\n",
"def MinimalCriticalCombinatorial(G, H, linf_constant = None):\n",
"\n",
" R = H.parent()\n",
" expanded_R.<t, lambda_, u> = PolynomialRing(R, 3)\n",
" vs = [expanded_R(z) for z in list(H.variables())] + [t, lambda_]\n",
" eH = expanded_R(H)\n",
" system = [expanded_R(z*H.derivative(z)) - lambda_ for z in H.variables()] + [expanded_R(H), expanded_R(H.subs({z: z*t for z in H.variables()}))]\n",
"\n",
" linf = None\n",
" if linf_constant:\n",
" linf = u - linf_constant\n",
" # Find the linear form if not specified (there is a bit of randomness here)\n",
" else:\n",
" flat = expanded_R.flattening_morphism()\n",
" flatsystem = [flat(f) for f in system]\n",
" maxcoeff = max([max(f.coefficients()) for f in flatsystem])\n",
" maxdegree = max([f.degree() for f in flatsystem])\n",
" linf = u - sum([randint(-maxcoeff*maxdegree-31, maxcoeff*maxdegree+31)*z for z in vs])\n",
"\n",
" P, Qs = Kronecker(system, u, lambda_, vs, linf)\n",
" U, N, kappa = NumericalKronecker(P, Qs)\n",
"\n",
" O = P.parent()\n",
" x0 = O.gens()[-1] # TODO: get x0 a different way\n",
" u = P.variables()[0] # ordered ring\n",
"\n",
" # Move P and Qs into the univariate polynomial ring Q[u]\n",
" A = PolynomialRing(QQ, u)\n",
" P = A(P.polynomial(u).map_coefficients(lambda x: x.constant_coefficient()))\n",
" Qs = [A(Q.polynomial(u).map_coefficients(lambda x: x.constant_coefficient())) for Q in Qs]\n",
"\n",
" Qt = Qs[-2] # Qs ordering is H.variables() + [t, lambda_]\n",
" Pd = P.derivative()\n",
"\n",
" # TODO Is using Ns in this way numerically valid? The precision needed depends on the polynomial 1-t\n",
" # exact_criticals = [U[root_i] for root_i in range(len(U)) if (N[root_i][-2]-1).contains_zero()]\n",
"\n",
" # Solutions to Pt are solutions to the system where t is not 1\n",
" one_minus_t = gcd(Pd-Qt, P)\n",
" Pt, _ = P.quo_rem(gcd(one_minus_t, P))\n",
" #Pt = Pt/Pt.content()\n",
" tol = 1/(1+max([abs(coeff) for (_, coeff) in Pt.dict().items()]))\n",
" # TODO replace tol\n",
"\n",
" # We are interested if those solutions have t in (0,1), to detect if solutions with t = 1 are minimal\n",
" R_t = []\n",
" R_exact = []\n",
" for root_i in range(len(U)):\n",
" if abs(Pt(U[root_i])) < tol:\n",
" R_t.append((root_i, N[root_i]))\n",
" else:\n",
" R_exact.append((root_i, N[root_i]))\n",
"\n",
" # Must change to support complex solutions\n",
" pos_R_exact = [(i, row) for i, row in R_exact if all([not coord < 0 for coord in row[:-2]])]\n",
" minimals = []\n",
" for i, candidate in pos_R_exact:\n",
" ts = [row[-2] for row in R_t if all([row[i] == candidate[i] for i in range(len(row)-2)])]\n",
" if any([0 < t and t < 1 for t in ts]):\n",
" continue\n",
" minimals.append((i, candidate))\n",
"\n",
" # Change the equations to only deal with t=1 solutions\n",
" newP = one_minus_t\n",
" newPd = one_minus_t.derivative()\n",
" _, invPd, _ = xgcd(Pd, newP)\n",
" Qs = [(Q*newPd*invPd).quo_rem(newP)[1] for Q in Qs]\n",
" P = newP\n",
" Pd = newPd\n",
"\n",
" # Find annihilating polynomials of coordinates\n",
" # See Melczer and Salvy paper Corollary 53 for the derivation of this bound\n",
" modBoundPerVariable = []\n",
" precisionPerVariable = []\n",
" annPolys = []\n",
"\n",
" # Temporarily move back to the ring O so that x0 is there\n",
" AnnRing = PolynomialRing(QQ, x0) # May need to be QQ(i) eventually\n",
" for var_k in range(len(Qs)):\n",
" Q = Qs[var_k]\n",
" resultant = O(P).resultant(O(Pd)*x0-O(Q), u) # TODO: primitive part?\n",
" resultant = AnnRing(resultant.polynomial(x0))\n",
" r_squarefree, _ = resultant.quo_rem(gcd(resultant, resultant.derivative()))\n",
" res = r_squarefree #/r_squarefree.content()\n",
" annPolys.append(res)\n",
" d = res.degree()\n",
" res_norm = max(res.dict().values())\n",
" # Need to know M_k(T) to this accuracy:\n",
" mod_bound = ((d^2+1)*res_norm^2)^(d^2-1)* (res_norm*d)^(2*d^2)\n",
" modBoundPerVariable.append(mod_bound)\n",
"\n",
" # For this we need to know T to this accuracy\n",
" max_abs_coord_value = max([abs(row[var_k]) for row in N])\n",
" res_max = sum([abs(coeff)*(max_abs_coord_value+1)^k for k, coeff in res.dict().items()])\n",
" precisionPerVariable.append(mod_bound * res_max)\n",
"\n",
" kappa2 = max(kappa, ceil(log(max(precisionPerVariable),2).upper()))\n",
"\n",
" # Start with the minimal positive real critical point\n",
" min_i, min_p = minimals[0]\n",
" torus = [U[min_i]]\n",
"\n",
" n_vars = len(min_p) - 2\n",
" # Form a list of candidate roots that might have the same modulus (up to the precision we know already)\n",
" cand = [(i,row) for i, row in R_exact if i != min_i and all([(abs(row[k])-abs(min_p[k])).contains_zero() for k in range(n_vars)])]\n",
"\n",
" # Refine each candidate point to the new precision\n",
" for i, row in cand:\n",
" precise_coords = [NewtonRefine(annPolys[k], row[k].center(), kappa2) for k in range(n_vars)]\n",
" # Verify annPolys[k](|x|) * annPolys[k](-|x|) is 0\n",
" if all([(annPolys[k](abs(precise_coords[k]))*annPolys[k](-abs(precise_coords[k]))).contains_zero() for k in range(n_vars)]):\n",
" torus.append(U[i])\n",
"\n",
" return P, Qs, torus, kappa2"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Example 1\n",
"$$F(x, y) = \\frac{1}{1 - x - y}$$"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"%display latex"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"R.<x, y> = PolynomialRing(QQ)\n",
"F = 1/(1 - x - y)\n",
"\n",
"P, Qs, P_zeros, _ = MinimalCriticalCombinatorial(F.numerator(), F.denominator())\n",
"\n",
"minimal_critical_points_num = [\n",
" tuple(q_poly(u=v)/P.diff()(u=v) for q_poly in Qs[:-2]) for v in P_zeros\n",
"]\n",
"minimal_critical_points_num"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"minimal_critical_points_alg = []\n",
"for v in P_zeros:\n",
" NF.<U> = NumberField(P, embedding=v.center())\n",
" mcp = tuple(q_poly(u=U)/P.diff()(u=U) for q_poly in Qs[:-2])\n",
" minimal_critical_points_alg.append(mcp)\n",
"\n",
"minimal_critical_points_alg"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing\n",
"\n",
"FFPD = FractionWithFactoredDenominatorRing(R, SR)\n",
"F_as_fraction = FFPD(F.numerator() / F.denominator().factor().unit(), F.denominator().factor())\n",
"F_as_fraction"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"asy, exp_growth, subexp = F_as_fraction.asymptotics(p={x: 1/2, y: 1/2}, alpha=[1, 1, 1], N=1)\n",
"asy"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Example 2\n",
"$$F(x, y) = \\frac{1}{1 - x - y - z(1-x)(1-y)}$$"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"R.<x,y,z> = PolynomialRing(QQ,3)\n",
"F = 1/(1 - x - y - z*(1-x)*(1-y))\n",
"\n",
"P, Qs, P_zeros, _ = MinimalCriticalCombinatorial(F.numerator(), F.denominator())\n",
"\n",
"minimal_critical_points_num = [\n",
" tuple(q_poly(u=v)/P.diff()(u=v) for q_poly in Qs[:-2]) for v in P_zeros\n",
"]\n",
"minimal_critical_points_num"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"minimal_critical_points_alg = []\n",
"for v in P_zeros:\n",
" NF.<U> = NumberField(P, embedding=v.center())\n",
" mcp = tuple(q_poly(u=U)/P.diff()(u=U) for q_poly in Qs[:-2])\n",
" minimal_critical_points_alg.append(mcp)\n",
"\n",
"minimal_critical_points_alg"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"tuple(SR(coord) for coord in minimal_critical_points_alg[0])"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"p = {\n",
" x: minimal_critical_points_alg[0][0],\n",
" y: minimal_critical_points_alg[0][1],\n",
" z: minimal_critical_points_alg[0][2]\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing\n",
"\n",
"FFPD = FractionWithFactoredDenominatorRing(R, SR)\n",
"F_as_fraction = FFPD(F.numerator() / F.denominator().factor().unit(), F.denominator().factor())\n",
"F_as_fraction"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"asy, exp_growth, subexp = F_as_fraction.asymptotics(p=p, alpha=[1, 1, 1], N=1)\n",
"asy.simplify_rational()"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Example 3\n",
"$$F(x, y) = \\frac{1}{1 - (x^2 + xy + y^2 + z^2)}$$"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"R.<x,y,z> = PolynomialRing(QQ,3)\n",
"F = 1/(1-(x^2+x*y+y^2+z^2))\n",
"\n",
"P, Qs, P_zeros, _ = MinimalCriticalCombinatorial(F.numerator(), F.denominator())\n",
"\n",
"minimal_critical_points_num = [\n",
" tuple(q_poly(u=v)/P.diff()(u=v) for q_poly in Qs[:-2]) for v in P_zeros\n",
"]\n",
"minimal_critical_points_num"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"P.factor()"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"(P1, P2) = P.factor()\n",
"P1 = P1[0]\n",
"P2 = P2[0]\n",
"P1, P2"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"[P1(u=v) for v in P_zeros]"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"[P2(u=v) for v in P_zeros]"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"minimal_critical_points_alg = []\n",
"for v in P_zeros:\n",
" NF.<U> = NumberField(P2, embedding=v.center())\n",
" mcp = tuple(q_poly(u=U)/P.diff()(u=U) for q_poly in Qs[:-2])\n",
" minimal_critical_points_alg.append(mcp)\n",
"\n",
"minimal_critical_points_alg"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"tuple(SR(coord) for coord in minimal_critical_points_alg[0])"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"p_dicts = [{\n",
" x: mcp[0],\n",
" y: mcp[1],\n",
" z: mcp[2]\n",
"} for mcp in minimal_critical_points_alg]"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing\n",
"\n",
"FFPD = FractionWithFactoredDenominatorRing(R, SR)\n",
"F_as_fraction = FFPD(F.numerator() / F.denominator().factor().unit(), F.denominator().factor())\n",
"F_as_fraction"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"asy, exp_growth, subexp = F_as_fraction.asymptotics(p=p_dicts[0], alpha=[1, 1, 1], N=1)\n",
"asy.simplify_rational()"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"asy, exp_growth, subexp = F_as_fraction.asymptotics(p=p_dicts[1], alpha=[1, 1, 1], N=1)\n",
"asy.simplify_rational()"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"asy, exp_growth, subexp = F_as_fraction.asymptotics(p=p_dicts[2], alpha=[1, 1, 1], N=1)\n",
"asy.simplify_rational()"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
"asy, exp_growth, subexp = F_as_fraction.asymptotics(p=p_dicts[3], alpha=[1, 1, 1], N=1)\n",
"asy.simplify_rational()"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
]
}
],
"metadata": {
"kernelspec": {
"name": "sagemath",
"display_name": "SageMath 9.1",
"language": "sage"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Raw
Dockerfile
FROM sagemath/sagemath:9.3.beta8
COPY --chown=sage:sage . ${HOME}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment