rburing's profile - activity

Write R.ideal( instead of Ideal( which tries to guess the ring that you want (and gets it wrong here).

Can you give a reference to a (mathematical) definition?

A simple alternative to 'setting it to be real' is to make a substitution such as this one:

whatever.subs({conjugate(H) : H, real_part(H) : H, imag_part(H) : 0})

It depends on what kind of expressions you want to get rid of.

Nonzero prime ideals in a Dedekind domain – such as $\mathfrak{p} = (1-\sqrt{-2})$ in $O_K$ – are maximal, so the quotient is a field, called the residue field in Sage:

K.<a> = QuadraticField(-2)
OK = K.ring_of_integers()
P = OK.ideal(1-a)
OK_mod_P = P.residue_field()

Then you can factor polynomials over $O_K/\mathfrak{p}$, e.g. as follows:

sage: R.<x> = PolynomialRing(ZZ)
sage: f = x^3 + 5
sage: f_mod_P = f.change_ring(OK_mod_P); f_mod_P
x^3 + 2
sage: f_mod_P.factor()
(x + 2)^3

You can also lift the factorization and double-check the result (e.g. for fun):

sage: S.<y> = PolynomialRing(OK)
sage: factorization_lift = prod(S(g)^e for (g,e) in f_mod_P.factor()); factorization_lift
y^3 + 6*y^2 + 12*y + 8
sage: S(f) - factorization_lift
-6*y^2 - 12*y - 3
sage: S(f) - factorization_lift in S.ideal(P)
True

You're welcome! :)

This works:

S = LieAlgebra(QQ, 10, step=3, names='X', naming='linear')

and this works:

S = LieAlgebra(QQ, 11, step=2, names='X', naming='linear')

but this doesn't work:

S = LieAlgebra(QQ, 11, step=3, names='X', naming='linear')

though it should?

Traceback:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-79-8672234ea66b> in <module>()
----> 1 S = LieAlgebra(QQ, Integer(11), step=Integer(3), names='X', naming='linear')

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/misc/classcall_metaclass.pyx in sage.misc.classcall_metaclass.ClasscallMetaclass.__call__ (build/cythonized/sage/misc/classcall_metaclass.c:1640)()
    327         """
    328         if cls.classcall is not None:
--> 329             return cls.classcall(cls, *args, **kwds)
    330         else:
    331             # Fast version of type.__call__(cls, *args, **kwds)

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/algebras/lie_algebras/lie_algebra.pyc in __classcall_private__(cls, R, arg0, arg1, names, index_set, abelian, nilpotent, category, **kwds)
    439                 from sage.algebras.lie_algebras.nilpotent_lie_algebra import FreeNilpotentLieAlgebra
    440                 del kwds["step"]
--> 441                 return FreeNilpotentLieAlgebra(R, arg1, step, names=names, **kwds)
    442             elif nilpotent:
    443                 raise ValueError("free nilpotent Lie algebras must have a"

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/misc/classcall_metaclass.pyx in sage.misc.classcall_metaclass.ClasscallMetaclass.__call__ (build/cythonized/sage/misc/classcall_metaclass.c:1640)()
    327         """
    328         if cls.classcall is not None:
--> 329             return cls.classcall(cls, *args, **kwds)
    330         else:
    331             # Fast version of type.__call__(cls, *args, **kwds)

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/algebras/lie_algebras/nilpotent_lie_algebra.pyc in __classcall_private__(cls, R, r, s, names, naming, category, **kwds)
    325         return super(FreeNilpotentLieAlgebra, cls).__classcall__(
    326             cls, R,r, s, names=tuple(names), naming=naming,
--> 327             category=category, **kwds)
    328 
    329     def __init__(self, R, r, s, names, naming, category, **kwds):

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/misc/cachefunc.pyx in sage.misc.cachefunc.CachedFunction.__call__ (build/cythonized/sage/misc/cachefunc.c:6301)()
   1003                 return self.cache[k]
   1004         except KeyError:
-> 1005             w = self.f(*args, **kwds)
   1006             self.cache[k] = w
   1007             return w

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/structure/unique_representation.py in __classcall__(cls, *args, **options)
   1025             True
   1026         """
-> 1027         instance = typecall(cls, *args, **options)
   1028         assert isinstance( instance, cls )
   1029         if instance.__class__.__reduce__ == CachedRepresentation.__reduce__:

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/misc/classcall_metaclass.pyx in sage.misc.classcall_metaclass.typecall (build/cythonized/sage/misc/classcall_metaclass.c:2090)()
    494             TypeError: Argument 'cls' has incorrect type (expected type, got classobj)
    495     """
--> 496     return (<PyTypeObject*>type).tp_call(cls, args, kwds)
    497 
    498 # Class for timing::

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/algebras/lie_algebras/nilpotent_lie_algebra.pyc in __init__(self, R, r, s, names, naming, category, **kwds)
    398                     for X_ind, X in basis_by_deg[dx]:
    399                         for Y_ind, Y in basis_by_deg[dy]:
--> 400                             Z = L[X, Y]
    401                             if not Z.is_zero():
    402                                 s_coeff[(X_ind, Y_ind)] = {W_ind: Z[W.leading_support()]

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/algebras/lie_algebras/lie_algebra.pyc in __getitem__(self, x)
    573                 return x[1].ideal(x[0])
    574             # Otherwise it is the bracket of two elements
--> 575             return self(x[0])._bracket_(self(x[1]))
    576         return super(LieAlgebra, self).__getitem__(x)
    577 

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/algebras/lie_algebras/lie_algebra_element.pyx in sage.algebras.lie_algebras.lie_algebra_element.FreeLieAlgebraElement._bracket_ (build/cythonized/sage/algebras/lie_algebras/lie_algebra_element.c:19132)()
   1452                     a, b = mr, ml
   1453                     cr = -cr
-> 1454                 for b_elt, coeff in self.parent()._rewrite_bracket(a, b).iteritems():
   1455                     d[b_elt] = d.get(b_elt, zero) + cl * cr * coeff
   1456                     if d[b_elt] == zero:

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/misc/cachefunc.pyx in sage.misc.cachefunc.CachedMethodCaller.__call__ (build/cythonized/sage/misc/cachefunc.c:10824)()
   1951                 return cache[k]
   1952         except KeyError:
-> 1953             w = self._instance_call(*args, **kwds)
   1954             cache[k] = w
   1955             return w

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/misc/cachefunc.pyx in sage.misc.cachefunc.CachedMethodCaller._instance_call (build/cythonized/sage/misc/cachefunc.c:10280)()
   1827             True
   1828         """
-> 1829         return self.f(self._instance, *args, **kwds)
   1830 
   1831     cdef fix_args_kwds(self, tuple args, dict kwds):

/opt/sagemath-8.4/local/lib/python2.7/site-packages/sage/algebras/lie_algebras/free_lie_algebra.pyc in _rewrite_bracket(self, l, r)
    690             # For a similar reason, we have b >= c.
    691             # Compute the left summand
--> 692             for m, inner_coeff in iteritems(self._rewrite_bracket(l._right, r)):
    693                 if l._left == m:
    694                     continue

AttributeError: 'sage.algebras.lie_algebras.lie_algebra_element.Lie' object has no attribute '_right'

This has been submitted as trac ticket #27018 and subsequently it was fixed.

Yes, this is a bug. Typing Q.multiplicity?? we see that it calls X.multiplicity(Q) and X.multiplicity?? shows that it changes the monomial ordering to negdegrevlex and then tries to use the interface to Singular. However the interface to Singular is broken for univariate polynomial rings with "local" monomial ordering (it always uses the global ordering instead); I reported this as trac ticket #27479.

As a temporary workaround I guess you can add a variable and set it to zero:

sage: A2.<x,y> = AffineSpace(QQ,2)
sage: X = A2.subscheme([x^1789+x,y])
sage: Q = X([0,0])
sage: Q.multiplicity()
1

You can do it like this (also simplifying the notation a bit):

sage: R.<x> = PolynomialRing(GF(3))
sage: S.<i> = R.quotient(x^2 + 1)
sage: sum(p^2 for p in S)
0

We can explain the result (for this choice of $f$) as follows: $p(i) = a+bi$ implies $p(i)^2 = a^2 - b^2 + 2abi$ and the sum runs over all $a$ and $b$, so the "real part" of the sum is $$\sum_{a,b} a^2 - b^2 = \sum_a a^2 - \sum_b b^2 = 0,$$ and the "imaginary part" is $$\sum_{a,b} 2ab = \sum_{a,b} ab + \sum_{a,b} ab = \sum_{a,b} ab + \sum_{a,b} (-a)b = \sum_{a,b} ab - ab = 0.$$

This does not quite do what you're asking, but is close:

sage: from sage.symbolic.constants import Constant
sage: l = Constant('l').expression()
sage: g.diff(z)(l)
D[0](g)(l)

The quotient ring is not so nice because it has zero divisors, but see my answer below.

We can set up the problem as follows:

n = 3
jordan_type = (2,1)

assert sum(jordan_type) == n
import itertools
A_coeff_names = {(i,j) : 'a_{}_{}'.format(i,j) for (i,j) in itertools.combinations(range(n),2)}
R = PolynomialRing(QQ, names=A_coeff_names.values())
A_coeff = {idx : R(A_coeff_names[idx]) for idx in A_coeff_names}
A = Matrix(R, n, A_coeff)
J = block_diagonal_matrix([jordan_block(0, s) for s in jordan_type], subdivide=False)

A necessary condition for $gAg^{-1} = J$ is that the rank of $A^k$ equals the rank of $J^k$ for $1 \leq k < n$.

(And since $A$ has all zeros as eigenvalues, this will give the right Jordan type.)

In particular it is necessary that the $(\operatorname{rank}(J^k)+1)$-minors of $A^k$ vanish; these are polynomial equations upon the coefficients of $A$:

A_coeff_eqns = []
J_power = identity_matrix(QQ, n)
A_power = identity_matrix(R, n)
for k in range(n):
    r = J_power.rank()
    # need (r+1) x (r+1) minors to vanish
    A_coeff_eqns.extend([eqn for eqn in A_power.minors(r+1) if eqn != 0])
    J_power *= J
    A_power *= A

We can solve these equations symbolically, substitute a solution into $A$, change to an exact ring (the fraction field of the polynomial ring in the new variables), calculate the Jordan form (to be sure), and change the names of the variables back to the originals where possible:

A_coeff_symb = map(SR, A_coeff_names.values())
A_coeff_eqns_symb = map(SR, A_coeff_eqns)
for sol in solve(A_coeff_eqns_symb, A_coeff_symb):
    # substitute symbolic solution
    A_new = A.change_ring(SR).apply_map(lambda x: x.subs(sol))
    A_new_vars = list(set(sum([list(x.variables()) for x in A_new.list()], [])))
    # change ring to fraction field of polynomial ring in the new variables
    A_new_ring = PolynomialRing(QQ, names=A_new_vars).fraction_field()
    A_new_rat = A_new.change_ring(A_new_ring)
    if A_new_rat.jordan_form(subdivide=False) != J:
        continue
    # change names of variables back to the originals where possible
    A_final_subs = {A_new_rat[i,j] : A_coeff[(i,j)] for (i,j) in A_coeff if A_new_rat[i,j] in A_new_vars}
    A_final = A_new_rat.apply_map(lambda x: x.subs(A_final_subs))
    print A_final
    print 
    print A_final.jordan_form(subdivide=False, transformation=True)[1]
    print
    print '---'
    print

Output:

[    0 a_0_1 a_0_2]
[    0     0     0]
[    0     0     0]

[         a_0_1              0              0]
[             0              1              1]
[             0              0 a_0_1/(-a_0_2)]

---

[    0     0 a_0_2]
[    0     0 a_1_2]
[    0     0     0]

[a_0_2     0     1]
[a_1_2     0     0]
[    0     1     0]

---

Beware that this assumes Sage can solve the polynomial equations symbolically, which may not be true for large n. To give Sage an easier time (though the solutions may become less pretty), you can replace the list of equations by a Groebner basis for the same ideal:

A_coeff_eqns = R.ideal(A_coeff_eqns).groebner_basis()

I don't know how to fix that code, but as an alternative you can use desolve_system_rk4:

var('y0 y1 t')
ode_rhs = [y1, -2/t*y1 - y0^(3/2)]
points = desolve_system_rk4(ode_rhs,[y0,y1],ics=[0.1,1,0],ivar=t,end_points=10,step=0.1)

To plot y0 (i.e. $x$) against t:

ty0_points = [ [i,j] for i,j,k in points]
list_plot(ty0_points, plotjoined=True)

To plot the curve in the (y0,y1)-plane (together with the vector field):

y0y1_points = [ [j,k] for i,j,k in points]
list_plot(y0y1_points, plotjoined=True, color='red') + sum(arrow2d(p[1:], vector(p[1:]) + vector([eqn_rhs.subs({t : p[0], y0 : p[1], y1: p[2]}) for eqn_rhs in ode_rhs]).normalized()*0.03,arrowsize=1,color='blue') for p in points if p[1] >= 0)

The vector field only makes sense when y0 >= 0 (same for the ODE and its "solution" plotted above). I don't know what the solver does past that point; I wouldn't trust it.

You're welcome. You can vote and accept using the buttons to the left of the answer.

@Kuldeep you're right, I got the definition of algebraic connectivity wrong because Wikipedia was not clear. I updated both Wikipedia and my answer.

This is not my field, but it seems you can do

sage: map(lambda x: x.element, p)
[1, 4, 2, 3]
sage: sigma = Permutation(map(lambda x: x.element, p))
sage: sigma
[1, 4, 2, 3]
sage: sigma(2)
4

As part of the symbolic summation you are trying to access e.g. alpha[i] where i is an indeterminate variable; that doesn't work. Instead you should use the ordinary summation. Also, indices start from 0 in Python and SageMath.

So instead of e.g.

sum(area[i] * sin(2*alpha[i])^2, i, 1, 5)

it should be e.g.

sum(area[i] * sin(2*alpha[i])^2 for i in range(5))

Note also that you can do things like

n = 5
alpha = vector([var('alpha_%d' % (k+1)) for k in range(n)])

to save yourself some typing.

Why do you say it's not correct? When I run G = max((g for g in graphs(8) if g.is_connected() and g.girth() == 5), key=algebraic_connectivity) (or more efficiently G = max((g for g in graphs.nauty_geng('8 -c') if g.girth() == 5), key=algebraic_connectivity)) (or the longer code with analogous modification) then it seems correct.

The command

LinearCode(M)

works well on a Finite Field F, but doesn't work for a finite ring. Is there any way to construct a linear code over a finite ring using sage? In other words, could we construct a module over ZZ/nZZ? I try to understand the trac tickets #6452 and #19345 for constructing a linear code over a finite ring. Unfortunately, I failed to understand.

For convenience we can define

algebraic_connectivity = lambda g: sorted(g.laplacian_matrix().eigenvalues())[1]

Then we can obtain at once

G = max((g for g in graphs(7) if g.girth() == 4), key=algebraic_connectivity)

But this assumes there is only one maximum. If we want to be more on the safe side:

my_graphs = [g for g in graphs(7) if g.girth() == 4]
from collections import defaultdict
acs = defaultdict(list)
for g in my_graphs:
    ac = algebraic_connectivity(g)
    acs[ac].append(g)
max_ac = max(acs.keys())
print 'There are {} graphs of (maximal) algebraic connectivity {}'.format(len(acs[max_ac]), max_ac)
for g in acs[max_ac]:
    show(g)

Output:

There are 1 graphs of (maximal) algebraic connectivity 3

From the documentation of solve:

algorithm - string (default: ‘maxima’);

Sage uses Maxima by default.

Also, more precisely, the source code of solve includes:

if len(s) == 0: # if Maxima's solve gave no solutions, try its to_poly_solve
    try:
        s = m.to_poly_solve(variables)
    except Exception: # if that gives an error, stick with no solutions
        s = []

if len(s) == 0: # if to_poly_solve gave no solutions, try use_grobner
    try:
        s = m.to_poly_solve(variables,'use_grobner=true')
    except Exception: # if that gives an error, stick with no solutions
        s = []

So try Maxima's solve and to_poly_solve, and try the latter with use_grobner=true.

An alternative approach is to write $A = g^{-1}Jg$ with $g$ symbolic, and solve for $g$ such that $A$ is strictly upper triangular ($n(n+1)/2$ polynomial equations for the $n^2$ entries of $g$). Would that suffice, or not? What quotient ring did you want to work in?

The generator is defined correctly, but you're not using it in the right way.

Expression count() creates a new instance of the generator; so the first time you call next() it will give zero.

Instead, you should create it once, and then call next() several times on the resulting object, e.g.:

sage: mycount = count()
sage: [mycount.next() for i in range(18)]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
sage: mycount = count()
sage: [next(mycount) for i in range(18)]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
sage: next(mycount)
18
sage: next(mycount)
19

Can you give a more explicit example?

This forum is for questions about the SageMath mathematics software, not the Sage accounting software.

For questions about the accounting software, see e.g. their support page.

You can change the base ring of the matrix to the ring of rational functions (as you did), factor the characteristic polynomial, and set the factors equal to zero symbolically again (and solve).

R = PolynomialRing(QQ, names='a,b,c,d')
M_rat = M.change_ring(R.fraction_field())
S.<x> = R[]
f = M_rat.characteristic_polynomial(x)
f_factors = f.factor()
linear_factors = [(g,e) for (g,e) in f_factors if g.degree() == 1]
rational_roots = [(solve(SR(g)==0,SR(x),solution_dict=True)[0][SR(x)],e) for (g,e) in linear_factors]
nonlinear_factors = [(g,e) for (g,e) in f_factors if g.degree() > 1]
nonrational_roots = sum([map(lambda eqn: (eqn.rhs(), e), solve(SR(g)==0,SR(x))) for (g,e) in nonlinear_factors], [])
eigenvalues = rational_roots + nonrational_roots

This runs in 33 seconds on my machine. The result is

[(3*(a^2 - b^2 + 2*b*c - c^2)*d + 6*a + 1, 1),
 (-3*(a^2 - b^2 - 2*a*c + c^2)*d + 6*b + 1, 1),
 (-3*(a^2 - 2*a*b + b^2 - c^2)*d + 6*c + 1, 1),
 ((3*a^2 + 2*a*b + 3*b^2 + 2*(a + b)*c + 3*c^2)*d - 2*a - 2*b - 2*c + 1, 2),
 ((3*a^2 - 3*b^2 - 2*b*c - 3*c^2)*d - 2*a + 1, 3),
 (-(3*a^2 - 3*b^2 + 2*a*c + 3*c^2)*d - 2*b + 1, 3),
 (-(3*a^2 + 2*a*b + 3*b^2 - 3*c^2)*d - 2*c + 1, 3),
 ((3*a^2 - 2*a*b + 3*b^2 - 2*(a + b)*c + 3*c^2)*d + 2*a + 2*b + 2*c - 4*sqrt((a^2*b^2 + (a^2 - a*b + b^2)*c^2 - (a^2*b + a*b^2)*c)*d^2 + a^2 - a*b + b^2 - (a + b)*c + c^2 + (a^2*b + a*b^2 + (a + b)*c^2 + (a^2 - 6*a*b + b^2)*c)*d) + 1,
  1),
 ((3*a^2 - 2*a*b + 3*b^2 - 2*(a + b)*c + 3*c^2)*d + 2*a + 2*b + 2*c + 4*sqrt((a^2*b^2 + (a^2 - a*b + b^2)*c^2 - (a^2*b + a*b^2)*c)*d^2 + a^2 - a*b + b^2 - (a + b)*c + c^2 + (a^2*b + a*b^2 + (a + b)*c^2 + (a^2 - 6*a*b + b^2)*c)*d) + 1,
  1)]

Maybe an algorithm like this could be included into Sage in the future.

(The hard part is to detect when this can be applied.)