Python sympy 模块，Symbol() 实例源码

def test_clambdify():
    x = Symbol('x')
    y = Symbol('y')
    z = Symbol('z')
    f1 = sqrt(x*y)
    pf1 = lambdify((x, y), f1, 'math')
    cf1 = clambdify((x, y), f1)
    for i in xrange(10):
        if cf1(i, 10 - i) != pf1(i, 10 - i):
            raise ValueError("Values should be equal")
    f2 = (x - y) / z * pi
    pf2 = lambdify((x, y, z), f2, 'math')
    cf2 = clambdify((x, y, z), f2)
    if round(pf2(1, 2, 3), 14) != round(cf2(1, 2, 3), 14):
        raise ValueError("Values should be equal")
    # FIXME: slight difference in precision

def error_bound(n):
    e = Symbol('e')
    y = Symbol('y')
    growth_function_bound = generate_growth_function_bound(50)
    a = original_vc_bound(n, 0.05, growth_function_bound)
    b = nsolve(Eq(rademacher_penalty_bound(n, 0.05, growth_function_bound), y), 1)
    c = nsolve(Eq(parrondo_van_den_broek_right(e, n, 0.05, growth_function_bound), e), 1)
    d = nsolve(Eq(devroye(e, n, 0.05, growth_function_bound), e), 1)
    return a, b, c, d


# def test_three():
#     e = Symbol('e')
#     y = Symbol('y')
#     n = Symbol('n')
#     growth_function_bound = generate_growth_function_bound(50)
#     a = original_vc_bound(5, 0.05, growth_function_bound)
#     b = nsolve(Eq(rademacher_penalty_bound(5, 0.05, growth_function_bound), y), 5)
#     c = nsolve(Eq(parrondo_van_den_broek_right(e, 5, 0.05, growth_function_bound), e), 1)
#     d = nsolve(Eq(devroye(e, 5, 0.05, growth_function_bound), e), 1)
#     return a, b, c, d

def __init__(self, index):
        self.points = numpy.linspace(-1.0, 1.0, index+1)
        self.degree = index + 1 if index % 2 == 0 else index

        # Formula (26) from
        # <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
        # Note that Sympy carries out all operations in rationals, i.e.,
        # _exactly_. Only at the end, the rational is converted into a float.
        n = index
        self.weights = numpy.empty(n+1)
        t = sympy.Symbol('t')
        for r in range(n+1):
            # Compare with get_weights().
            f = sympy.prod([(t - i) for i in range(n+1) if i != r])
            alpha = 2 * \
                (-1)**(n-r) * sympy.integrate(f, (t, 0, n)) \
                / (math.factorial(r) * math.factorial(n-r)) \
                / index
            self.weights[r] = alpha
        return

def __init__(self, index):
        self.points = numpy.linspace(-1.0, 1.0, index+2)[1:-1]
        self.degree = index if (index+1) % 2 == 0 else index - 1
        #
        n = index+1
        self.weights = numpy.empty(n-1)
        t = sympy.Symbol('t')
        for r in range(1, n):
            # Compare with get_weights().
            f = sympy.prod([(t - i) for i in range(1, n) if i != r])
            alpha = 2 * \
                (-1)**(n-r+1) * sympy.integrate(f, (t, 0, n)) \
                / (math.factorial(r-1) * math.factorial(n-1-r)) \
                / n
            self.weights[r-1] = alpha
        return

def obtain_atomic_formulas(file):
    propositions = set()
    for line in file:
        if line.startswith("("):
            prop_char = set()
            for char in line:
                #print(str(char))
                if(str(char).isalpha()):
                    prop_char.add(str(char))
            for item in prop_char:
                new = Symbol(item)
                propositions.add(new)
    return propositions                     #returns the set of propositions involved in the set of rules


# Parses each line of the rule file to create a a dictionary of rules, distinguishing the item, body and head. The key is the name of the rule
# while the value is the Rule object itself

def obtain_atomic_formulas(file):
    propositions = set()
    for line in file:
        if line.startswith("("):
            prop_char = set()
            for char in line:
                #print(str(char))
                if(str(char).isalpha()):
                    prop_char.add(str(char))
            for item in prop_char:
                new = Symbol(item)
                propositions.add(new)
    return propositions


# Parses each line of the rule file to create a a dictionary of rules, distinguishing the item, body and head. The key is the name of the rule
# while the value is the Rule object itself

def add_proposition(propositions, p):
    _p = p.strip()
    _p = re.sub(r'\s+', '', _p)
    _p = _p.replace("~", "")
    _p = _p.replace("&", ",")
    _p = _p.replace("|", ",")
    _p = _p.replace("(", "")
    _p = _p.replace(")", "")
    _p = _p.replace("->", ",")
    _p = _p.replace("!", "")
    new_props = _p.split(",")
    for prop in new_props:
        if prop == "":
            continue 
        new = Symbol(prop)
        propositions.add(new)

#Used to reconstruct worlds in light of any constraints that might be included in the file

def test_sqrt_symbolic_denest():
    x = Symbol('x')
    z = sqrt(((1 + sqrt(sqrt(2 + x) + 3))**2).expand())
    assert sqrtdenest(z) == sqrt((1 + sqrt(sqrt(2 + x) + 3))**2)
    z = sqrt(((1 + sqrt(sqrt(2 + cos(1)) + 3))**2).expand())
    assert sqrtdenest(z) == 1 + sqrt(sqrt(2 + cos(1)) + 3)
    z = ((1 + cos(2))**4 + 1).expand()
    assert sqrtdenest(z) == z
    z = sqrt(((1 + sqrt(sqrt(2 + cos(3*x)) + 3))**2 + 1).expand())
    assert sqrtdenest(z) == z
    c = cos(3)
    c2 = c**2
    assert sqrtdenest(sqrt(2*sqrt(1 + r3)*c + c2 + 1 + r3*c2)) == \
        -1 - sqrt(1 + r3)*c
    ra = sqrt(1 + r3)
    z = sqrt(20*ra*sqrt(3 + 3*r3) + 12*r3*ra*sqrt(3 + 3*r3) + 64*r3 + 112)
    assert sqrtdenest(z) == z

def dotest(s):
    x = Symbol("x")
    y = Symbol("y")
    l = [
        Rational(2),
        Float("1.3"),
        x,
        y,
        pow(x, y)*y,
        5,
        5.5
    ]
    for x in l:
        for y in l:
            s(x, y)
    return True

def test_doit():
    from sympy import Symbol
    p = Symbol('p', positive=True)
    n = Symbol('n', negative=True)
    np = Symbol('np', nonpositive=True)
    nn = Symbol('nn', nonnegative=True)

    assert Gt(p, 0).doit() is S.true
    assert Gt(p, 1).doit() == Gt(p, 1)
    assert Ge(p, 0).doit() is S.true
    assert Le(p, 0).doit() is S.false
    assert Lt(n, 0).doit() is S.true
    assert Le(np, 0).doit() is S.true
    assert Gt(nn, 0).doit() == Gt(nn, 0)
    assert Lt(nn, 0).doit() is S.false

    assert Eq(x, 0).doit() == Eq(x, 0)

def canonical_variables(self):
        """Return a dictionary mapping any variable defined in
        ``self.variables`` as underscore-suffixed numbers
        corresponding to their position in ``self.variables``. Enough
        underscores are added to ensure that there will be no clash with
        existing free symbols.

        Examples
        ========

        >>> from sympy import Lambda
        >>> from sympy.abc import x
        >>> Lambda(x, 2*x).canonical_variables
        {x: 0_}
        """
        if not hasattr(self, 'variables'):
            return {}
        u = "_"
        while any(s.name.endswith(u) for s in self.free_symbols):
            u += "_"
        name = '%%i%s' % u
        V = self.variables
        return dict(list(zip(V, [C.Symbol(name % i, **v.assumptions0)
            for i, v in enumerate(V)])))

def _recursive_call(expr_to_call, on_args):
        def the_call_method_is_overridden(expr):
            for cls in getmro(type(expr)):
                if '__call__' in cls.__dict__:
                    return cls != Basic

        if callable(expr_to_call) and the_call_method_is_overridden(expr_to_call):
            if isinstance(expr_to_call, C.Symbol):  # XXX When you call a Symbol it is
                return expr_to_call               # transformed into an UndefFunction
            else:
                return expr_to_call(*on_args)
        elif expr_to_call.args:
            args = [Basic._recursive_call(
                sub, on_args) for sub in expr_to_call.args]
            return type(expr_to_call)(*args)
        else:
            return expr_to_call

def _coeff_isneg(a):
    """Return True if the leading Number is negative.

    Examples
    ========

    >>> from sympy.core.function import _coeff_isneg
    >>> from sympy import S, Symbol, oo, pi
    >>> _coeff_isneg(-3*pi)
    True
    >>> _coeff_isneg(S(3))
    False
    >>> _coeff_isneg(-oo)
    True
    >>> _coeff_isneg(Symbol('n', negative=True)) # coeff is 1
    False

    """

    if a.is_Mul:
        a = a.args[0]
    return a.is_Number and a.is_negative

def _diff_wrt(self):
        """Allow derivatives wrt Derivatives if it contains a function.

        Examples
        ========

            >>> from sympy import Function, Symbol, Derivative
            >>> f = Function('f')
            >>> x = Symbol('x')
            >>> Derivative(f(x),x)._diff_wrt
            True
            >>> Derivative(x**2,x)._diff_wrt
            False
        """
        if self.expr.is_Function:
            return True
        else:
            return False

def __new__(cls, variables, expr):
        from sympy.sets.sets import FiniteSet
        try:
            for v in variables if iterable(variables) else [variables]:
                if not v.is_Symbol:
                    raise TypeError("v is not a symbol")
        except (AssertionError, AttributeError):
            raise ValueError('variable is not a Symbol: %s' % v)
        try:
            variables = Tuple(*variables)
        except TypeError:
            variables = Tuple(variables)
        if len(variables) == 1 and variables[0] == expr:
            return S.IdentityFunction

        obj = Expr.__new__(cls, Tuple(*variables), S(expr))
        obj.nargs = FiniteSet(len(variables))
        return obj

def test_unify_iter():
    expr = Add(1, 2, 3, evaluate=False)
    a, b, c = map(Symbol, 'abc')
    pattern = Add(a, c, evaluate=False)
    assert is_associative(deconstruct(pattern))
    assert is_commutative(deconstruct(pattern))

    result   = list(unify(expr, pattern, {}, (a, c)))
    expected = [{a: 1, c: Add(2, 3, evaluate=False)},
                {a: 1, c: Add(3, 2, evaluate=False)},
                {a: 2, c: Add(1, 3, evaluate=False)},
                {a: 2, c: Add(3, 1, evaluate=False)},
                {a: 3, c: Add(1, 2, evaluate=False)},
                {a: 3, c: Add(2, 1, evaluate=False)},
                {a: Add(1, 2, evaluate=False), c: 3},
                {a: Add(2, 1, evaluate=False), c: 3},
                {a: Add(1, 3, evaluate=False), c: 2},
                {a: Add(3, 1, evaluate=False), c: 2},
                {a: Add(2, 3, evaluate=False), c: 1},
                {a: Add(3, 2, evaluate=False), c: 1}]

    assert iterdicteq(result, expected)

def test_field_assumptions():
    X = MatrixSymbol('X', 4, 4)
    Y = MatrixSymbol('Y', 4, 4)
    assert ask(Q.real_elements(X), Q.real_elements(X))
    assert not ask(Q.integer_elements(X), Q.real_elements(X))
    assert ask(Q.complex_elements(X), Q.real_elements(X))
    assert ask(Q.real_elements(X+Y), Q.real_elements(X)) is None
    assert ask(Q.real_elements(X+Y), Q.real_elements(X) & Q.real_elements(Y))
    from sympy.matrices.expressions.hadamard import HadamardProduct
    assert ask(Q.real_elements(HadamardProduct(X, Y)),
                    Q.real_elements(X) & Q.real_elements(Y))
    assert ask(Q.complex_elements(X+Y), Q.real_elements(X) & Q.complex_elements(Y))

    assert ask(Q.real_elements(X.T), Q.real_elements(X))
    assert ask(Q.real_elements(X.I), Q.real_elements(X) & Q.invertible(X))
    assert ask(Q.real_elements(Trace(X)), Q.real_elements(X))
    assert ask(Q.integer_elements(Determinant(X)), Q.integer_elements(X))
    assert not ask(Q.integer_elements(X.I), Q.integer_elements(X))
    alpha = Symbol('alpha')
    assert ask(Q.real_elements(alpha*X), Q.real_elements(X) & Q.real(alpha))
    assert ask(Q.real_elements(LofLU(X)), Q.real_elements(X))

def _pop_symbol_list(self, lists):
        symbol_lists = []
        for l in lists:
            mark = True
            for s in l:
                if s is not None and not isinstance(s, Symbol):
                    mark = False
                    break
            if mark:
                lists.remove(l)
                symbol_lists.append(l)
        if len(symbol_lists) == 1:
            return symbol_lists[0]
        elif len(symbol_lists) == 0:
            return []
        else:
            raise ValueError("Only one list of Symbols "
                             "can be given for a color scheme.")

def equation(self, type='cosine'):
        """
        Returns equation of the wave.

        Examples
        ========

        >>> from sympy import symbols
        >>> from sympy.physics.optics import TWave
        >>> A, phi, f = symbols('A, phi, f')
        >>> w = TWave(A, f, phi)
        >>> w.equation('cosine')
        A*cos(2*pi*f*t + phi)
        """
        if not isinstance(type, str):
            raise TypeError("type can only be a string.")
        if type == 'cosine':
            return self._amplitude*cos(self.angular_velocity*Symbol('t') + self._phase)

def test_twave():
    A1, phi1, A2, phi2, f = symbols('A1, phi1, A2, phi2, f')
    c = Symbol('c')  # Speed of light in vacuum
    n = Symbol('n')  # Refractive index
    t = Symbol('t')  # Time
    w1 = TWave(A1, f, phi1)
    w2 = TWave(A2, f, phi2)
    assert w1.amplitude == A1
    assert w1.frequency == f
    assert w1.phase == phi1
    assert w1.wavelength == c/(f*n)
    assert w1.time_period == 1/f
    w3 = w1 + w2
    assert w3.amplitude == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)
    assert w3.frequency == f
    assert w3.wavelength == c/(f*n)
    assert w3.time_period == 1/f
    assert w3.angular_velocity == 2*pi*f
    assert w3.equation() == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)*cos(2*pi*f*t + phi1 + phi2)

def test_dagger():
    x = symbols('x')
    expr = Dagger(x)
    assert str(expr) == 'Dagger(x)'
    ascii_str = \
"""\
 +\n\
x \
"""
    ucode_str = \
u("""\
 †\n\
x \
""")
    assert pretty(expr) == ascii_str
    assert upretty(expr) == ucode_str
    assert latex(expr) == r'x^{\dag}'
    sT(expr, "Dagger(Symbol('x'))")

def test_tensor_product():
    n = Symbol('n')
    hs1 = ComplexSpace(2)
    hs2 = ComplexSpace(n)

    h = hs1*hs2
    assert isinstance(h, TensorProductHilbertSpace)
    assert h.dimension == 2*n
    assert h.spaces == (hs1, hs2)

    h = hs2*hs2
    assert isinstance(h, TensorPowerHilbertSpace)
    assert h.base == hs2
    assert h.exp == 2
    assert h.dimension == n**2

    f = FockSpace()
    h = hs1*hs2*f
    assert h.dimension == oo

def test_tensor_power():
    n = Symbol('n')
    hs1 = ComplexSpace(2)
    hs2 = ComplexSpace(n)

    h = hs1**2
    assert isinstance(h, TensorPowerHilbertSpace)
    assert h.base == hs1
    assert h.exp == 2
    assert h.dimension == 4

    h = hs2**3
    assert isinstance(h, TensorPowerHilbertSpace)
    assert h.base == hs2
    assert h.exp == 3
    assert h.dimension == n**3

def test_units():
    assert (5*m/s * day) / km == 432
    assert foot / meter == Rational('0.3048')
    # amu is a pure mass so mass/mass gives a number, not an amount (mol)
    assert str(grams/(amu).n(2)) == '6.0e+23'

    # Light from the sun needs about 8.3 minutes to reach earth
    t = (1*au / speed_of_light).evalf() / minute
    assert abs(t - 8.31) < 0.1

    assert sqrt(m**2) == m
    assert (sqrt(m))**2 == m

    t = Symbol('t')
    assert integrate(t*m/s, (t, 1*s, 5*s)) == 12*m*s
    assert (t * m/s).integrate((t, 1*s, 5*s)) == 12*m*s

def _mat_inv_mul(A, B):
    """
    Computes A^-1 * B symbolically w/ substitution, where B is not
    necessarily a vector, but can be a matrix.

    """

    r1, c1 = A.shape
    r2, c2 = B.shape
    temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j) + str(r1 * i)))
    temp2 = Matrix(r2, c2, lambda i, j: Symbol('y' + str(j) + str(r2 * i)))
    for i in range(len(temp1)):
        if A[i] == 0:
            temp1[i] = 0
    for i in range(len(temp2)):
        if B[i] == 0:
            temp2[i] = 0
    temp3 = []
    for i in range(c2):
        temp3.append(temp1.LDLsolve(temp2[:, i]))
    temp3 = Matrix([i.T for i in temp3]).T
    return temp3.subs(dict(list(zip(temp1, A)))).subs(dict(list(zip(temp2, B))))

def _print_MatrixElement(self, expr):
        from sympy.matrices import MatrixSymbol
        from sympy import Symbol
        if (isinstance(expr.parent, MatrixSymbol)
                and expr.i.is_number and expr.j.is_number):
            return self._print(
                    Symbol(expr.parent.name + '_%d%d'%(expr.i, expr.j)))
        else:
            prettyFunc = self._print(expr.parent)
            prettyIndices = self._print_seq((expr.i, expr.j), delimiter=', '
                    ).parens(left='[', right=']')[0]
            pform = prettyForm(binding=prettyForm.FUNC,
                    *stringPict.next(prettyFunc, prettyIndices))

            # store pform parts so it can be reassembled e.g. when powered
            pform.prettyFunc = prettyFunc
            pform.prettyArgs = prettyIndices

            return pform

def _print_Function(self, e, sort=False):
        # XXX works only for applied functions
        func = e.func
        args = e.args
        if sort:
            args = sorted(args, key=default_sort_key)

        func_name = func.__name__

        prettyFunc = self._print(C.Symbol(func_name))
        prettyArgs = prettyForm(*self._print_seq(args).parens())
        #postioning func_name
        mid = prettyArgs.height()//2
        if mid > 2:
            prettyFunc.baseline = -mid + 1

        pform = prettyForm(
            binding=prettyForm.FUNC, *stringPict.next(prettyFunc, prettyArgs))

        # store pform parts so it can be reassembled e.g. when powered
        pform.prettyFunc = prettyFunc
        pform.prettyArgs = prettyArgs

        return pform

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 sympy.ask(~sympy.Q.zero(sympy.log(base))):
            return rule
        elif sympy.ask(sympy.Q.zero(sympy.log(base))):
            return ConstantRule(1, 1, symbol)

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

def test_numerically(f, g, z=None, tol=1.0e-6, a=2, b=-1, c=3, d=1):
    """
    Test numerically that f and g agree when evaluated in the argument z.

    If z is None, all symbols will be tested. This routine does not test
    whether there are Floats present with precision higher than 15 digits
    so if there are, your results may not be what you expect due to round-
    off errors.

    Examples
    ========

    >>> from sympy import sin, cos
    >>> from sympy.abc import x
    >>> from sympy.utilities.randtest import test_numerically as tn
    >>> tn(sin(x)**2 + cos(x)**2, 1, x)
    True
    """
    f, g, z = Tuple(f, g, z)
    z = [z] if isinstance(z, Symbol) else (f.free_symbols | g.free_symbols)
    reps = list(zip(z, [random_complex_number(a, b, c, d) for zi in z]))
    z1 = f.subs(reps).n()
    z2 = g.subs(reps).n()
    return comp(z1, z2, tol)

def rationalize(tokens, local_dict, global_dict):
    """Converts floats into ``Rational``. Run AFTER ``auto_number``."""
    result = []
    passed_float = False
    for toknum, tokval in tokens:
        if toknum == NAME:
            if tokval == 'Float':
                passed_float = True
                tokval = 'Rational'
            result.append((toknum, tokval))
        elif passed_float == True and toknum == NUMBER:
            passed_float = False
            result.append((STRING, tokval))
        else:
            result.append((toknum, tokval))

    return result


#: Standard transformations for :func:`parse_expr`.
#: Inserts calls to :class:`Symbol`, :class:`Integer`, and other SymPy
#: datatypes and allows the use of standard factorial notation (e.g. ``x!``).

def test_function_exponentiation():
    cases = {
        'sin**2(x)': 'sin(x)**2',
        'exp^y(z)': 'exp(z)^y',
        'sin**2(E^(x))': 'sin(E^(x))**2'
    }
    transformations = standard_transformations + (convert_xor,)
    transformations2 = transformations + (function_exponentiation,)
    for case in cases:
        implicit = parse_expr(case, transformations=transformations2)
        normal = parse_expr(cases[case], transformations=transformations)
        assert(implicit == normal)

    other_implicit = ['x y', 'x sin x', '2x', 'sin x',
                      'cos 2*x', 'sin cos x']
    for case in other_implicit:
        raises(SyntaxError,
               lambda: parse_expr(case, transformations=transformations2))

    assert parse_expr('x**2', local_dict={ 'x': sympy.Symbol('x') },
                      transformations=transformations2) == parse_expr('x**2')

def test_maxima_functions():
    assert parse_maxima('expand( (x+1)^2)') == x**2 + 2*x + 1
    assert parse_maxima('factor( x**2 + 2*x + 1)') == (x + 1)**2
    assert parse_maxima('2*cos(x)^2 + sin(x)^2') == 2*cos(x)**2 + sin(x)**2
    assert parse_maxima('trigexpand(sin(2*x)+cos(2*x))') == \
        -1 + 2*cos(x)**2 + 2*cos(x)*sin(x)
    assert parse_maxima('solve(x^2-4,x)') == [-2, 2]
    assert parse_maxima('limit((1+1/x)^x,x,inf)') == E
    assert parse_maxima('limit(sqrt(-x)/x,x,0,minus)') == -oo
    assert parse_maxima('diff(x^x, x)') == x**x*(1 + log(x))
    assert parse_maxima('sum(k, k, 1, n)', name_dict=dict(
        n=Symbol('n', integer=True),
        k=Symbol('k', integer=True)
    )) == (n**2 + n)/2
    assert parse_maxima('product(k, k, 1, n)', name_dict=dict(
        n=Symbol('n', integer=True),
        k=Symbol('k', integer=True)
    )) == factorial(n)
    assert parse_maxima('ratsimp((x^2-1)/(x+1))') == x - 1
    assert Abs( parse_maxima(
        'float(sec(%pi/3) + csc(%pi/3))') - 3.154700538379252) < 10**(-5)

def main():
    a = sympy.Symbol('a')
    b = sympy.Symbol('b')
    e = (a + 2*b)**5

    print("\nExpression : ")
    print()
    pprint(e)
    print("\n\nDifferentiating w.r.t. a:")
    print()
    pprint(e.diff(a))
    print("\n\nDifferentiating w.r.t. b:")
    print()
    pprint(e.diff(b))
    print("\n\nSecond derivative of the above result w.r.t. a:")
    print()
    pprint(e.diff(b).diff(a, 2))
    print("\n\nExpanding the above result:")
    print()
    pprint(e.expand().diff(b).diff(a, 2))
    print()

def main():
    x = Symbol("x")
    a = Symbol("a")
    h = Symbol("h")

    show( limit(sqrt(x**2 - 5*x + 6) - x, x, oo), -Rational(5)/2 )

    show( limit(x*(sqrt(x**2 + 1) - x), x, oo), Rational(1)/2 )

    show( limit(x - sqrt3(x**3 - 1), x, oo), Rational(0) )

    show( limit(log(1 + exp(x))/x, x, -oo), Rational(0) )

    show( limit(log(1 + exp(x))/x, x, oo), Rational(1) )

    show( limit(sin(3*x)/x, x, 0), Rational(3) )

    show( limit(sin(5*x)/sin(2*x), x, 0), Rational(5)/2 )

    show( limit(((x - 1)/(x + 1))**x, x, oo), exp(-2))

def main():
    x = sympy.Symbol('x')
    y = sympy.Symbol('y')

    e = 1/sympy.cos(x)
    print()
    pprint(e)
    print('\n')
    pprint(e.subs(sympy.cos(x), y))
    print('\n')
    pprint(e.subs(sympy.cos(x), y).subs(y, x**2))

    e = 1/sympy.log(x)
    e = e.subs(x, sympy.Float("2.71828"))
    print('\n')
    pprint(e)
    print('\n')
    pprint(e.evalf())
    print()

def state_update_equation(self, state_update_equation):
        if state_update_equation is None:
            if self.dim_state > 0:
                state_update_equation = self.state
            else:
                state_update_equation = self.input
        assert state_update_equation.atoms(sp.Symbol) <= set(
            self.constants_values.keys()) | set([dynamicsymbols._t])
        self._state_update_equation = state_update_equation
        if self.dim_state:
            assert find_dynamicsymbols(state_update_equation) <= \
                set(self.state)
            self.state_update_equation_function = self.code_generator(
                    [dynamicsymbols._t] + sp.flatten(self.state),
                    self._state_update_equation.subs(self.constants_values),
                    **self.code_generator_args
            )
        else:
            assert find_dynamicsymbols(state_update_equation) <= \
                set(self.input)
            self.state_update_equation_function = self.code_generator(
                    [dynamicsymbols._t] + sp.flatten(self.input),
                    self._state_update_equation.subs(self.constants_values),
                    **self.code_generator_args
            )

def event_variable_equation(self, event_variable_equation):
        assert event_variable_equation.atoms(sp.Symbol) <= set(
            self.constants_values.keys()) | set([dynamicsymbols._t])
        self._event_variable_equation = event_variable_equation
        if self.dim_state:
            assert find_dynamicsymbols(event_variable_equation) <= \
                set(self.state)
            self.event_variable_equation_function = self.code_generator(
                [dynamicsymbols._t] + sp.flatten(self.state),
                self._event_variable_equation.subs(self.constants_values),
                **self.code_generator_args
            )
        else:
            assert find_dynamicsymbols(event_variable_equation) <= \
                set(self.input)
            self.event_variable_equation_function = self.code_generator(
                [dynamicsymbols._t] + sp.flatten(self.input),
                self._event_variable_equation.subs(self.constants_values),
                **self.code_generator_args
            )

def state_equation(self, state_equation):
        if state_equation is None:  # or other checks?
            state_equation = empty_array()
        else:
            assert len(state_equation) == len(self.state)
            assert find_dynamicsymbols(state_equation) <= (
                    set(self.state) | set(self.input)
                )
            assert state_equation.atoms(sp.Symbol) <= (
                    set(self.constants_values.keys())
                    | set([dynamicsymbols._t])
                )

        self._state_equation = state_equation
        self.update_state_equation_function()

        self.state_jacobian_equation = grad(self.state_equation, self.state)
        self.update_state_jacobian_function()

        self.input_jacobian_equation = grad(self.state_equation, self.input)
        self.update_input_jacobian_function()

def output_equation(self, output_equation):
        if isinstance(output_equation, sp.Expr):
            output_equation = Array([output_equation])
        if output_equation is None and self.dim_state == 0:
            output_equation = empty_array()
        else:
            if output_equation is None:
                output_equation = self.state

            assert output_equation.atoms(sp.Symbol) <= (
                    set(self.constants_values.keys())
                    | set([dynamicsymbols._t])
                   )

            if self.dim_state:
                assert find_dynamicsymbols(output_equation) <= set(self.state)
            else:
                assert find_dynamicsymbols(output_equation) <= set(self.input)

        self.dim_output = len(output_equation)

        self._output_equation = output_equation
        self.update_output_equation_function()

def _print_Symbol(self, sym):
        """Print a Symbol object.

        :param sym: The expression.
        :rtype : bce.dom.mathml.all.Base
        :return: The printed MathML object.
        """

        assert isinstance(sym, _sympy.Symbol)

        if self.__ph_enabled and sym.name.startswith(self.__ph_prefix):
            return _mathml.SubComponent(
                _mathml.TextComponent(self.__ph_prefix.lower()),
                _mathml.TextComponent(sym.name[len(self.__ph_prefix):])
            )
        else:
            return _mathml.TextComponent(sym.name)

def free_parameters(expr):
    """
    Returns the free parameters of a given expression. In general, this
    corresponds to length of a For loop.

    expr: sympy.Expr
        any sympy expression or pyccel.ast.core object
    """
    args = []
    if isinstance(expr, For):
        b = expr.iterable.args[0]
        e = expr.iterable.args[1]
        if isinstance(b, Symbol):
            args.append(str(b))
        if isinstance(e, Symbol):
            args.append(str(e))
        for i in expr.body:
            args += free_parameters(i)
        args = set(args)
        args = list(args)

    return args
# ...

# ...

def do_arg(a):
    if isinstance(a, str):
        arg = Symbol(a, integer=True)
    elif isinstance(a, (Integer, Float)):
        arg = a
    elif isinstance(a, ArithmeticExpression):
        arg = a.expr
        if isinstance(arg, (Symbol, Variable)):
            arg = Symbol(arg.name, integer=True)
        else:
            arg = convert_to_integer_expression(arg)
    else:
        raise Exception('Wrong instance in do_arg')

    return arg
# ...

# ... TODO improve. this version is not working with function calls

def expr(self):
        args = sp_symbols(self.args)

        # ... we update the namespace
        ls = args
        if isinstance(args, Symbol):
            ls = [args]

        for i in ls:
            namespace[str(i)] = i
        # ...

        # ... we do it here after the namespace has been updated
        e = self.rhs.expr
        # ...

        # ... we clean the namespace
        for i in ls:
            namespace.pop(str(i))
        # ...

        return Lambda(args, e)

def expr(self):
        """
        Process a Trailer by returning the approriate objects from
        pyccel.ast.core
        """
        self.update()
        args = []
        for a in self.args:
            if isinstance(a, ArithmeticExpression):
                arg = do_arg(a)

                # TODO treat n correctly
                n = Symbol('n', integer=True)
                i = Idx(arg, n)
                args.append(i)
            elif isinstance(a, BasicSlice):
                arg = a.expr
                args.append(arg)
            else:
                raise Exception('Wrong instance')
        return args

def __new__(cls, lhs, rhs, strict=False, status=None, like=None):
        cls._strict = strict
        if strict:
            lhs = sympify(lhs)
            rhs = sympify(rhs)
            # Tuple of things that can be on the lhs of an assignment
            assignable = (Symbol, MatrixSymbol, MatrixElement, Indexed, Idx)
            #if not isinstance(lhs, assignable):
            #    raise TypeError("Cannot assign to lhs of type %s." % type(lhs))
            # Indexed types implement shape, but don't define it until later. This
            # causes issues in assignment validation. For now, matrices are defined
            # as anything with a shape that is not an Indexed
            lhs_is_mat = hasattr(lhs, 'shape') and not isinstance(lhs, Indexed)
            rhs_is_mat = hasattr(rhs, 'shape') and not isinstance(rhs, Indexed)
            # If lhs and rhs have same structure, then this assignment is ok
            if lhs_is_mat:
                if not rhs_is_mat:
                    raise ValueError("Cannot assign a scalar to a matrix.")
                elif lhs.shape != rhs.shape:
                    raise ValueError("Dimensions of lhs and rhs don't align.")
            elif rhs_is_mat and not lhs_is_mat:
                raise ValueError("Cannot assign a matrix to a scalar.")
        return Basic.__new__(cls, lhs, rhs, status, like)

def wiener_attack(e, n):
    pq = partial_quotiens(e, n)
    c = convergents(pq)
    x = Symbol('x')
    print 'done'
    for (k, d) in c:
        if d != 196176397553781094149364022413702698618416570027142099596910222792686220199799634430113040703408219640551905602804050764120518037301734891772260502596755642106767631991390703026953707346708504142195119698789998222566587930259938572562702895638995035040020739921966671123993267592460124879940448923:
            continue
        if k != 0:
            y = n - phiN(e, d, k) + 1
            roots = solve(x**2 - y*x + n, x)
            if len(roots) == 2:
                p = roots[0]
                q = roots[1]
                if p * q == n:
                    break
    return p, q

def __init__(self, n, e):
        self.d = None
        self.p = None
        self.q = None
        sys.setrecursionlimit(100000)
        frac = self.rational_to_contfrac(e, n)
        convergents = self.convergents_from_contfrac(frac)

        for (k,d) in convergents:
            if k!=0 and (e*d-1)%k == 0:
                phi = (e*d-1)//k
                s = n - phi + 1
                discr = s*s - 4*n
                if(discr>=0):
                    t = self.is_perfect_square(discr)
                    if t!=-1 and (s+t)%2==0:
                        self.d = d
                        x = Symbol('x')
                        roots = solve(x**2 - s*x + n, x)
                        if len(roots) == 2:
                            self.p = roots[0]
                            self.q = roots[1]
                        break

def __init__(self, n, e):
        self.p = None
        self.q = None
        pq = self.partial_quotiens(e, n)
        c = self.convergents(pq)
        x = Symbol('x')
        for (k, d) in c:
            if k != 0:
                y = n - self.phiN(e, d, k) + 1
                roots = solve(x**2 - y*x + n, x)
                if len(roots) == 2:
                    p = roots[0]
                    q = roots[1]
                    if p * q == n:
                        self.p = p
                        self.q = q
                        break

def dof_unknowns(control_volume_dimensions):
    unknown_list = []

    for path in control_volume_dimensions:
        for equation_key, equation_value in path.items():
            for variable in sp.sympify(equation_value).atoms(sp.Symbol):
                if variable not in unknown_list:
                    unknown_list.append(variable)

    unknowns = len(unknown_list)

    return unknowns


# This function uses the dof_unknowns, along with the setup matrix from the first function
# to determine how many degrees of freedom exist in a given control volume

def get_S95(b0, sigma):

    S95 = sy.Symbol('S95', positive = True, real = True)
    b = sy.Symbol('b', positive = True)
    chi21 = sy.Symbol('chi21')
    chi22 = sy.Symbol('chi22')

    chi2 = 3.84
    N = b0 

    replacements = [(b, (b0 - S95 - sigma**2)/2 + 1./2*((b0 - S95 - sigma**2)**2 + 4*(sigma**2*N - S95*sigma**2 + S95*b0))**0.5)]

    replacements2 = [(S95, 0.)]

    chi21 = -2*( N*sy.log(S95 + b) - (S95 + b) - ((b-b0)/sigma)**2)

    chi21 = chi21.subs(replacements)
    chi22 = chi21.subs(replacements2)

    eq = chi2 - chi21 + chi22

    S95_new = sy.nsolve(eq, S95, 1)

    return float(S95_new)