"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.

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 = m.degree(n) - denom_monom.degree(n) + t*(m.degree(z) - denom_monom.degree(z))
    t_degree = m.degree(t) - denom_monom.degree(t)
    print(c/denom_coeff, SR(t)^t_degree * SR(n)^SR(n_degree))

Output:

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.

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())

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)
[898]

Alternatively:

sage: R.ideal([Phi, X-33]).variety()
[{Y: 898, X: 33}]

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

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$.
</div>

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

Over what field?

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

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()) + \
                             TermOrder('wdegrevlex',list(range(1,m+1))))
    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 g.lt().variables())]
    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

Example:

sage: R.<x> = QQ[]
sage: [min(f.degree() 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.

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.

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()
3^3
sage: B.is_field()
True

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)

How about this?

var('x', latex_name=r"\bbm{x}")

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.

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

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.

See https://en.wikipedia.org/wiki/Sign_fu...

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:

plot(sgn(x)*abs(x)^(1/3),(x,-13,13))

What do you expect to get instead?

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).

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

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.

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)

etc.

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

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

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.

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?

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

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 K.degree().

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

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)

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.

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.

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.