minpoly raises 'NotAlgebraic' for tan(13*pi/45)

[hyadav2k]

This issue was encountered when we tried to construct a DomainMatrix from a matrix setting extension=True:
Try creating DomainMatrix from M = Matrix([[1, -1, 8], [-sqrt(3), tan(13*pi/45), 0]]) as DOM = DomainMatrix.from_Matrix(M, extension=True)

Since tan(13*pi/45) is an algebraic element(13/45 is rational), minpoly shouldn't be raising an error, rather returning a Poly/ Expr.
The error in terminal:

In [1]: minpoly(tan(13*pi/45))
---------------------------------------------------------------------------
NotAlgebraic                              Traceback (most recent call last)
<ipython-input-1-c699111c7c38> in <module>
----> 1 minpoly(tan(13*pi/45))

~\sympy\sympy\polys\numberfields.py in minimal_polynomial(ex, x, compose, polys, domain)
    670
    671     if compose:
--> 672         result = _minpoly_compose(ex, x, domain)
    673         result = result.primitive()[1]
    674         c = result.coeff(x**degree(result, x))

~\sympy\sympy\polys\numberfields.py in _minpoly_compose(ex, x, dom)
    584         res = _minpoly_rootof(ex, x)
    585     else:
--> 586         raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)
    587     return res
    588

NotAlgebraic: tan(13*pi/45) doesn't seem to be an algebraic element

On looking closely, it is found that this is caused as there is no computational method for tan in _minpoly_compose as there are for sin and cos. Refer here

So, on adding condition for tan,

diff --git a/sympy/polys/numberfields.py b/sympy/polys/numberfields.py
index 79f60d3e89..52e28b191c 100644
--- a/sympy/polys/numberfields.py
+++ b/sympy/polys/numberfields.py
@@ -11,7 +11,7 @@
 from sympy.core.exprtools import Factors
 from sympy.core.function import _mexpand
 from sympy.functions.elementary.exponential import exp
-from sympy.functions.elementary.trigonometric import cos, sin
+from sympy.functions.elementary.trigonometric import cos, sin, tan
 from sympy.ntheory import sieve
 from sympy.ntheory.factor_ import divisors
 from sympy.polys.densetools import dup_eval
@@ -578,6 +578,8 @@ def _minpoly_compose(ex, x, dom):
         res = _minpoly_sin(ex, x)
     elif ex.__class__ is cos:
         res = _minpoly_cos(ex, x)
+    elif ex.__class__ is tan:
+        res = _minpoly_sin(sin(ex.args[0]), x) * minimal_polynomial(1/cos(ex.args[0]))
     elif ex.__class__ is exp:
         res = _minpoly_exp(ex, x)
     elif ex.__class__ is CRootOf:

This is gets evaluted though as:

>>> minpoly(tan(13*pi/45))
(_x**12 + 24*_x**11 - 144*_x**10 - 248*_x**9 + 1680*_x**8 + 864*_x**7 - 7168*_x**6 - 1152*_x**5 + 13824*_x**4 + 512*_x**3 - 12288*_x**2 + 4096)*(16777216*_x**24 - 100663296*_x**22 + 264241152*_x**20 - 398196736*_x**18 + 379846656*_x**16 - 238436352*_x**14 + 99147776*_x**12 - 26797056*_x**10 + 4485888*_x**8 - 423872*_x**6 + 18912*_x**4 - 288*_x**2 + 1)

Even after the diffs, our main issue of constructing the DomainMatrix is unsolved.

[oscarbenjamin]

We can't just multiply the minimal polynomials like that to get the minimal polynomial of the product. I think the correct answer is

In [12]: arg = 13*pi/45

In [13]: minpoly(sin(arg)/cos(arg))
Out[13]: 
 24        22          20           18            16            14            12            10            8           6          4        2    
x   - 852⋅x   + 26562⋅x   - 288452⋅x   + 1465839⋅x   - 3852456⋅x   + 5407388⋅x   - 3994152⋅x   + 1475055⋅x  - 250372⋅x  + 15810⋅x  - 276⋅x  + 1

although minpoly is very slow at computing this.

We can check this numerically:

In [17]: s = minpoly(sin(arg)/cos(arg))

In [18]: Poly(s).nroots()
Out[18]: 
[-28.6362532829156, -4.01078093353584, -2.4750868534163, -2.0503038415793, -1.48256096851274, -1.27994163219308, -0.965688774807074, -0.624869351909327, -0.
531709431661479, -0.286745385758808, -0.140540834702391, -0.0699268119435104, 0.0699268119435104, 0.140540834702391, 0.286745385758808, 0.531709431661479, 0
.624869351909327, 0.965688774807074, 1.27994163219308, 1.48256096851274, 2.0503038415793, 2.4750868534163, 4.01078093353584, 28.6362532829156]

In [19]: tan(arg).n()
Out[19]: 1.27994163219308

In [22]: Poly(s).is_irreducible
Out[22]: True

[jksuom]

The minimal polynomial of tan(x), for x a rational multiple of pi, can be computed in the same way as for sin(x) and cos(x). Writing Euler's formula in the form

    cos(x)*(1 + I*tan(x)) = exp(I*x)

and raising to the n'th power, we get

    cos(x)**n*(P(tan(x)) + I*Q(tan(x))) = cos(n*x) + I*sin(n*x),

where P(tan(x)) (Q(tan(x)) is a polynomial in even (odd) powers of tan(x). Then

   P(tan(x)) = 0  or Q(tan(x)) = 0

if n*x is an odd resp. even multiple of pi/2.

[oscarbenjamin]

So like this:

In [12]: arg = 13*pi/45

In [13]: arg / (pi/2)
Out[13]: 
26
──
45

In [14]: x = Symbol('x', real=True)

In [16]: factor(im((1+I*x)**45))
Out[16]: 
  ⎛ 2    ⎞ ⎛ 4       2    ⎞ ⎛ 6       4       2    ⎞ ⎛ 8       6        4     
x⋅⎝x  - 3⎠⋅⎝x  - 10⋅x  + 5⎠⋅⎝x  - 33⋅x  + 27⋅x  - 3⎠⋅⎝x  - 92⋅x  + 134⋅x  - 28

  2    ⎞ ⎛ 24        22          20           18            16            14  
⋅x  + 1⎠⋅⎝x   - 852⋅x   + 26562⋅x   - 288452⋅x   + 1465839⋅x   - 3852456⋅x   +

          12            10            8           6          4        2    ⎞
 5407388⋅x   - 3994152⋅x   + 1475055⋅x  - 250372⋅x  + 15810⋅x  - 276⋅x  + 1⎠

[hyadav2k]

@oscarbenjamin, @jksuom According to the theory described above, tried implementing _minpoly_tan as:

def _minpoly_tan(ex):
    """
    Returns the minimal polynomial of ``tan(ex)``
    see https://github.com/sympy/sympy/issues/21430
    """
    x = Symbol('x', real=True)
    c, a = ex.args[0].as_coeff_Mul()
    if a is pi:
        if c.is_rational:
            c = c * 2
            n = c.q
            if c.p % 2 == 0:
                res = im((1+I*x)**n)
            else:
                res = re((1+I*x)**n)

            return res
    raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)

This seems to work for our case:

>>> p = minpoly(tan(13*pi/45))
>>> p
x**45 - 990*x**43 + 148995*x**41 - 8145060*x**39 + 215553195*x**37 - 3190187286*x**35 + 28760021745*x**33 - 166871334960*x**31 + 646626422970*x**29 - 1715884494940*x**27 + 3169870830126*x**25 - 4116715363800*x**23 + 3773655750150*x**21 - 2438362177020*x**19 + 1103068603890*x**17 - 344867425584*x**15 + 73006209045*x**13 - 10150595910*x**11 + 886163135*x**9 - 45379620*x**7 + 1221759*x**5 - 14190*x**3 + 45*x
>>> s = Poly(p)
>>> a = s.nroots()
>>> a
[-28.6362532829156, -9.51436445422259, -5.67128181961771, -4.01078093353584, -3.07768353717525, -2.47508685341630, -2.05030384157930, -1.73205080756888, -1.48256096851274, -1.27994163219308, -1.11061251482919, -0.965688774807074, -0.839099631177280, -0.726542528005361, -0.624869351909327, -0.531709431661479, -0.445228685308536, -0.363970234266202, -0.286745385758808, -0.212556561670022, -0.140540834702391, -0.0699268119435104, 0, 0.0699268119435104, 0.140540834702391, 0.212556561670022, 0.286745385758808, 0.363970234266202, 0.445228685308536, 0.531709431661479, 0.624869351909327, 0.726542528005361, 0.839099631177280, 0.965688774807074, 1.11061251482919, 1.27994163219308, 1.48256096851274, 1.73205080756888, 2.05030384157930, 2.47508685341630, 3.07768353717525, 4.01078093353584, 5.67128181961771, 9.51436445422259, 28.6362532829156]
>>> b = tan(13*pi/45).n()
>>> b
1.27994163219308
>>> b in a
True

Also tried for the other case when P(tan(x)) = 0 :

>>> p = minpoly(tan(5*pi/46))
>>> p
-23*x**22 + 1771*x**20 - 33649*x**18 + 245157*x**16 - 817190*x**14 + 1352078*x**12 - 1144066*x**10 + 490314*x**8 - 100947*x**6 + 8855*x**4 - 253*x**2 + 1
>>> s = Poly(p)
>>> a = s.nroots()
>>> a
[-7.27554032218333, -3.56904674291811, -2.30223090468814, -1.64442994147187, -1.22916512298902, -0.933934804224966, -0.705877079496608, -0.518158235091440, -0.355400085379703, -0.207802389613187, -0.0684018739205795, 0.0684018739205795, 0.207802389613187, 0.355400085379703, 0.518158235091440, 0.705877079496608, 0.933934804224966, 1.22916512298902, 1.64442994147187, 2.30223090468814, 3.56904674291811, 7.27554032218333]
>>> b = tan(5*pi/46).n()
>>> b
0.355400085379703
>>> b in a
True

Kindly share if this will work.
Thanks.

[jksuom]

res = im((1+I*x)**n) is a succinct way of writing the polynomial but not very efficient for production code. I would build the polynomial in a for loop adding one monomial at a time and computing the binomial coefficients (with alternating signs) recursively. For the minimal polynomial, one has to factor this polynomial and choose the appropriate factor.

[hyadav2k]

res = im((1+I*x)**n) is a succinct way of writing the polynomial but not very efficient for production code.

Sure, thanks.

I would build the polynomial in a for loop adding one monomial at a time and computing the binomial coefficients (with alternating signs) recursively. For the minimal polynomial, one has to factor this polynomial and choose the appropriate factor.

@jksuom, the below implementation is not entirely as you instructed, we are using binomial_coefficients_list instead of computing the binomial coefficients recursively.

def _minpoly_tan(ex):
    """
    Returns the minimal polynomial of ``tan(ex)``
    see https://github.com/sympy/sympy/issues/21430
    """
    x = Symbol('x', real=True)
    c, a = ex.args[0].as_coeff_Mul()
    if a is pi:
        if c.is_rational:
            c = c * 2
            n = c.q
            a = binomial_coefficients_list(n)
            if c.p % 2 == 0:
                coeff = {j : coeff*I**j for j, coeff in enumerate(a) if j % 2 != 0}
            else:
                coeff = {j : coeff*I**j for j, coeff in enumerate(a) if j % 2 == 0}

            a = [x**i*coeff[i] for i in coeff.keys()]
            r = sum(a)            
            _, factors = factor_list(r)
            res = _choose_factor(factors, x, ex)
            return res

    raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)

Kindly review.
Thanks.

[jksuom]

I would rather avoid introducing I in the computations. Otherwise, binomial_coefficient_list could be fairly usable. I would use looping over range(0, n+1, 2) or range(1, n+1, 2).

[hyadav2k]

I would rather avoid introducing I in the computations. Otherwise,binomial_coefficient_list` could be fairly usable.

Can we introduce an array b = [1, I, -1, -I] , now we can regulate as:

if c.p % 2 == 0:
    coeff = {j : a[j]*b[j%4] for j in range(1, n+1, 2)}
else:
    coeff = {j : a[j]*b[j%4] for j in range(0, n+1, 2)}

Else we can add an option complex=True to binomial_coefficient_list which can return for (1 + I*x)**n, but I am not sure, if this would be helpful in other places.

One more thing to ask, do you prefer using list over dict?
Thanks.

[jksuom]

The coefficients should all be real. They could be computed recursively in the same way as binomial_coefficient_list does. I would also build a list of monomials with proper coefficients. Then the polynomial would be Add(*list).

[hyadav2k]

@jksuom, as instructed, we know compute the coefficients analogous to binomial_coefficient_list. d is the list of monomials with real coefficients with signs accordingly for (1 + I*x)**n.

def _minpoly_tan(ex):
    """
    Returns the minimal polynomial of ``tan(ex)``
    see https://github.com/sympy/sympy/issues/21430
    """
    x = Symbol('x', real=True)
    c, a = ex.args[0].as_coeff_Mul()
    if a is pi:
        if c.is_rational:
            c = c * 2
            n = int(c.q)
            d = [1] * (n + 1)
            a = 1
            for k in range(1, n + 1):
                a = (a * (n - k + 1))//k
                if k % 4 == 2 or k % 4 == 3:
                    d[k] = -a
                else:
                    d[k] = a

            if c.p % 2 == 0:
                coeff = {j : d[j] for j in range(1, n+1, 2)}
            else:
                coeff = {j : d[j] for j in range(0, n+1, 2)}

            a = [x**i*coeff[i] for i in coeff.keys()]
            r = Add(*a)            
            _, factors = factor_list(r)
            res = _choose_factor(factors, x, ex)
            return res

    raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)

Knidly review.
Thanks,