def test_add_1():
    assert Order(x + x) == Order(x)
    assert Order(3 * x - 2 * x ** 2) == Order(x)
    assert Order(1 + x) == Order(1, x)
    assert Order(1 + 1 / x) == Order(1 / x)
    assert Order(ln(x) + 1 / ln(x)) == Order(ln(x))
    assert Order(exp(1 / x) + x) == Order(exp(1 / x))
    assert Order(exp(1 / x) + 1 / x ** 20) == Order(exp(1 / x))

def test_ansi_math1_codegen():
    # not included: log10
    from sympy import acos, asin, atan, ceiling, cos, cosh, floor, log, ln, \
        sin, sinh, sqrt, tan, tanh, N
    x = symbols('x')
    name_expr = [
        ("test_fabs", abs(x)),
        ("test_acos", acos(x)),
        ("test_asin", asin(x)),
        ("test_atan", atan(x)),
        ("test_ceil", ceiling(x)),
        ("test_cos", cos(x)),
        ("test_cosh", cosh(x)),
        ("test_floor", floor(x)),
        ("test_log", log(x)),
        ("test_ln", ln(x)),
        ("test_sin", sin(x)),
        ("test_sinh", sinh(x)),
        ("test_sqrt", sqrt(x)),
        ("test_tan", tan(x)),
        ("test_tanh", tanh(x)),
    ]
    numerical_tests = []
    for name, expr in name_expr:
        for xval in 0.2, 0.5, 0.8:
            expected = N(expr.subs(x, xval))
            numerical_tests.append((name, (xval,), expected, 1e-14))
    run_cc_test("ansi_math1", name_expr, numerical_tests)

def test_intrinsic_math1_codegen():
    # not included: log10
    from sympy import acos, asin, atan, ceiling, cos, cosh, floor, log, ln, \
        sin, sinh, sqrt, tan, tanh, N
    x = symbols('x')
    name_expr = [
        ("test_fabs", abs(x)),
        ("test_acos", acos(x)),
        ("test_asin", asin(x)),
        ("test_atan", atan(x)),
        ("test_cos", cos(x)),
        ("test_cosh", cosh(x)),
        ("test_log", log(x)),
        ("test_ln", ln(x)),
        ("test_sin", sin(x)),
        ("test_sinh", sinh(x)),
        ("test_sqrt", sqrt(x)),
        ("test_tan", tan(x)),
        ("test_tanh", tanh(x)),
    ]
    numerical_tests = []
    for name, expr in name_expr:
        for xval in 0.2, 0.5, 0.8:
            expected = N(expr.subs(x, xval))
            numerical_tests.append((name, (xval,), expected, 1e-14))
    for lang, commands in valid_lang_commands:
        if lang == "C":
            name_expr_C = [("test_floor", floor(x)), ("test_ceil", ceiling(x))]
        else:
            name_expr_C = []
        run_test("intrinsic_math1", name_expr + name_expr_C, numerical_tests, lang, commands)

    def eta_fil(self, x, V_app, apprx=(0, 0, 0, 0)):
        m_eff = self.m_r * const.electron_mass

        mpmath.mp.dps = 20
        x0 = Symbol('x0')  # eta_fil
        x1 = Symbol('x1')  # eta_ac
        x2 = Symbol('x2')  # eta_hop
        x3 = Symbol('x3')  # V_tunnel

        f0 = const.Boltzmann * self.T / (1 - self.alpha) / const.elementary_charge / self.z * \
             ln(self.A_fil/self.A_ac*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1) + 1) - x1# eta_ac = f(eta_fil) x1 = f(x0)
        f1 = x*2*const.Boltzmann*self.T/self.a/self.z/const.elementary_charge*\
             asinh(self.j_0et/self.j_0hop*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1)) - x2# eta_hop = f(eta_fil)
        f2 = x1 - x0 + x2 - x3

        f3 = -V_app + ((self.C * 3 * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
             exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge))) * self.A_fil*x3)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*x0/const.Boltzmann/self.T) - 1))) * (self.R_el + self.R_S + self.rho_fil*(self.L - x) / self.A_fil) \
             + x3

        eta_fil, eta_ac, eta_hop, V_tunnel = nsolve((f0, f1, f2, f3), [x0, x1, x2, x3], apprx)
        eta_fil = np.real(np.complex128(eta_fil))
        eta_ac = np.real(np.complex128(eta_ac))
        eta_hop = np.real(np.complex128(eta_hop))
        V_tunnel = np.real(np.complex128(V_tunnel))
        current = ((self.C * 3 * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
            exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge))) * self.A_fil*V_tunnel)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*eta_fil/const.Boltzmann/self.T) - 1)))
        print(eta_fil, eta_ac, eta_hop, V_tunnel)
        # print(eta_ac - eta_fil + eta_hop - V_tunnel)
        return eta_fil, eta_ac, eta_hop, V_tunnel, current

def test_solve_relation_for_derivatives():
    from sympy import symbols, Function, ln

    x, y = symbols("x, y")
    f = symbols("f", cls=Function)(x, y)
    relation = ln(x * y + f) - f
    subsd = {x: 2, y: 1}
    diff_wrt = {x: 2, y: 1}
    df, df_err = solve_relation_for_derivatives(
        relation, {x: 2, y: 1}, f, {((x, 0), (y, 0)): 1.0}, diff_wrt, verbose=True, abstol=1e-8
    )
    assert abs(df[((x, 1), (y, 0))] - 0.46594127238) < 1e-8
    assert abs(df[((x, 0), (y, 1))] - 0.9318825447699826) < 1e-8
    assert abs(df[((x, 1), (y, 1))] + 0.17057414955751987) < 1e-8
    assert abs(df[((x, 2), (y, 1))] - 0.37120710736327495) < 1e-8

    x, y = symbols("x, y")
    # x_arg, y_arg = symbols('x_arg, y_arg')
    g = symbols("g", cls=Function)(x, y)  # (x_arg, y_arg)
    relation = x ** 2 + y ** 2 - g  # Circle
    subsd = {x: 1.0, y: 0.0}
    initial_guess_func_vals = {((x, 0),): 0.5}
    diff_wrt = {x: 1}
    dg, dg_err = solve_relation_for_derivatives(
        relation, subsd, g, initial_guess_func_vals, diff_wrt, verbose=True, abstol=1e-7
    )
    assert abs(dg[((x, 1),)] - 2.0) < 1e-7

 def test_general_symbolic(self):
     a, b, c, d, e, x = syp.symbols('a,b,c,d,e,x')
     # natural log function
     y =  a*syp.ln(b*x**c + d)+e
     phys = self.phys
     phys['pore.item1'] = 0.16e-14
     phys['pore.item2'] = 4
     phys['pore.item3'] = 1.4
     phys['pore.item4'] = 0.133
     phys['pore.item5'] = -5.1e-14
     phys.add_model(propname='pore.source1',
                    model=pm.generic_source_term.natural_logarithm,
                    A1='pore.item1',
                    A2='pore.item2',
                    A3='pore.item3',
                    A4='pore.item4',
                    A5='pore.item5',
                    X='pore.mole_fraction',
                    regen_mode='normal')
     arg_map=collections.OrderedDict([('a', 'pore.item1'),
                                      ('b', 'pore.item2'),
                                      ('c', 'pore.item3'),
                                      ('d', 'pore.item4'),
                                      ('e', 'pore.item5'),
                                      ('x', 'pore.mole_fraction')])
     phys.add_model(propname='pore.general',
                    model=op.models.physics.generic_source_term.general_symbolic,
                    eqn=y, arg_map=arg_map,
                    regen_mode='normal')
     assert np.allclose(phys['pore.source1.rate'], phys['pore.general.rate'])
     assert np.allclose(phys['pore.source1.S1'], phys['pore.general.S1'])
     assert np.allclose(phys['pore.source1.S2'], phys['pore.general.S2'])

def euro_call_exact(S, K, T, r, sigma):

    # S: spot price
    # K: strike price
    # T: time to maturity
    # r: interest rate
    # sigma: volatility of underlying asset

    N = systats.Normal(0.0, 1.0)

    d1 = (sy.ln(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
    d2 = (sy.ln(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

    call = S * N.cdf(d1) - K * sy.exp(-r * T) * N.cdf(d2)

    return call

def exact_put_option_price(S, K, T, r, sigma):
    
    #S: spot price
    #K: strike price
    #T: time to maturity
    #r: interest rate
    #sigma: volatility of underlying asset
    
    N = systats.Normal(0.0, 1.0)
    
    d1 = (sy.ln(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
    d2 = (sy.ln(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
    
    put_price = K * sy.exp(-r * T) * N.cdf(-d2) - S * N.cdf(-d1)
    
    return put_price

def exact_euro_put_div(S, K, T, r, q, sigma):

    # S: spot price
    # K: strike price
    # T: time to maturity
    # r: interest rate
    # q: rate of continuous dividend paying asset
    # sigma: volatility of underlying asset

    N = systats.Normal(0.0, 1.0)

    d1 = (sy.ln(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
    d2 = (sy.ln(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

    put = K * sy.exp(-r * T) * N.cdf(-d2) - S * sy.exp(-q * T) * N.cdf(-d1)

    return put

def getJeffreysPrior(rv):
    """
    Uses SymPy to determine the Jeffreys prior of a random variable analytically.

    Args:
        rv: SymPy RandomSymbol, corresponding to a probability distribution

    Returns:
        list: List containing Jeffreys prior in symbolic form and corresponding lambda function

    Example:
        rate = Symbol('rate', positive=True)
        rv = stats.Exponential('exponential', rate)
        print getJeffreysPrior(rv)

        >>> (1/rate, <function <lambda> at 0x0000000007F79AC8>)
    """

    # get support of random variable
    support = rv._sorted_args[0].distribution.set

    # get list of free parameters
    parameters = list(rv._sorted_args[0].distribution.free_symbols)
    x = abc.x

    # symbolic probability density function
    symPDF = density(rv)(x)

    # compute Fisher information matrix
    dim = len(parameters)
    G = Matrix.zeros(dim, dim)

    func = summation if support.is_iterable else integrate
    for i in range(0, dim):
        for j in range(0, dim):
            G[i, j] = func(simplify(symPDF *
                                    diff(ln(symPDF), parameters[i]) *
                                    diff(ln(symPDF), parameters[j])),
                           (x, support.inf, support.sup))

    # symbolic Jeffreys prior
    symJeff = simplify(sqrt(G.det()))

    # return symbolic Jeffreys prior and corresponding lambda function
    return symJeff, lambdify(parameters, symJeff, 'numpy')

def get_functions(var):
    # Denne generator laver alle spørgsmålene.
    yield exp(a*var)
    yield ln(a*var)
    yield a*(var**4)
    yield a/(var**4)
    yield sqrt(1+(var**4))
    yield (1+var)**a
    yield exp(a*var**2)

def euro_vanilla_exact(S, K, T, r, sigma, option="call"):

    # S: spot price
    # K: strike price
    # T: time to maturity
    # r: interest rate
    # sigma: volatility of underlying asset

    N = systats.Normal(0.0, 1.0)

    d1 = (sy.ln(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
    d2 = (sy.ln(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

    if option == "call":
        result = S * N.cdf(d1) - K * sy.exp(-r * T) * N.cdf(d2)
    if option == "put":
        result = K * sy.exp(-r * T) * N.cdf(-d2) - S * N.cdf(-d1)

    return result

def test_ansi_math1_codegen():
    # not included: log10
    from sympy import (acos, asin, atan, ceiling, cos, cosh, floor, log, ln,
        sin, sinh, sqrt, tan, tanh, Abs)
    x = symbols('x')
    name_expr = [
        ("test_fabs", Abs(x)),
        ("test_acos", acos(x)),
        ("test_asin", asin(x)),
        ("test_atan", atan(x)),
        ("test_ceil", ceiling(x)),
        ("test_cos", cos(x)),
        ("test_cosh", cosh(x)),
        ("test_floor", floor(x)),
        ("test_log", log(x)),
        ("test_ln", ln(x)),
        ("test_sin", sin(x)),
        ("test_sinh", sinh(x)),
        ("test_sqrt", sqrt(x)),
        ("test_tan", tan(x)),
        ("test_tanh", tanh(x)),
    ]
    result = codegen(name_expr, "C89", "file", header=False, empty=False)
    assert result[0][0] == "file.c"
    assert result[0][1] == (
        '#include "file.h"\n#include <math.h>\n'
        'double test_fabs(double x) {\n   double test_fabs_result;\n   test_fabs_result = fabs(x);\n   return test_fabs_result;\n}\n'
        'double test_acos(double x) {\n   double test_acos_result;\n   test_acos_result = acos(x);\n   return test_acos_result;\n}\n'
        'double test_asin(double x) {\n   double test_asin_result;\n   test_asin_result = asin(x);\n   return test_asin_result;\n}\n'
        'double test_atan(double x) {\n   double test_atan_result;\n   test_atan_result = atan(x);\n   return test_atan_result;\n}\n'
        'double test_ceil(double x) {\n   double test_ceil_result;\n   test_ceil_result = ceil(x);\n   return test_ceil_result;\n}\n'
        'double test_cos(double x) {\n   double test_cos_result;\n   test_cos_result = cos(x);\n   return test_cos_result;\n}\n'
        'double test_cosh(double x) {\n   double test_cosh_result;\n   test_cosh_result = cosh(x);\n   return test_cosh_result;\n}\n'
        'double test_floor(double x) {\n   double test_floor_result;\n   test_floor_result = floor(x);\n   return test_floor_result;\n}\n'
        'double test_log(double x) {\n   double test_log_result;\n   test_log_result = log(x);\n   return test_log_result;\n}\n'
        'double test_ln(double x) {\n   double test_ln_result;\n   test_ln_result = log(x);\n   return test_ln_result;\n}\n'
        'double test_sin(double x) {\n   double test_sin_result;\n   test_sin_result = sin(x);\n   return test_sin_result;\n}\n'
        'double test_sinh(double x) {\n   double test_sinh_result;\n   test_sinh_result = sinh(x);\n   return test_sinh_result;\n}\n'
        'double test_sqrt(double x) {\n   double test_sqrt_result;\n   test_sqrt_result = sqrt(x);\n   return test_sqrt_result;\n}\n'
        'double test_tan(double x) {\n   double test_tan_result;\n   test_tan_result = tan(x);\n   return test_tan_result;\n}\n'
        'double test_tanh(double x) {\n   double test_tanh_result;\n   test_tanh_result = tanh(x);\n   return test_tanh_result;\n}\n'
    )
    assert result[1][0] == "file.h"
    assert result[1][1] == (
        '#ifndef PROJECT__FILE__H\n#define PROJECT__FILE__H\n'
        'double test_fabs(double x);\ndouble test_acos(double x);\n'
        'double test_asin(double x);\ndouble test_atan(double x);\n'
        'double test_ceil(double x);\ndouble test_cos(double x);\n'
        'double test_cosh(double x);\ndouble test_floor(double x);\n'
        'double test_log(double x);\ndouble test_ln(double x);\n'
        'double test_sin(double x);\ndouble test_sinh(double x);\n'
        'double test_sqrt(double x);\ndouble test_tan(double x);\n'
        'double test_tanh(double x);\n#endif\n'
    )

def test_ansi_math1_codegen():
    # not included: log10
    from sympy import acos, asin, atan, ceiling, cos, cosh, floor, log, ln, sin, sinh, sqrt, tan, tanh, N, Abs

    x = symbols("x")
    name_expr = [
        ("test_fabs", Abs(x)),
        ("test_acos", acos(x)),
        ("test_asin", asin(x)),
        ("test_atan", atan(x)),
        ("test_ceil", ceiling(x)),
        ("test_cos", cos(x)),
        ("test_cosh", cosh(x)),
        ("test_floor", floor(x)),
        ("test_log", log(x)),
        ("test_ln", ln(x)),
        ("test_sin", sin(x)),
        ("test_sinh", sinh(x)),
        ("test_sqrt", sqrt(x)),
        ("test_tan", tan(x)),
        ("test_tanh", tanh(x)),
    ]
    result = codegen(name_expr, "C", "file", header=False, empty=False)
    assert result[0][0] == "file.c"
    assert result[0][1] == (
        '#include "file.h"\n#include <math.h>\n'
        "double test_fabs(double x) {\n   return fabs(x);\n}\n"
        "double test_acos(double x) {\n   return acos(x);\n}\n"
        "double test_asin(double x) {\n   return asin(x);\n}\n"
        "double test_atan(double x) {\n   return atan(x);\n}\n"
        "double test_ceil(double x) {\n   return ceil(x);\n}\n"
        "double test_cos(double x) {\n   return cos(x);\n}\n"
        "double test_cosh(double x) {\n   return cosh(x);\n}\n"
        "double test_floor(double x) {\n   return floor(x);\n}\n"
        "double test_log(double x) {\n   return log(x);\n}\n"
        "double test_ln(double x) {\n   return log(x);\n}\n"
        "double test_sin(double x) {\n   return sin(x);\n}\n"
        "double test_sinh(double x) {\n   return sinh(x);\n}\n"
        "double test_sqrt(double x) {\n   return sqrt(x);\n}\n"
        "double test_tan(double x) {\n   return tan(x);\n}\n"
        "double test_tanh(double x) {\n   return tanh(x);\n}\n"
    )
    assert result[1][0] == "file.h"
    assert result[1][1] == (
        "#ifndef PROJECT__FILE__H\n#define PROJECT__FILE__H\n"
        "double test_fabs(double x);\ndouble test_acos(double x);\n"
        "double test_asin(double x);\ndouble test_atan(double x);\n"
        "double test_ceil(double x);\ndouble test_cos(double x);\n"
        "double test_cosh(double x);\ndouble test_floor(double x);\n"
        "double test_log(double x);\ndouble test_ln(double x);\n"
        "double test_sin(double x);\ndouble test_sinh(double x);\n"
        "double test_sqrt(double x);\ndouble test_tan(double x);\n"
        "double test_tanh(double x);\n#endif\n"
    )

def natural_logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
    r"""
    Calculates the rate, as well as slope and intercept of the following
    function at the given value of *x*:

        .. math::
            rate =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}

    Parameters
    ----------
    A1 -> A5 : string
        The dictionary keys on the target object containing the coefficients
        values to be used in the source term model

    X : string or float/int or array/list
        The dictionary key on the target objecxt containing the the quantity
        of interest

    Returns
    -------
    A dictionary containing the following three items:

        **'rate'** - The value of the source term function at the given X.

        **'S1'** - The slope of the source term function at the given X.

        **'S2'** - The intercept of the source term function at the given X.

    The slope and intercept provide a linearized source term equation about the
    current value of X as follow:

        .. math::
            rate = S_{1}   X  +  S_{2}

    """
    A = _parse_args(target=target, key=A1, default=1.0)
    B = _parse_args(target=target, key=A2, default=1.0)
    C = _parse_args(target=target, key=A3, default=1.0)
    D = _parse_args(target=target, key=A4, default=1.0)
    E = _parse_args(target=target, key=A5, default=0.0)
    X = target[X]
    # Symbols used in symbolic function
    a, b, c, d, e, x = _syp.symbols('a,b,c,d,e,x')
    # Equation
    y = a*_syp.ln(b*x**c + d) + e
    # Callable functions
    r, s1, s2 = _build_func(eq=y, a=a, b=b, c=c, d=d, e=e, x=x)
    # Values
    r_val = r(A, B, C, D, E, X)
    s1_val = s1(A, B, C, D, E, X)
    s2_val = s2(A, B, C, D, E, X)
    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
    return values

def test_order_at_infinity():
    assert Order(1 + x, (x, oo)) == Order(x, (x, oo))
    assert Order(3*x, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo))*3 == Order(x, (x, oo))
    assert -28*Order(x, (x, oo)) == Order(x, (x, oo))
    assert Order(Order(x, (x, oo)), (x, oo)) == Order(x, (x, oo))
    assert Order(Order(x, (x, oo)), (y, oo)) == Order(x, (x, oo), (y, oo))
    assert Order(3, (x, oo)) == Order(1, (x, oo))
    assert Order(x**2 + x + y, (x, oo)) == O(x**2, (x, oo))
    assert Order(x**2 + x + y, (y, oo)) == O(y, (y, oo))

    assert Order(2*x, (x, oo))*x == Order(x**2, (x, oo))
    assert Order(2*x, (x, oo))/x == Order(1, (x, oo))
    assert Order(2*x, (x, oo))*x*exp(1/x) == Order(x**2*exp(1/x), (x, oo))
    assert Order(2*x, (x, oo))*x*exp(1/x)/ln(x)**3 == Order(x**2*exp(1/x)*ln(x)**-3, (x, oo))

    assert Order(x, (x, oo)) + 1/x == 1/x + Order(x, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo)) + 1 == 1 + Order(x, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo)) + x == x + Order(x, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo)) + x**2 == x**2 + Order(x, (x, oo))
    assert Order(1/x, (x, oo)) + 1/x**2 == 1/x**2 + Order(1/x, (x, oo)) == Order(1/x, (x, oo))
    assert Order(x, (x, oo)) + exp(1/x) == exp(1/x) + Order(x, (x, oo))

    assert Order(x, (x, oo))**2 == Order(x**2, (x, oo))

    assert Order(x, (x, oo)) + Order(x**2, (x, oo)) == Order(x**2, (x, oo))
    assert Order(x, (x, oo)) + Order(x**-2, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo)) + Order(1/x, (x, oo)) == Order(x, (x, oo))

    assert Order(x, (x, oo)) - Order(x, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo)) + Order(1, (x, oo)) == Order(x, (x, oo))
    assert Order(x, (x, oo)) + Order(x**2, (x, oo)) == Order(x**2, (x, oo))
    assert Order(1/x, (x, oo)) + Order(1, (x, oo)) == Order(1, (x, oo))
    assert Order(x, (x, oo)) + Order(exp(1/x), (x, oo)) == Order(x, (x, oo))
    assert Order(x**3, (x, oo)) + Order(exp(2/x), (x, oo)) == Order(x**3, (x, oo))
    assert Order(x**-3, (x, oo)) + Order(exp(2/x), (x, oo)) == Order(exp(2/x), (x, oo))

    # issue 7207
    assert Order(exp(x), (x, oo)).expr == Order(2*exp(x), (x, oo)).expr == exp(x)
    assert Order(y**x, (x, oo)).expr == Order(2*y**x, (x, oo)).expr == exp(log(y)*x)

def test_order_at_infinity():
    assert Order(1 + x, x, oo) == Order(x, x, oo)
    assert Order(3*x, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo)*3 == Order(x, x, oo)
    assert -28*Order(x, x, oo) == Order(x, x, oo)
    assert Order(Order(x, x, oo)) == Order(x, x, oo)
    assert Order(Order(x, x, oo), y) == Order(Order(x, x, oo), x, y)
    assert Order(3, x, oo) == Order(1, x, oo)
    assert Order(x**2 + x + y, x, oo) == O(x**2, x, oo)
    assert Order(x**2 + x + y, y, oo) == O(y, y, oo)

    assert Order(2*x, x, oo)*x == Order(x**2, x, oo)
    assert Order(2*x, x, oo)/x == Order(1, x, oo)
    assert Order(2*x, x, oo)*x*exp(1/x) == Order(x**2*exp(1/x), x, oo)
    assert Order(2*x, x, oo)*x*exp(1/x)/ln(x)**3 == Order(x**2*exp(1/x)*ln(x)**-3, x, oo)

    assert Order(x, x, oo) + 1/x == 1/x + Order(x, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo) + 1 == 1 + Order(x, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo) + x == x + Order(x, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo) + x**2 == x**2 + Order(x, x, oo)
    assert Order(1/x, x, oo) + 1/x**2 == 1/x**2 + Order(1/x, x, oo) == Order(1/x, x, oo)
    assert Order(x, x, oo) + exp(1/x) == exp(1/x) + Order(x, x, oo)

    assert Order(x, x, oo)**2 == Order(x**2, x, oo)
    assert Order(x**3, x, oo)**-2 == Order(x**-6, x, oo)

    assert Order(x, x, oo) + Order(x**2, x, oo) == Order(x**2, x, oo)
    assert Order(x, x, oo) + Order(x**-2, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo) + Order(1/x, x, oo) == Order(x, x, oo)

    assert Order(x, x, oo) - Order(x, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo) + Order(1, x, oo) == Order(x, x, oo)
    assert Order(x, x, oo) + Order(x**2, x, oo) == Order(x**2, x, oo)
    assert Order(1/x, x, oo) + Order(1, x, oo) == Order(1, x, oo)
    assert Order(x, x, oo) + Order(exp(1/x), x, oo) == Order(x, x, oo)
    assert Order(x**3, x, oo) + Order(exp(2/x), x, oo) == Order(x**3, x, oo)
    assert Order(x**-3, x, oo) + Order(exp(2/x), x, oo) == Order(exp(2/x), x, oo)

def test_get_cb_from_rel__1():
    from sympy import symbols, Function, ln

    x, y = symbols("x,y")
    f = symbols("f", cls=Function)(x, y)
    relation = ln(x * y + f) - f
    subsd = {x: 2, y: 1}
    subs_fs = frozenset(subsd.items())  # Hashable
    varied = f
    initial_guess = 1.0
    cb = get_cb_from_rel(relation, subs_fs, f)
    assert cb(1.14619322062) < 1e-12

    cb_numexpr = get_cb_from_rel(relation, subs_fs, f, use_numexpr=True)
    assert cb_numexpr(1.14619322062) < 1e-12

    cb_autowrap = get_cb_from_rel(relation, subs_fs, f, (), bincompile=True, unitless=True)
    assert cb_autowrap(1.14619322062) < 1e-12

 def print_Log(self, rule):
     with self.new_step():
         if rule.base == sympy.E:
             self.append("The derivative of {} is {}.".format(
                 self.format_math(rule.context),
                 self.format_math(diff(rule))
             ))
         else:
             # This case shouldn't come up often, seeing as SymPy
             # automatically applies the change-of-base identity
             self.append("The derivative of {} is {}.".format(
                 self.format_math(sympy.log(rule.symbol, rule.base,
                                            evaluate=False)),
                 self.format_math(1/(rule.arg * sympy.ln(rule.base)))))
             self.append("So {}".format(
                 self.format_math(sympy.Eq(
                     sympy.Derivative(rule.context, rule.symbol),
                     diff(rule)))))

def test_solve_relation_num(use_numexpr=False, bincompile=False):
    from sympy import symbols, Function, ln

    x, y = symbols("x,y")
    f = symbols("f", cls=Function)(x, y)
    relation = ln(x * y + f) - f
    subsd = {x: 2, y: 1}
    varied = f
    initial_guess = 1.0
    unitless = True
    x, y = solve_relation_num(
        relation, subsd, varied, initial_guess, (), use_numexpr, unitless, bincompile, verbose=True, yabstol=1e-9
    )

    assert abs(x - 1.14619322062) < 1e-8

    x, y = symbols("x,y")
    g = symbols("g", cls=Function)(x, y)
    relation = x ** 2 + y ** 2 - g  # Circle
    g_val, delta_rel = solve_relation_num(relation, {x: 1, y: 0}, g, 0.1, (), use_numexpr, verbose=True, abstol=1e-8)
    assert abs(g_val - 1) < 1e-7