def test_trig():
    arg = FiniteSet(0, *(pi/x for x in (12, 6, 4, 3, 2, 1)))
    arg += FiniteSet(*sum(((pi-x, pi+x, 2*pi-x) for x in arg), ()))
    arg += FiniteSet(*(-x for x in arg))

    def emit(name, func):
        for x in sorted(arg):
            if func(x) != sp.zoo:
                test(name, x, sp.simplify(sp.sqrtdenest(func(x))))
        test(name, S.NaN, S.NaN)

    def aemit(aname, func, domain):
        for x, y in sorted((sp.simplify(sp.sqrtdenest(func(x))), x)
                           for x in arg & domain):
            test(aname, x, y)
        test(aname, S.NaN, S.NaN)

    emit("sin", sp.sin)
    emit("cos", sp.cos)
    emit("tan", sp.tan)
    aemit("asin", sp.sin, Interval(-pi/2, pi/2))
    aemit("acos", sp.cos, Interval(0, pi))
    aemit("atan", sp.tan, Interval(-pi/2, pi/2, True, True))

    for x in sorted(arg & Interval(-pi, pi, True, True)):
        a = sp.simplify(sp.sqrtdenest(sp.sin(x)))
        b = sp.simplify(sp.sqrtdenest(sp.cos(x)))
        f = sp.sqrt(3) if \
            sp.sqrt(3) == abs(sp.simplify(a/b).as_numer_denom()[0]) else 1
        if b:
            test("atan2", sp.simplify(a/abs(b))*f, sp.sign(b)*f, x)
        if a and a != b:
            test("atan2", sp.sign(a)*f, sp.simplify(b/abs(a))*f, x)
    test("atan2", S.NaN, 1, S.NaN)
    test("atan2", 1, S.NaN, S.NaN)

def test_rewrite():
    x, y = Symbol('x', real=True), Symbol('y')
    assert Heaviside(x).rewrite(Piecewise) == (
        Piecewise((0, x < 0), (Heaviside(0), Eq(x, 0)), (1, x > 0)))
    assert Heaviside(y).rewrite(Piecewise) == (
        Piecewise((0, y < 0), (Heaviside(0), Eq(y, 0)), (1, y > 0)))
    assert Heaviside(x, y).rewrite(Piecewise) == (
        Piecewise((0, x < 0), (y, Eq(x, 0)), (1, x > 0)))
    assert Heaviside(x, 0).rewrite(Piecewise) == (
        Piecewise((0, x <= 0), (1, x > 0)))
    assert Heaviside(x, 1).rewrite(Piecewise) == (
        Piecewise((0, x < 0), (1, x >= 0)))

    assert Heaviside(x).rewrite(sign) == (sign(x)+1)/2
    assert Heaviside(y).rewrite(sign) == Heaviside(y)
    assert Heaviside(x, S.Half).rewrite(sign) == (sign(x)+1)/2
    assert Heaviside(x, y).rewrite(sign) == Heaviside(x, y)

    assert DiracDelta(y).rewrite(Piecewise) == Piecewise((DiracDelta(0), Eq(y, 0)), (0, True))
    assert DiracDelta(y, 1).rewrite(Piecewise) == DiracDelta(y, 1)
    assert DiracDelta(x - 5).rewrite(Piecewise) == (
        Piecewise((DiracDelta(0), Eq(x - 5, 0)), (0, True)))

    assert (x*DiracDelta(x - 10)).rewrite(SingularityFunction) == x*SingularityFunction(x, 10, -1)
    assert 5*x*y*DiracDelta(y, 1).rewrite(SingularityFunction) == 5*x*y*SingularityFunction(y, 0, -2)
    assert DiracDelta(0).rewrite(SingularityFunction) == SingularityFunction(0, 0, -1)
    assert DiracDelta(0, 1).rewrite(SingularityFunction) == SingularityFunction(0, 0, -2)

    assert Heaviside(x).rewrite(SingularityFunction) == SingularityFunction(x, 0, 0)
    assert 5*x*y*Heaviside(y + 1).rewrite(SingularityFunction) == 5*x*y*SingularityFunction(y, -1, 0)
    assert ((x - 3)**3*Heaviside(x - 3)).rewrite(SingularityFunction) == (x - 3)**3*SingularityFunction(x, 3, 0)
    assert Heaviside(0).rewrite(SingularityFunction) == SingularityFunction(0, 0, 0)

def sp_derive():

    import sympy as sp

    vars = 'G_s, s_n, s_p_n, w_n, dw_n, ds_n, G_s, G_w, c, phi'

    syms = sp.symbols(vars)

    for var, sym in zip(vars.split(','), syms):
        globals()[var.strip()] = sym

    s_n1 = s_n + ds_n
    w_n1 = w_n + dw_n

    tau_trial = G_s * (s_n1 - s_p_n)

    print 'diff', sp.diff(tau_trial, ds_n)

    print tau_trial

    sig_n1 = G_w * w_n1

    print sig_n1

    tau_fr = (c + sig_n1 * sp.tan(phi)) * sp.Heaviside(sig_n1 - c / sp.tan(phi))

    print tau_fr

    d_tau_fr = sp.diff(tau_fr, dw_n)

    print d_tau_fr

    f_trial = sp.abs(tau_trial) - tau_fr

    print f_trial

    d_gamma = f_trial / G_s

    print 'd_gamma'
    sp.pretty_print(d_gamma)

    print 'd_gamma_s'
    sp.pretty_print(sp.diff(d_gamma, ds_n))

    print 'tau_n1'
    tau_n1 = sp.simplify(tau_trial - d_gamma * G_s * sp.sign(tau_trial))
    sp.pretty_print(tau_n1)

    print 'dtau_n1_w'
    dtau_n1_w = sp.diff(tau_n1, dw_n)
    sp.pretty_print(dtau_n1_w)

    print 'dtau_n1_s'
    dtau_n1_s = sp.diff(d_gamma * sp.sign(tau_trial), ds_n)
    print dtau_n1_s

    s_p_n1 = s_p_n + d_gamma * sp.sign(tau_trial)

    print s_p_n1

def test_fcode_sign():  #issue 12267
    x=symbols('x')
    y=symbols('y', integer=True)
    z=symbols('z', complex=True)
    assert fcode(sign(x), standard=95, source_format='free') == "merge(0d0, dsign(1d0, x), x == 0d0)"
    assert fcode(sign(y), standard=95, source_format='free') == "merge(0, isign(1, y), y == 0)"
    assert fcode(sign(z), standard=95, source_format='free') == "merge(cmplx(0d0, 0d0), z/abs(z), abs(z) == 0d0)"
    raises(NotImplementedError, lambda: fcode(sign(x)))

def test_issue_3068X():
    eps = pi**-1500
    big = pi**1000
    one = cos(x)**2 + sin(x)**2
    e = big*one - big + eps
    assert sign(simplify(e)) == 1
    for xi in (111, 11, 1, S(1)/10):
        assert sign(e.subs(x, xi)) == 1

def test_basic1():
    assert limit(x, x, oo) == oo
    assert limit(x, x, -oo) == -oo
    assert limit(-x, x, oo) == -oo
    assert limit(x**2, x, -oo) == oo
    assert limit(-x**2, x, oo) == -oo
    assert limit(x*log(x), x, 0, dir="+") == 0
    assert limit(1/x, x, oo) == 0
    assert limit(exp(x), x, oo) == oo
    assert limit(-exp(x), x, oo) == -oo
    assert limit(exp(x)/x, x, oo) == oo
    assert limit(1/x - exp(-x), x, oo) == 0
    assert limit(x + 1/x, x, oo) == oo
    assert limit(x - x**2, x, oo) == -oo
    assert limit((1 + x)**(1 + sqrt(2)), x, 0) == 1
    assert limit((1 + x)**oo, x, 0) == oo
    assert limit((1 + x)**oo, x, 0, dir='-') == 0
    assert limit((1 + x + y)**oo, x, 0, dir='-') == (1 + y)**(oo)
    assert limit(y/x/log(x), x, 0) == -oo*sign(y)
    assert limit(cos(x + y)/x, x, 0) == sign(cos(y))*oo
    assert limit(gamma(1/x + 3), x, oo) == 2
    assert limit(S.NaN, x, -oo) == S.NaN
    assert limit(Order(2)*x, x, S.NaN) == S.NaN
    assert limit(1/(x - 1), x, 1, dir="+") == oo
    assert limit(1/(x - 1), x, 1, dir="-") == -oo
    assert limit(1/(5 - x)**3, x, 5, dir="+") == -oo
    assert limit(1/(5 - x)**3, x, 5, dir="-") == oo
    assert limit(1/sin(x), x, pi, dir="+") == -oo
    assert limit(1/sin(x), x, pi, dir="-") == oo
    assert limit(1/cos(x), x, pi/2, dir="+") == -oo
    assert limit(1/cos(x), x, pi/2, dir="-") == oo
    assert limit(1/tan(x**3), x, (2*pi)**(S(1)/3), dir="+") == oo
    assert limit(1/tan(x**3), x, (2*pi)**(S(1)/3), dir="-") == -oo
    assert limit(1/cot(x)**3, x, (3*pi/2), dir="+") == -oo
    assert limit(1/cot(x)**3, x, (3*pi/2), dir="-") == oo

    # test bi-directional limits
    assert limit(sin(x)/x, x, 0, dir="+-") == 1
    assert limit(x**2, x, 0, dir="+-") == 0
    assert limit(1/x**2, x, 0, dir="+-") == oo

    # test failing bi-directional limits
    raises(ValueError, lambda: limit(1/x, x, 0, dir="+-"))
    # approaching 0
    # from dir="+"
    assert limit(1 + 1/x, x, 0) == oo
    # from dir='-'
    # Add
    assert limit(1 + 1/x, x, 0, dir='-') == -oo
    # Pow
    assert limit(x**(-2), x, 0, dir='-') == oo
    assert limit(x**(-3), x, 0, dir='-') == -oo
    assert limit(1/sqrt(x), x, 0, dir='-') == (-oo)*I
    assert limit(x**2, x, 0, dir='-') == 0
    assert limit(sqrt(x), x, 0, dir='-') == 0
    assert limit(x**-pi, x, 0, dir='-') == oo*sign((-1)**(-pi))
    assert limit((1 + cos(x))**oo, x, 0) == oo

def base_solution_linear(c, a, b, t=None):
    """
    Return the base solution for a linear diophantine equation with two
    variables. Called repeatedly by diop_linear().

    Usage
    =====

        base_solution_linear(c, a, b, t) -> a, b, c are Integers as in a*x + b*y = c
        and t is the parameter to be used in the solution.

    Details
    =======

        ``c`` is the constant term in a*x + b*y = c
        ``a`` is the integer coefficient of x in a*x + b*y = c
        ``b`` is the integer coefficient of y in a*x + b*y = c
        ``t`` is the parameter to be used in the solution


    Examples
    ========

    >>> from sympy.solvers.diophantine import base_solution_linear
    >>> from sympy.abc import t
    >>> base_solution_linear(5, 2, 3) # equation 2*x + 3*y = 5
    (-5, 5)
    >>> base_solution_linear(0, 5, 7) # equation 5*x + 7*y = 0
    (0, 0)
    >>> base_solution_linear(5, 2, 3, t) # equation 2*x + 3*y = 5
    (3*t - 5, -2*t + 5)
    >>> base_solution_linear(0, 5, 7, t) # equation 5*x + 7*y = 0
    (7*t, -5*t)
    """
    d = igcd(a, igcd(b, c))
    a = a // d
    b = b // d
    c = c // d

    if c == 0:
        if t != None:
            return (b*t , -a*t)
        else:
            return (S.Zero, S.Zero)
    else:
        x0, y0, d = extended_euclid(int(abs(a)), int(abs(b)))

        x0 = x0 * sign(a)
        y0 = y0 * sign(b)

        if divisible(c, d):
            if t != None:
                return (c*x0 + b*t, c*y0 - a*t)
            else:
                return (Integer(c*x0), Integer(c*y0))
        else:
            return (None, None)

def test_basic1():
    assert limit(x, x, oo) == oo
    assert limit(x, x, -oo) == -oo
    assert limit(-x, x, oo) == -oo
    assert limit(x**2, x, -oo) == oo
    assert limit(-x**2, x, oo) == -oo
    assert limit(x*log(x), x, 0, dir="+") == 0
    assert limit(1/x, x, oo) == 0
    assert limit(exp(x), x, oo) == oo
    assert limit(-exp(x), x, oo) == -oo
    assert limit(exp(x)/x, x, oo) == oo
    assert limit(1/x - exp(-x), x, oo) == 0
    assert limit(x + 1/x, x, oo) == oo
    assert limit(x - x**2, x, oo) == -oo
    assert limit((1 + x)**(1 + sqrt(2)), x, 0) == 1
    assert limit((1 + x)**oo, x, 0) == oo
    assert limit((1 + x)**oo, x, 0, dir='-') == 0
    assert limit((1 + x + y)**oo, x, 0, dir='-') == (1 + y)**(oo)
    assert limit(y/x/log(x), x, 0) == -oo*sign(y)
    assert limit(cos(x + y)/x, x, 0) == sign(cos(y))*oo
    raises(NotImplementedError, lambda: limit(Sum(1/x, (x, 1, y)) -
           log(y), y, oo))
    raises(NotImplementedError, lambda: limit(Sum(1/x, (x, 1, y)) - 1/y, y, oo))
    assert limit(gamma(1/x + 3), x, oo) == 2
    assert limit(S.NaN, x, -oo) == S.NaN
    assert limit(Order(2)*x, x, S.NaN) == S.NaN
    assert limit(gamma(1/x + 3), x, oo) == 2
    assert limit(S.NaN, x, -oo) == S.NaN
    assert limit(Order(2)*x, x, S.NaN) == S.NaN
    assert limit(1/(x - 1), x, 1, dir="+") == oo
    assert limit(1/(x - 1), x, 1, dir="-") == -oo
    assert limit(1/(5 - x)**3, x, 5, dir="+") == -oo
    assert limit(1/(5 - x)**3, x, 5, dir="-") == oo
    assert limit(1/sin(x), x, pi, dir="+") == -oo
    assert limit(1/sin(x), x, pi, dir="-") == oo
    assert limit(1/cos(x), x, pi/2, dir="+") == -oo
    assert limit(1/cos(x), x, pi/2, dir="-") == oo
    assert limit(1/tan(x**3), x, (2*pi)**(S(1)/3), dir="+") == oo
    assert limit(1/tan(x**3), x, (2*pi)**(S(1)/3), dir="-") == -oo
    assert limit(1/cot(x)**3, x, (3*pi/2), dir="+") == -oo
    assert limit(1/cot(x)**3, x, (3*pi/2), dir="-") == oo

    # approaching 0
    # from dir="+"
    assert limit(1 + 1/x, x, 0) == oo
    # from dir='-'
    # Add
    assert limit(1 + 1/x, x, 0, dir='-') == -oo
    # Pow
    assert limit(x**(-2), x, 0, dir='-') == oo
    assert limit(x**(-3), x, 0, dir='-') == -oo
    assert limit(1/sqrt(x), x, 0, dir='-') == (-oo)*I
    assert limit(x**2, x, 0, dir='-') == 0
    assert limit(sqrt(x), x, 0, dir='-') == 0
    assert limit((1 + cos(x))**oo, x, 0) == oo

def test_issue_6167_6151():
    n = pi**1000
    i = int(n)
    assert sign(n - i) == 1
    assert abs(n - i) == n - i
    eps = pi**-1500
    big = pi**1000
    one = cos(x)**2 + sin(x)**2
    e = big*one - big + eps
    assert sign(simplify(e)) == 1
    for xi in (111, 11, 1, S(1)/10):
        assert sign(e.subs(x, xi)) == 1

 def eval(cls, arg):
     if not arg.is_Atom:
         c, arg_ = factor_terms(arg).as_coeff_Mul()
         if arg_.is_Mul:
             arg_ = Mul(*[a if (sign(a) not in (-1, 1)) else
                 sign(a) for a in arg_.args])
         arg_ = sign(c)*arg_
     else:
         arg_ = arg
     x, y = re(arg_), im(arg_)
     rv = atan2(y, x)
     if rv.is_number and not rv.atoms(AppliedUndef):
         return rv
     if arg_ != arg:
         return cls(arg_, evaluate=False)

 def _updated_range(r, first):
     st = sign(r.step)*step
     if r.start.is_finite:
         rv = Range(first, r.stop, st)
     else:
         rv = Range(r.start, first + st, st)
     return rv

def frictionforce(rbtdef, ifunc=None):
    """Generate friction forces (Coulomb/viscouse model plus offset)."""

    if not ifunc:
        ifunc = identity

    fric = zeros((rbtdef.dof, 1))

    if rbtdef.frictionmodel is None or len(rbtdef.frictionmodel) == 0:
        pass
    else:
        askedterms = set(rbtdef.frictionmodel)
        if askedterms.issubset(_frictionterms):
            if "viscous" in askedterms:
                for i in range(rbtdef.dof):
                    fric[i] += rbtdef.fv[i] * rbtdef.dq[i]
            if "Coulomb" in askedterms:
                for i in range(rbtdef.dof):
                    fric[i] += rbtdef.fc[i] * sign(rbtdef.dq[i])
            if "offset" in askedterms:
                for i in range(rbtdef.dof):
                    fric[i] += rbtdef.fo[i]
            fric[i] = ifunc(fric[i])
        else:
            raise Exception(
                "Friction model terms '%s' not understanded. Use None or a"
                " combination of %s." % (str(askedterms - _frictionterms), _frictionterms)
            )

    return fric

def frictionforce(rbtdef, ifunc=None):
    '''Generate friction forces (Coulomb/viscouse model plus offset).'''

    if not ifunc:
        ifunc = identity

    fric = zeros(rbtdef.dof, 1)

    if rbtdef.frictionmodel is None or len(rbtdef.frictionmodel) == 0:
        pass
    else:
        askedterms = set(rbtdef.frictionmodel)
        if askedterms.issubset(_frictionterms):
            if 'viscous' in askedterms:
                for i in range(rbtdef.dof):
                    fric[i] += rbtdef.fv[i] * rbtdef.dq[i]
            if 'Coulomb' in askedterms:
                for i in range(rbtdef.dof):
                    fric[i] += rbtdef.fc[i] * sign(rbtdef.dq[i])
            if 'offset' in askedterms:
                for i in range(rbtdef.dof):
                    fric[i] += rbtdef.fo[i]
            fric[i] = ifunc(fric[i])
        else:
            raise Exception(
                'Friction model terms \'%s\' not understanded. Use None or a'
                ' combination of %s.' %
                (str(askedterms - _frictionterms), _frictionterms))

    return fric

    def eval(cls, t, tprime, d_i, d_j, l):
        if t.is_Number and tprime.is_Number and d_i.is_Number and d_j.is_Number and l.is_Number:

            diff_t = t - tprime
            l2 = l * l
            h = h(t, tprime, d_i, d_j, l)
            half_l_di = 0.5 * l * d_i
            arg_1 = half_l_di + tprime / l
            arg_2 = half_l_di - (t - tprime) / l
            ln_part_1 = ln_diff_erf(arg_1, arg_2)
            arg_1 = half_l_di
            arg_2 = half_l_di - t / l
            sign_val = sign(t / l)
            ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t / l)

            base = (
                (0.5 * d_i * l2 * (d_i + d_j) - 1) * h
                + (
                    -diff_t * sign_val * exp(half_l_di * half_l_di - d_i * diff_t + ln_part_1)
                    + t * sign_val * exp(half_l_di * half_l_di - d_i * t - d_j * tprime + ln_part_2)
                )
                + l
                / sqrt(pi)
                * (
                    -exp(-diff_t * diff_t / l2)
                    + exp(-tprime * tprime / l2 - d_i * t)
                    + exp(-t * t / l2 - d_j * tprime)
                    - exp(-(d_i * t + d_j * tprime))
                )
            )
            return base / (d_i + d_j)

    def __new__(cls, lhs, rhs, rop=None, **assumptions):
        lhs = _sympify(lhs)
        rhs = _sympify(rhs)
        if cls is not Relational:
            rop_cls = cls
        else:
            try:
                rop_cls = Relational.ValidRelationOperator[rop]
            except KeyError:
                msg = "Invalid relational operator symbol: '%r'"
                raise ValueError(msg % repr(rop))
        if lhs.is_number and rhs.is_number and (rop_cls in (Equality, Unequality) or lhs.is_real and rhs.is_real):
            diff = lhs - rhs
            know = diff.equals(0, failing_expression=True)
            if know is True:  # exclude failing expression case
                Nlhs = S.Zero
            elif know is False:
                from sympy import sign

                Nlhs = sign(diff.n(1))
            else:
                Nlhs = None
                lhs = know
                rhs = S.Zero
            if Nlhs is not None:
                return rop_cls._eval_relation(Nlhs, S.Zero)

        obj = Expr.__new__(rop_cls, lhs, rhs, **assumptions)
        return obj

def _compute_tau(model, robo, j):
    """
    Compute the joint torque by subtracting the effect of friction.

    Args:
        model: An instance of DynModel
        robo: An instance of Robot
        j: joint number

    Returns:
        An instance of DynModel that contains all the new values.
    """
    tau_j = 0
    # local variables
    qdot_j = robo.qdots[j]
    gamma_j = robo.torques[j]
    f_cj = robo.dyns[j].frc
    f_vj = robo.dyns[j].frv
    # actual computation
    coriolis_friction_term = f_cj * sign(qdot_j)
    viscous_friction_term = f_vj * qdot_j
    tau_j = gamma_j - coriolis_friction_term - viscous_friction_term
    # store in model
    model.taus[j] = tau_j
    return model

 def _eval_derivative(self, x):
     if self.args[0].is_real or self.args[0].is_imaginary:
         return Derivative(self.args[0], x, evaluate=True) \
             * sign(conjugate(self.args[0]))
     return (re(self.args[0]) * Derivative(re(self.args[0]), x,
         evaluate=True) + im(self.args[0]) * Derivative(im(self.args[0]),
             x, evaluate=True)) / Abs(self.args[0])

def _compute_joint_torque(model, robo, j):
    """
    Compute the joint torque for joint j.

    Args:
        model: An instance of DynModel
        robo: An instance of Robot
        j: joint number

    Returns:
        An instance of DynModel that contains all the new values.
    """
    # local variables
    qdot_j = robo.qdots[j]
    qddot_j = robo.qddots[j]
    j_a_j = robo.geos[j].axisa
    ia_j = robo.dyns[j].ia
    f_cj = robo.dyns[j].frc
    f_vj = robo.dyns[j].frv
    j_f_j = model.wrenchs[j].val
    # actual computation
    wrench_term = j_f_j.transpose() * j_a_j
    actuator_inertia_term = Matrix([ia_j * qddot_j])
    coriolis_friction_term = Matrix([f_cj * sign(qdot_j)])
    viscous_friction_term = Matrix([f_vj * qdot_j])
    gamma_j = wrench_term + actuator_inertia_term + \
        viscous_friction_term + coriolis_friction_term
    # store computed torque in model
    model.torques[j] = gamma_j[0, 0]
    return model

 def _eval_nseries(self, x, n, logx):
     direction = self.args[0].leadterm(x)[0]
     s = self.args[0]._eval_nseries(x, n=n, logx=logx)
     when = Eq(direction, 0)
     return Piecewise(
         ((s.subs(direction, 0)), when),
         (sign(direction)*s, True),
     )

def test_Abs():
    x, y = symbols('x,y')
    assert Abs(0) == 0
    assert Abs(1) == 1
    assert Abs(-1)== 1
    assert Abs(nan) == nan
    x = Symbol('x',real=True)
    n = Symbol('n',integer=True)
    assert x**(2*n) == Abs(x)**(2*n)
    assert Abs(x).diff(x) == sign(x)
    assert abs(x) == Abs(x) # Python built-in
    assert Abs(x)**3 == x**2*Abs(x)
    assert (Abs(x)**(3*n)).args == (Abs(x), 3*n) # leave symbolic odd unchanged
    assert (1/Abs(x)).args == (Abs(x), -1)
    assert 1/Abs(x)**3 == 1/(x**2*Abs(x))
    assert Abs(x).diff(x) == sign(x)
    assert conjugate(Abs(x)) == Abs(x)