Hello, @Grim! I am not quite sure about Schrödinger's Equation, but I know how to obtain $k=\frac{n\pi}{a}$. I believe this will result a helpful version of your code:

k, a = var('k, a', domain='positive')
psi = function('psi')(x)
eqs = diff(psi, x, 2) + k^(2)*psi == 0
EQS = desolve(eqs, psi, ivar=x)
A, B = var('A, B')
EQS = EQS.subs(_K1=B, _K2=A)
show(LatexExpr(r'\psi(x) = '), EQS)
eq1 = EQS.subs(x=0) == 0  # this is psi(0) = 0
eq2 = EQS.subs(x=a) == 0  # this is psi(a) = 0
show(eq1)
show(eq2)

The first line defines the symbolic variables k and a to be positive, so you can avoid using the assume command. The second line defines psi as the unknown function; the third line defines the ordinary differential equation; and the fourth one solves the differential equation.

Now, things get a little complicated. If you write print(EQS), you will notice that there are a couple of new variables called _K1 and _K2. They are integration variables (byproducts of the resolution of the ODE). However, they are not symbolic variables that you have defined, so you won't be able to use them outside of EQS (for example, to use the solve command or any other). So, I replaced them with A and B in the sixth line. The seventh line (the first show command) is there only to check that I have obtained the same definition of $\psi$ as you in your code (it will show you $\psi(x)=A\cos(kx)+B\sin(kx)$). Then, eq1 and eq2 represent your boundary conditions $\psi(0)=0$ and $\psi(a)=0$, respectively.

At this point, there are many ways in which you can obtain the values of $A$, $B$ and $k$. The first method is manually solving equations eq1 and eq2, which are shown by the last two lines of my code. You will see that

$$A=0$$ and $$A\cos(ak)+B\sin(ak)=0$$

If you replace the first equation in the second, you obtain $B\sin(ax)=0$. Now, since $B>0$ (according to your code), you must have that $\sin(ak)=0$, which implies that $k=\frac{n\pi}{a}$. Finally, in this case, $B$ can be any positive number (regardless of its value, your ODE and boundary conditions are satisfied.)

The second way to solve this is using Sage itself. You have two equations, namely eq1 and eq2, so you will need two variables. You can write

sol = solve([eq1, eq2], A, k)
show(sol)

just after the previous code. You will get something like

$$\left[\left[A = 0, k = \frac{\pi + 2\pi z_{1879}}{a}\right], \left[A = 0, k = \frac{2 \pi z_{1900}}{a}\right]\right]$$

The variables z1879 and z1900 are similar to the _K1 and _K2 that I mentioned above. More precisely, they are integer parameters; just mentally replace them with $m$ and $n$. Then you have two possible solutions. After further analysis, you will notice that these two solutions can be written more simply as $A=0$ and $k=\frac{n\pi}{a}$. Once again, $B$ has no restrictions, besides being positive.

You can also use the Sympy package to solve these equations:

sol = solve([eq1, eq2], A, k, algorithm='sympy')
show(sol)

You will get

$$\left[\left\{A : 0, k : 0\right\}, \left\{A : 0, k : \frac{\pi}{a}\right\}\right]$$

These are just two particular solutions that you can easily generalize using reduction formulas and properties of the trigonometric functions.

What happens if you use solve([eq1, eq2], A, B) or solve([eq1, eq2], B, k)? You will get $[[A=0, B=0]]$ and $[]$, respectively. The first is a correct solution to both equations, but not one that you are looking for. The second is Sage telling you that could not find a solution. If you use SymPy, you will get $\left[\left\{B: -\frac{A}{\tan(ak)}\right\}\right]$, which is a way of admission of defeat, but it is more useful, because it hints you that it is Sage's (actually SymPy's) fault, while the $[]$ before could lead you to think that there is no solution.

I would like to point out a few interesting teachings derived from your question:

1. Do not depend on any software to solve a problem; it's just a tool.
2. SageMath is excellent in solving problems, but it is not perfect. Try different approaches (e.g., as we did when we used SymPy), especially if the problem is complex.
3. You can use a software to help you, but the more of the resolution you can do manually (using ingenuity), the more sure you will be of your result. In the end, if you can do it manually, why bother using Sage? In the end, there are things that are easy for a human mind, but very difficult for a software.

@dsejas -- this might be due to a bug in the way LaTeX is processed on Ask Sage.

It seems some of the backslashes need escaping, ie we need \\ instead of \.

In my experience, not all backslashes need this, but if you do them all you're safe.

It seems to depend on the following character, and seems

• required for \{, \}, \p
• not required for \c, \s, \t

and one could of course prepare a more complete list.

The preview when composing a question or answer behaves differently, making this all the more puzzling.

I'm trying to solve the Schrödinger equation through Sage, for the potential well.

I wrote:

k, a = var('k, a')
psi = function('psi')(x)
eqs = diff(psi, x, 2) + k^(2)*psi == 0
assume(k > 0, a > 0)
EQS = function('EQS')(x)
EQS = desolve(eqs, psi, ivar=x)
A, B, c = var('A, B, c')
result0 = A*cos(k*0) + B*sin(k*0) == 0
result1 = B*sin(k*a) == 0
assume(B > 0)

How can i find the constants A and B through boundary conditions?

My conditions are: Psi(0)=0 and Psi(a)=0.

If i put it in desolve function it return me zero.

So i tried to work with the equations separated but i cannot figure out how to get k=n*pi/a (the expected result).

Duplicate of Ask Sage question 55455.

The difference comes from different methods being available depending on the underlying finite field implementation.

The givaro implementation, available only for prime powers less than 2^16, provides a fetch_int method.

Below we illustrate this and propose a fetch_int function to provide the functionality in other cases.

Illustration

sage: p = 13
sage: F.<a> = GF(p^2)
sage: b = F.fetch_int(150)
sage: c = F.fetch_int(97)
sage: b, c, b + c
(11*a + 7, 7*a + 6, 5*a)


sage: pp = 4091
sage: FF.<aa> = GF(pp^2)
sage: FF.fetch_int(150)

sage: type(F)
<class 'sage.rings.finite_rings.finite_field_givaro.FiniteField_givaro_with_category'>
sage: type(FF)
<class 'sage.rings.finite_rings.finite_field_pari_ffelt.FiniteField_pari_ffelt_with_category'>

sage: FF.<aa> = GF(pp^2, impl='givaro')
Traceback (most recent call last)
...
ValueError: q must be < 2^16



sage: F.fetch_int??
...
    Given an integer `n` return a finite field element in ``self``
    which equals `n` under the condition that :meth:`gen()` is set to
    :meth:`characteristic()`.
...

Define fetch_int as a function:

def fetch_int(F, n):
    return FF(n.digits(base=F.characteristic()))

Then use it as follows:

sage: fetch_int(FF, 12275)
3*a + 2
sage: 3*pp + 2
sage: n
12275

Making all finite fields provide fetch_int regardless of the implementation is the object of

• Sage Trac ticket 31294: Implement fetch_int method for non-givaro finite fields

Here is a piece of code trying to compute the rank of all elliptic curves $E(a,b)$ of the shape $y^2=x^3+ax+b$ for $a,b\in[2010, 2021]$. Sometimes, there it will be "harder" to compute the rank, so the code will not deliver a computed rank. We fill in a dictionary with keys $(a,b)$ and values the corresponding rank for the key when it could be computed, and None otherwise.

R, dic = [2010 .. 2021], {}
for a, b in cartesian_product([R, R]):
    try:
        E = EllipticCurve(QQ, [a, b])
        E.two_descent(second_limit=13, verbose=False)
        r = E.rank(only_use_mwrank=False)
        dic[(a, b)] = r
        print(f'({a}, {b}) -> {r}')
    except Exception:
        dic[(a, b)] = None

To see which cases could not be computed...

[key for key, val in dic.items() if val is None]

... and get an empty list, so in this case there was no None value. But well, in other cases we may have the situation. (Also maybe consider first commenting out the two_descent line, as it may be time consuming; but it may allow finding more curve ranks.)

I see you added an answer. You can accept that answer to mark the question solved.

Another way (less complete than @slelievre's answer but more efficient for just a quick computation) is to compute directly in the free algebra, using the following instead of the product:

def fast_product(a, b):
    res = F.zero()
    a = F(a)
    b = F(b)
    data_a = [(w.to_word(), cf) for w, cf in a] 
    data_b = [(w.to_word(), cf) for w, cf in b]
    for wa, cfa  in data_a:
        for wb, cfb in data_b:
            if len(wb) + len(wa) <= k:
                res += cfa * cfb * F.monomial(wa) * F.monomial(wb)  
    return res

I assume that you mean elliptic curves over ℚ.

I would do two or three things. First, wrap your call to E.rank() in a try/except block to catch cases where an error is raised. Second, create an E.mwrank_curve() object, which you can ask its rank and ask whether the result is proved correct (this is delivered by mwrank):

sage: E = EllipticCurve([0,0,1,-7,6])
sage: Em = E.mwrank_curve()
sage: Em.rank()
3
sage: Em.certain()
True

sage: E = EllipticCurve([0, -1, 1, -929, -10595])
sage: Em = E.mwrank_curve()
sage: Em.rank()
0
sage: Em.certain()
False

In the second example, the rank really is 0 but E has Sha of order 4 and mwrank is not able to tell the difference.

I hope this helps.

Maxima 5.21.1 gives

-1/x-gamma\_incomplete(0,-x)-gamma\_incomplete(-1,-x)

for

integrate(diff((exp(x) - 1)/x, x), x)

which is correct from what I can tell (agrees numerically with the original expression and has the same derivative).

The result isn't as simple as it could be because the integration algorithm is phrased in more general terms, such that the integrand you specified is a special case of some general form. Often that's the most effective way to calculate integrals, since you can cover a lot of special cases with one general form.

The command you want requires the GAP package HAPprime, no longer being developed: its repository at GitHub is archived, and some of the code has been migrated to HAP.

That said, you can still install it by downloading the package's compressed tarball and unpacking it in your ~/.gap/pkg/ folder (after creating it if necessary).

To do that on CoCalc, open a CoCalc terminal and run the following commands (skip the initial $ on each line which represents the prompt):

$ GAP_PKG='~/.gap/pkg'
$ mkdir -p $GAP_PKG
$ cd $GAP_PKG
$ URL=https://files.gap-system.org/gap48/tar.gz/packages
$ TGZ=happrime-0.6.tar.gz
$ wget $URL/$TGZ
$ tar xf $TGZ

Then you can use the desired GAP command, either directly in GAP or in Sage (via gap or, better, libgap) after loading the package HAPprime.

In Sage:

sage: libgap.LoadPackage("happrime")
#I  HAP warning: Set BROWSER_PATH manually if needed.
true
sage: G = libgap.SmallGroup(64, 134)
sage: R = libgap.ResolutionPrimePowerGroupRadical(G, 10)
sage: R
Resolution of length 10 in characteristic 2 for <pc group of size 64 with 6 generators> .
No contracting homotopy available.
A partial contracting homotopy is available.

Note about zero and one:

• Fx.zero() directly provides the zero element of Fx while Fx(0) uses the integer zero and then turns it into an element of Fx.
• Likewise, Fx.one() directly provides the one element of Fx while Fx(1) turns the integer one into an element of Fx.

So one can use E.encode(Fx.one()) and E.encode(Fx.zero()).

One could name the polynomials 0, 1 and a for easy access, e.g.

p0 = Fx.zero()
p1 = Fx.one()
pa = Fx(a)

and then use e.g. E.encode(pa) instead of E.encode(Fx(a)).

Welcome to Ask Sage! Thank you for your question!

Please provide self-contained code that can be copied and pasted.

Currently, we get name errors, as neither A, B, C, Rien are defined.

Calling A the table or matrix when one of the duelists is also A seems error-prone.

Thanks for reporting. This appears to reveal a bug. Fixing it is now tracked at

• Sage Trac ticket 31266: Signed tensor products does not allow tensor products with non-signed modules

Thanks for reporting. This appears to be a bug.

Fixing it is now tracked at

• Sage Trac ticket 31266: Signed tensor products does not allow tensor products with non-signed modules

Note that the code

a = SteenrodAlgebra(2).an_element()
M = CombinatorialFreeModule(GF(2), ['s', 't', 'u'])
s = M.basis()['s']
T = tensor([a, s])

works fine in Sage 8.8 but fails in Sage 8.9.

This change might have to do with the changes in Sage Trac ticket 25603.

Indentation is fixed now.

I have defined and generated a system like below:

def p(*t):
    return SR.var(('p' + '_{}' * len(t)).format(*t))

f = [[[sum(p(i, a, k, b) for b in (1 .. 2)) for a in (1 .. 2)]
      for k in (1 .. 2)] for i in (1 .. 2)]
b = [p(a, b) for a in (1 .. 2) for b in (1 .. 2) for k in (1 .. 2)]
v = vector([f[i][j][k] for i in (0 .. 1) for j in (0 .. 1) for k in (0 .. 1)])
e = [bi - vi ==0 for bi, vi in zip(b, v)]

I would like to know if it is possible to extract the matrix coefficients of the equations in 'e'. The difficulty comes from the definition of the variables 'p()'.

For example: in 'e', I have 8 expressions as below:

e:
[p_1_1 - p_1_1_1_1==0, p_1_1 - p_1_2_1_1==0, p_1_2 - p_1_1_2_1==0, p_1_2 - p_1_2_2_1==0, 
p_2_1 - p_2_1_1_1==0, p_2_1 - p_2_2_1_1==0,  p_2_2 - p_2_1_2_1==0, p_2_2 - p_2_2_2_1==0]

So, my question is if it is possible to extract the coefficients of all the variables p_ _ and p_ _ _ _ considering the lexicography order and taking the p's with 2 entrees as the firsts variables, i.e., with the order: p_1_1, p_1_2, p_2_1, p_2_2, p_1_1_1_1, p_1_1_1_2, ... and so on.

It means, that I have a system with 4 variables with 2 entrees p_ _ and more 16 with 4 entrees p_ _ _ _. I would like to write these coefficients using matrices in this case with 20 variables and 8 equations given in e. Precisely, I want to know the degrees of freedom of the system.

The function individual_duel_scores below returns a dictionary of dictionaries with the duel scores.

One can use it to build a matrix for all the duel scores.

The function:

def individual_duel_scores(duel_winners, duel_scores):
    r"""
    Return a dictionary of dictionaries giving the scores
    of each players against each other player.
    """
    from collections import defaultdict
    players = sorted(set(flatten(duel_winners)))
    scores = defaultdict(lambda: defaultdict(int))
    for ((a, b), _), (x, y) in zip(duel_winners, duel_scores):
        scores[a][b] = x
        scores[b][a] = y
    return scores

Define the players:

sage: A, B, C = 'ABC'

The duel winners and the scores (called l1 and l2 in the question):

sage: duel_winners = [([𝙰, 𝙱], 𝙰), ([𝙰, 𝙲], 𝙲), ([𝙱, 𝙲], 𝙱)]
sage: duel_scores = [[33, 27], [25, 35], [42, 18]]

Individual scores:

sage: scores = individual_duel_scores(duel_winners, duel_scores)

Extract the players from the scores:

sage: players = sorted(scores)
sage: players
['A', 'B', 'C']

Some scores:

sage: scores[A][C]
25

sage: scores[A][A]
0

The scores as a matrix:

sage: matrix(ZZ, 3, [scores[a][b] for a in players for b in players])
[ 0 33 25]
[27  0 42]
[35 18  0]

The code

a = SteenrodAlgebra(2).an_element()
M = CombinatorialFreeModule(GF(2), ['s', 't', 'u'])
s = M.basis()['s']
T = tensor([a, s])

works fine in Sage 8.8 but fails in Sage 8.9.

This might have to do with the changes in Sage Trac ticket 25603.

To format code blocks in questions, answers, or comments, separate them by a blank line above and below, and indent them by four spaces. Cay you edit your comment to do that?

Or you could turn your comment into an answer, since it actually answers the question.

Setting zmin and zmax to get a chosen z-range

Missing functionality

Inspecting the plot method of hyperplanes arrangement reveals this is not implemented in Sage yet.

Searching Ask Sage and sage-support, we find that many 3d plotting functions similarly lack the ability to limit the z-range. Providing this missing functionality is now tracked at

• Sage Trac ticket 31264: Allow setting zmin, zmax in plot3d and other 3d plots

Workaround

Some of the questions on this theme were answered with a suggestion to use implicit_plot3d. For hyperplanes, this works very well.

We can therefore mimic the plot method of hyperplane arrangements and cook up a function which, if given only an x-range and a y-range, will simply call the plot method (which adapts the z-range to the full range of z values corresponding to the x-range and y-range), but if additionally given a z-range, will use implicit plots to limit the z-range.

Here is such a function:

def ha_plot(A, ranges, **opt):
    r"""
    Return a plot of this hyperplane arrangement.

    INPUT:

    - ``A`` -- a hyperplane arrangement

    - ``ranges`` -- a list consisting of an x-range, a y-range,
      and optionally a z-range
    """
    from colorsys import hsv_to_rgb
    if len(ranges) == 2:
        return A.plot(ranges=ranges, **opt)
    elif len(ranges) != 3:
        raise ValueError("ranges should consist in two or three tuples")
    hh = [h.coefficients() for h in A.hyperplanes()]
    N = len(hh)
    HSV_tuples = [(i*1.0/N, 0.8, 0.9) for i in range(N)]
    cols = [hsv_to_rgb(*x) for x in HSV_tuples]
    ff = lambda a, b, c, d: lambda x, y, z: d + a*x + b*y + c*z
    ip = lambda f, col: implicit_plot3d(f, *ranges, color=col, **opt)
    hp = lambda a, b, c, d, col: ip(ff(a, b, c, d), col)
    return sum(hp(a, b, c, d, col) for (d, a, b, c), col in zip(hh, cols))

Usage is as follows.

One can provide only an x-range and a y-range:

sage: ha_plot(A, [(-1, 1), (-1, 1)]).show(aspect_ratio=(1, 1, 0.25))

and get the exact same picture as in the question.

Or one can additionally provide a z-range:

sage: ha_plot(A, [(-1, 1), (-1, 1), (-1, 1)]).show(aspect_ratio=1)

Optional extra arguments are passed to implicit_plot3d.

This for instance allows to set the opacity:

sage: ha_plot(A, [(-1, 1), (-1, 1), (-1, 1)], opacity=0.5).show(aspect_ratio=1)

The function could be improved to accept a choice of colours. Currently it uses the same default as the plot method of hyperplane arrangements.

I am plotting arrangements of hyperplanes. Like this one:

H3.<x, y, z> = HyperplaneArrangements(QQ)
A = H3([(1, 2, 1), 0],
       [(-4, -3, -2), 0],
       [(3, 3, -2), 0],
       [(1, -4, 3), 0],
       [(-2, 2, 1), 0])
A.plot(ranges=[(-1, 1), (-1, 1)], aspect_ratio=(1, 1, 0.25))

The problem is I need the z axis to range from -1 to 1 rather than from -4 to 4.

Is there a way to change that?

Duplicate of Ask Sage question 55355.

Having defined

• the finite field K of size 4 in a
• a polynomial ring R1 over K

one can ask for a random element using the random_element method.

The random polynomial I obtained that way was of degree two:

sage: K.<a> = FiniteField(4)
sage: R1 = PolynomialRing(K, ['z%s' % p for p in range(1, 3)])
sage: R1.random_element()
(a)*z1*z2 + (a + 1)*z2^2 + z1 + (a)*z2

If I only want a random linear polynomial from this ring, how can I get one?

Suppose I take a polynomial from $K[x]$, say $x^2 + 5x$, and another polynomial from $K[y]$, say $y^3$. I want to formally multiply them and get $x^2y^3 + 5xy^3$ as the output.

How can I do that?

Note: K is the finite field of size 4 in a.

My efforts:

sage: K.<a> = FiniteField(4)

sage: R_1 = PolynomialRing(K, ['z%s' % p for p in range(1, 3)])
sage: R_2 = PolynomialRing(K, ['x%s' % p for p in range(1, 3)])

sage: pp = R_1.random_element()
sage: pp
(a + 1)*z1^2 + (a)*z1*z2 + (a + 1)*z2 + (a)

sage: qq = R_2.random_element()
sage: qq
x1*x2 + (a + 1)*x2^2 + x1 + 1

When I do pp*qq I get the following output

unsupported operand parent(s) for *: 'Multivariate Polynomial Ring
in z1, z2 over Finite Field in a of size 2^2' and 'Multivariate Polynomial
Ring in x1, x2 over Finite Field in a of size 2^2'

The .mws file is really a text file, you can open it with a text editor and sort of make sense of it.

Does Sage have a general implementation of the "affine diagonalization algorithm" for n-dimensional vector spaces?

I found some pseudo-code, see page 15 (in section 3 "affine diagonalization" which begins on page 12) of

• Lecture 21: Surfaces in Robust geometric computation

Searching online led me also to this answer from 2008 where someone does it with Maple:

• Maple primes: Normal forms for quadratic functions

but unfortunately, I cannot open the .mws file (I get an error message concerning the version number).

Thank you very much for the help.

Regarding the "\tiny doesn't work" part, you should add a $ before the LaTeX formula and you should either escape the \ with another \ or use a raw string:

sage: cond1 = ['$\\tiny{}{}{}{}$'.format(*w) for w in Ap]
sage: cond1 = [r'$\tiny{}{}{}{}$'.format(*w) for w in Ap]

Regarding the table, it needs a list of lists, not just a list, so use extra brackets; and you want headers for the columns, so they form a row of headers or "header row".

sage: t3 = table([total_score_de_Borda], header_row=cond3)
sage: t3
  $A$   $B$   $C$   $D$
+-----+-----+-----+-----+
  717   691   669   613

I change my question to clarify. In the following code (I have integrated the solution for \tiny):

First define a function:

def All_pref(cand=["A", "B", "C", "D"], code=1) :
    ncand = len(cand)
    Scand = sorted(Set(cand))
    all_pref = Arrangements(Scand, ncand).list()
    all_pref1 = [str(Word(x)) for x in all_pref]
    if code == 1:
        return ncand
    if code == 2:
        return all_pref1

Then use it:

cand = ["A", "B", "C", "D"]
ne1 = [18, 16, 14, 12, 11, 20, 19, 14, 16, 12, 2, 1, 0, 0, 20, 16, 13, 15, 11, 10, 9, 8, 7, 5]
ne = ne1
Ap = All_pref(cand, 2) 
cond1 = ['$\\tiny${}{}{}{}$'.format(*w) for w in Ap]
cond2 = [""] + [r'${}$'.format(*w) for w in cand]
cond3 = [r'${}$'.format(*w) for w in cand]
t1 = table([ne], header_row=cond1)
rank = [[x.find(l) for x in Ap] for l in cand] 
score_de_Borda = [[abs(x.find(l) - len(cand)) for x in Ap] for l in cand]
t2 = table(rank, header_row=cond1, header_column=cond2)
total_score_de_Borda = [add([score_de_Borda[j][i]*ne[i] for i in range(24)]) for j in range(len(cand))]
# t3 = table(total_score_de_Borda, header_column=cond3)
show(cond3,total_score_de_Borda)

everything works fine apart from the header_column=cond3 (without it I have no problem) and as show(cond3, total_score_de_Borda) shows, the two elements have the same length.

Why not open a new question:

• click the "Edit" button under the question
• copy the "follow-up" part
• open a new question
• paste the "follow-up" part there
• remove it from here

The new question can say "this is a follow up to Ask Sage question 55323" and you can add a comment here saying "See a follow-up question at Ask Sage question 55xxx".

I have problems to implement the following set up:

I want to have a field $K = \mathbb{Q}_3$ and an extension $L = \mathbb{Q}_3(\alpha)$ over $K$ where $f:=\min_K(\alpha) = x^4 - 3x^2 + 18$. This extension has degree $4$ and ramification index $2$. Furthermore, let $F/K$ be the unique unramified extension of $K$ of degree $4$ which is generated by a primitive $5$-th root of unity $\zeta_5$.

Then one can show that $\varphi: \alpha \mapsto \frac{(2 \alpha^2 - 3)\sqrt{-\frac{2}{7}}}{\alpha},$ $\zeta_5 \mapsto \zeta_5^3$ is an element of the Galois group of $LF/K$. Now let $L' = (LF)^{\langle \varphi \rangle}$. This must be a quadratic and totally ramified extension of $K$. There are only two possibilities for that: $K(\sqrt{3})$ or $K(\sqrt{-3})$.

Question: How to determine whether $L' = K(\sqrt{3})$ or $L' = K(\sqrt{-3})$ (or equivalently, $\varphi(\sqrt{3}) = \sqrt{3}$ or $\varphi(\sqrt{-3}) = \sqrt{-3}$)?

Since I only have only superficial knowledge about Sage, I was not even able to set up the easy things like the field $L$ properly. When I use

K = Qp(3)
R.<x> = ZZ[]
f = x^4 - 3*x^2 + 18
L.<alpha> = K.extension(f)

I get an error that my polynomial $f$ must be either unramified or Eisenstein (which of course does not exist since $L/K$ is neither unramified nor totally ramified). Furthermore, I have no idea how to approach with my problem with Sage otherwise. And since computation by hand is pretty hard in this case (I already tried!), it would be nice to solve with problem here, so I can use it for similar computations in the future.

Could you please help me with this problem?

Welcome to Ask Sage! Thank you for your question!

Edited answer to partially address multicharacter generator names.

Sadly they are still not latexing properly.