Ask Your Question

rburing's profile - activity

2020-05-30 03:58:26 -0500 commented question Are there fast(er) ways to compute inverses of Hermitian matrices?

"How symbolic" are they? What do the entries look like? Polynomials? Rational functions? Worse? You might benefit from changing the ring. Also, inverting is usually best avoided. Are you sure it's necessary? In any case, an example would help.

2020-05-27 06:00:58 -0500 answered a question extract coefficients from products and sums of descending products

If you don't mind I will just focus on the example of desc_prod_recip as in your previous question. It should be possible to generalize this approach. You want to work with a symbolic expression like

sage: n,t = var('n t')
sage: E1 = desc_prod_recip(n,t).expand().simplify(); E1
1/384*n^(-t - 4)*t^8 + 1/96*n^(-t - 4)*t^7 + 1/48*n^(-t - 3)*t^6 - 1/576*n^(-t - 4)*t^6 + 1/48*n^(-t - 3)*t^5 - 1/30*n^(-t - 4)*t^5 + 1/8*n^(-t - 2)*t^4 - 1/16*n^(-t - 3)*t^4 - 5/1152*n^(-t - 4)*t^4 - 1/12*n^(-t - 2)*t^3 - 1/48*n^(-t - 3)*t^3 + 1/32*n^(-t - 4)*t^3 + 1/2*n^(-t - 1)*t^2 - 1/8*n^(-t - 2)*t^2 + 1/24*n^(-t - 3)*t^2 + 1/288*n^(-t - 4)*t^2 - 1/2*n^(-t - 1)*t + 1/12*n^(-t - 2)*t - 1/120*n^(-t - 4)*t + 1/n^t

This is almost a rational function in n (which would be easy to manipulate), except that you have this pesky symbolic exponent t. Suggestion: introduce a new variable z which represents n^t, so that the expression becomes a rational function in n and z.

def desc_prod_recip(x,y,z=None):
    if z is None: z = x^y
    t0 = 1
    t1 = y*(y-1)/(2*x)
    t2 = (y*(y+1)*(y-1)*(3*y-2))/(24*x^2)
    t3 = (y^2*(y-1)^2*(y+1)*(y+2))/(48*x^3)
    t4 = y*(y-1)*(y+1)*(y+2)*(y+3)*(15*y^3-15*y^2-10*y+8)/(5760*x^4)
    return (z^(-1))*(t0 + t1 + t2 + t3 + t4)

To simplify our job, use the fraction field of a polynomial ring rather than the symbolic ring:

sage: R.<n,t,z> = PolynomialRing(QQ)
sage: E1 = desc_prod_recip(n,t,z); E1
(34560*t^8 + 276480*n*t^6 + 138240*t^7 + 1658880*n^2*t^4 + 276480*n*t^5 - 23040*t^6 + 6635520*n^3*t^2 - 1105920*n^2*t^3 - 829440*n*t^4 - 442368*t^5 + 13271040*n^4 - 6635520*n^3*t - 1658880*n^2*t^2 - 276480*n*t^3 - 57600*t^4 + 1105920*n^2*t + 552960*n*t^2 + 414720*t^3 + 46080*t^2 - 110592*t)/(13271040*n^4*z)

To get the (actual) exponents of n, get the exponents of monomials in n and z and weigh them appropriately:

denom = E1.denominator()
denom_coeff = denom.coefficients()[0]
denom_monom = denom.monomials()[0]
for (c,m) in E1.numerator():
    n_degree = - + t*( -
    t_degree = -
    print(c/denom_coeff, SR(t)^t_degree * SR(n)^SR(n_degree))


1/384 n^(-t - 4)*t^8
1/48 n^(-t - 3)*t^6
1/96 n^(-t - 4)*t^7
1/8 n^(-t - 2)*t^4
1/48 n^(-t - 3)*t^5
-1/576 n^(-t - 4)*t^6
1/2 n^(-t - 1)*t^2
-1/12 n^(-t - 2)*t^3
-1/16 n^(-t - 3)*t^4
-1/30 n^(-t - 4)*t^5
1 1/(n^t)
-1/2 n^(-t - 1)*t
-1/8 n^(-t - 2)*t^2
-1/48 n^(-t - 3)*t^3
-5/1152 n^(-t - 4)*t^4
1/12 n^(-t - 2)*t
1/24 n^(-t - 3)*t^2
1/32 n^(-t - 4)*t^3
1/288 n^(-t - 4)*t^2
-1/120 n^(-t - 4)*t

You can also collect terms according to n_degree, etc.

2020-05-27 04:42:53 -0500 answered a question Warning against a deprecated function

You are implicitly treating the expression (x-(sum(vector(M)))/len(vector(M)))^2 as a function.

You should be more explicit about it, e.g. replacing it by lambda x: (x-(sum(vector(M)))/len(vector(M)))^2.

Note that it would be more efficient to compute (sum(vector(M)))/len(vector(M)) once instead of doing it anew every time.

To be explicit:

def varelmat(M) :
    return sum((vector(M)).apply_map(lambda x: (x-(sum(vector(M)))/len(vector(M)))^2))/len(vector(M))

or, alternatively e.g.

def varelmat2(M):
    avg = sum(M.list())/len(M.list())
    return vector(M - avg*ones_matrix(M.base_ring(),M.nrows(), M.ncols())).norm(2)^2 / len(M.list())
2020-05-26 11:35:47 -0500 answered a question How do i solve a 2 variable polynomial over 1 variable

The polynomial Phi can only be evaluated at two inputs (using the function call syntax). The method roots is only available for single variable polynomials. Your Phi(33,Y) is (in SageMath) still a multivariable polynomial (in which X does not appear). To convert it to a single variable polynomial in Y, use the univariate_polynomial method:

sage: f = Phi.subs(X=33).univariate_polynomial(); f
Y^3 + 6150*Y^2 + 5541*Y + 1175
sage: f.parent()
Univariate Polynomial Ring in Y over Finite Field of size 8009
sage: f.roots()
[(898, 1)]
sage: f.roots(multiplicities=False)


sage: R.ideal([Phi, X-33]).variety()
[{Y: 898, X: 33}]
2020-05-19 00:56:54 -0500 received badge  Nice Answer (source)
2020-05-18 12:57:15 -0500 commented question Elliptic Curve over Tower of Finite Field

Why do you want to use change_ring? Why not define E over the big field directly?

2020-05-18 06:43:30 -0500 answered a question Navigation inside a notebook

You can use Markdown and HTML.

For example, if you add a Markdown cell (Cell → Cell Type → Markdown) then you can add headings

# Chapter 1
Some text.
## Subsection 1.1
Some more text.

which will render approximately like this:

Chapter 1

Some text.

Subsection 1.1

Some more text.

but the Jupyter notebook also adds links which you can see by hovering over the headings.

Then you can link to these headings in a Markdown cell using the link syntax:

See [Subsection 1.1](#Subsection-1.1).

By using HTML you can also link to places which are not headings, e.g. if you make a div with an id

<div id="Example-1">
<b>Example 1.</b> $3^2 + 4^2 = 5^2$.

then you can link to it by [Example 1](#Example-1) in a Markdown cell.

2020-05-14 05:23:43 -0500 commented question How can I find the points on the hyperelliptic curve?

Over what field?

2020-05-12 01:05:25 -0500 commented question Latex and SageMath

Well, yes, but there are many ways to typeset matrices in LaTeX. Can you give some example(s) which you want to convert?

2020-05-09 12:34:31 -0500 answered a question Implicitization by symmetric polynomials

As suggested in my comment: you can eliminate $x_1,…,x_t$ from the ideal $$I = \langle y_i - e_i(p_1(x_1,\ldots,x_t),\ldots,p_m(x_1,\ldots,x_t)) : i = 1,\ldots, m \rangle$$ in the ring $\mathbb{Q}[x_1,\ldots,x_t,y_1,\ldots,y_m]$, where $e_i$ are the elementary symmetric polynomials in m variables.

For efficiency reasons we compute the elimination ideal $I \cap \mathbb{Q}[y_1,\ldots,y_m]$ manually, using a Groebner basis $G$ of $I$ with respect to a block ordering where the $x$'s are the greatest and the $y$'s have a weighted degrevlex ordering where $y_k$ has degree $k$ (simulating the degree of $e_k$), so that substituting $y_i = e_i(x_1,\ldots,x_t)$ into the polynomials in $G \cap \mathbb{Q}[y_1,\ldots,y_m]$ yields symmetric dependencies of minimal degree.

def symmetric_dependencies(p):
    R = p[0].parent()
    m = len(p)
    S = PolynomialRing(R.base_ring(),
                       names=list(R.gens()) + ['y{}'.format(k+1) for k in range(m)],
                       order=TermOrder('degrevlex', R.ngens()) + \
    p = [S(f) for f in p]
    X = S.gens()[0:R.ngens()]
    Y = S.gens()[R.ngens():]
    E = SymmetricFunctions(QQ).elementary()
    Es = [E[k+1].expand(m) for k in range(m)]
    I = S.ideal([Y[k] - Es[k](p) for k in range(m)])
    G = I.groebner_basis()
    GY = [g for g in G if all(v in Y for v in]
    E_subs = {Y[k] : Es[k] for k in range(m)}
    dependencies = [f.subs(E_subs) for f in GY]
    #assert all(f(p) == 0 for f in dependencies)
    return dependencies


sage: R.<x> = QQ[]
sage: [min( for f in symmetric_dependencies([x, x^2, x^(2*k)])) for k in range(1,8)]
[4, 5, 8, 11, 12, 12, 12]

It seems to be slightly faster than your function on this example.

2020-05-09 09:22:42 -0500 received badge  Good Answer (source)
2020-05-09 04:53:41 -0500 commented question Implicitization by symmetric polynomials

Can you add your code and a sample case? Maybe naive idea: eliminate $x_1,\ldots,x_t$ from $\langle y_i - e_i(p_1(x_1,\ldots,x_t),\ldots,p_m(x_1,\ldots,x_t)) \rangle$ where $e_i$ are the elementary symmetric polynomials in $m$ variables.

2020-05-08 09:39:34 -0500 answered a question Polynomial ring over the ring of integers modulo 3

The most literal interpretation is to build the quotient ring $(\mathbb{Z}/3\mathbb{Z})[x]/(x^p - x - 1)$:

sage: p = 3
sage: A.<x> = PolynomialRing(Zmod(3))
sage: B.<y> = A.quotient(x^p - x - 1)
sage: B.cardinality().factor()
sage: B.is_field()

Since $3$ is prime you can also replace Zmod(3) by GF(3).

If $x^p - x - 1$ is irreducible (for example for $p=3$) then it is a modulus for the field with $3^p$ elements:

sage: C.<z> = GF(3^p, modulus=x^p - x - 1)
2020-05-04 13:41:23 -0500 received badge  Nice Answer (source)
2020-05-04 11:35:16 -0500 answered a question Escape characters for LaTeX

How about this?

var('x', latex_name=r"\bbm{x}")
2020-04-30 13:55:13 -0500 received badge  Good Answer (source)
2020-04-30 09:41:47 -0500 commented question RuntimeError Groebner basis for a Boolean system

What is happening is: an ideal with generators $x_i^2 + x_i$ etc. plus your generators is created in the polynomial ring over GF(2) with lexicographic ordering, and a Groebner basis of this has to be computed, and for this Singular chooses an algorithm which doesn't work here. I think it's not easy to work around this on the SageMath level, so I reported it on the Singular trac.

2020-04-30 07:41:31 -0500 received badge  Nice Answer (source)
2020-04-30 06:45:05 -0500 answered a question Find expansion of polynomial in an ideal

You can do this with the lift method of p, passing the ideal I as the argument:

sage: R.<x,y> = QQ[]
sage: I = R.ideal([x-1,y-1])
sage: p = x^2 - y^2 - 2*x + 2*y
sage: c = p.lift(I); c
[x - 1, -y + 1]
sage: sum(c_k*p_k for (c_k,p_k) in zip(c,I.gens()))
x^2 - y^2 - 2*x + 2*y
2020-04-25 12:03:01 -0500 received badge  Good Answer (source)
2020-04-25 11:01:27 -0500 commented answer How to correctly plot x^(1/3)

To get the real cube root of $-x$ (where $x>0$) you just have to add a minus sign in front of the real cube root of $x$; that's what this does.

2020-04-25 10:57:42 -0500 commented answer How to correctly plot x^(1/3)
2020-04-25 10:04:02 -0500 received badge  Nice Answer (source)
2020-04-25 09:18:15 -0500 answered a question How to correctly plot x^(1/3)

SageMath sometimes chooses complex cube roots, which explains this behavior.

There's not much you can do about this internal choice. But you can do this:


cube root graph

2020-04-25 07:31:05 -0500 commented question How to correctly plot x^(1/3)

What do you expect to get instead?

2020-04-24 00:09:54 -0500 commented question Unexpected error in a notebook

It says you cannot subtract a graphics object from an integer. Look for subtraction (with -) in the code that gives the error, and check that the values of the variables involved are the ones you expect. Probably you defined a variable, say g, at some point, and you expect it to have this value, but somewhere in between you probably accidentally redefined g to be something else (maybe used it as an index in a loop or something).

2020-04-23 10:01:44 -0500 commented question graphics not save from one cell to the other

It should work. What do you mean by "doesn't work"? What happens exactly?

2020-04-23 05:49:45 -0500 commented answer Choosing the way a function is displayed

Yes, if you need irrational coefficients you should change QQ to something else. If the coefficients are numerical real or complex you could use e.g. RDF or CDF. If the coefficients are algebraic numbers then you could use e.g. QQbar or a NumberField. You could also replace an irrational by an extra variable in A and do the substitution later.

2020-04-23 05:29:11 -0500 answered a question Choosing the way a function is displayed

Yes. I guess you have it as a symbolic expression, like this:

sage: var('p,x,D,w0')
sage: z = -(D*p-p*x - D + w0)/(p-1)

You can consider it as an element of the ring $\mathbb{Q}(p,D,w_0)[x]$, and then get the coefficients:

sage: A = PolynomialRing(QQ, names='p,D,w0').fraction_field()
sage: B = PolynomialRing(A, names='x')
sage: f = B(z); f
(p/(p - 1))*x + (-p*D + D - w0)/(p - 1)
sage: f.coefficients()
[(-p*D + D - w0)/(p - 1), p/(p - 1)]
sage: f.monomial_coefficient(B(x))
p/(p - 1)
sage: f.monomial_coefficient(B(1))
(-p*D + D - w0)/(p - 1)

You might want to convert these back to symbolic expressions again:

sage: SR(f.monomial_coefficient(B(x)))
p/(p - 1)
sage: SR(f.monomial_coefficient(B(1)))
-(D*p - D + w0)/(p - 1)

Or avoid symbolic expressions entirely, by defining z as a polynomial directly:

sage: A.<p,D,w0> = PolynomialRing(QQ)
sage: B.<x> = PolynomialRing(A.fraction_field())
sage: z = -(D*p-p*x - D + w0)/(p-1); z
(p/(p - 1))*x + (-p*D + D - w0)/(p - 1)
sage: z.monomial_coefficient(x)
p/(p - 1)


2020-04-13 08:04:14 -0500 commented question product of functions

A symbolic product in SageMath is currently not implemented, but you could take a logarithm and use a symbolic sum.

2020-04-11 10:45:10 -0500 commented question I'm looking for a program in Sage which allows me to know if two seidel matrices from a graph are similar by a signed permutation matrix

There is already an answer to this question here: how to test if two matrices are similar by a signed permutation matrix

2020-04-11 08:55:57 -0500 commented question Definition of symbolic functions on path algebra

Do I understand correctly that you want a "generic" derivation $d$, not any particular one? (A particular one would be easier.) I'm not an expert on path algebras, but: for an element a, you can get its terms by doing list(a) (that takes care of linearity), and you should manipulate the QuiverPaths somehow to get the "factors", to implement the product rule.

2020-04-11 06:45:13 -0500 commented question Definition of symbolic functions on path algebra

It seems SageMath only has an implementation for derivations over commutative rings. I guess you could do something yourself though. Which derivations do you want to create, and what do you want to do with them?

2020-04-11 03:12:57 -0500 commented question Definition of symbolic functions on path algebra

What do you actually want to achieve? Symbolic functions are probably not the answer.

2020-04-09 15:28:52 -0500 received badge  Nice Answer (source)
2020-04-09 07:41:52 -0500 answered a question Compute dimension of vector subspace

Take the block matrix over $\mathbb{Q}$ defined by

m2 = block_matrix([[m[i,j].matrix() for j in range(M.ncols())] for i in range(M.nrows())])

and divide its rank by

P.S. To get random matrices of lower rank you can use e.g. M.random_element(density=0.1).

2020-04-07 12:46:26 -0500 commented question Compute dimension of vector subspace

How about this matrix? m2 = block_matrix([[m[i,j].matrix() for j in range(M.ncols())] for i in range(M.nrows())], subdivide=False)

2020-04-06 13:47:31 -0500 commented answer Random errors when using Singular via Sage

You're welcome. When using Sage, I find it easiest to use the Sage interface. Mostly for the ease of writing code (which should not be forgotten), but also because the runtime cost compared to using the programs/interfaces directly is often negligible. For example here the most time is probably spent on things like the GTZ algorithm and computing a Gröbner basis (as part of solve). I only use interfaces to other programs directly if it cannot be avoided.

2020-04-06 12:33:54 -0500 answered a question Random errors when using Singular via Sage

The answer is that the output of Singular's primdecGTZ is not deterministic; particularly the ordering can change. You implicitly assumed that it was deterministic, by using indices to access the components. Sometimes you get a component which does not intersect your sphere, so J becomes the unit ideal, and hence solve fails.

2020-04-06 11:58:34 -0500 commented question Random errors when using Singular via Sage

I'll have a look, but in the meantime note that everything you are doing can be done in ordinary PolynomialRings in SageMath (which will use Singular in the background), e.g. using the methods of an ideal primary_decomposition (with argument algorithm='gtz') and variety.