def shape(self):
        """Returns a list with dimensions of each index.

        Dimensions is a property of the array, not of the indices.  Still, if
        the ``IndexedBase`` does not define a shape attribute, it is assumed
        that the ranges of the indices correspond to the shape of the array.

        >>> from sympy import IndexedBase, Idx, symbols
        >>> n, m = symbols('n m', integer=True)
        >>> i = Idx('i', m)
        >>> j = Idx('j', m)
        >>> A = IndexedBase('A', shape=(n, n))
        >>> B = IndexedBase('B')
        >>> A[i, j].shape
        (n, n)
        >>> B[i, j].shape
        (m, m)
        """
        from sympy.utilities.misc import filldedent

        if self.base.shape:
            return self.base.shape
        try:
            return Tuple(*[i.upper - i.lower + 1 for i in self.indices])
        except AttributeError:
            raise IndexException(filldedent("""
                Range is not defined for all indices in: %s""" % self))
        except TypeError:
            raise IndexException(filldedent("""
                Shape cannot be inferred from Idx with
                undefined range: %s""" % self))

 def __getitem__(self, key):
     if not isinstance(key, tuple) and isinstance(key, slice):
         from sympy.matrices.expressions.slice import MatrixSlice
         return MatrixSlice(self, key, (0, None, 1))
     if isinstance(key, tuple) and len(key) == 2:
         i, j = key
         if isinstance(i, slice) or isinstance(j, slice):
             from sympy.matrices.expressions.slice import MatrixSlice
             return MatrixSlice(self, i, j)
         i, j = sympify(i), sympify(j)
         if self.valid_index(i, j) != False:
             return self._entry(i, j)
         else:
             raise IndexError("Invalid indices (%s, %s)" % (i, j))
     elif isinstance(key, (SYMPY_INTS, Integer)):
         # row-wise decomposition of matrix
         rows, cols = self.shape
         # allow single indexing if number of columns is known
         if not isinstance(cols, Integer):
             raise IndexError(filldedent('''
                 Single indexing is only supported when the number
                 of columns is known.'''))
         key = sympify(key)
         i = key // cols
         j = key % cols
         if self.valid_index(i, j) != False:
             return self._entry(i, j)
         else:
             raise IndexError("Invalid index %s" % key)
     elif isinstance(key, (Symbol, Expr)):
             raise IndexError(filldedent('''
                 Only integers may be used when addressing the matrix
                 with a single index.'''))
     raise IndexError("Invalid index, wanted %s[i,j]" % self)

    def __new__(cls, label, range=None, **kw_args):
        from sympy.utilities.misc import filldedent

        if isinstance(label, basestring):
            label = Symbol(label, integer=True)
        label, range = map(sympify, (label, range))

        if not label.is_integer:
            raise TypeError("Idx object requires an integer label.")

        elif is_sequence(range):
            if len(range) != 2:
                raise ValueError(filldedent("""
                    Idx range tuple must have length 2, but got %s""" % len(range)))
            for bound in range:
                if not (bound.is_integer or abs(bound) is S.Infinity):
                    raise TypeError("Idx object requires integer bounds.")
            args = label, Tuple(*range)
        elif isinstance(range, Expr):
            if not (range.is_integer or range is S.Infinity):
                raise TypeError("Idx object requires an integer dimension.")
            args = label, Tuple(0, range - 1)
        elif range:
            raise TypeError(filldedent("""
                The range must be an ordered iterable or
                integer SymPy expression."""))
        else:
            args = label,

        obj = Expr.__new__(cls, *args, **kw_args)
        return obj

 def _contains(self, other):
     from sympy.matrices import Matrix
     from sympy.solvers.solveset import solveset, linsolve
     from sympy.utilities.iterables import iterable, cartes
     L = self.lamda
     if self._is_multivariate():
         if not iterable(L.expr):
             if iterable(other):
                 return S.false
             return other.as_numer_denom() in self.func(
                 Lambda(L.variables, L.expr.as_numer_denom()), self.base_set)
         if len(L.expr) != len(self.lamda.variables):
             raise NotImplementedError(filldedent('''
 Dimensions of input and output of Lambda are different.'''))
         eqs = [expr - val for val, expr in zip(other, L.expr)]
         variables = L.variables
         free = set(variables)
         if all(i.is_number for i in list(Matrix(eqs).jacobian(variables))):
             solns = list(linsolve([e - val for e, val in
             zip(L.expr, other)], variables))
         else:
             syms = [e.free_symbols & free for e in eqs]
             solns = {}
             for i, (e, s, v) in enumerate(zip(eqs, syms, other)):
                 if not s:
                     if e != v:
                         return S.false
                     solns[vars[i]] = [v]
                     continue
                 elif len(s) == 1:
                     sy = s.pop()
                     sol = solveset(e, sy)
                     if sol is S.EmptySet:
                         return S.false
                     elif isinstance(sol, FiniteSet):
                         solns[sy] = list(sol)
                     else:
                         raise NotImplementedError
                 else:
                     raise NotImplementedError
             solns = cartes(*[solns[s] for s in variables])
     else:
         # assume scalar -> scalar mapping
         solnsSet = solveset(L.expr - other, L.variables[0])
         if solnsSet.is_FiniteSet:
             solns = list(solnsSet)
         else:
             raise NotImplementedError(filldedent('''
             Determining whether an ImageSet contains %s has not
             been implemented.''' % func_name(other)))
     for soln in solns:
         try:
             if soln in self.base_set:
                 return S.true
         except TypeError:
             return self.base_set.contains(soln.evalf())
     return S.false

    def _solve_reduced_system(system, gens, entry=False):
        """Recursively solves reduced polynomial systems. """
        if len(system) == len(gens) == 1:
            zeros = list(roots(system[0], gens[-1]).keys())
            return [(zero,) for zero in zeros]

        basis = groebner(system, gens, polys=True)

        if len(basis) == 1 and basis[0].is_ground:
            if not entry:
                return []
            else:
                return None

        univariate = list(filter(_is_univariate, basis))

        if len(univariate) == 1:
            f = univariate.pop()
        else:
            raise NotImplementedError(filldedent('''
                only zero-dimensional systems supported
                (finite number of solutions)
                '''))

        gens = f.gens
        gen = gens[-1]

        zeros = list(roots(f.ltrim(gen)).keys())

        if not zeros:
            return []

        if len(basis) == 1:
            return [(zero,) for zero in zeros]

        solutions = []

        for zero in zeros:
            new_system = []
            new_gens = gens[:-1]

            for b in basis[:-1]:
                eq = _subs_root(b, gen, zero)

                if eq is not S.Zero:
                    new_system.append(eq)

            for solution in _solve_reduced_system(new_system, new_gens):
                solutions.append(solution + (zero,))

        if solutions and len(solutions[0]) != len(gens):
            raise NotImplementedError(filldedent('''
                only zero-dimensional systems supported
                (finite number of solutions)
                '''))
        return solutions

def conjugate_gauss_beams(wavelen, waist_in, waist_out, **kwargs):
    """
    Find the optical setup conjugating the object/image waists.

    Parameters
    ==========

    wavelen : the wavelength of the beam
    waist_in and waist_out : the waists to be conjugated
    f : the focal distance of the element used in the conjugation

    Returns
    =======

    a tuple containing (s_in, s_out, f)
    s_in : the distance before the optical element
    s_out : the distance after the optical element
    f : the focal distance of the optical element

    Examples
    ========

    >>> from sympy.physics.optics import conjugate_gauss_beams
    >>> from sympy import symbols, factor
    >>> l, w_i, w_o, f = symbols('l w_i w_o f')

    >>> conjugate_gauss_beams(l, w_i, w_o, f=f)[0]
    f*(1 - sqrt(w_i**2/w_o**2 - pi**2*w_i**4/(f**2*l**2)))

    >>> factor(conjugate_gauss_beams(l, w_i, w_o, f=f)[1])
    f*w_o**2*(w_i**2/w_o**2 - sqrt(w_i**2/w_o**2 -
              pi**2*w_i**4/(f**2*l**2)))/w_i**2

    >>> conjugate_gauss_beams(l, w_i, w_o, f=f)[2]
    f
    """
    #TODO add the other possible arguments
    wavelen, waist_in, waist_out = map(sympify, (wavelen, waist_in, waist_out))
    m = waist_out / waist_in
    z = waist2rayleigh(waist_in, wavelen)
    if len(kwargs) != 1:
        raise ValueError("The function expects only one named argument")
    elif 'dist' in kwargs:
        raise NotImplementedError(filldedent('''
            Currently only focal length is supported as a parameter'''))
    elif 'f' in kwargs:
        f = sympify(kwargs['f'])
        s_in = f * (1 - sqrt(1/m**2 - z**2/f**2))
        s_out = gaussian_conj(s_in, z, f)[0]
    elif 's_in' in kwargs:
        raise NotImplementedError(filldedent('''
            Currently only focal length is supported as a parameter'''))
    else:
        raise ValueError(filldedent('''
            The functions expects the focal length as a named argument'''))
    return (s_in, s_out, f)

 def __new__(cls, *args, **kwargs):
     from sympy.matrices.immutable import ImmutableDenseMatrix
     from sympy.matrices import zeros
     from sympy.matrices.matrices import MatrixBase
     from sympy.utilities.iterables import is_sequence
     isMat = lambda i: getattr(i, 'is_Matrix', False)
     if len(args) != 1 or \
             not is_sequence(args[0]) or \
             len(set([isMat(r) for r in args[0]])) != 1:
         raise ValueError(filldedent('''
             expecting a sequence of 1 or more rows
             containing Matrices.'''))
     rows = args[0] if args else []
     if not isMat(rows):
         if rows and isMat(rows[0]):
             rows = [rows]  # rows is not list of lists or []
         # regularity check
         # same number of matrices in each row
         blocky = ok = len(set([len(r) for r in rows])) == 1
         if ok:
             # same number of rows for each matrix in a row
             for r in rows:
                 ok = len(set([i.rows for i in r])) == 1
                 if not ok:
                     break
             blocky = ok
             # same number of cols for each matrix in each col
             for c in range(len(rows[0])):
                 ok = len(set([rows[i][c].cols
                     for i in range(len(rows))])) == 1
                 if not ok:
                     break
         if not ok:
             # same total cols in each row
             ok = len(set([
                 sum([i.cols for i in r]) for r in rows])) == 1
             if blocky and ok:
                 raise ValueError(filldedent('''
                     Although this matrix is comprised of blocks,
                     the blocks do not fill the matrix in a
                     size-symmetric fashion. To create a full matrix
                     from these arguments, pass them directly to
                     Matrix.'''))
             raise ValueError(filldedent('''
                 When there are not the same number of rows in each
                 row's matrices or there are not the same number of
                 total columns in each row, the matrix is not a
                 block matrix. If this matrix is known to consist of
                 blocks fully filling a 2-D space then see
                 Matrix.irregular.'''))
     mat = ImmutableDenseMatrix(rows, evaluate=False)
     obj = Basic.__new__(cls, mat)
     return obj

def reduce_inequalities(inequalities, symbols=[]):
    """Reduce a system of inequalities with rational coefficients.

    Examples
    ========

    >>> from sympy import sympify as S, Symbol
    >>> from sympy.abc import x, y
    >>> from sympy.solvers.inequalities import reduce_inequalities

    >>> reduce_inequalities(0 <= x + 3, [])
    And(-3 <= x, x < oo)

    >>> reduce_inequalities(0 <= x + y*2 - 1, [x])
    x >= -2*y + 1
    """
    if not iterable(inequalities):
        inequalities = [inequalities]
    inequalities = [sympify(i) for i in inequalities]

    gens = set().union(*[i.free_symbols for i in inequalities])

    if not iterable(symbols):
        symbols = [symbols]
    symbols = (set(symbols) or gens) & gens
    if any(i.is_real is False for i in symbols):
        raise TypeError(
            filldedent(
                """
            inequalities cannot contain symbols that are not real."""
            )
        )

    # make vanilla symbol real
    recast = dict([(i, Dummy(i.name, real=True)) for i in gens if i.is_real is None])
    inequalities = [i.xreplace(recast) for i in inequalities]
    symbols = set([i.xreplace(recast) for i in symbols])

    # prefilter
    keep = []
    for i in inequalities:
        if isinstance(i, Relational):
            i = i.func(i.lhs.as_expr() - i.rhs.as_expr(), 0)
        elif i not in (True, False):
            i = Eq(i, 0)
        if i == True:
            continue
        elif i == False:
            return S.false
        if i.lhs.is_number:
            raise NotImplementedError("could not determine truth value of %s" % i)
        keep.append(i)
    inequalities = keep
    del keep

    # solve system
    rv = _reduce_inequalities(inequalities, symbols)

    # restore original symbols and return
    return rv.xreplace(dict([(v, k) for k, v in recast.items()]))

    def __new__(
        cls, center=None, hradius=None, vradius=None, eccentricity=None, **kwargs):
        hradius = sympify(hradius)
        vradius = sympify(vradius)

        eccentricity = sympify(eccentricity)

        if center is None:
            center = Point(0, 0)
        else:
            center = Point(center, dim=2)

        if len(center) != 2:
            raise ValueError('The center of "{0}" must be a two dimensional point'.format(cls))

        if len(list(filter(lambda x: x is not None, (hradius, vradius, eccentricity)))) != 2:
            raise ValueError(filldedent('''
                Exactly two arguments of "hradius", "vradius", and
                "eccentricity" must not be None.'''))

        if eccentricity is not None:
            if hradius is None:
                hradius = vradius / sqrt(1 - eccentricity**2)
            elif vradius is None:
                vradius = hradius * sqrt(1 - eccentricity**2)

        if hradius == vradius:
            return Circle(center, hradius, **kwargs)

        if hradius == 0 or vradius == 0:
            return Segment(Point(center[0] - hradius, center[1] - vradius), Point(center[0] + hradius, center[1] + vradius))

        return GeometryEntity.__new__(cls, center, hradius, vradius, **kwargs)

 def _solve_relational(r):
     if sym not in r.free_symbols:
         nonsymfail(r)
     rv = _solve_inequality(r, sym)
     if isinstance(rv, Relational):
         free = rv.args[1].free_symbols
         if rv.args[0] != sym or sym in free:
             raise NotImplementedError(filldedent('''
                 Unable to solve relational
                 %s for %s.''' % (r, sym)))
         if rv.rel_op == '==':
             # this equality has been affirmed to have the form
             # Eq(sym, rhs) where rhs is sym-free; it represents
             # a zero-width interval which will be ignored
             # whether it is an isolated condition or contained
             # within an And or an Or
             rv = S.false
         elif rv.rel_op == '!=':
             try:
                 rv = Or(sym < rv.rhs, sym > rv.rhs)
             except TypeError:
                 # e.g. x != I ==> all real x satisfy
                 rv = S.true
     elif rv == (S.NegativeInfinity < sym) & (sym < S.Infinity):
         rv = S.true
     return rv

 def __new__(cls, lhs, rhs, rop=None, **assumptions):
     # If called by a subclass, do nothing special and pass on to Expr.
     if cls is not Relational:
         return Expr.__new__(cls, lhs, rhs, **assumptions)
     # If called directly with an operator, look up the subclass
     # corresponding to that operator and delegate to it
     try:
         cls = cls.ValidRelationOperator[rop]
         rv = cls(lhs, rhs, **assumptions)
         # /// drop when Py2 is no longer supported
         # validate that Booleans are not being used in a relational
         # other than Eq/Ne;
         if isinstance(rv, (Eq, Ne)):
             pass
         elif isinstance(rv, Relational):  # could it be otherwise?
             from sympy.core.symbol import Symbol
             from sympy.logic.boolalg import Boolean
             for a in rv.args:
                 if isinstance(a, Symbol):
                     continue
                 if isinstance(a, Boolean):
                     from sympy.utilities.misc import filldedent
                     raise TypeError(filldedent('''
                         A Boolean argument can only be used in
                         Eq and Ne; all other relationals expect
                         real expressions.
                     '''))
         # \\\
         return rv
     except KeyError:
         raise ValueError(
             "Invalid relational operator symbol: %r" % rop)

def mechanics_printing(**kwargs):

    # mechanics_printing has slightly different functionality in 0.7.5 but
    # shouldn't fundamentally need a deprecation warning so we do this
    # little wrapper that gives the warning that things have changed.

    # TODO : Remove this warning in the release after SymPy 0.7.5

    # The message is only printed if this function is called with no args,
    # as was the previous only way to call it.

    def dict_is_empty(D):
        for k in D:
            return False
        return True

    if dict_is_empty(kwargs):
        msg = (
            "See the doc string for slight changes to this function: "
            "keyword args may be needed for the desired effect. "
            "Otherwise use sympy.physics.vector.init_vprinting directly."
        )
        SymPyDeprecationWarning(filldedent(msg)).warn()

    init_vprinting(**kwargs)

 def __new__(cls, dist):
     if not isinstance(dist, (ContinuousDistribution, DiscreteDistribution)):
         raise ValueError(filldedent('''CompoundDistribution can only be
          initialized from ContinuousDistribution or DiscreteDistribution
          '''))
     _args = dist.args
     if not any([isinstance(i, RandomSymbol) for i in _args]):
         return dist
     return Basic.__new__(cls, dist)

 def _eval_rewrite_as_ITE(self, *args):
     byfree = {}
     args = list(args)
     default = any(c == True for b, c in args)
     for i, (b, c) in enumerate(args):
         if not isinstance(b, Boolean) and b != True:
             raise TypeError(filldedent('''
                 Expecting Boolean or bool but got `%s`
                 ''' % func_name(b)))
         if c == True:
             break
         # loop over independent conditions for this b
         for c in c.args if isinstance(c, Or) else [c]:
             free = c.free_symbols
             x = free.pop()
             try:
                 byfree[x] = byfree.setdefault(
                     x, S.EmptySet).union(c.as_set())
             except NotImplementedError:
                 if not default:
                     raise NotImplementedError(filldedent('''
                         A method to determine whether a multivariate
                         conditional is consistent with a complete coverage
                         of all variables has not been implemented so the
                         rewrite is being stopped after encountering `%s`.
                         This error would not occur if a default expression
                         like `(foo, True)` were given.
                         ''' % c))
             if byfree[x] in (S.UniversalSet, S.Reals):
                 # collapse the ith condition to True and break
                 args[i] = list(args[i])
                 c = args[i][1] = True
                 break
         if c == True:
             break
     if c != True:
         raise ValueError(filldedent('''
             Conditions must cover all reals or a final default
             condition `(foo, True)` must be given.
             '''))
     last, _ = args[i]  # ignore all past ith arg
     for a, c in reversed(args[:i]):
         last = ITE(c, a, last)
     return _canonical(last)

 def _calc_limit(a, b):
     """
     replace d with a, using subs if possible, otherwise limit
     where sign of b is considered
     """
     avals = list(set([_calc_limit_1(Fi, a, b) for Fi in F]))
     if len(avals) > 1:
         raise ValueError(filldedent('''
         The mapping between F(x) and f(u) did not
         give a unique limit.'''))
     return avals[0]

 def __new__(cls, *args):
     if len(args) == 1 and isinstance(args[0], Matrix) \
             and args[0].shape == (2, 1):
         temp = args[0]
     elif len(args) == 2:
         temp = ((args[0],), (args[1],))
     else:
         raise ValueError(filldedent('''
             Expecting 2x1 Matrix or the 2 elements of
             the Matrix but got %s''' % str(args)))
     return Matrix.__new__(cls, temp)

    def reflect(self, line):
        """Override GeometryEntity.reflect since the radius
        is not a GeometryEntity.

        Examples
        ========

        >>> from sympy import Circle, Line
        >>> Circle((0, 1), 1).reflect(Line((0, 0), (1, 1)))
        Circle(Point(1, 0), -1)
        >>> from sympy import Ellipse, Line, Point
        >>> Ellipse(Point(3, 4), 1, 3).reflect(Line(Point(0, -4), Point(5, 0)))
        Traceback (most recent call last):
        ...
        NotImplementedError:
        General Ellipse is not supported but the equation of the reflected
        Ellipse is given by the zeros of: f(x, y) = (9*x/41 + 40*y/41 +
        37/41)**2 + (40*x/123 - 3*y/41 - 364/123)**2 - 1

        Notes
        =====

        Until the general ellipse (with no axis parallel to the x-axis) is
        supported a NotImplemented error is raised and the equation whose
        zeros define the rotated ellipse is given.

        """
        def _uniquely_named_symbol(xname, *exprs):
            """Return a symbol which, when printed, will have a name unique
            from any other already in the expressions given. The name is made
            unique by prepending underscores.
            """
            prefix = '%s'
            x = prefix % xname
            syms = set.union(*[e.free_symbols for e in exprs])
            while any(x == str(s) for s in syms):
                prefix = '_' + prefix
                x = prefix % xname
            return _symbol(x)

        if line.slope in (0, oo):
            c = self.center
            c = c.reflect(line)
            return self.func(c, -self.hradius, self.vradius)
        else:
            x, y = [_uniquely_named_symbol(name, self, line) for name in 'xy']
            expr = self.equation(x, y)
            p = Point(x, y).reflect(line)
            result = expr.subs(zip((x, y), p.args
                               ), simultaneous=True)
            raise NotImplementedError(filldedent(
                'General Ellipse is not supported but the equation '
                'of the reflected Ellipse is given by the zeros of: ' +
                "f(%s, %s) = %s" % (str(x), str(y), str(result))))

    def __new__(cls, base, *args, **kw_args):
        from sympy.utilities.misc import filldedent

        if not args:
            raise IndexException("Indexed needs at least one index.")
        if isinstance(base, (basestring, Symbol)):
            base = IndexedBase(base)
        elif not isinstance(base, IndexedBase):
            raise TypeError(filldedent("""
                Indexed expects string, Symbol or IndexedBase as base."""))
        return Expr.__new__(cls, base, *args, **kw_args)

 def __init__(self, *args):
     if len(args) == 4:
         temp = ((args[0], args[1]), (args[2], args[3]))
     elif len(args) == 1 \
          and isinstance(args[0], Matrix) \
          and args[0].shape == (2, 2):
         temp = args[0]
     else:
         raise ValueError(filldedent('''
             Expecting 2x2 Matrix or the 4 elements of
             the Matrix but got %s''' % str(args)))
     Matrix.__init__(self, temp)

def exprcondpair_new(cls, expr, cond):
    expr = as_Basic(expr)
    if cond == True:
        return Tuple.__new__(cls, expr, true)
    elif cond == False:
        return Tuple.__new__(cls, expr, false)

    if not isinstance(cond, Boolean):
        raise TypeError(filldedent('''
            Second argument must be a Boolean,
            not `%s`''' % func_name(cond)))
    return Tuple.__new__(cls, expr, cond)