rburing's profile - activity

You're welcome. You can accept an answer using the checkmark button on the left of the answer.

This method is currently available only for the Givaro and NTL implementations of finite fields.

You can specify the implementation using the impl keyword to GF, e.g.:

sage: F.<x> = GF(5^1, impl='givaro')
sage: F.fetch_int(4)
4

In theory, this method could easily be added to all implementations. You might open a trac ticket for it and/or bring it up on the sage-devel mailing list.

1. The method calls GAP's ConjugacyClassesSubgroups to get the conjugacy classes and then GAP's Representative to get a representative of each. This does not have the property you want. For example SymmetricGroup(4).conjugacy_classes_subgroups() includes both $K=\langle (1,3)(2,4)\rangle$ and $L=\langle (3,4), (1,2)(3,4) \rangle$ where a conjugate of $K$ is a subgroup of $L$ but $K$ is not a subgroup of $L$.

   See symmetric group: get back conjugacy class from its generators if you want the actual conjugacy classes as GAP objects; you will be able to call the methods Representative() and AsList() on them, and on each subgroup you can also call IsSubgroup(H).

2. I am not a group theorist but your algorithm seems fine to me.

Edit: Maybe I don't have to tell you this but, caveat: being conjugate in a subgroup is a stronger requirement than being conjugate in the full group (because there are fewer elements to conjugate by).

I'm pretty sure it can't be done with a built-in function, but I have written one (rather, a class) as an answer to a previous (slightly different) question: is it possible to round numbers in symbolic expression.

from sage.symbolic.expression_conversions import ExpressionTreeWalker

class SubstituteNumericalApprox(ExpressionTreeWalker):
    def __init__(self, **kwds):
        self.kwds = kwds

    def pyobject(self, ex, obj):
        if not isinstance(obj, Integer) and hasattr(obj, 'numerical_approx'):
            return obj.numerical_approx(**self.kwds)
        else:
            return obj

Here I added an exception for integers, so it works in your use case:

sage: var('t')
sage: f(t) = (4/27*t^9*log(t)^2 - 32/243*t^9*log(t) + 59/2187*t^9 - 2/21*t^7*log(t)^2)
sage: SubstituteNumericalApprox()(f(t))
0.148148148148148*t^9*log(t)^2 - 0.131687242798354*t^9*log(t) + 0.0269775948788294*t^9 - 0.0952380952380952*t^7*log(t)^2

Or in the definition of f:

sage: f(t) = SubstituteNumericalApprox()(4/27*t^9*log(t)^2 - 32/243*t^9*log(t) + 59/2187*t^9 - 2/21*t^7*log(t)^2)
sage: f
t |--> 0.148148148148148*t^9*log(t)^2 - 0.131687242798354*t^9*log(t) + 0.0269775948788294*t^9 - 0.0952380952380952*t^7*log(t)^2

It could be a nice idea to have such a class included in SageMath.

You want to access the generators of R as a tuple:

sage: R = PolynomialRing(QQ, 3, names='x'); R
Multivariate Polynomial Ring in x0, x1, x2 over Rational Field
sage: x = R.gens(); x
(x0, x1, x2)
sage: x[0]
x0

If you want to use some strange alternative indexing, then you can achieve it with a function.

Are the inequalities always polynomial? Probably you can use cylindrical algebraic decomposition.

Numerically or symbolically?

You're welcome. You can accept the answer by pressing the check mark button under the voting arrows (also on all of your other answered questions).

By the way, J = R.ideal(jacobian(f,R.gens()).list()) and J = R.ideal(f.gradient()) and J = f.jacobian_ideal().

If you do explicit division of polynomials with /, you (reasonably) end up in the fraction field of the polynomial ring. What you want to do instead is division with //, which is division (ignoring the remainder) in the polynomial ring. Since you know the remainder is zero, this will be the actual quotient, as a member of the polynomial ring.

(If you have an element q in the fraction field of R and it has denominator 1, then you can also convert it to an element of the polynomial ring by R(q).)

Polynomial rings have monomial orderings (in this case it is degrevlex by default), and the terms of polynomials are automatically ordered from greatest to smallest (with respect to the monomial ordering), as you can see e.g. with print(f).

You can get the exponents you want like this:

sage: [e for e in f.exponents() if all(e_i % 5 == 0 for e_i in e)]
[(10, 5, 0)]

Or the indices of the respective terms (again, with respect to the monomial ordering):

sage: idc = [k for (k,e) in enumerate(f.exponents()) if all(e_i % 5 == 0 for e_i in e)]; idc
[0]

Or the terms themselves:

sage: sum(prod(t) for (k,t) in enumerate(list(f)) if k in idc)
x^10*y^5

A bit more efficiently:

sage: terms = list(f)
sage: sum(prod(terms[k]) for k in idc)
x^10*y^5

There is no inconsistency or bug, but a misunderstanding: the linear equation for v1 has a single solution only if the coefficient of v1 is nonzero. If the coefficient of v1 is zero, then the equation puts no restriction on v1. (So the linear equation is not always enough to determine v1.)

The coefficient of v1 in the linear equation is di + 3/5*a^3 + 3/5*a^2 + 4/5. Indeed, it vanishes for some of the solutions (to check: do Id.variety(QQbar), substitute each solution into the coefficient, and call the radical_expression() method).

See https://ask.sagemath.org/question/493...

I don't know what's hard about it (can't answer that part), but you can use some of the options that solve offers:

sage: solve((x-1)^.5-(2*x-3)^.5-(3*x-4)^.5==0,x,to_poly_solve=True)
[x == (3/2)]
sage: solve((x-1)^(1/2)-(2*x-3)^(1/2)-(3*x-4)^(1/2)==0,x,algorithm='sympy') # note the 1/2
[x == (3/2)]

The first thing this method does is compute a.minpoly()...

You're welcome :) This approach is not exact but if for an algebraic number z you manage to find an f in this way, then you can check whether it is actually exact e.g. with f(z).is_zero(). (If this is the case, it is not guaranteed that f will be of minimal degree, but it is likely, and you could check it in another way, e.g. by f.is_irreducible().)

In the example in the answer, the vanishing locus $V(I)$ of $I$ is the intersection of the cylinder $x^2+y^2=1$ and the plane $z=0$, so it is just a circle in the plane, i.e. a 1-dimensional algebraic variety (and this dimension is the same as the Krull dimension of the quotient ring). The quotient ring can be identified with the ring of polynomial functions on $V(I)$, i.e. the set of polynomial functions on the circle. On the one hand the circle is 1-dimensional (Krull dimension of the quotient), but on the other hand the functions $1, y, y^2, y^3 \ldots$ on the circle are all linearly independent because there is no equation that could relate them (they are part of an infinite normal basis, so the vector space dimension of the quotient is infinite).

By the dimension of an ideal $J$, I (and Sage and Singular) mean the Krull dimension of the ring modulo $J$. Yes, the Krull dimension of the quotient is 2, but (hence) the vector space dimension of the quotient is $\infty$. The normal basis is a basis of the quotient as a vector space, so it is infinite. A "witness" of the fact that the Krull dimension of the quotient is (at least) 2 would be: a chain of prime ideals of length 2 in the quotient.

You have indeed found a bug in Singular and/or SageMath. According to the documentation of Singular's kbase it should return -1 when the quotient is not finite-dimensional, but here it returns [0]. In SageMath's normal_basis() there was code added in trac ticket #29543 that should filter out the zero, so that ticket could be a natural place to bring this up. In any case, $0$ can never be an element of any vector space basis, and if an ideal is positive-dimensional then the quotient ring has an infinite basis as a vector space, which can be accessed degree-wise as described in the answer.

If you change R to have the lex monomial ordering by PolynomialRing(QQ, order='lex'), and then compute a Groebner basis by list(Id.groebner_basis(algorithm='singular:slimgb')) (avoiding a Singular bug), you get a "triangular" system of equations for your variety. If it were 0-dimensional, then you could start with the last equation and get a finite number of solutions, then use those to solve the second-to-last equation and get finitely many solutions, etc., all the way down. Since the ideal is not 0-dimensional, you will instead reach some point where you have infinitely many solutions, e.g. here you get c12 + c14 where neither value has been fixed before, so you can fix one of them and continue, e.g. adding c14-1, c11-1, c4-1 makes it 0-dimensional. (To check this you have to change the monomial ordering back, to avoid the same Singular bug.) Normally one would then call .variety() at this point, but an alternative that happens to work here is to solve symbolically:

sage: I = Id + R.ideal(c14-1, c11-1, c4-1)
sage: solve(list(map(SR,list(I.groebner_basis(algorithm='singular:slimgb')))),list(map(SR,R.gens())))[0]
[a0 == (1/2),
 a1 == (-1/2),
 a10 == (1/3),
 a11 == (1/3),
 a12 == (-1/6),
 a13 == (-1/6),
 a14 == (-1/6),
 a15 == (-1/6),
 a2 == (1/3),
 a3 == (1/3),
 a4 == 0,
 a5 == -1/6*sqrt(3),
 a6 == 1/6*sqrt(3),
 a7 == 1/6*sqrt(3),
 a8 == -1/6*sqrt(3),
 a9 == (1/3),
 b0 == -1/2*(-1)^(1/4),
 b1 == 1/2*(-1)^(1/4),
 b10 == (1/2*I),
 b11 == (1/12),
 b12 == (-1/36),
 b13 == (1/6*I),
 b14 == (1/36),
 b15 == (1/6*I),
 b2 == (-1/3),
 b3 == (2/9*I),
 b4 == (1/9),
 b5 == 0,
 b6 == (1/12),
 b7 == (1/6),
 b8 == (1/2*I),
 b9 == (1/12*I),
 c0 == 1/3*I*(-1)^(1/4),
 c1 == -1/3*I*(-1)^(1/4),
 c10 == (-1/6*I),
 c11 == 1,
 c12 == -1,
 c13 == (-1/6*I),
 c14 == 1,
 c15 == (-1/6*I),
 c2 == (-1/3),
 c3 == (-1/2*I),
 c4 == 1,
 c5 == 0,
 c6 == 1,
 c7 == (1/6),
 c8 == (-1/6*I),
 c9 == (-1/3*I),
 s2 == -sqrt(2),
 s3 == -sqrt(3)]

I'm not sure how minpoly() works but I'm pretty sure the following numerical approach (using PSLQ) is different:

z = AA(sqrt(3)+sqrt(2))
max_degree = 4
from mpmath import mp, findpoly
mp.dps = 15 # decimal places
f = ZZ['x'](findpoly(z.numerical_approx(digits=mp.dps), max_degree))
print(f)
print(f(z.numerical_approx(digits=mp.dps)))
print(f(z).is_zero())
print(f.roots(SR))

Output:

x^4 - 10*x^2 + 1
1.31006316905768e-14
True
[(-sqrt(2*sqrt(6) + 5), 1), (sqrt(2*sqrt(6) + 5), 1), (-sqrt(-2*sqrt(6) + 5), 1), (sqrt(-2*sqrt(6) + 5), 1)]

It can be used to find minimal polynomials with some degree of uncertainty, but likely at a higher speed. I hope it can be useful.

Edit: The uncertainty of z being a root can be eliminated by checking f(z).is_zero(), and the uncertainty of minimality could be eliminated e.g. by checking f.is_irreducible().

The hyperplane trick works, but .variety() runs into this Singular bug because there are more than 30-something variables, and unfortunately .variety() currently does not allow control over the Groebner basis algorithm that it uses internally.

This output should never occur, because Id.normal_basis() should return a list/sequence of monomials which form a vector space basis of the quotient ring, and 0 is not a monomial (it is the "empty sum of monomials").

If Id.dimension() > 0 then the quotient is not finite-dimensional (i.e. the normal basis is infinite), and the convention seems to be for normal_basis() to return the empty list [] in this case (it does not mean that the normal basis is actually empty, it just means we cannot really return the whole thing in a reasonable way). In recent versions of SageMath you can still obtain part of the (infinite) basis, by passing a degree argument to normal_basis; this will return all monomials of the specified degree in the basis:

sage: R.<x,y,z>=PolynomialRing(QQ)
sage: I=R.ideal([x^2+y^2-1, z])
sage: I.dimension()
1
sage: I.normal_basis(1)
[y, x]
sage: I.normal_basis(2)
[y^2, x*y]

In the above output z is absent because $z=0$ in the quotient, and x^2 is absent because $x^2 = 1 - y^2$ in the quotient, etc.

You're welcome! (I updated the answer to make the substitution method more readable.)

You can set all variables to zero:

sage: var("x,y,z")
sage: expr = x*y+z^2+4
sage: expr.subs({v : 0 for v in expr.variables()})
4

If you are working with polynomials, consider using polynomial ring instead:

sage: R.<x,y,z> = PolynomialRing(QQ)
sage: expr = x*y+z^2+4
sage: expr.constant_coefficient()
4

Also you can convert back and forth between symbolic expressions and polynomials:

sage: var("x,y,z")
sage: expr = x*y+z^2+4
sage: R = PolynomialRing(QQ, names='x,y,z')
sage: SR(R(expr).constant_coefficient())
4

It seems Maxima can't handle the symbolic wrapper around AA elements. Try SR(AA(sqrt(3))).numerical_approx() for numerics, or AA(sqrt(3)).radical_expression() for an exact expression. Not all algebraics are expressible in terms of radicals, so this is not a good approach in general. Also solve may return only approximate solutions in more complicated cases. I would instead create an ideal I in a polynomial ring and call I.variety(AA) or I.variety(QQbar).

You can convert algebraic numbers to symbolic expressions using SR(...). Probably you would rather want to define an ideal in a polynomial ring, and compute a Gröbner basis and/or the associated variety (if the system has finitely many solutions). Can you add the system you actually want to solve?

Excess variables are interpreted as parameters, and only solutions that are valid for all values of the parameters are sought. This explains the empty list in the second case.

The value of 0 == 0 is True because 0 is an Integer and equality of Integers is defined that way. To construct the symbolic equation $0=0$, you need the value 0 as a symbolic expression, so that == will construct a symbolic equation:

sage: SR.zero() == SR.zero()
0 == 0

Here SR is the name of the Symbolic Ring, and the zero() method of any ring returns its zero. You can also convert the integer 0 to a symbolic expression:

sage: SR(0) == 0
0 == 0