Solve polynomial--documentation guide

[bertiewooster]

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

• other
  • A new, user-friendly page on solving (finding the roots of) polynomial equations algebraically was added to the Guides section.

[sympy-bot]

✅

Hi, I am the SymPy bot (v167). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• other
  • A new, user-friendly page on solving (finding the roots of) polynomial equations algebraically was added to the Guides section. (#24013 by @bertiewooster)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.12.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* other
  * A new, user-friendly page on solving (finding the roots of) polynomial equations algebraically was added to the Guides section.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[sympy-bot]

🟠

Hi, I am the SymPy bot (v167). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

• 164ccb8:
  • doc/src/guides/solving/solve-polynomial.md

If these files were added/deleted on purpose, you can ignore this message.

[bertiewooster]

I just started this page and seek guidance on what it should include:

• Poly() function? (Is there a good discussion of that function? That API link doesn't give much detail on it, for example why you should use it instead of a bare expression such as x^2 + x.)
• solve_poly_system(), in addition to solve?
  • @asmeurer is it this function that will give the multiplicity of a root (how many times it appears as a solution), or some other function?

[oscarbenjamin]

There are lots of functions:

• factor
• roots
• nroots
• RootOf
• real_roots (aka Poly().real_roots)
• Poly().all_roots (there should be an all_roots function)

The factor function can factorise a polynomial in a given polynomial ring which can reveal roots lie in the coefficient ring e.g. if the polynomial has rational coefficients then factor will reveal any rational roots. If the coefficients are polynomials involving say a symbol a with rational coefficients then any roots that are polynomial functions of a with rational coefficients will be revealed:

In [12]: p = expand((x - a**2)*(x + a + a**3))

In [13]: p
Out[13]: 
   5    3      3    2            2
- a  + a ⋅x - a  - a ⋅x + a⋅x + x 

In [14]: factor(p)
Out[14]: 
⎛   2    ⎞ ⎛ 3        ⎞
⎝- a  + x⎠⋅⎝a  + a + x⎠

The roots function uses a combination of techniques (factorisation, decomposition, radical formulae) to find expressions in radicals if possible for the roots. When it can find some radical expressions for the roots it returns them along with their multiplicity. This function will fail for most high degree polynomials because of the Abel-Ruffini theorem:
https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem

The roots function is the only method that can give explicit expressions for the roots of polynomials that have symbolic coefficients (i.e. there are symbols in the coefficients) if factor does not reveal them.

The nroots function can compute numerical approximations of the roots of any polynomial whose coefficients can be numerically evaluated with evalf (i.e. they do not have free symbols):

In [15]: nroots(x**2 + sqrt(2)*x + pi)
Out[15]: [-0.707106781186548 - 1.62529771229452⋅ⅈ, -0.707106781186548 + 1.62529771229452⋅ⅈ]

There is no restriction on the degree for which nroots can be used but it can fail sometimes for polynomials that are numerically ill conditioned e.g. Wilkinson's polynomial:
https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial

The nroots function is analogous to numpy's roots function. Usually the difference between these two is just that nroots is more accurate but slower.

The RootOf, real_roots and all_roots functions can compute all the roots of any polynomial with purely rational coefficients but of any degree. RootOf represent those roots semi-implicitly when the irreducible factors are degree three or more. These are a general representation of algebraic numbers:
https://en.wikipedia.org/wiki/Algebraic_number

Symbolic algorithms that need to work consistently with the roots of potentially high degree polynomials generally do so without needing explicit representations of those roots (as would be returned by solve or roots) and instead use things like discriminants, resultants and Groebner bases. There are numerical problems with the quadratic, cubic and quartic and other formulas and radicals are difficult to manipulate symbolically. It's important to consider whether you really explicit solutions at all or what you want to do with them because potentially the radical formulae returned by solve or roots are actually not useful in practice.

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [41d90958]       [d44f27e1]
     <sympy-1.11.1^0>                 
-         960±3μs          621±4μs     0.65  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-     2.79±0.01ms         1.17±0ms     0.42  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.64±0.01ms         1.71±0ms     0.30  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[oscarbenjamin]

The examples here are too simple. I think it's important to illustrate that RootOf, real_roots and all_roots can find all the roots exactly of a polynomial of arbitrarily large degree. The examples here are only for quadratic in which case RootOf uses radical formulae but its real power is in being able to avoid radicals and not being limited by the Abel-Ruffini theorem. To demonstrate this it is necessary to give examples that show the actual RootOf:

In [7]: p = x**5 - x + 1

In [8]: roots(p)
Out[8]: {}

In [9]: real_roots(p)
Out[9]: 
⎡       ⎛ 5           ⎞⎤
⎣CRootOf⎝x  - x + 1, 0⎠⎦

In [10]: Poly(p, x).all_roots()
Out[10]: 
⎡       ⎛ 5           ⎞         ⎛ 5           ⎞         ⎛ 5           ⎞         ⎛ 5           ⎞    
⎣CRootOf⎝x  - x + 1, 0⎠, CRootOf⎝x  - x + 1, 1⎠, CRootOf⎝x  - x + 1, 2⎠, CRootOf⎝x  - x + 1, 3⎠, CR

     ⎛ 5           ⎞⎤
ootOf⎝x  - x + 1, 4⎠⎦

In [11]: r1, r2, r3, r4, r5 = Poly(p, x).all_roots()

In [12]: r1
Out[12]: 
       ⎛ 5           ⎞
CRootOf⎝x  - x + 1, 0⎠

In [13]: r1.n(100)
Out[13]: 
-1.167303978261418684256045899854842180720560371525489039140082449275651903429527053180685205049728
673

In [14]: r1.is_real
Out[14]: True

In [15]: r2.n()
Out[15]: -0.181232444469875 - 1.08395410131771⋅ⅈ

In [16]: r3.n()
Out[16]: -0.181232444469875 + 1.08395410131771⋅ⅈ

In [17]: r2
Out[17]: 
       ⎛ 5           ⎞
CRootOf⎝x  - x + 1, 1⎠

In [18]: r2.conjugate()
Out[18]: 
       ⎛ 5           ⎞
CRootOf⎝x  - x + 1, 2⎠

In [19]: r2.is_real
Out[19]: False

Using solve will also give the complex roots where possible but it is less efficient than using all_roots directly:

In [20]: solve(p, x)
Out[20]: 
⎡       ⎛ 5           ⎞         ⎛ 5           ⎞         ⎛ 5           ⎞         ⎛ 5           ⎞    
⎣CRootOf⎝x  - x + 1, 0⎠, CRootOf⎝x  - x + 1, 1⎠, CRootOf⎝x  - x + 1, 2⎠, CRootOf⎝x  - x + 1, 3⎠, CR

     ⎛ 5           ⎞⎤
ootOf⎝x  - x + 1, 4⎠⎦

Note that although RootOf is in some sense numeric it does exactly represent the root in a way that can be manipulated symbolically as well as computed to arbitrary precision. The RootOf representation makes it possible to precisely:

1. Compute all roots of a polynomial with exact rational coefficients.
2. Decide exactly the multiplicity of every root.
3. Determine exactly whether roots are real or not.
4. Order the real and complex roots precisely.
5. Know which roots are complex conjugate pairs of each other.
6. Say precisely which roots are rational vs irrational.
7. Represent every possible algebraic number exactly.

The other numerical methods like numpy.roots, nroots, nsolve etc cannot do any of these things robustly or simply at all. Similarly the radical expressions returned by solve or roots cannot do these things robustly when numerically evaluated because of issues such as casus irreducibilus:

In [22]: r1, r2, r3 = roots(2*x**3 - 9*x**2 - 6*x + 3)

In [23]: r1
Out[23]: 
                            ___________
3           13             ╱ 39   13⋅ⅈ 
─ + ───────────────── + 3 ╱  ── + ──── 
2         ___________   ╲╱   8     4   
         ╱ 39   13⋅ⅈ                   
    4⋅3 ╱  ── + ────                   
      ╲╱   8     4                     

In [24]: r1.n()
Out[24]: 5.03651666318515 - 0.e-21⋅ⅈ

In [25]: r2
Out[25]: 
                                                      ___________
3                 13                 ⎛  1   √3⋅ⅈ⎞    ╱ 39   13⋅ⅈ 
─ + ────────────────────────────── + ⎜- ─ + ────⎟⋅3 ╱  ── + ──── 
2                      ___________   ⎝  2    2  ⎠ ╲╱   8     4   
      ⎛  1   √3⋅ⅈ⎞    ╱ 39   13⋅ⅈ                                
    4⋅⎜- ─ + ────⎟⋅3 ╱  ── + ────                                
      ⎝  2    2  ⎠ ╲╱   8     4                                  

In [26]: r2.n()
Out[26]: -0.876359818402395 + 0.e-20⋅ⅈ

https://en.wikipedia.org/wiki/Casus_irreducibilis

I think it's important to show the complicated expressions that are generated for cubics and quartics by roots and solve explicitly:

In [27]: r1, r2, r3, r4 = roots(x**4 + 3*x**2 + 2*x + 1)

In [28]: r1
Out[28]: 
            _______________________________________________________________________________________
           ╱            ______________                                                             
          ╱            ╱   1   √237⋅ⅈ                                4                             
         ╱   -4 - 2⋅3 ╱  - ─ + ──────  + ──────────────────────────────────────────────────────────
        ╱           ╲╱     8     36              __________________________________________________
       ╱                                        ╱                                   ______________ 
      ╱                                        ╱                7                  ╱   1   √237⋅ⅈ  
     ╱                                        ╱   -2 + ──────────────────── + 2⋅3 ╱  - ─ + ──────  
    ╱                                        ╱               ______________     ╲╱     8     36    
   ╱                                        ╱               ╱   1   √237⋅ⅈ                         
  ╱                                        ╱           6⋅3 ╱  - ─ + ──────                         
╲╱                                       ╲╱              ╲╱     8     36                           
───────────────────────────────────────────────────────────────────────────────────────────────────
                                                             2                                     

________________________                                                             
                                                                                     
            7                                                                        
 - ────────────────────                                                              
         ______________            __________________________________________________
        ╱   1   √237⋅ⅈ            ╱                                   ______________ 
   6⋅3 ╱  - ─ + ──────           ╱                7                  ╱   1   √237⋅ⅈ  
     ╲╱     8     36            ╱   -2 + ──────────────────── + 2⋅3 ╱  - ─ + ──────  
                               ╱               ______________     ╲╱     8     36    
                              ╱               ╱   1   √237⋅ⅈ                         
                             ╱           6⋅3 ╱  - ─ + ──────                         
                           ╲╱              ╲╱     8     36                           
──────────────────────── - ──────────────────────────────────────────────────────────
                                                       2                             

In [29]: r1.n()
Out[29]: -0.349745826211722 - 0.438990337475312⋅ⅈ

It's worth considering that if all you can do to make sense of the expression is evaluate it numerically then perhaps it would have been better just to compute the roots numerically:

In [30]: nroots(x**4 + 3*x**2 + 2*x + 1)
Out[30]: 
[-0.349745826211722 - 0.438990337475312⋅ⅈ, -0.349745826211722 + 0.438990337475312⋅ⅈ, 0.349745826211
722 - 1.74697789611327⋅ⅈ, 0.349745826211722 + 1.74697789611327⋅ⅈ]

If you wanted an exact representation then in fact the best exact representation is RootOf which can also be evaluated numerically:

In [31]: Poly(x**4 + 3*x**2 + 2*x + 1).all_roots()
Out[31]: 
⎡       ⎛ 4      2             ⎞         ⎛ 4      2             ⎞         ⎛ 4      2             ⎞ 
⎣CRootOf⎝x  + 3⋅x  + 2⋅x + 1, 0⎠, CRootOf⎝x  + 3⋅x  + 2⋅x + 1, 1⎠, CRootOf⎝x  + 3⋅x  + 2⋅x + 1, 2⎠,

        ⎛ 4      2             ⎞⎤
 CRootOf⎝x  + 3⋅x  + 2⋅x + 1, 3⎠⎦

If you want numeric approximations of the real roots but you want to know exactly which roots are real then the best method is real_roots with evalf:

In [34]: [r.n(2) for r in real_roots(x**4 - 3*x**2 - 2*x + 1)]
Out[34]: [0.34, 1.9]

Note that none of the approximate numeric methods can distinguish precisely which roots are real. Consider Wilkinson's polynomial:

In [35]: p = prod((x - i) for i in range(1, 31))

In [36]: p
Out[36]: 
(x - 30)⋅(x - 29)⋅(x - 28)⋅(x - 27)⋅(x - 26)⋅(x - 25)⋅(x - 24)⋅(x - 23)⋅(x - 22)⋅(x - 21)⋅(x - 20)⋅
(x - 19)⋅(x - 18)⋅(x - 17)⋅(x - 16)⋅(x - 15)⋅(x - 14)⋅(x - 13)⋅(x - 12)⋅(x - 11)⋅(x - 10)⋅(x - 9)⋅(
x - 8)⋅(x - 7)⋅(x - 6)⋅(x - 5)⋅(x - 4)⋅(x - 3)⋅(x - 2)⋅(x - 1)

In [37]: p = expand(p)

In [38]: p
Out[38]: 
 30        29           28             27               26                 25                  24  
x   - 465⋅x   + 103385⋅x   - 14631225⋅x   + 1480321269⋅x   - 114009431445⋅x   + 6949189247325⋅x   -

                  23                      22                       21                         20   
 344092707928125⋅x   + 14097793282984515⋅x   - 484338676679532675⋅x   + 14090257524223082475⋅x   - 

                       19                           18                             17              
349600545868057540875⋅x   + 7435941626111727234855⋅x   - 136055808711963322871175⋅x   + 21458832493

                16                               15                                14              
34501452139775⋅x   - 29197210605623737977801375⋅x   + 342563613932937660652700640⋅x   - 34602661104

                   13                                  12                                   11     
93898677911394000⋅x   + 30006513636556697864066736800⋅x   - 222457423246962063058403076000⋅x   + 14

                               10                                    9                             
01937624086807501691142239744⋅x   - 7454161471690660700139655157760⋅x  + 33114629767614997850763390

        8                                      7                                      6            
570240⋅x  - 121365366674745136523074652102400⋅x  + 360930788158836812805614538878976⋅x  - 851899888

                          5                                       4                                
505423112503184251412480⋅x  + 1547794975254719737111781253120000⋅x  - 20707922020245946836608666419

       3                                       2                                                   
20000⋅x  + 1902893785240928209998216560640000⋅x  - 1059681761389533859949327155200000⋅x + 265252859

                        
812191058636308480000000

In [39]: factor(p)
Out[39]: 
(x - 30)⋅(x - 29)⋅(x - 28)⋅(x - 27)⋅(x - 26)⋅(x - 25)⋅(x - 24)⋅(x - 23)⋅(x - 22)⋅(x - 21)⋅(x - 20)⋅
(x - 19)⋅(x - 18)⋅(x - 17)⋅(x - 16)⋅(x - 15)⋅(x - 14)⋅(x - 13)⋅(x - 12)⋅(x - 11)⋅(x - 10)⋅(x - 9)⋅(
x - 8)⋅(x - 7)⋅(x - 6)⋅(x - 5)⋅(x - 4)⋅(x - 3)⋅(x - 2)⋅(x - 1)

In [42]: roots(p, multiple=True)
Out[42]: 
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
 28, 29, 30]

In [40]: nroots(p)
---------------------------------------------------------------------------
NoConvergence

In [41]: np.roots(Poly(p).all_coeffs())
Out[41]: 
array([32.07636983+0.j        , 31.29103418+2.57752104j,
       31.29103418-2.57752104j, 29.17520185+4.8010735j ,
       29.17520185-4.8010735j , 26.19348536+6.28087068j,
       26.19348536-6.28087068j, 22.87013876+6.83360078j,
       22.87013876-6.83360078j, 19.66310309+6.55232555j,
       19.66310309-6.55232555j, 16.84020824+5.6911237j ,
       16.84020824-5.6911237j , 14.47784031+4.52035773j,
       14.47784031-4.52035773j, 12.50458382+3.23884483j,
       12.50458382-3.23884483j, 11.53631012+0.j        ,
       10.72355387+1.9306279j , 10.72355387-1.9306279j ,
        9.02521141+0.80837533j,  9.02521141-0.80837533j,
        8.83820369+0.j        ,  7.02160259+0.j        ,
        5.99874634+0.j        ,  5.00004638+0.j        ,
        3.99999929+0.j        ,  2.99999999+0.j        ,
        2.        +0.j        ,  1.        +0.j        ])

Here the roots should all be real but np.roots thinks that they are complex. There are significant errors in the roots which should all be integers. By contrast at least nroots lets you know that this polynomial is numerically problematic rather than just returning incorrect results. Note that this is not really a deficiency of np.roots but rather a limitation of fixed precision floating point arithmetic. Only arbitrary precision arithmetic (e.g. nroots or nsolve) can get numerically close approximations for all polynomials. Here nroots just needs permission to try harder:

In [43]: nroots(p, maxsteps=100)
Out[43]: 
[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0]

It looks as if nroots can distinguish real and complex roots but this is heuristic and will sometimes fail. Likewise it often looks like nroots or np.roots can identify multiple roots but it's not precise:

In [46]: p = expand((x - S(1)/3)**2)

In [47]: p
Out[47]: 
 2   2⋅x   1
x  - ─── + ─
      3    9

In [48]: coeffs = Poly(p).all_coeffs()

In [49]: coeffs
Out[49]: [1, -2/3, 1/9]

In [50]: np.roots(coeffs)
Out[50]: array([0.33333333+4.39029672e-09j, 0.33333333-4.39029672e-09j])

In [51]: np.roots(coeffs)[0]
Out[51]: (0.33333333333333337+4.390296719884795e-09j)

In [52]: np.roots(coeffs)[1]
Out[52]: (0.33333333333333337-4.390296719884795e-09j)

In [53]: nroots(p)
Out[53]: [0.333333333333333 - 3.71750825013552e-12⋅ⅈ, 0.333333333333333 + 4.02352310590228e-12⋅ⅈ]

Only exact arithmetic can make the distinction between roots that are close and multiple roots via the square free factorisation (sqf_factor):
https://en.wikipedia.org/wiki/Square-free_polynomial
Both factor, RootOf, solve, roots etc can find this:

In [55]: factor(p)
Out[55]: 
         2
(3⋅x - 1) 
──────────
    9     

In [56]: solve(p)
Out[56]: [1/3]

In [57]: roots(p)
Out[57]: {1/3: 2}

In [58]: real_roots(p)
Out[58]: [1/3, 1/3]

Another thing that I think is worth mentioning is that sometimes for a very high degree polynomial it can be more efficient to use nsolve to find a single root. For example if you have a high-degree polynomial and know that it has a unique positive root which is all that you want to find then it might be more efficient to find that using nsolve with the bisection method rather than computing e.g. all the real roots and then separating out the only one that you are interested in at the end.

[oscarbenjamin]

I think it would better to have a clear list of the possible functions to use near the top of this page with a short (e.g. one sentence) summary of when and why each might be used.

[oscarbenjamin]

There is an important point that users should be made aware of but isn't really covered here because you don't show any examples that have big horrible messy output. The point is this:

The cubic and quartic formulae give horrible expressions that are often not useful in practice.

It is important to understand this point in order to understand why it is better to use RootOf even if radical expressions are possible. This is the reason that both roots and solve have flags that can be used to disable the use of these formulae.

The followup point which is mentioned is this:

There is no general formula for quintics or higher degree polynomials.

The problem with emphasising this second point about the nonexistence of formulae in some cases is that it overlooks the fact that you probably don't want to use the formulae that do exist either.

When RootOf isn't possible e.g. because of symbolic coefficients then it is better not to try and compute the roots at all for anything more complicated than a quadratic equation. It is certainly too much to explain how to handle that situation in this guide but in fact the real solution when RootOf and factor don't work is that it is better just to stick with the polynomial itself and consider it as a representation of its own roots rather than demanding any explicit representation for the roots.

Techniques that use the polynomial as an implicit representation of its roots can be much more powerful and not limited by the issues discussed in this article. I'm not sure how best to convey this point to users but what it basically means is that the whole premise that you want to "solve a polynomial" is often misguided: usually the (factorised) polynomial already is the best way to represent its roots.

[bertiewooster]

Unfortunately, the SciPy documentation site is down, leading SymPy's sphinx tests and documentation builds to fail because of the intersphinx reference.

[asmeurer]

There is an example of a quartic here, but it's printed in str form which doesn't line wrap, so perhaps the point isn't brought home very well. Maybe we should include the LaTeX form of the expression just to show how complex it really it is.

[asmeurer]

We could also pprint it, but I'd rather just show it in LaTeX. pprint has issues rendering on some machines (at least until we address #15700).

[bertiewooster]

While I appreciate the point of emphasizing a very complex quartic solution by making all of it appear on the screen without scrolling, I don't think it's worth the effort--using either dollar math or {math}`...{math}` produces even worse formatting, namely the long formula escaping the horizontal bounds of the body and appearing to the right of the right nav pane:

I just added a comment to the code block that the radical formula returns very complex terms.

This was something we discussed before for another example that produced very complex terms, and I believe @asmeurer suggested changing the CSS to wrap within code blocks. If we want to emphasize long (many-term) results, I think we should make an issue to do that general solution.

[asmeurer]

You should use latex(formula) to get the latex. The screenshow you show there doesn't have valid LaTeX math. The correct version should be taller and less wide, due to the nested square roots. It should be overflowing under the sidebar either way. That's something that needs to fixed in the CSS (it's possible we are messing it up somewhere with our custom CSS, or else it's an upstream issue with Furo).

You can also hard wrap the math, I believe using \\ in the formula (unfortunately, MathJax 3 does not support automatically wrapping formulas, see #22467).

[bertiewooster]

You should use latex(formula) to get the latex. The screenshow you show there doesn't have valid LaTeX math. The correct version should be taller and less wide, due to the nested square roots. It should be overflowing under the sidebar either way. That's something that needs to fixed in the CSS (it's possible we are messing it up somewhere with our custom CSS, or else it's an upstream issue with Furo).

You can also hard wrap the math, I believe using \\ in the formula (unfortunately, MathJax 3 does not support automatically wrapping formulas, see #22467).

Got it, thanks. I converted the formula to LaTeX using SymPy's latex command. The only trick was I had to then replace all double backslashes \\ with single backslashes \. Perhaps \\ works in reStructuredText and \ works in Markdown? Is there a flag in the latex() to output \ instead of \\?

Regarding web page formatting, that expression fits in the body section of the page. I experimented and, if the expression were longer, it would break out of the body section, in which case a line break would be helpful. (I played around and a a line break works if you use double dollar signs, e.g. $$a\\b$$). Once CircleCI builds the docs, I'll check that the LaTeX expression fits on my phone.

[bertiewooster]

Once CircleCI builds the docs, I'll check that the LaTeX expression fits on my phone.

The expression did extend beyond the body width on my phone, so I pushed a commit with a line break. I'm not thrilled with how it looks on a computer screen because the horizontal line of the radical on the second line could be visually confused with a fraction line at first glance:

[bertiewooster]

The formatting is better on a phone with the hard line break: The expression is in a scrollable container, so it doesn't extend beyond the body width. (It's unclear to me why the expression isn't in a scrollable container when it's on one line.)

iPhone screenshot

[asmeurer]

The formatting is better on a phone with the hard line break: The expression is in a scrollable container, so it doesn't extend beyond the body width. (It's unclear to me why the expression isn't in a scrollable container when it's on one line.)

Inline math is not expected to be that big. It's supposed to be "inline" with the body text. Large equations like this should definitely use the non-inline version.

[bertiewooster]

The formatting is better on a phone with the hard line break: The expression is in a scrollable container, so it doesn't extend beyond the body width. (It's unclear to me why the expression isn't in a scrollable container when it's on one line.)

Inline math is not expected to be that big. It's supposed to be "inline" with the body text. Large equations like this should definitely use the non-inline version.

I removed the hard line break but kept the non-inline formatting by using two dollar signs before and after the LaTeX, and it produces a scrollable container as desired.

[asmeurer]

The cubic formula looks better now. Are there any other objections to merging this? As @bertiewooster noted, the link to the issues page will be removed from all the guides in a separate pull request.