Repository: sympy/sympy

Issue:
the first element of `sympy.Array` object won't show in latex output
when the array length is greater than 10,  the latex environment will be `\begin{array} \end{array}`, and because of the missing mandatory arguments of `array` environment , the first element 0 is taken for the arguments.
![image](https://user-images.githubusercontent.com/66580331/173239078-44806fc5-980f-4403-8a11-ea3988431b0c.png)

it may be due to the method `sympy.printing.latex.LatexPrinter._print_NDimArray` 
![image](https://user-images.githubusercontent.com/66580331/173239113-4a62fad1-590e-4fe6-bcc9-5dc4bf7aa924.png)

one alternative solution is using `array` environment like `\beg{array}{} ... \end{array}`

and the same problem in matrix printing is handled properly
![image](https://user-images.githubusercontent.com/66580331/173239697-cd5bcc13-f69f-4ba1-8712-27d1e28a8b39.png)



Discussion:
sympy-bot: 

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

:x: There was an issue with the release notes. **Please do not close this pull request;** instead edit the description after reading the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

* No release notes were found. Please edit the PR description and write the release notes under `<!-- BEGIN RELEASE NOTES -->`.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- 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. -->
    Fixes #23620

    #### Brief description of what is fixed or changed
    after adding `{}` for the required argumets, `Array` object can be normally displayed
    in latex output
    ![屏幕截图 2022-06-13 092124](https://user-images.githubusercontent.com/66580331/173263951-a489c868-fa02-49a5-a693-32cc3fb63837.png)


    #### 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 -->

    <!-- END RELEASE NOTES -->


</details><p>

eagleoflqj: Could you please
1. Add an entry to mailmap according to https://github.com/sympy/sympy/wiki/Development-workflow#add-your-name-and-email-address-to-the-mailmap-file
2. Add a test at https://github.com/sympy/sympy/blob/master/sympy/printing/tests/test_latex.py. It should be a new `test_Array` function, containing the 2 cases in #23620.
eagleoflqj: Changes look good, thank you!
Ideally there doesn’t need to be merge commits. If you are familiar with git please squash commits, or I can do it for you.
viocha: > Changes look good, thank you! Ideally there doesn’t need to be merge commits. If you are familiar with git please squash commits, or I can do it for you.

i'm new to git and have no idea how to do it. thanks for your trouble.
github-actions[bot]: 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)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [77f1d79c]       [cbd00f11]
     <sympy-1.10.1^0>                 
+         115±1ms          214±8ms     1.86  sum.TimeSum.time_doit

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
eagleoflqj: Just a reminder: next time please avoid directly committing to `master` branch of your fork. Instead create a new branch to work on.

Patch:
diff --git a/.mailmap b/.mailmap
index 7e7b6f987b..9af816f66c 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1462,6 +1462,7 @@ tools4origins <tools4origins@gmail.com>
 user202729 <25191436+user202729@users.noreply.github.com>
 vedantc98 <vedantc98@gmail.com>
 vezeli <37907135+vezeli@users.noreply.github.com>
+viocha <66580331+viocha@users.noreply.github.com>
 vishal <vishal.panjwani15@gmail.com>
 w495 <w495@yandex-team.ru>
 wuyudi <wuyudi119@163.com>
diff --git a/sympy/functions/elementary/hyperbolic.py b/sympy/functions/elementary/hyperbolic.py
index 72a2955152..c2bebb3623 100644
--- a/sympy/functions/elementary/hyperbolic.py
+++ b/sympy/functions/elementary/hyperbolic.py
@@ -12,8 +12,7 @@
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.integers import floor
 from sympy.functions.elementary.trigonometric import (acot, asin, atan, cos,
-                                              cot, sin, tan,
-                                              _imaginary_unit_as_coefficient)
+                                                      cot, sin, tan)
 from sympy.polys.specialpolys import symmetric_poly
 
 
@@ -129,7 +128,7 @@ def eval(cls, arg):
             if arg is S.ComplexInfinity:
                 return S.NaN
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 return S.ImaginaryUnit * sin(i_coeff)
@@ -326,7 +325,7 @@ def eval(cls, arg):
             if arg is S.ComplexInfinity:
                 return S.NaN
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 return cos(i_coeff)
@@ -570,7 +569,7 @@ def eval(cls, arg):
             if arg is S.ComplexInfinity:
                 return S.NaN
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 if i_coeff.could_extract_minus_sign():
@@ -783,7 +782,7 @@ def eval(cls, arg):
             if arg is S.ComplexInfinity:
                 return S.NaN
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 if i_coeff.could_extract_minus_sign():
@@ -1152,7 +1151,7 @@ def eval(cls, arg):
             if arg.is_zero:
                 return S.Zero
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 return S.ImaginaryUnit * asin(i_coeff)
@@ -1240,32 +1239,6 @@ def fdiff(self, argindex=1):
         else:
             raise ArgumentIndexError(self, argindex)
 
-    @staticmethod
-    @cacheit
-    def _acosh_table():
-        return {
-            S.ImaginaryUnit: log(S.ImaginaryUnit*(1 + sqrt(2))),
-            -S.ImaginaryUnit: log(-S.ImaginaryUnit*(1 + sqrt(2))),
-            S.Half: S.Pi/3,
-            Rational(-1, 2): S.Pi*Rational(2, 3),
-            sqrt(2)/2: S.Pi/4,
-            -sqrt(2)/2: S.Pi*Rational(3, 4),
-            1/sqrt(2): S.Pi/4,
-            -1/sqrt(2): S.Pi*Rational(3, 4),
-            sqrt(3)/2: S.Pi/6,
-            -sqrt(3)/2: S.Pi*Rational(5, 6),
-            (sqrt(3) - 1)/sqrt(2**3): S.Pi*Rational(5, 12),
-            -(sqrt(3) - 1)/sqrt(2**3): S.Pi*Rational(7, 12),
-            sqrt(2 + sqrt(2))/2: S.Pi/8,
-            -sqrt(2 + sqrt(2))/2: S.Pi*Rational(7, 8),
-            sqrt(2 - sqrt(2))/2: S.Pi*Rational(3, 8),
-            -sqrt(2 - sqrt(2))/2: S.Pi*Rational(5, 8),
-            (1 + sqrt(3))/(2*sqrt(2)): S.Pi/12,
-            -(1 + sqrt(3))/(2*sqrt(2)): S.Pi*Rational(11, 12),
-            (sqrt(5) + 1)/4: S.Pi/5,
-            -(sqrt(5) + 1)/4: S.Pi*Rational(4, 5)
-            }
-
     @classmethod
     def eval(cls, arg):
         arg = sympify(arg)
@@ -1285,7 +1258,28 @@ def eval(cls, arg):
                 return S.Pi*S.ImaginaryUnit
 
         if arg.is_number:
-            cst_table = cls._acosh_table()
+            cst_table = {
+                S.ImaginaryUnit: log(S.ImaginaryUnit*(1 + sqrt(2))),
+                -S.ImaginaryUnit: log(-S.ImaginaryUnit*(1 + sqrt(2))),
+                S.Half: S.Pi/3,
+                Rational(-1, 2): S.Pi*Rational(2, 3),
+                sqrt(2)/2: S.Pi/4,
+                -sqrt(2)/2: S.Pi*Rational(3, 4),
+                1/sqrt(2): S.Pi/4,
+                -1/sqrt(2): S.Pi*Rational(3, 4),
+                sqrt(3)/2: S.Pi/6,
+                -sqrt(3)/2: S.Pi*Rational(5, 6),
+                (sqrt(3) - 1)/sqrt(2**3): S.Pi*Rational(5, 12),
+                -(sqrt(3) - 1)/sqrt(2**3): S.Pi*Rational(7, 12),
+                sqrt(2 + sqrt(2))/2: S.Pi/8,
+                -sqrt(2 + sqrt(2))/2: S.Pi*Rational(7, 8),
+                sqrt(2 - sqrt(2))/2: S.Pi*Rational(3, 8),
+                -sqrt(2 - sqrt(2))/2: S.Pi*Rational(5, 8),
+                (1 + sqrt(3))/(2*sqrt(2)): S.Pi/12,
+                -(1 + sqrt(3))/(2*sqrt(2)): S.Pi*Rational(11, 12),
+                (sqrt(5) + 1)/4: S.Pi/5,
+                -(sqrt(5) + 1)/4: S.Pi*Rational(4, 5)
+            }
 
             if arg in cst_table:
                 if arg.is_extended_real:
@@ -1414,7 +1408,7 @@ def eval(cls, arg):
                 from sympy.calculus.accumulationbounds import AccumBounds
                 return S.ImaginaryUnit*AccumBounds(-S.Pi/2, S.Pi/2)
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 return S.ImaginaryUnit * atan(i_coeff)
@@ -1523,7 +1517,7 @@ def eval(cls, arg):
             if arg is S.ComplexInfinity:
                 return S.Zero
 
-            i_coeff = _imaginary_unit_as_coefficient(arg)
+            i_coeff = arg.as_coefficient(S.ImaginaryUnit)
 
             if i_coeff is not None:
                 return -S.ImaginaryUnit * acot(i_coeff)
@@ -1609,10 +1603,26 @@ def fdiff(self, argindex=1):
         else:
             raise ArgumentIndexError(self, argindex)
 
-    @staticmethod
-    @cacheit
-    def _asech_table():
-            return {
+    @classmethod
+    def eval(cls, arg):
+        arg = sympify(arg)
+
+        if arg.is_Number:
+            if arg is S.NaN:
+                return S.NaN
+            elif arg is S.Infinity:
+                return S.Pi*S.ImaginaryUnit / 2
+            elif arg is S.NegativeInfinity:
+                return S.Pi*S.ImaginaryUnit / 2
+            elif arg.is_zero:
+                return S.Infinity
+            elif arg is S.One:
+                return S.Zero
+            elif arg is S.NegativeOne:
+                return S.Pi*S.ImaginaryUnit
+
+        if arg.is_number:
+            cst_table = {
                 S.ImaginaryUnit: - (S.Pi*S.ImaginaryUnit / 2) + log(1 + sqrt(2)),
                 -S.ImaginaryUnit: (S.Pi*S.ImaginaryUnit / 2) + log(1 + sqrt(2)),
                 (sqrt(6) - sqrt(2)): S.Pi / 12,
@@ -1639,27 +1649,6 @@ def _asech_table():
                 (-sqrt(6) - sqrt(2)): 7*S.Pi / 12,
             }
 
-    @classmethod
-    def eval(cls, arg):
-        arg = sympify(arg)
-
-        if arg.is_Number:
-            if arg is S.NaN:
-                return S.NaN
-            elif arg is S.Infinity:
-                return S.Pi*S.ImaginaryUnit / 2
-            elif arg is S.NegativeInfinity:
-                return S.Pi*S.ImaginaryUnit / 2
-            elif arg.is_zero:
-                return S.Infinity
-            elif arg is S.One:
-                return S.Zero
-            elif arg is S.NegativeOne:
-                return S.Pi*S.ImaginaryUnit
-
-        if arg.is_number:
-            cst_table = cls._asech_table()
-
             if arg in cst_table:
                 if arg.is_extended_real:
                     return cst_table[arg]*S.ImaginaryUnit
@@ -1745,25 +1734,6 @@ def fdiff(self, argindex=1):
         else:
             raise ArgumentIndexError(self, argindex)
 
-    @staticmethod
-    @cacheit
-    def _acsch_table():
-        return {
-                S.ImaginaryUnit: -S.Pi / 2,
-                S.ImaginaryUnit*(sqrt(2) + sqrt(6)): -S.Pi / 12,
-                S.ImaginaryUnit*(1 + sqrt(5)): -S.Pi / 10,
-                S.ImaginaryUnit*2 / sqrt(2 - sqrt(2)): -S.Pi / 8,
-                S.ImaginaryUnit*2: -S.Pi / 6,
-                S.ImaginaryUnit*sqrt(2 + 2/sqrt(5)): -S.Pi / 5,
-                S.ImaginaryUnit*sqrt(2): -S.Pi / 4,
-                S.ImaginaryUnit*(sqrt(5)-1): -3*S.Pi / 10,
-                S.ImaginaryUnit*2 / sqrt(3): -S.Pi / 3,
-                S.ImaginaryUnit*2 / sqrt(2 + sqrt(2)): -3*S.Pi / 8,
-                S.ImaginaryUnit*sqrt(2 - 2/sqrt(5)): -2*S.Pi / 5,
-                S.ImaginaryUnit*(sqrt(6) - sqrt(2)): -5*S.Pi / 12,
-                S(2): -S.ImaginaryUnit*log((1+sqrt(5))/2),
-            }
-
     @classmethod
     def eval(cls, arg):
         arg = sympify(arg)
@@ -1783,7 +1753,21 @@ def eval(cls, arg):
                 return - log(1 + sqrt(2))
 
         if arg.is_number:
-            cst_table = cls._acsch_table()
+            cst_table = {
+                S.ImaginaryUnit: -S.Pi / 2,
+                S.ImaginaryUnit*(sqrt(2) + sqrt(6)): -S.Pi / 12,
+                S.ImaginaryUnit*(1 + sqrt(5)): -S.Pi / 10,
+                S.ImaginaryUnit*2 / sqrt(2 - sqrt(2)): -S.Pi / 8,
+                S.ImaginaryUnit*2: -S.Pi / 6,
+                S.ImaginaryUnit*sqrt(2 + 2/sqrt(5)): -S.Pi / 5,
+                S.ImaginaryUnit*sqrt(2): -S.Pi / 4,
+                S.ImaginaryUnit*(sqrt(5)-1): -3*S.Pi / 10,
+                S.ImaginaryUnit*2 / sqrt(3): -S.Pi / 3,
+                S.ImaginaryUnit*2 / sqrt(2 + sqrt(2)): -3*S.Pi / 8,
+                S.ImaginaryUnit*sqrt(2 - 2/sqrt(5)): -2*S.Pi / 5,
+                S.ImaginaryUnit*(sqrt(6) - sqrt(2)): -5*S.Pi / 12,
+                S(2): -S.ImaginaryUnit*log((1+sqrt(5))/2),
+            }
 
             if arg in cst_table:
                 return cst_table[arg]*S.ImaginaryUnit
diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py
index 15a482f159..73df754ca2 100644
--- a/sympy/functions/elementary/trigonometric.py
+++ b/sympy/functions/elementary/trigonometric.py
@@ -5,7 +5,7 @@
 from sympy.core.expr import Expr
 from sympy.core.function import Function, ArgumentIndexError, PoleError, expand_mul
 from sympy.core.logic import fuzzy_not, fuzzy_or, FuzzyBool, fuzzy_and
-from sympy.core.numbers import igcdex, Rational, pi, Integer, Float
+from sympy.core.numbers import igcdex, Rational, pi, Integer
 from sympy.core.relational import Ne
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol, Dummy
@@ -20,19 +20,6 @@
 from sympy.polys.specialpolys import symmetric_poly
 from sympy.utilities.iterables import numbered_symbols
 
-
-###############################################################################
-########################## UTILITIES ##########################################
-###############################################################################
-
-
-def _imaginary_unit_as_coefficient(arg):
-    """ Helper to extract symbolic coefficient for imaginary unit """
-    if isinstance(arg, Float):
-        return None
-    else:
-        return arg.as_coefficient(S.ImaginaryUnit)
-
 ###############################################################################
 ########################## TRIGONOMETRIC FUNCTIONS ############################
 ###############################################################################
@@ -127,7 +114,7 @@ def _peeloff_pi(arg):
     pi_coeff = S.Zero
     rest_terms = []
     for a in Add.make_args(arg):
-        K = a.coeff(pi)
+        K = a.coeff(S.Pi)
         if K and K.is_rational:
             pi_coeff += K
         else:
@@ -181,12 +168,12 @@ def _pi_coeff(arg, cycles=1):
 
     """
     arg = sympify(arg)
-    if arg is pi:
+    if arg is S.Pi:
         return S.One
     elif not arg:
         return S.Zero
     elif arg.is_Mul:
-        cx = arg.coeff(pi)
+        cx = arg.coeff(S.Pi)
         if cx:
             c, x = cx.as_coeff_Mul()  # pi is not included as coeff
             if c.is_Float:
@@ -293,20 +280,20 @@ def eval(cls, arg):
         if isinstance(arg, AccumBounds):
             from sympy.sets.sets import FiniteSet
             min, max = arg.min, arg.max
-            d = floor(min/(2*pi))
+            d = floor(min/(2*S.Pi))
             if min is not S.NegativeInfinity:
-                min = min - d*2*pi
+                min = min - d*2*S.Pi
             if max is not S.Infinity:
-                max = max - d*2*pi
-            if AccumBounds(min, max).intersection(FiniteSet(pi/2, pi*Rational(5, 2))) \
+                max = max - d*2*S.Pi
+            if AccumBounds(min, max).intersection(FiniteSet(S.Pi/2, S.Pi*Rational(5, 2))) \
                     is not S.EmptySet and \
-                    AccumBounds(min, max).intersection(FiniteSet(pi*Rational(3, 2),
-                        pi*Rational(7, 2))) is not S.EmptySet:
+                    AccumBounds(min, max).intersection(FiniteSet(S.Pi*Rational(3, 2),
+                        S.Pi*Rational(7, 2))) is not S.EmptySet:
                 return AccumBounds(-1, 1)
-            elif AccumBounds(min, max).intersection(FiniteSet(pi/2, pi*Rational(5, 2))) \
+            elif AccumBounds(min, max).intersection(FiniteSet(S.Pi/2, S.Pi*Rational(5, 2))) \
                     is not S.EmptySet:
                 return AccumBounds(Min(sin(min), sin(max)), 1)
-            elif AccumBounds(min, max).intersection(FiniteSet(pi*Rational(3, 2), pi*Rational(8, 2))) \
+            elif AccumBounds(min, max).intersection(FiniteSet(S.Pi*Rational(3, 2), S.Pi*Rational(8, 2))) \
                         is not S.EmptySet:
                 return AccumBounds(-1, Max(sin(min), sin(max)))
             else:
@@ -318,7 +305,7 @@ def eval(cls, arg):
         if arg.could_extract_minus_sign():
             return -cls(-arg)
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import sinh
             return S.ImaginaryUnit*sinh(i_coeff)
@@ -335,7 +322,7 @@ def eval(cls, arg):
                     return S.NegativeOne**(pi_coeff - S.Half)
 
             if not pi_coeff.is_Rational:
-                narg = pi_coeff*pi
+                narg = pi_coeff*S.Pi
                 if narg != arg:
                     return cls(narg)
                 return None
@@ -345,21 +332,21 @@ def eval(cls, arg):
             if pi_coeff.is_Rational:
                 x = pi_coeff % 2
                 if x > 1:
-                    return -cls((x % 1)*pi)
+                    return -cls((x % 1)*S.Pi)
                 if 2*x > 1:
-                    return cls((1 - x)*pi)
-                narg = ((pi_coeff + Rational(3, 2)) % 2)*pi
+                    return cls((1 - x)*S.Pi)
+                narg = ((pi_coeff + Rational(3, 2)) % 2)*S.Pi
                 result = cos(narg)
                 if not isinstance(result, cos):
                     return result
-                if pi_coeff*pi != arg:
-                    return cls(pi_coeff*pi)
+                if pi_coeff*S.Pi != arg:
+                    return cls(pi_coeff*S.Pi)
                 return None
 
         if arg.is_Add:
             x, m = _peeloff_pi(arg)
             if m:
-                m = m*pi
+                m = m*S.Pi
                 return sin(m)*cos(x) + cos(m)*sin(x)
 
         if arg.is_zero:
@@ -428,7 +415,7 @@ def _eval_rewrite_as_Pow(self, arg, **kwargs):
             return I*x**-I/2 - I*x**I /2
 
     def _eval_rewrite_as_cos(self, arg, **kwargs):
-        return cos(arg - pi/2, evaluate=False)
+        return cos(arg - S.Pi/2, evaluate=False)
 
     def _eval_rewrite_as_tan(self, arg, **kwargs):
         tan_half = tan(S.Half*arg)
@@ -451,7 +438,7 @@ def _eval_rewrite_as_csc(self, arg, **kwargs):
         return 1/csc(arg)
 
     def _eval_rewrite_as_sec(self, arg, **kwargs):
-        return 1/sec(arg - pi/2, evaluate=False)
+        return 1/sec(arg - S.Pi/2, evaluate=False)
 
     def _eval_rewrite_as_sinc(self, arg, **kwargs):
         return arg*sinc(arg)
@@ -497,9 +484,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         from sympy.calculus.accumulationbounds import AccumBounds
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        n = x0/pi
+        n = x0/S.Pi
         if n.is_integer:
-            lt = (arg - n*pi).as_leading_term(x)
+            lt = (arg - n*S.Pi).as_leading_term(x)
             return (S.NegativeOne**n)*lt
         if x0 is S.ComplexInfinity:
             x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
@@ -601,7 +588,7 @@ def eval(cls, arg):
             return S.NaN
 
         if isinstance(arg, AccumBounds):
-            return sin(arg + pi/2)
+            return sin(arg + S.Pi/2)
         elif isinstance(arg, SetExpr):
             return arg._eval_func(cls)
 
@@ -611,7 +598,7 @@ def eval(cls, arg):
         if arg.could_extract_minus_sign():
             return cls(-arg)
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import cosh
             return cosh(i_coeff)
@@ -628,7 +615,7 @@ def eval(cls, arg):
                     return S.Zero
 
             if not pi_coeff.is_Rational:
-                narg = pi_coeff*pi
+                narg = pi_coeff*S.Pi
                 if narg != arg:
                     return cls(narg)
                 return None
@@ -648,10 +635,10 @@ def eval(cls, arg):
                 q = pi_coeff.q
                 p = pi_coeff.p % (2*q)
                 if p > q:
-                    narg = (pi_coeff - 1)*pi
+                    narg = (pi_coeff - 1)*S.Pi
                     return -cls(narg)
                 if 2*p > q:
-                    narg = (1 - pi_coeff)*pi
+                    narg = (1 - pi_coeff)*S.Pi
                     return -cls(narg)
 
                 # If nested sqrt's are worse than un-evaluation
@@ -669,11 +656,11 @@ def eval(cls, arg):
                     120: (40, 60)
                     }
                 if q in table2:
-                    a, b = p*pi/table2[q][0], p*pi/table2[q][1]
+                    a, b = p*S.Pi/table2[q][0], p*S.Pi/table2[q][1]
                     nvala, nvalb = cls(a), cls(b)
                     if None in (nvala, nvalb):
                         return None
-                    return nvala*nvalb + cls(pi/2 - a)*cls(pi/2 - b)
+                    return nvala*nvalb + cls(S.Pi/2 - a)*cls(S.Pi/2 - b)
 
                 if q > 12:
                     return None
@@ -683,7 +670,7 @@ def eval(cls, arg):
                     return chebyshevt(pi_coeff.p, cts).expand()
 
                 if 0 == q % 2:
-                    narg = (pi_coeff*2)*pi
+                    narg = (pi_coeff*2)*S.Pi
                     nval = cls(narg)
                     if None == nval:
                         return None
@@ -695,7 +682,7 @@ def eval(cls, arg):
         if arg.is_Add:
             x, m = _peeloff_pi(arg)
             if m:
-                m = m*pi
+                m = m*S.Pi
                 return cos(m)*cos(x) - sin(m)*sin(x)
 
         if arg.is_zero:
@@ -764,7 +751,7 @@ def _eval_rewrite_as_Pow(self, arg, **kwargs):
             return x**I/2 + x**-I/2
 
     def _eval_rewrite_as_sin(self, arg, **kwargs):
-        return sin(arg + pi/2, evaluate=False)
+        return sin(arg + S.Pi/2, evaluate=False)
 
     def _eval_rewrite_as_tan(self, arg, **kwargs):
         tan_half = tan(S.Half*arg)**2
@@ -823,7 +810,7 @@ def ipartfrac(r, factors=None):
 
         if pi_coeff.is_integer:
             # it was unevaluated
-            return self.func(pi_coeff*pi)
+            return self.func(pi_coeff*S.Pi)
 
         if not pi_coeff.is_Rational:
             return None
@@ -900,7 +887,7 @@ def _fermatCoords(n):
 
         if not pi_coeff.q % 2:  # recursively remove factors of 2
             pico2 = pi_coeff*2
-            nval = cos(pico2*pi).rewrite(sqrt)
+            nval = cos(pico2*S.Pi).rewrite(sqrt)
             x = (pico2 + 1)/2
             sign_cos = -1 if int(x) % 2 else 1
             return sign_cos*sqrt( (1 + nval)/2 )
@@ -908,12 +895,12 @@ def _fermatCoords(n):
         FC = _fermatCoords(pi_coeff.q)
         if FC:
             decomp = ipartfrac(pi_coeff, FC)
-            X = [(x[1], x[0]*pi) for x in zip(decomp, numbered_symbols('z'))]
+            X = [(x[1], x[0]*S.Pi) for x in zip(decomp, numbered_symbols('z'))]
             pcls = cos(sum([x[0] for x in X]))._eval_expand_trig().subs(X)
             return pcls.rewrite(sqrt)
         else:
             decomp = ipartfrac(pi_coeff)
-            X = [(x[1], x[0]*pi) for x in zip(decomp, numbered_symbols('z'))]
+            X = [(x[1], x[0]*S.Pi) for x in zip(decomp, numbered_symbols('z'))]
             pcls = cos(sum([x[0] for x in X]))._eval_expand_trig().subs(X)
             return pcls
 
@@ -956,9 +943,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         from sympy.calculus.accumulationbounds import AccumBounds
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        n = (x0 + pi/2)/pi
+        n = (x0 + S.Pi/2)/S.Pi
         if n.is_integer:
-            lt = (arg - n*pi + pi/2).as_leading_term(x)
+            lt = (arg - n*S.Pi + S.Pi/2).as_leading_term(x)
             return (S.NegativeOne**n)*lt
         if x0 is S.ComplexInfinity:
             x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
@@ -1058,13 +1045,13 @@ def eval(cls, arg):
 
         if isinstance(arg, AccumBounds):
             min, max = arg.min, arg.max
-            d = floor(min/pi)
+            d = floor(min/S.Pi)
             if min is not S.NegativeInfinity:
-                min = min - d*pi
+                min = min - d*S.Pi
             if max is not S.Infinity:
-                max = max - d*pi
+                max = max - d*S.Pi
             from sympy.sets.sets import FiniteSet
-            if AccumBounds(min, max).intersection(FiniteSet(pi/2, pi*Rational(3, 2))):
+            if AccumBounds(min, max).intersection(FiniteSet(S.Pi/2, S.Pi*Rational(3, 2))):
                 return AccumBounds(S.NegativeInfinity, S.Infinity)
             else:
                 return AccumBounds(tan(min), tan(max))
@@ -1072,7 +1059,7 @@ def eval(cls, arg):
         if arg.could_extract_minus_sign():
             return -cls(-arg)
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import tanh
             return S.ImaginaryUnit*tanh(i_coeff)
@@ -1083,7 +1070,7 @@ def eval(cls, arg):
                 return S.Zero
 
             if not pi_coeff.is_Rational:
-                narg = pi_coeff*pi
+                narg = pi_coeff*S.Pi
                 if narg != arg:
                     return cls(narg)
                 return None
@@ -1106,8 +1093,8 @@ def eval(cls, arg):
                     else:
                         return table10[n]
                 if not pi_coeff.q % 2:
-                    narg = pi_coeff*pi*2
-                    cresult, sresult = cos(narg), cos(narg - pi/2)
+                    narg = pi_coeff*S.Pi*2
+                    cresult, sresult = cos(narg), cos(narg - S.Pi/2)
                     if not isinstance(cresult, cos) \
                             and not isinstance(sresult, cos):
                         if sresult == 0:
@@ -1124,14 +1111,14 @@ def eval(cls, arg):
                     120: (40, 60)
                     }
                 if q in table2:
-                    nvala, nvalb = cls(p*pi/table2[q][0]), cls(p*pi/table2[q][1])
+                    nvala, nvalb = cls(p*S.Pi/table2[q][0]), cls(p*S.Pi/table2[q][1])
                     if None in (nvala, nvalb):
                         return None
                     return (nvala - nvalb)/(1 + nvala*nvalb)
-                narg = ((pi_coeff + S.Half) % 1 - S.Half)*pi
+                narg = ((pi_coeff + S.Half) % 1 - S.Half)*S.Pi
                 # see cos() to specify which expressions should  be
                 # expanded automatically in terms of radicals
-                cresult, sresult = cos(narg), cos(narg - pi/2)
+                cresult, sresult = cos(narg), cos(narg - S.Pi/2)
                 if not isinstance(cresult, cos) \
                         and not isinstance(sresult, cos):
                     if cresult == 0:
@@ -1143,7 +1130,7 @@ def eval(cls, arg):
         if arg.is_Add:
             x, m = _peeloff_pi(arg)
             if m:
-                tanm = tan(m*pi)
+                tanm = tan(m*S.Pi)
                 if tanm is S.ComplexInfinity:
                     return -cot(x)
                 else: # tanm == 0
@@ -1195,7 +1182,7 @@ def taylor_term(n, x, *previous_terms):
             return S.NegativeOne**a*b*(b - 1)*B/F*x**n
 
     def _eval_nseries(self, x, n, logx, cdir=0):
-        i = self.args[0].limit(x, 0)*2/pi
+        i = self.args[0].limit(x, 0)*2/S.Pi
         if i and i.is_Integer:
             return self.rewrite(cos)._eval_nseries(x, n=n, logx=logx)
         return Function._eval_nseries(self, x, n=n, logx=logx)
@@ -1257,7 +1244,7 @@ def _eval_rewrite_as_sin(self, x, **kwargs):
         return 2*sin(x)**2/sin(2*x)
 
     def _eval_rewrite_as_cos(self, x, **kwargs):
-        return cos(x - pi/2, evaluate=False)/cos(x)
+        return cos(x - S.Pi/2, evaluate=False)/cos(x)
 
     def _eval_rewrite_as_sincos(self, arg, **kwargs):
         return sin(arg)/cos(arg)
@@ -1292,9 +1279,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         from sympy.functions.elementary.complexes import re
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        n = 2*x0/pi
+        n = 2*x0/S.Pi
         if n.is_integer:
-            lt = (arg - n*pi/2).as_leading_term(x)
+            lt = (arg - n*S.Pi/2).as_leading_term(x)
             return lt if n.is_even else -1/lt
         if x0 is S.ComplexInfinity:
             x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
@@ -1400,12 +1387,12 @@ def eval(cls, arg):
             return S.NaN
 
         if isinstance(arg, AccumBounds):
-            return -tan(arg + pi/2)
+            return -tan(arg + S.Pi/2)
 
         if arg.could_extract_minus_sign():
             return -cls(-arg)
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import coth
             return -S.ImaginaryUnit*coth(i_coeff)
@@ -1416,17 +1403,17 @@ def eval(cls, arg):
                 return S.ComplexInfinity
 
             if not pi_coeff.is_Rational:
-                narg = pi_coeff*pi
+                narg = pi_coeff*S.Pi
                 if narg != arg:
                     return cls(narg)
                 return None
 
             if pi_coeff.is_Rational:
                 if pi_coeff.q in (5, 10):
-                    return tan(pi/2 - arg)
+                    return tan(S.Pi/2 - arg)
                 if pi_coeff.q > 2 and not pi_coeff.q % 2:
-                    narg = pi_coeff*pi*2
-                    cresult, sresult = cos(narg), cos(narg - pi/2)
+                    narg = pi_coeff*S.Pi*2
+                    cresult, sresult = cos(narg), cos(narg - S.Pi/2)
                     if not isinstance(cresult, cos) \
                             and not isinstance(sresult, cos):
                         return 1/sresult + cresult/sresult
@@ -1443,14 +1430,14 @@ def eval(cls, arg):
                 q = pi_coeff.q
                 p = pi_coeff.p % q
                 if q in table2:
-                    nvala, nvalb = cls(p*pi/table2[q][0]), cls(p*pi/table2[q][1])
+                    nvala, nvalb = cls(p*S.Pi/table2[q][0]), cls(p*S.Pi/table2[q][1])
                     if None in (nvala, nvalb):
                         return None
                     return (1 + nvala*nvalb)/(nvalb - nvala)
-                narg = (((pi_coeff + S.Half) % 1) - S.Half)*pi
+                narg = (((pi_coeff + S.Half) % 1) - S.Half)*S.Pi
                 # see cos() to specify which expressions should be
                 # expanded automatically in terms of radicals
-                cresult, sresult = cos(narg), cos(narg - pi/2)
+                cresult, sresult = cos(narg), cos(narg - S.Pi/2)
                 if not isinstance(cresult, cos) \
                         and not isinstance(sresult, cos):
                     if sresult == 0:
@@ -1462,7 +1449,7 @@ def eval(cls, arg):
         if arg.is_Add:
             x, m = _peeloff_pi(arg)
             if m:
-                cotm = cot(m*pi)
+                cotm = cot(m*S.Pi)
                 if cotm is S.ComplexInfinity:
                     return cot(x)
                 else: # cotm == 0
@@ -1514,7 +1501,7 @@ def taylor_term(n, x, *previous_terms):
             return S.NegativeOne**((n + 1)//2)*2**(n + 1)*B/F*x**n
 
     def _eval_nseries(self, x, n, logx, cdir=0):
-        i = self.args[0].limit(x, 0)/pi
+        i = self.args[0].limit(x, 0)/S.Pi
         if i and i.is_Integer:
             return self.rewrite(cos)._eval_nseries(x, n=n, logx=logx)
         return self.rewrite(tan)._eval_nseries(x, n=n, logx=logx)
@@ -1549,7 +1536,7 @@ def _eval_rewrite_as_sin(self, x, **kwargs):
         return sin(2*x)/(2*(sin(x)**2))
 
     def _eval_rewrite_as_cos(self, x, **kwargs):
-        return cos(x)/cos(x - pi/2, evaluate=False)
+        return cos(x)/cos(x - S.Pi/2, evaluate=False)
 
     def _eval_rewrite_as_sincos(self, arg, **kwargs):
         return cos(arg)/sin(arg)
@@ -1584,9 +1571,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         from sympy.functions.elementary.complexes import re
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        n = 2*x0/pi
+        n = 2*x0/S.Pi
         if n.is_integer:
-            lt = (arg - n*pi/2).as_leading_term(x)
+            lt = (arg - n*S.Pi/2).as_leading_term(x)
             return 1/lt if n.is_even else -lt
         if x0 is S.ComplexInfinity:
             x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
@@ -1648,7 +1635,7 @@ def _eval_is_zero(self):
     def _eval_subs(self, old, new):
         arg = self.args[0]
         argnew = arg.subs(old, new)
-        if arg != argnew and (argnew/pi).is_integer:
+        if arg != argnew and (argnew/S.Pi).is_integer:
             return S.ComplexInfinity
         return cot(argnew)
 
@@ -1682,10 +1669,10 @@ def eval(cls, arg):
                 q = pi_coeff.q
                 p = pi_coeff.p % (2*q)
                 if p > q:
-                    narg = (pi_coeff - 1)*pi
+                    narg = (pi_coeff - 1)*S.Pi
                     return -cls(narg)
                 if 2*p > q:
-                    narg = (1 - pi_coeff)*pi
+                    narg = (1 - pi_coeff)*S.Pi
                     if cls._is_odd:
                         return cls(narg)
                     elif cls._is_even:
@@ -1863,9 +1850,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         from sympy.functions.elementary.complexes import re
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        n = (x0 + pi/2)/pi
+        n = (x0 + S.Pi/2)/S.Pi
         if n.is_integer:
-            lt = (arg - n*pi + pi/2).as_leading_term(x)
+            lt = (arg - n*S.Pi + S.Pi/2).as_leading_term(x)
             return (S.NegativeOne**n)/lt
         if x0 is S.ComplexInfinity:
             x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
@@ -1964,9 +1951,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         from sympy.functions.elementary.complexes import re
         arg = self.args[0]
         x0 = arg.subs(x, 0).cancel()
-        n = x0/pi
+        n = x0/S.Pi
         if n.is_integer:
-            lt = (arg - n*pi).as_leading_term(x)
+            lt = (arg - n*S.Pi).as_leading_term(x)
             return (S.NegativeOne**n)/lt
         if x0 is S.ComplexInfinity:
             x0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
@@ -2104,79 +2091,75 @@ class InverseTrigonometricFunction(Function):
     _singularities = (S.One, S.NegativeOne, S.Zero, S.ComplexInfinity)  # type: tTuple[Expr, ...]
 
     @staticmethod
-    @cacheit
     def _asin_table():
         # Only keys with could_extract_minus_sign() == False
         # are actually needed.
         return {
-            sqrt(3)/2: pi/3,
-            sqrt(2)/2: pi/4,
-            1/sqrt(2): pi/4,
-            sqrt((5 - sqrt(5))/8): pi/5,
-            sqrt(2)*sqrt(5 - sqrt(5))/4: pi/5,
-            sqrt((5 + sqrt(5))/8): pi*Rational(2, 5),
-            sqrt(2)*sqrt(5 + sqrt(5))/4: pi*Rational(2, 5),
-            S.Half: pi/6,
-            sqrt(2 - sqrt(2))/2: pi/8,
-            sqrt(S.Half - sqrt(2)/4): pi/8,
-            sqrt(2 + sqrt(2))/2: pi*Rational(3, 8),
-            sqrt(S.Half + sqrt(2)/4): pi*Rational(3, 8),
-            (sqrt(5) - 1)/4: pi/10,
-            (1 - sqrt(5))/4: -pi/10,
-            (sqrt(5) + 1)/4: pi*Rational(3, 10),
-            sqrt(6)/4 - sqrt(2)/4: pi/12,
-            -sqrt(6)/4 + sqrt(2)/4: -pi/12,
-            (sqrt(3) - 1)/sqrt(8): pi/12,
-            (1 - sqrt(3))/sqrt(8): -pi/12,
-            sqrt(6)/4 + sqrt(2)/4: pi*Rational(5, 12),
-            (1 + sqrt(3))/sqrt(8): pi*Rational(5, 12)
+            sqrt(3)/2: S.Pi/3,
+            sqrt(2)/2: S.Pi/4,
+            1/sqrt(2): S.Pi/4,
+            sqrt((5 - sqrt(5))/8): S.Pi/5,
+            sqrt(2)*sqrt(5 - sqrt(5))/4: S.Pi/5,
+            sqrt((5 + sqrt(5))/8): S.Pi*Rational(2, 5),
+            sqrt(2)*sqrt(5 + sqrt(5))/4: S.Pi*Rational(2, 5),
+            S.Half: S.Pi/6,
+            sqrt(2 - sqrt(2))/2: S.Pi/8,
+            sqrt(S.Half - sqrt(2)/4): S.Pi/8,
+            sqrt(2 + sqrt(2))/2: S.Pi*Rational(3, 8),
+            sqrt(S.Half + sqrt(2)/4): S.Pi*Rational(3, 8),
+            (sqrt(5) - 1)/4: S.Pi/10,
+            (1 - sqrt(5))/4: -S.Pi/10,
+            (sqrt(5) + 1)/4: S.Pi*Rational(3, 10),
+            sqrt(6)/4 - sqrt(2)/4: S.Pi/12,
+            -sqrt(6)/4 + sqrt(2)/4: -S.Pi/12,
+            (sqrt(3) - 1)/sqrt(8): S.Pi/12,
+            (1 - sqrt(3))/sqrt(8): -S.Pi/12,
+            sqrt(6)/4 + sqrt(2)/4: S.Pi*Rational(5, 12),
+            (1 + sqrt(3))/sqrt(8): S.Pi*Rational(5, 12)
         }
 
-
     @staticmethod
-    @cacheit
     def _atan_table():
         # Only keys with could_extract_minus_sign() == False
         # are actually needed.
         return {
-            sqrt(3)/3: pi/6,
-            1/sqrt(3): pi/6,
-            sqrt(3): pi/3,
-            sqrt(2) - 1: pi/8,
-            1 - sqrt(2): -pi/8,
-            1 + sqrt(2): pi*Rational(3, 8),
-            sqrt(5 - 2*sqrt(5)): pi/5,
-            sqrt(5 + 2*sqrt(5)): pi*Rational(2, 5),
-            sqrt(1 - 2*sqrt(5)/5): pi/10,
-            sqrt(1 + 2*sqrt(5)/5): pi*Rational(3, 10),
-            2 - sqrt(3): pi/12,
-            -2 + sqrt(3): -pi/12,
-            2 + sqrt(3): pi*Rational(5, 12)
+            sqrt(3)/3: S.Pi/6,
+            1/sqrt(3): S.Pi/6,
+            sqrt(3): S.Pi/3,
+            sqrt(2) - 1: S.Pi/8,
+            1 - sqrt(2): -S.Pi/8,
+            1 + sqrt(2): S.Pi*Rational(3, 8),
+            sqrt(5 - 2*sqrt(5)): S.Pi/5,
+            sqrt(5 + 2*sqrt(5)): S.Pi*Rational(2, 5),
+            sqrt(1 - 2*sqrt(5)/5): S.Pi/10,
+            sqrt(1 + 2*sqrt(5)/5): S.Pi*Rational(3, 10),
+            2 - sqrt(3): S.Pi/12,
+            -2 + sqrt(3): -S.Pi/12,
+            2 + sqrt(3): S.Pi*Rational(5, 12)
         }
 
     @staticmethod
-    @cacheit
     def _acsc_table():
         # Keys for which could_extract_minus_sign()
         # will obviously return True are omitted.
         return {
-            2*sqrt(3)/3: pi/3,
-            sqrt(2): pi/4,
-            sqrt(2 + 2*sqrt(5)/5): pi/5,
-            1/sqrt(Rational(5, 8) - sqrt(5)/8): pi/5,
-            sqrt(2 - 2*sqrt(5)/5): pi*Rational(2, 5),
-            1/sqrt(Rational(5, 8) + sqrt(5)/8): pi*Rational(2, 5),
-            2: pi/6,
-            sqrt(4 + 2*sqrt(2)): pi/8,
-            2/sqrt(2 - sqrt(2)): pi/8,
-            sqrt(4 - 2*sqrt(2)): pi*Rational(3, 8),
-            2/sqrt(2 + sqrt(2)): pi*Rational(3, 8),
-            1 + sqrt(5): pi/10,
-            sqrt(5) - 1: pi*Rational(3, 10),
-            -(sqrt(5) - 1): pi*Rational(-3, 10),
-            sqrt(6) + sqrt(2): pi/12,
-            sqrt(6) - sqrt(2): pi*Rational(5, 12),
-            -(sqrt(6) - sqrt(2)): pi*Rational(-5, 12)
+            2*sqrt(3)/3: S.Pi/3,
+            sqrt(2): S.Pi/4,
+            sqrt(2 + 2*sqrt(5)/5): S.Pi/5,
+            1/sqrt(Rational(5, 8) - sqrt(5)/8): S.Pi/5,
+            sqrt(2 - 2*sqrt(5)/5): S.Pi*Rational(2, 5),
+            1/sqrt(Rational(5, 8) + sqrt(5)/8): S.Pi*Rational(2, 5),
+            2: S.Pi/6,
+            sqrt(4 + 2*sqrt(2)): S.Pi/8,
+            2/sqrt(2 - sqrt(2)): S.Pi/8,
+            sqrt(4 - 2*sqrt(2)): S.Pi*Rational(3, 8),
+            2/sqrt(2 + sqrt(2)): S.Pi*Rational(3, 8),
+            1 + sqrt(5): S.Pi/10,
+            sqrt(5) - 1: S.Pi*Rational(3, 10),
+            -(sqrt(5) - 1): S.Pi*Rational(-3, 10),
+            sqrt(6) + sqrt(2): S.Pi/12,
+            sqrt(6) - sqrt(2): S.Pi*Rational(5, 12),
+            -(sqrt(6) - sqrt(2)): S.Pi*Rational(-5, 12)
         }
 
 
@@ -2255,9 +2238,9 @@ def eval(cls, arg):
             elif arg.is_zero:
                 return S.Zero
             elif arg is S.One:
-                return pi/2
+                return S.Pi/2
             elif arg is S.NegativeOne:
-                return -pi/2
+                return -S.Pi/2
 
         if arg is S.ComplexInfinity:
             return S.ComplexInfinity
@@ -2270,7 +2253,7 @@ def eval(cls, arg):
             if arg in asin_table:
                 return asin_table[arg]
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import asinh
             return S.ImaginaryUnit*asinh(i_coeff)
@@ -2324,9 +2307,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 < S.NegativeOne:
-            return -pi - self.func(x0)
+            return -S.Pi - self.func(x0)
         elif im(cdir) > 0 and x0.is_real and x0 > S.One:
-            return pi - self.func(x0)
+            return S.Pi - self.func(x0)
         return self.func(x0)
 
     def _eval_nseries(self, x, n, logx, cdir=0):  # asin
@@ -2339,7 +2322,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # asin
             f = arg1.as_leading_term(x)
             g = (arg1 - f)/ f
             if not g.is_meromorphic(x, 0):   # cannot be expanded
-                return O(1) if n == 0 else pi/2 + O(sqrt(x))
+                return O(1) if n == 0 else S.Pi/2 + O(sqrt(x))
             res1 = sqrt(S.One + g)._eval_nseries(x, n=n, logx=logx)
             res = (res1.removeO()*sqrt(f)).expand()
             return ser.removeO().subs(t, res).expand().powsimp() + O(x**n, x)
@@ -2351,7 +2334,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # asin
             f = arg1.as_leading_term(x)
             g = (arg1 - f)/ f
             if not g.is_meromorphic(x, 0):   # cannot be expanded
-                return O(1) if n == 0 else -pi/2 + O(sqrt(x))
+                return O(1) if n == 0 else -S.Pi/2 + O(sqrt(x))
             res1 = sqrt(S.One + g)._eval_nseries(x, n=n, logx=logx)
             res = (res1.removeO()*sqrt(f)).expand()
             return ser.removeO().subs(t, res).expand().powsimp() + O(x**n, x)
@@ -2362,13 +2345,13 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # asin
         if cdir != 0:
             cdir = self.args[0].dir(x, cdir)
         if im(cdir) < 0 and arg0.is_real and arg0 < S.NegativeOne:
-            return -pi - res
+            return -S.Pi - res
         elif im(cdir) > 0 and arg0.is_real and arg0 > S.One:
-            return pi - res
+            return S.Pi - res
         return res
 
     def _eval_rewrite_as_acos(self, x, **kwargs):
-        return pi/2 - acos(x)
+        return S.Pi/2 - acos(x)
 
     def _eval_rewrite_as_atan(self, x, **kwargs):
         return 2*atan(x/(1 + sqrt(1 - x**2)))
@@ -2380,7 +2363,7 @@ def _eval_rewrite_as_acot(self, arg, **kwargs):
         return 2*acot((1 + sqrt(1 - arg**2))/arg)
 
     def _eval_rewrite_as_asec(self, arg, **kwargs):
-        return pi/2 - asec(1/arg)
+        return S.Pi/2 - asec(1/arg)
 
     def _eval_rewrite_as_acsc(self, arg, **kwargs):
         return acsc(1/arg)
@@ -2464,11 +2447,11 @@ def eval(cls, arg):
             elif arg is S.NegativeInfinity:
                 return S.NegativeInfinity*S.ImaginaryUnit
             elif arg.is_zero:
-                return pi/2
+                return S.Pi/2
             elif arg is S.One:
                 return S.Zero
             elif arg is S.NegativeOne:
-                return pi
+                return S.Pi
 
         if arg is S.ComplexInfinity:
             return S.ComplexInfinity
@@ -2480,7 +2463,7 @@ def eval(cls, arg):
             elif -arg in asin_table:
                 return pi/2 + asin_table[-arg]
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             return pi/2 - asin(arg)
 
@@ -2502,7 +2485,7 @@ def eval(cls, arg):
     @cacheit
     def taylor_term(n, x, *previous_terms):
         if n == 0:
-            return pi/2
+            return S.Pi/2
         elif n < 0 or n % 2 == 0:
             return S.Zero
         else:
@@ -2526,7 +2509,7 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 < S.NegativeOne:
-            return 2*pi - self.func(x0)
+            return 2*S.Pi - self.func(x0)
         elif im(cdir) > 0 and x0.is_real and x0 > S.One:
             return -self.func(x0)
         return self.func(x0)
@@ -2560,7 +2543,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # acos
             f = arg1.as_leading_term(x)
             g = (arg1 - f)/ f
             if not g.is_meromorphic(x, 0):   # cannot be expanded
-                return O(1) if n == 0 else pi + O(sqrt(x))
+                return O(1) if n == 0 else S.Pi + O(sqrt(x))
             res1 = sqrt(S.One + g)._eval_nseries(x, n=n, logx=logx)
             res = (res1.removeO()*sqrt(f)).expand()
             return ser.removeO().subs(t, res).expand().powsimp() + O(x**n, x)
@@ -2571,20 +2554,20 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # acos
         if cdir != 0:
             cdir = self.args[0].dir(x, cdir)
         if im(cdir) < 0 and arg0.is_real and arg0 < S.NegativeOne:
-            return 2*pi - res
+            return 2*S.Pi - res
         elif im(cdir) > 0 and arg0.is_real and arg0 > S.One:
             return -res
         return res
 
     def _eval_rewrite_as_log(self, x, **kwargs):
-        return pi/2 + S.ImaginaryUnit*\
+        return S.Pi/2 + S.ImaginaryUnit*\
             log(S.ImaginaryUnit*x + sqrt(1 - x**2))
 
     def _eval_rewrite_as_asin(self, x, **kwargs):
-        return pi/2 - asin(x)
+        return S.Pi/2 - asin(x)
 
     def _eval_rewrite_as_atan(self, x, **kwargs):
-        return atan(sqrt(1 - x**2)/x) + (pi/2)*(1 - x*sqrt(1/x**2))
+        return atan(sqrt(1 - x**2)/x) + (S.Pi/2)*(1 - x*sqrt(1/x**2))
 
     def inverse(self, argindex=1):
         """
@@ -2593,13 +2576,13 @@ def inverse(self, argindex=1):
         return cos
 
     def _eval_rewrite_as_acot(self, arg, **kwargs):
-        return pi/2 - 2*acot((1 + sqrt(1 - arg**2))/arg)
+        return S.Pi/2 - 2*acot((1 + sqrt(1 - arg**2))/arg)
 
     def _eval_rewrite_as_asec(self, arg, **kwargs):
         return asec(1/arg)
 
     def _eval_rewrite_as_acsc(self, arg, **kwargs):
-        return pi/2 - acsc(1/arg)
+        return S.Pi/2 - acsc(1/arg)
 
     def _eval_conjugate(self):
         z = self.args[0]
@@ -2685,19 +2668,19 @@ def eval(cls, arg):
             if arg is S.NaN:
                 return S.NaN
             elif arg is S.Infinity:
-                return pi/2
+                return S.Pi/2
             elif arg is S.NegativeInfinity:
-                return -pi/2
+                return -S.Pi/2
             elif arg.is_zero:
                 return S.Zero
             elif arg is S.One:
-                return pi/4
+                return S.Pi/4
             elif arg is S.NegativeOne:
-                return -pi/4
+                return -S.Pi/4
 
         if arg is S.ComplexInfinity:
             from sympy.calculus.accumulationbounds import AccumBounds
-            return AccumBounds(-pi/2, pi/2)
+            return AccumBounds(-S.Pi/2, S.Pi/2)
 
         if arg.could_extract_minus_sign():
             return -cls(-arg)
@@ -2707,7 +2690,7 @@ def eval(cls, arg):
             if arg in atan_table:
                 return atan_table[arg]
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import atanh
             return S.ImaginaryUnit*atanh(i_coeff)
@@ -2751,9 +2734,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if re(cdir) < 0 and re(x0).is_zero and im(x0) > S.One:
-            return self.func(x0) - pi
+            return self.func(x0) - S.Pi
         elif re(cdir) > 0 and re(x0).is_zero and im(x0) < S.NegativeOne:
-            return self.func(x0) + pi
+            return self.func(x0) + S.Pi
         return self.func(x0)
 
     def _eval_nseries(self, x, n, logx, cdir=0):  # atan
@@ -2763,12 +2746,12 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # atan
             cdir = self.args[0].dir(x, cdir)
         if arg0 is S.ComplexInfinity:
             if re(cdir) > 0:
-                return res - pi
+                return res - S.Pi
             return res
         if re(cdir) < 0 and re(arg0).is_zero and im(arg0) > S.One:
-            return res - pi
+            return res - S.Pi
         elif re(cdir) > 0 and re(arg0).is_zero and im(arg0) < S.NegativeOne:
-            return res + pi
+            return res + S.Pi
         return res
 
     def _eval_rewrite_as_log(self, x, **kwargs):
@@ -2777,9 +2760,9 @@ def _eval_rewrite_as_log(self, x, **kwargs):
 
     def _eval_aseries(self, n, args0, x, logx):
         if args0[0] is S.Infinity:
-            return (pi/2 - atan(1/self.args[0]))._eval_nseries(x, n, logx)
+            return (S.Pi/2 - atan(1/self.args[0]))._eval_nseries(x, n, logx)
         elif args0[0] is S.NegativeInfinity:
-            return (-pi/2 - atan(1/self.args[0]))._eval_nseries(x, n, logx)
+            return (-S.Pi/2 - atan(1/self.args[0]))._eval_nseries(x, n, logx)
         else:
             return super()._eval_aseries(n, args0, x, logx)
 
@@ -2790,7 +2773,7 @@ def inverse(self, argindex=1):
         return tan
 
     def _eval_rewrite_as_asin(self, arg, **kwargs):
-        return sqrt(arg**2)/arg*(pi/2 - asin(1/sqrt(1 + arg**2)))
+        return sqrt(arg**2)/arg*(S.Pi/2 - asin(1/sqrt(1 + arg**2)))
 
     def _eval_rewrite_as_acos(self, arg, **kwargs):
         return sqrt(arg**2)/arg*acos(1/sqrt(1 + arg**2))
@@ -2802,7 +2785,7 @@ def _eval_rewrite_as_asec(self, arg, **kwargs):
         return sqrt(arg**2)/arg*asec(sqrt(1 + arg**2))
 
     def _eval_rewrite_as_acsc(self, arg, **kwargs):
-        return sqrt(arg**2)/arg*(pi/2 - acsc(sqrt(1 + arg**2)))
+        return sqrt(arg**2)/arg*(S.Pi/2 - acsc(sqrt(1 + arg**2)))
 
 
 class acot(InverseTrigonometricFunction):
@@ -2883,11 +2866,11 @@ def eval(cls, arg):
             elif arg is S.NegativeInfinity:
                 return S.Zero
             elif arg.is_zero:
-                return pi/ 2
+                return S.Pi/ 2
             elif arg is S.One:
-                return pi/4
+                return S.Pi/4
             elif arg is S.NegativeOne:
-                return -pi/4
+                return -S.Pi/4
 
         if arg is S.ComplexInfinity:
             return S.Zero
@@ -2903,13 +2886,13 @@ def eval(cls, arg):
                     ang -= pi
                 return ang
 
-        i_coeff = _imaginary_unit_as_coefficient(arg)
+        i_coeff = arg.as_coefficient(S.ImaginaryUnit)
         if i_coeff is not None:
             from sympy.functions.elementary.hyperbolic import acoth
             return -S.ImaginaryUnit*acoth(i_coeff)
 
         if arg.is_zero:
-            return pi*S.Half
+            return S.Pi*S.Half
 
         if isinstance(arg, cot):
             ang = arg.args[0]
@@ -2931,7 +2914,7 @@ def eval(cls, arg):
     @cacheit
     def taylor_term(n, x, *previous_terms):
         if n == 0:
-            return pi/2  # FIX THIS
+            return S.Pi/2  # FIX THIS
         elif n < 0 or n % 2 == 0:
             return S.Zero
         else:
@@ -2947,12 +2930,12 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             cdir = arg.dir(x, cdir)
         if x0.is_zero:
             if re(cdir) < 0:
-                return self.func(x0) - pi
+                return self.func(x0) - S.Pi
             return self.func(x0)
         if re(cdir) > 0 and re(x0).is_zero and im(x0) > S.Zero and im(x0) < S.One:
-            return self.func(x0) + pi
+            return self.func(x0) + S.Pi
         if re(cdir) < 0 and re(x0).is_zero and im(x0) < S.Zero and im(x0) > S.NegativeOne:
-            return self.func(x0) - pi
+            return self.func(x0) - S.Pi
         return self.func(x0)
 
     def _eval_nseries(self, x, n, logx, cdir=0):  # acot
@@ -2964,19 +2947,19 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # acot
             cdir = self.args[0].dir(x, cdir)
         if arg0.is_zero:
             if re(cdir) < 0:
-                return res - pi
+                return res - S.Pi
             return res
         if re(cdir) > 0 and re(arg0).is_zero and im(arg0) > S.Zero and im(arg0) < S.One:
-            return res + pi
+            return res + S.Pi
         if re(cdir) < 0 and re(arg0).is_zero and im(arg0) < S.Zero and im(arg0) > S.NegativeOne:
-            return res - pi
+            return res - S.Pi
         return res
 
     def _eval_aseries(self, n, args0, x, logx):
         if args0[0] is S.Infinity:
-            return (pi/2 - acot(1/self.args[0]))._eval_nseries(x, n, logx)
+            return (S.Pi/2 - acot(1/self.args[0]))._eval_nseries(x, n, logx)
         elif args0[0] is S.NegativeInfinity:
-            return (pi*Rational(3, 2) - acot(1/self.args[0]))._eval_nseries(x, n, logx)
+            return (S.Pi*Rational(3, 2) - acot(1/self.args[0]))._eval_nseries(x, n, logx)
         else:
             return super(atan, self)._eval_aseries(n, args0, x, logx)
 
@@ -2992,7 +2975,7 @@ def inverse(self, argindex=1):
 
     def _eval_rewrite_as_asin(self, arg, **kwargs):
         return (arg*sqrt(1/arg**2)*
-                (pi/2 - asin(sqrt(-arg**2)/sqrt(-arg**2 - 1))))
+                (S.Pi/2 - asin(sqrt(-arg**2)/sqrt(-arg**2 - 1))))
 
     def _eval_rewrite_as_acos(self, arg, **kwargs):
         return arg*sqrt(1/arg**2)*acos(sqrt(-arg**2)/sqrt(-arg**2 - 1))
@@ -3004,7 +2987,7 @@ def _eval_rewrite_as_asec(self, arg, **kwargs):
         return arg*sqrt(1/arg**2)*asec(sqrt((1 + arg**2)/arg**2))
 
     def _eval_rewrite_as_acsc(self, arg, **kwargs):
-        return arg*sqrt(1/arg**2)*(pi/2 - acsc(sqrt((1 + arg**2)/arg**2)))
+        return arg*sqrt(1/arg**2)*(S.Pi/2 - acsc(sqrt((1 + arg**2)/arg**2)))
 
 
 class asec(InverseTrigonometricFunction):
@@ -3077,9 +3060,9 @@ def eval(cls, arg):
             elif arg is S.One:
                 return S.Zero
             elif arg is S.NegativeOne:
-                return pi
+                return S.Pi
         if arg in [S.Infinity, S.NegativeInfinity, S.ComplexInfinity]:
-            return pi/2
+            return S.Pi/2
 
         if arg.is_number:
             acsc_table = cls._acsc_table()
@@ -3126,7 +3109,7 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if im(cdir) < 0 and x0.is_real and x0 > S.Zero and x0 < S.One:
             return -self.func(x0)
         elif im(cdir) > 0 and x0.is_real and x0 < S.Zero and x0 > S.NegativeOne:
-            return 2*pi - self.func(x0)
+            return 2*S.Pi - self.func(x0)
         return self.func(x0)
 
     def _eval_nseries(self, x, n, logx, cdir=0):  # asec
@@ -3160,7 +3143,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # asec
         if im(cdir) < 0 and arg0.is_real and arg0 > S.Zero and arg0 < S.One:
             return -res
         elif im(cdir) > 0 and arg0.is_real and arg0 < S.Zero and arg0 > S.NegativeOne:
-            return 2*pi - res
+            return 2*S.Pi - res
         return res
 
     def _eval_is_extended_real(self):
@@ -3170,22 +3153,22 @@ def _eval_is_extended_real(self):
         return fuzzy_or(((x - 1).is_nonnegative, (-x - 1).is_nonnegative))
 
     def _eval_rewrite_as_log(self, arg, **kwargs):
-        return pi/2 + S.ImaginaryUnit*log(S.ImaginaryUnit/arg + sqrt(1 - 1/arg**2))
+        return S.Pi/2 + S.ImaginaryUnit*log(S.ImaginaryUnit/arg + sqrt(1 - 1/arg**2))
 
     def _eval_rewrite_as_asin(self, arg, **kwargs):
-        return pi/2 - asin(1/arg)
+        return S.Pi/2 - asin(1/arg)
 
     def _eval_rewrite_as_acos(self, arg, **kwargs):
         return acos(1/arg)
 
     def _eval_rewrite_as_atan(self, arg, **kwargs):
-        return sqrt(arg**2)/arg*(-pi/2 + 2*atan(arg + sqrt(arg**2 - 1)))
+        return sqrt(arg**2)/arg*(-S.Pi/2 + 2*atan(arg + sqrt(arg**2 - 1)))
 
     def _eval_rewrite_as_acot(self, arg, **kwargs):
-        return sqrt(arg**2)/arg*(-pi/2 + 2*acot(arg - sqrt(arg**2 - 1)))
+        return sqrt(arg**2)/arg*(-S.Pi/2 + 2*acot(arg - sqrt(arg**2 - 1)))
 
     def _eval_rewrite_as_acsc(self, arg, **kwargs):
-        return pi/2 - acsc(arg)
+        return S.Pi/2 - acsc(arg)
 
 
 class acsc(InverseTrigonometricFunction):
@@ -3239,9 +3222,9 @@ def eval(cls, arg):
             if arg is S.NaN:
                 return S.NaN
             elif arg is S.One:
-                return pi/2
+                return S.Pi/2
             elif arg is S.NegativeOne:
-                return -pi/2
+                return -S.Pi/2
         if arg in [S.Infinity, S.NegativeInfinity, S.ComplexInfinity]:
             return S.Zero
 
@@ -3295,9 +3278,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if cdir != 0:
             cdir = arg.dir(x, cdir)
         if im(cdir) < 0 and x0.is_real and x0 > S.Zero and x0 < S.One:
-            return pi - self.func(x0)
+            return S.Pi - self.func(x0)
         elif im(cdir) > 0 and x0.is_real and x0 < S.Zero and x0 > S.NegativeOne:
-            return -pi - self.func(x0)
+            return -S.Pi - self.func(x0)
         return self.func(x0)
 
     def _eval_nseries(self, x, n, logx, cdir=0):  # acsc
@@ -3329,9 +3312,9 @@ def _eval_nseries(self, x, n, logx, cdir=0):  # acsc
         if cdir != 0:
             cdir = self.args[0].dir(x, cdir)
         if im(cdir) < 0 and arg0.is_real and arg0 > S.Zero and arg0 < S.One:
-            return pi - res
+            return S.Pi - res
         elif im(cdir) > 0 and arg0.is_real and arg0 < S.Zero and arg0 > S.NegativeOne:
-            return -pi - res
+            return -S.Pi - res
         return res
 
     def _eval_rewrite_as_log(self, arg, **kwargs):
@@ -3341,16 +3324,16 @@ def _eval_rewrite_as_asin(self, arg, **kwargs):
         return asin(1/arg)
 
     def _eval_rewrite_as_acos(self, arg, **kwargs):
-        return pi/2 - acos(1/arg)
+        return S.Pi/2 - acos(1/arg)
 
     def _eval_rewrite_as_atan(self, arg, **kwargs):
-        return sqrt(arg**2)/arg*(pi/2 - atan(sqrt(arg**2 - 1)))
+        return sqrt(arg**2)/arg*(S.Pi/2 - atan(sqrt(arg**2 - 1)))
 
     def _eval_rewrite_as_acot(self, arg, **kwargs):
-        return sqrt(arg**2)/arg*(pi/2 - acot(1/sqrt(arg**2 - 1)))
+        return sqrt(arg**2)/arg*(S.Pi/2 - acot(1/sqrt(arg**2 - 1)))
 
     def _eval_rewrite_as_asec(self, arg, **kwargs):
-        return pi/2 - asec(arg)
+        return S.Pi/2 - asec(arg)
 
 
 class atan2(InverseTrigonometricFunction):
@@ -3462,8 +3445,8 @@ def eval(cls, y, x):
         if x is S.NegativeInfinity:
             if y.is_zero:
                 # Special case y = 0 because we define Heaviside(0) = 1/2
-                return pi
-            return 2*pi*(Heaviside(re(y))) - pi
+                return S.Pi
+            return 2*S.Pi*(Heaviside(re(y))) - S.Pi
         elif x is S.Infinity:
             return S.Zero
         elif x.is_imaginary and y.is_imaginary and x.is_number and y.is_number:
@@ -3475,21 +3458,21 @@ def eval(cls, y, x):
                 return atan(y/x)
             elif x.is_negative:
                 if y.is_negative:
-                    return atan(y/x) - pi
+                    return atan(y/x) - S.Pi
                 elif y.is_nonnegative:
-                    return atan(y/x) + pi
+                    return atan(y/x) + S.Pi
             elif x.is_zero:
                 if y.is_positive:
-                    return pi/2
+                    return S.Pi/2
                 elif y.is_negative:
-                    return -pi/2
+                    return -S.Pi/2
                 elif y.is_zero:
                     return S.NaN
         if y.is_zero:
             if x.is_extended_nonzero:
-                return pi*(S.One - Heaviside(x))
+                return S.Pi*(S.One - Heaviside(x))
             if x.is_number:
-                return Piecewise((pi, re(x) < 0),
+                return Piecewise((S.Pi, re(x) < 0),
                                  (0, Ne(x, 0)),
                                  (S.NaN, True))
         if x.is_number and y.is_number:
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index d2e25ab893..82bfdad842 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -1874,6 +1874,8 @@ def _print_NDimArray(self, expr: NDimArray):
                     mat_str = 'array'
         block_str = r'\begin{%MATSTR%}%s\end{%MATSTR%}'
         block_str = block_str.replace('%MATSTR%', mat_str)
+        if mat_str == 'array':
+            block_str= block_str.replace('%s','{}%s')
         if self._settings['mat_delim']:
             left_delim: str = self._settings['mat_delim']
             right_delim = self._delim_dict[left_delim]


Test patch:
diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index 29de18b4c4..ea92baa72a 100644
--- a/sympy/printing/tests/test_latex.py
+++ b/sympy/printing/tests/test_latex.py
@@ -1,4 +1,4 @@
-from sympy import MatAdd, MatMul
+from sympy import MatAdd, MatMul, Array
 from sympy.algebras.quaternion import Quaternion
 from sympy.calculus.accumulationbounds import AccumBounds
 from sympy.combinatorics.permutations import Cycle, Permutation, AppliedPermutation
@@ -3101,3 +3101,11 @@ def test_printing_latex_array_expressions():
     M = MatrixSymbol("M", 3, 3)
     N = MatrixSymbol("N", 3, 3)
     assert latex(ArrayElement(M*N, [x, 0])) == "{{\\left(M N\\right)}_{x, 0}}"
+
+def test_Array():
+    arr = Array(range(10))
+    assert latex(arr) == r'\left[\begin{matrix}0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\end{matrix}\right]'
+
+    arr = Array(range(11))
+    # added empty arguments {}
+    assert latex(arr) == r'\left[\begin{array}{}0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10\end{array}\right]'