def test_zoo():
    b = Symbol('b', finite=True)
    nz = Symbol('nz', nonzero=True)
    p = Symbol('p', positive=True)
    n = Symbol('n', negative=True)
    im = Symbol('i', imaginary=True)
    c = Symbol('c', complex=True)
    pb = Symbol('pb', positive=True, finite=True)
    nb = Symbol('nb', negative=True, finite=True)
    imb = Symbol('ib', imaginary=True, finite=True)
    for i in [I, S.Infinity, S.NegativeInfinity, S.Zero, S.One, S.Pi, S.Half, S(3), log(3),
              b, nz, p, n, im, pb, nb, imb, c]:
        if i.is_finite and (i.is_real or i.is_imaginary):
            assert i + zoo is zoo
            assert i - zoo is zoo
            assert zoo + i is zoo
            assert zoo - i is zoo
        elif i.is_finite is not False:
            assert (i + zoo).is_Add
            assert (i - zoo).is_Add
            assert (zoo + i).is_Add
            assert (zoo - i).is_Add
        else:
            assert (i + zoo) is S.NaN
            assert (i - zoo) is S.NaN
            assert (zoo + i) is S.NaN
            assert (zoo - i) is S.NaN

        if fuzzy_not(i.is_zero) and (i.is_real or i.is_imaginary):
            assert i*zoo is zoo
            assert zoo*i is zoo
        elif i.is_zero:
            assert i*zoo is S.NaN
            assert zoo*i is S.NaN
        else:
            assert (i*zoo).is_Mul
            assert (zoo*i).is_Mul

        if fuzzy_not((1/i).is_zero) and (i.is_real or i.is_imaginary):
            assert zoo/i is zoo
        elif (1/i).is_zero:
            assert zoo/i is S.NaN
        elif i.is_zero:
            assert zoo/i is zoo
        else:
            assert (zoo/i).is_Mul

    assert (I*oo).is_Mul  # allow directed infinity
    assert zoo + zoo is S.NaN
    assert zoo * zoo is zoo
    assert zoo - zoo is S.NaN
    assert zoo/zoo is S.NaN
    assert zoo**zoo is S.NaN
    assert zoo**0 is S.One
    assert zoo**2 is zoo
    assert 1/zoo is S.Zero

    assert Mul.flatten([S(-1), oo, S(0)]) == ([S.NaN], [], None)

    def __new__(cls, *args):
        n = args[BinomialDistribution._argnames.index('n')]
        p = args[BinomialDistribution._argnames.index('p')]
        n_sym = sympify(n)
        p_sym = sympify(p)

        if fuzzy_not(fuzzy_and((n_sym.is_integer, n_sym.is_nonnegative))):
            raise ValueError("'n' must be positive integer. n = %s." % str(n))
        elif fuzzy_not(fuzzy_and((p_sym.is_nonnegative, (p_sym - 1).is_nonpositive))):
            raise ValueError("'p' must be: 0 <= p <= 1 . p = %s" % str(p))
        else:
            return super(BinomialDistribution, cls).__new__(cls, *args)

 def Basic(expr, assumptions):
     if expr.is_number:
         notpositive = fuzzy_not(AskPositiveHandler._number(expr, assumptions))
         if notpositive:
             return ask(Q.real(expr), assumptions)
         else:
             return notpositive

 def _eval_is_algebraic(self):
     s = self.func(*self.args)
     if s.func == self.func:
         if fuzzy_not(self.args[0].is_zero) and self.args[0].is_algebraic:
             return False
     else:
         return s.is_algebraic

 def _eval_power(self, other):
     if (
         fuzzy_not(self.args[0].is_zero) and
         other.is_integer and
         other.is_even
     ):
         return S.One

    def eval(cls, x, k=None):
        if k is S.Zero:
            return cls(x)
        elif k is None:
            k = S.Zero

        if k is S.Zero:
            if x is S.Zero:
                return S.Zero
            if x is S.Exp1:
                return S.One
            if x == -1/S.Exp1:
                return S.NegativeOne
            if x == -log(2)/2:
                return -log(2)
            if x is S.Infinity:
                return S.Infinity

        if fuzzy_not(k.is_zero):
            if x is S.Zero:
                return S.NegativeInfinity
        if k is S.NegativeOne:
            if x == -S.Pi/2:
                return -S.ImaginaryUnit*S.Pi/2
            elif x == -1/S.Exp1:
                return S.NegativeOne
            elif x == -2*exp(-2):
                return -Integer(2)

def refine_abs(expr, assumptions):
    """
    Handler for the absolute value.

    Examples
    ========

    >>> from sympy import Symbol, Q, refine, Abs
    >>> from sympy.assumptions.refine import refine_abs
    >>> from sympy.abc import x
    >>> refine_abs(Abs(x), Q.real(x))
    >>> refine_abs(Abs(x), Q.positive(x))
    x
    >>> refine_abs(Abs(x), Q.negative(x))
    -x

    """
    from sympy.core.logic import fuzzy_not

    arg = expr.args[0]
    if ask(Q.real(arg), assumptions) and fuzzy_not(ask(Q.negative(arg), assumptions)):
        # if it's nonnegative
        return arg
    if ask(Q.negative(arg), assumptions):
        return -arg

    def eval(cls, arg, k=0):
        """
        Returns a simplified form or a value of DiracDelta depending on the
        argument passed by the DiracDelta object.

        The ``eval()`` method is automatically called when the ``DiracDelta`` class
        is about to be instantiated and it returns either some simplified instance
        or the unevaluated instance depending on the argument passed. In other words,
        ``eval()`` method is not needed to be called explicitly, it is being called
        and evaluated once the object is called.

        Examples
        ========

        >>> from sympy import DiracDelta, S, Subs
        >>> from sympy.abc import x

        >>> DiracDelta(x)
        DiracDelta(x)

        >>> DiracDelta(x,1)
        DiracDelta(x, 1)

        >>> DiracDelta(1)
        0

        >>> DiracDelta(5,1)
        0

        >>> DiracDelta(0)
        DiracDelta(0)

        >>> DiracDelta(-1)
        0

        >>> DiracDelta(S.NaN)
        nan

        >>> DiracDelta(x).eval(1)
        0

        >>> DiracDelta(x - 100).subs(x, 5)
        0

        >>> DiracDelta(x - 100).subs(x, 100)
        DiracDelta(0)

        """
        k = sympify(k)
        if not k.is_Integer or k.is_negative:
            raise ValueError("Error: the second argument of DiracDelta must be \
            a non-negative integer, %s given instead." % (k,))
        arg = sympify(arg)
        if arg is S.NaN:
            return S.NaN
        if arg.is_nonzero:
            return S.Zero
        if fuzzy_not(im(arg).is_zero):
            raise ValueError("Function defined only for Real Values. Complex part: %s  found in %s ." % (repr(im(arg)), repr(arg)) )

    def _eval_is_even(self):
        is_integer = self.is_integer

        if is_integer:
            return fuzzy_not(self._eval_is_odd())

        elif is_integer is False:
            return False

def substitution_rule(integral):
    integrand, symbol = integral

    u_var = sympy.Dummy("u")
    substitutions = find_substitutions(integrand, symbol, u_var)
    if substitutions:
        ways = []
        for u_func, c, substituted in substitutions:
            subrule = integral_steps(substituted, u_var)
            if contains_dont_know(subrule):
                continue

            if sympy.simplify(c - 1) != 0:
                _, denom = c.as_numer_denom()
                if subrule:
                    subrule = ConstantTimesRule(c, substituted, subrule, substituted, u_var)

                if denom.free_symbols:
                    piecewise = []
                    could_be_zero = []

                    if isinstance(denom, sympy.Mul):
                        could_be_zero = denom.args
                    else:
                        could_be_zero.append(denom)

                    for expr in could_be_zero:
                        if not fuzzy_not(expr.is_zero):
                            substep = integral_steps(integrand.subs(expr, 0), symbol)

                            if substep:
                                piecewise.append((
                                    substep,
                                    sympy.Eq(expr, 0)
                                ))
                    piecewise.append((subrule, True))
                    subrule = PiecewiseRule(piecewise, substituted, symbol)

            ways.append(URule(u_var, u_func, c,
                              subrule,
                              integrand, symbol))

        if len(ways) > 1:
            return AlternativeRule(ways, integrand, symbol)
        elif ways:
            return ways[0]

    elif integrand.has(sympy.exp):
        u_func = sympy.exp(symbol)
        c = 1
        substituted = integrand / u_func.diff(symbol)
        substituted = substituted.subs(u_func, u_var)

        if symbol not in substituted.free_symbols:
            return URule(u_var, u_func, c,
                         integral_steps(substituted, u_var),
                         integrand, symbol)

 def _eval_is_rational(self):
     s = self.func(*self.args)
     if s.func == self.func:
         if s.exp is S.Zero:
             return True
         elif s.exp.is_rational and fuzzy_not(s.exp.is_zero):
             return False
     else:
         return s.is_rational

 def _eval_is_real(self):
     x = self.args[0]
     if len(self.args) == 1:
         k = S.Zero
     else:
         k = self.args[1]
     if k.is_zero:
         if (x + 1/S.Exp1).is_positive:
             return True
         elif (x + 1/S.Exp1).is_nonpositive:
             return False
     elif (k + 1).is_zero:
         if x.is_negative and (x + 1/S.Exp1).is_positive:
             return True
         elif x.is_nonpositive or (x + 1/S.Exp1).is_nonnegative:
             return False
     elif fuzzy_not(k.is_zero) and fuzzy_not((k + 1).is_zero):
         if x.is_real:
             return False

 def eval(cls, arg):
     arg = sympify(arg)
     if arg is S.NaN:
         return S.NaN
     elif fuzzy_not(im(arg).is_zero):
         raise ValueError("Function defined only for Real Values. Complex part: %s  found in %s ." % (repr(im(arg)), repr(arg)) )
     elif arg.is_negative:
         return S.Zero
     elif arg.is_positive:
         return S.One

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return ReciprocalRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if fuzzy_not(sympy.log(base).is_zero):
            return rule
        elif sympy.log(base).is_zero:
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule(
            [(ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)), (rule, True)], integrand, symbol
        )

    def eval(cls, i, j):
        """
        Evaluates the discrete delta function.

        Examples
        ========

        >>> from sympy.functions.special.tensor_functions import KroneckerDelta
        >>> from sympy.abc import i, j, k

        >>> KroneckerDelta(i, j)
        KroneckerDelta(i, j)
        >>> KroneckerDelta(i, i)
        1
        >>> KroneckerDelta(i, i + 1)
        0
        >>> KroneckerDelta(i, i + 1 + k)
        KroneckerDelta(i, i + k + 1)

        # indirect doctest

        """
        diff = i - j
        if diff.is_zero:
            return S.One
        elif fuzzy_not(diff.is_zero):
            return S.Zero

        if i.assumptions0.get("below_fermi") and \
                j.assumptions0.get("above_fermi"):
            return S.Zero
        if j.assumptions0.get("below_fermi") and \
                i.assumptions0.get("above_fermi"):
            return S.Zero
        # to make KroneckerDelta canonical
        # following lines will check if inputs are in order
        # if not, will return KroneckerDelta with correct order
        if i is not min(i, j, key=default_sort_key):
            return cls(j, i)

def refine_abs(expr, assumptions):
    """
    Handler for the absolute value.

    Examples
    ========

    >>> from sympy import Symbol, Q, refine, Abs
    >>> from sympy.assumptions.refine import refine_abs
    >>> from sympy.abc import x
    >>> refine_abs(Abs(x), Q.real(x))
    >>> refine_abs(Abs(x), Q.positive(x))
    x
    >>> refine_abs(Abs(x), Q.negative(x))
    -x

    """
    from sympy.core.logic import fuzzy_not
    from sympy import Abs
    arg = expr.args[0]
    if ask(Q.real(arg), assumptions) and \
            fuzzy_not(ask(Q.negative(arg), assumptions)):
        # if it's nonnegative
        return arg
    if ask(Q.negative(arg), assumptions):
        return -arg
    # arg is Mul
    if isinstance(arg, Mul):
        r = [refine(abs(a), assumptions) for a in arg.args]
        non_abs = []
        in_abs = []
        for i in r:
            if isinstance(i, Abs):
                in_abs.append(i.args[0])
            else:
                non_abs.append(i)
        return Mul(*non_abs) * Abs(Mul(*in_abs))

    def eval(cls, variable, offset, exponent):
        """
        Returns a simplified form or a value of Singularity Function depending on the
        argument passed by the object.

        The ``eval()`` method is automatically called when the ``SingularityFunction`` class
        is about to be instantiated and it returns either some simplified instance
        or the unevaluated instance depending on the argument passed. In other words,
        ``eval()`` method is not needed to be called explicitly, it is being called
        and evaluated once the object is called.

        Examples
        ========
        >>> from sympy import SingularityFunction, Symbol, nan
        >>> from sympy.abc import x, a, n
        >>> SingularityFunction(x, a, n)
        SingularityFunction(x, a, n)
        >>> SingularityFunction(5, 3, 2)
        4
        >>> SingularityFunction(x, a, nan)
        nan
        >>> SingularityFunction(x, 3, 0).subs(x, 3)
        1
        >>> SingularityFunction(x, a, n).eval(3, 5, 1)
        0
        >>> SingularityFunction(x, a, n).eval(4, 1, 5)
        243
        >>> x = Symbol('x', positive = True)
        >>> a = Symbol('a', negative = True)
        >>> n = Symbol('n', nonnegative = True)
        >>> SingularityFunction(x, a, n)
        (-a + x)**n
        >>> x = Symbol('x', negative = True)
        >>> a = Symbol('a', positive = True)
        >>> SingularityFunction(x, a, n)
        0

        """

        x = sympify(variable)
        a = sympify(offset)
        n = sympify(exponent)
        shift = (x - a)

        if fuzzy_not(im(shift).is_zero):
            raise ValueError("Singularity Functions are defined only for Real Numbers.")
        if fuzzy_not(im(n).is_zero):
            raise ValueError("Singularity Functions are not defined for imaginary exponents.")
        if shift is S.NaN or n is S.NaN:
            return S.NaN
        if (n + 2).is_negative:
            raise ValueError("Singularity Functions are not defined for exponents less than -2.")
        if shift.is_negative:
            return S.Zero
        if n.is_nonnegative and shift.is_nonnegative:
            return (x - a)**n
        if n == -1 or n == -2:
            if shift.is_negative or shift.is_positive:
                return S.Zero
            if shift.is_zero:
                return S.Infinity

 def _eval_is_nonzero(self):
     return fuzzy_not(self._args[0].is_zero)

 def _eval_Abs(self):
     if fuzzy_not(self.args[0].is_zero):
         return S.One

 def _eval_is_integer(self):
     from sympy.core.logic import fuzzy_and, fuzzy_not
     p, q = self.args
     if fuzzy_and([p.is_integer, q.is_integer, fuzzy_not(q.is_zero)]):
         return True