ASKSAGE: Sage Q&A Forum - Latest question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 02 Aug 2019 20:27:39 -0500Find solution to Polynomial Sequences without going through varietyhttp://ask.sagemath.org/question/47357/find-solution-to-polynomial-sequences-without-going-through-variety/Hi everyone,
I'm computing the Groebner basis of an ideal defined over the QQ ring. Once I have this Groebner basis, I would like to obtain a set of values that satisfy the equations in the Groebner basis. I know that the full set is going to be the variety of the ideal, but since this object is huge, I might not be interested in finding those values.
In Maple, after computing the Groebner basis, I'm able to call the solve() method on it and even set a maximum number of solutions I want to obtain. As a small example:
P.<x,y,z,t>=PolynomialRing(QQ,4)
I = P.ideal(x*(x-1), y*(y-1), z*(z-1), t - x*y*z)
gb = I.groebner_basis()
Here I could have called `I.variety()` or `gb.variety()` and obtained the same set of solutions:
sage: gbI.variety()
[{y: 0, z: 0, t: 0, x: 0},
{y: 0, z: 0, t: 0, x: 1},
{y: 1, z: 0, t: 0, x: 0},
{y: 1, z: 0, t: 0, x: 1},
{y: 0, z: 1, t: 0, x: 0},
{y: 0, z: 1, t: 0, x: 1},
{y: 1, z: 1, t: 0, x: 0},
{y: 1, z: 1, t: 1, x: 1}]
But I would like to know if it is possible to call a method like (I can do this in Maple):
solve(gb,[max_sol=2])
Such that I can obtain a subset of the variety instead of the whole set. My motivation is that the size of the initial system of polynomials that I have is considerably larger than this example, and finding the feasible solutions on the reduced Groebner basis is more manageable. I might also not be interested in all the elements in the variety. Finally, if I transform the Groebner basis in an ideal itself and try to compute the variety on that object
gbI = ideal(gb)
gbI.variety()
I find the following error
RuntimeError: error in Singular function call 'groebner':
int overflow in hilb 3
error occurred in or before standard.lib::stdhilb line 299: ` intvec hi = hilb( Id(1),1,W );`
expected intvec-expression. type 'help intvec;'
leaving standard.lib::stdhilb
leaving standard.lib::groebner
bernaldeFri, 02 Aug 2019 20:27:39 -0500http://ask.sagemath.org/question/47357/coercion between polynomial rings with block term ordershttp://ask.sagemath.org/question/45553/coercion-between-polynomial-rings-with-block-term-orders/ Hello,
I have the following problem:
I define a polynomial ring
R.<x,y> = PolynomialRing(QQ,2,order='deglex')
then I define a new polynomial ring, with an addition auxiliary variable and a polynomial in this ring
S = PolynomialRing(R.base(),['w1']+[str(v) for v in R.gens() ],order=TermOrder('lex',1)+R.term_order())
S.inject_variables(verbose=False)
f = w1*x^3 + w1*y + x^2 + y^2
then I define another polynomial ring, with another addition auxiliary variable. Coercing to this new ring works fine
T = PolynomialRing(S.base(),['w2']+[str(v) for v in S.gens() ],order=TermOrder('lex',1)+S.term_order())
f2 = T(f)
However the way backwards fails
S(f2)
> (...)
> /home/phil/Applications/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.pyc in _coerce_map_from_(self, P)
791 from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing
792 if is_MPolynomialRing(P) and self.variable_name() in P.variable_names():
--> 793 P_ = P.remove_var(self.variable_name())
794 return self.base_ring()!=P_ and self.base_ring().has_coerce_map_from(P_)
795
> /home/phil/Applications/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/multi_polynomial_ring_base.pyx in sage.rings.polynomial.multi_polynomial_ring_base.MPolynomialRing_base.remove_var (build/cythonized/sage/rings/polynomial/multi_polynomial_ring_base.c:5831)()
299 return PolynomialRing(self.base_ring(), vars, order=self.term_order())
300 except ValueError:
--> 301 raise ValueError("impossible to use the original term order (most likely because it was a block order). Please specify the term order for the subring")
302 else:
303 return PolynomialRing(self.base_ring(), vars, order=order)
> ValueError: impossible to use the original term order (most likely because it was a block order). Please specify the term order for the subring
Any ideas what I could do to prevent this problem? How would I specify the term order (as it was suggested)?
Kind regards,
PhilippPhilipp SchneiderTue, 26 Feb 2019 05:15:45 -0600http://ask.sagemath.org/question/45553/F4 and F5 implementation in Sagemathhttp://ask.sagemath.org/question/44638/f4-and-f5-implementation-in-sagemath/ Hi all,
Is there an implementation in Sagemath for the F4 and F5 algorithms that relates to Grobner bases? I found something about libsingular:slimgb, however I am not 100% sure if it is the right algorithm. Could you please let me know?
Thanks
LujminaTue, 11 Dec 2018 10:36:53 -0600http://ask.sagemath.org/question/44638/How to compute syzygy module of an ideal in a quotient ring?http://ask.sagemath.org/question/37320/how-to-compute-syzygy-module-of-an-ideal-in-a-quotient-ring/I am trying to compute the syzygy module of an ideal generated by two polynomials `<p,q>` modulo `I`, where `I` is another ideal. This means to compute a generating set `[(p1,q1),...,(ps,qs)]` of the module `{(g,h): gp+hq is in I}`. I know that in Sage, we can use singular command to compute syzygy module:
R.<x,y> = PolynomialRing(QQ, order='lex')
f=2*x^2+y
g=y
h=2*f+g
I=ideal(f,g,h)
M = I.syzygy_module();M
[ -2 -1 1]
[ -y 2*x^2 + y 0]
But this does not work with modulo `I`:
R.<x,y> = PolynomialRing(QQ, order='lex')
S.<a,b>=R.quo(x^2+y^2)
I=ideal(a^2,b^2);I
M = I.syzygy_module();M
Ideal (-b^2, b^2) of Quotient of Multivariate Polynomial Ring in x, y over Rational Field by the ideal (x^2 + y^2)
Error in lines 4-4
Traceback (most recent call last):
Is there a way to do that?KittyLTue, 18 Apr 2017 04:34:27 -0500http://ask.sagemath.org/question/37320/GCD of multivariable polynomials and conversion of Laurent polynomials to ordinary polynomialshttp://ask.sagemath.org/question/33260/gcd-of-multivariable-polynomials-and-conversion-of-laurent-polynomials-to-ordinary-polynomials/Let's assume that I am working with some set of Laurent polynomials in $\mathbb{C}[t_1^{\pm1}, \ldots, t_n^{\pm 1}]$.
R = LaurentPolynomialRing(CC, 't', n)
My first question: is there any method which would multiply elements of $R$ by a big enough monomial $t_1^{k_1}\cdot \ldots \cdot t_n^{k_n}$ to get rid of the negative powers and after that changed them to ordinary polynomials? I will consider resulting polynomials up to $t_1^{k_1}\cdot \ldots \cdot t_n^{k_n}$, so $k_i$'s don't matter that much as long as multiplication will result in a Laurent polynomial with non-negative powers (however I would like to avoid shifting exponents by some huge constant). At the moment I am using following workaround which I hope is redundant to some sage method.
# Since I am working over C I would like to get rid of summands coming from numeric errors
def clean_poly(laurent_poly, precision):
decomp = list(laurent_poly)
ret = 0
for (coeff, monom) in decomp:
if(coeff * coeff.conjugate() >= precision):
ret += coeff * monom
return ret
# Here I am finding lowest exponents of laurent_poly corresponding to particular t_i's
def find_shift(laurent_poly, n):
shift = [0 for i in range(0, n)]
for exp in laurent_poly.exponents():
for i in range(0, n):
shift[i] = max(shift[i], -exp[i])
return shift
def change_laurent_poly_to_ordinary(laurent_poly, R_ordinary, n):
laurent_poly = clean_poly(laurent_poly, 0.00000000001)
decomp = list(laurent_poly)
shift = find_shift(laurent_poly, n)
ret = 0
for (coeff, monom) in decomp:
exp = monom.exponents()
exp = exp[0]
summand = 1
for i in range(0, n):
summand *= (R_ordinary.gen(i)) ^ (exp[i] + shift[i])
ret += coeff * summand
return ret
Where $R_{\text{ordinary}} = \mathbb{C}[t_1, \ldots, t_n]$ is passed as argument in the last method.
R_ordinary = PolynomialRing(CC, 's', n)
My second and the most important question is how to fix the last method to return a polynomial which posses a gcd() method just like in this snippet which I've found in the documentation (I can't post a link due to insufficient karma):
sage: R, (x, y) = PolynomialRing(RationalField(), 2, 'xy').objgens()
sage: f = (x^3 + 2*y^2*x)^2
sage: g = x^2*y^2
sage: f.gcd(g)
x^2
Unfortunately at the moment method change_laurent_poly_to_ordinary returns class of type MPolynomial_polydict which doesn't have gcd() method, what results in the following error after attempting to use them:
Traceback (most recent call last):
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 905, in execute
exec compile(block+'\n', '', 'single') in namespace, locals
File "", line 28, in <module>
File "sage/structure/element.pyx", line 420, in sage.structure.element.Element.__getattr__ (/projects/sage/sage-6.10/src/build/cythonized/sage/structure/element.c:4675)
return getattr_from_other_class(self, P._abstract_element_class, name)
File "sage/structure/misc.pyx", line 259, in sage.structure.misc.getattr_from_other_class (/projects/sage/sage-6.10/src/build/cythonized/sage/structure/misc.c:1771)
raise dummy_attribute_error
AttributeError: 'MPolynomial_polydict' object has no attribute 'gcd'
I am working with sage version offered by cloud sagemath.EilenbergFri, 29 Apr 2016 17:21:39 -0500http://ask.sagemath.org/question/33260/Cannot mulyiply polynomial by matrix when ordering is explicitly specifiedhttp://ask.sagemath.org/question/31323/cannot-mulyiply-polynomial-by-matrix-when-ordering-is-explicitly-specified/ Consider the following code:
sage: F = GF(17)
sage: R.<x, y> = PolynomialRing(F)
sage: MS = MatrixSpace(F, 5, 4)
sage: x, y = R.gens()
sage: MS.random_element() * x # Good.
sage: MS.random_element() * y # Good.
# So far, so good. Now watch.
sage: R.<x, y> = PolynomialRing(F, order='lex')
sage: x, y = R.gens()
sage: MS.random_element() * x
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-12-928a7b216caf> in <module>()
----> 1 MS.random_element() * x
sage/structure/element.pyx in sage.structure.element.Matrix.__mul__ (/build/sagemath/src/sage-6.9/src/build/cythonized/sage/structure/element.c:23250)()
sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/build/sagemath/src/sage-6.9/src/build/cythonized/sage/structure/coerce.c:9739)()
TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 5 by 4 dense matrices over Finite Field of size 17' and 'Multivariate Polynomial Ring in x, y over Finite Field of size 17'
Does anybody know if this is intentional? In my opinion it shouldn't happen because lex ordering is already the implicit default. The bug does not seem to occur with degrevlex ordering, which is also weird. The bug also does not occur with square matrices, which is even weirder.
SageMath Version 6.9, Release Date: 2015-10-10d125qFri, 04 Dec 2015 09:43:59 -0600http://ask.sagemath.org/question/31323/derivative of MPolynomial_polydicthttp://ask.sagemath.org/question/11364/derivative-of-mpolynomial_polydict/In a for loop I generate a list of functions, then I take their derivatives with respect to different variables. For example, in the follwoing, `Y[0,0]` is a 'MPolynomial_polydict', and is equal to
x_0_1*x_0_1 + x_0_0*x_1_1 + x_0_2*x_0_2 + x_1_2*x_1_2 + x_0_0*x_2_2 + x_1_1*x_2_2 - 1.00000000000000
and `list_of_variables[0] = x_0_1`. When I take the derivative of the mentioned polynomial with respect to the mentioned variable I should get `2*x_0_1`, but instead I get:
sage: derivative(Y[0,0],list_of_variables[0])
x_0_1
When I try this manually, it gives me the correct answer though:
sage: derivative(x_0_1*x_0_1 + x_0_0*x_1_1 + x_0_2*x_0_2 + x_1_2*x_1_2 + x_0_0*x_2_2 + x_1_1*x_2_2 - 1.00000000000000,x_0_1)
2*x_0_1
So, for some reason it considers the two x_0_1's in the first term of my polynomial as independent variables when I construct them automatically. Any ideas on why this problem arises, and how to fix this?
**Edit:**
This is how I build my vaiables:
names = [ [[] for i in range(n)] for j in range(n) ]
for i in range(n):
for j in range(i+1):
names[i][j] = (SR('x_' + str(j) + '_' + str(i)))
names[j][i] = (SR('x_' + str(j) + '_' + str(i)))
R = PolynomialRing(CC,names[i][j] for i in range(n) for j in range(n)])
R.gens()
R.inject_variables()
**More edit:**
Here is a the complete code:
#########################################################################
# Build variables, and the matrix corresponding to it
#########################################################################
def build_variables(n):
names = [ [[] for i in range(n)] for j in range(n) ]
for i in range(n):
for j in range(i+1):
names[i][j] = (SR('x_' + str(j) + '_' + str(i)))
names[j][i] = (SR('x_' + str(j) + '_' + str(i)))
return(names)
#########################################################################
# Define the function f that maps a matrix to the coefficients of its characteristic polynomial
#########################################################################
def CharPoly(Mat):
X = matrix(Mat)
n = X.ncols()
C_X = X.characteristic_polynomial()
Y = []
for i in range(n):
Y.append(C_X[i])
return(Y)
############################################################################
# This solves that lambda problem
############################################################################
def lambda_siep(G,L,iter=100,epsilon = .1):
# G is any graph on n vertices
# L is the list of n desired distinct eigenvalues
# m is the number of itterations of the Newton's method
# epsilon: the off-diagonal entries will be equal to epsilon
n = G.order()
my_variables = build_variables(n)
R = PolynomialRing(CC,[my_variables[i][j] for i in range(n) for j in range(n)])
R.gens()
R.inject_variables()
X = [ [ R.gens()[n*i+j] for i in range(n) ] for j in range(n) ]
Y = matrix(CharPoly(X)) - matrix(CharPoly(diagonal_matrix(L)))
J = matrix(R,n)
for i in range(n):
for j in range(n):
J[i,j] = derivative(Y[0][i],my_variables[j][j])
B = diagonal_matrix(L) + epsilon * G.adjacency_matrix()
count = 0
while count < iter:
T = [ B[i,j] for i in range(n) for j in range(n)]
C = (J(T)).solve_right(Y(T).transpose())
LC = list(C)
B = B - diagonal_matrix([LC[i][0] for i in range(n)])
count = count + 1
return(B)k1Mon, 28 Apr 2014 19:27:25 -0500http://ask.sagemath.org/question/11364/