diff --git a/symengine/lib/symengine_wrapper.in.pyx b/symengine/lib/symengine_wrapper.in.pyx
index b37892a5..0db21878 100644
--- a/symengine/lib/symengine_wrapper.in.pyx
+++ b/symengine/lib/symengine_wrapper.in.pyx
@@ -1218,6 +1218,10 @@ cdef class Basic(object):
     def is_Relational(self):
         return False
 
+    @property
+    def is_Piecewise(self):
+        return False
+
     @property
     def is_Equality(self):
         return False
@@ -1610,6 +1614,10 @@ eulergamma = EulerGamma()
 
 cdef class Boolean(Expr):
 
+    @property
+    def func(self):
+        return self.__class__
+
     def logical_not(self):
         return c2py(<rcp_const_basic>(deref(symengine.rcp_static_cast_Boolean(self.thisptr)).logical_not()))
 
@@ -1710,11 +1718,27 @@ class Xor(Boolean):
 
 
 class Relational(Boolean):
+    """A relation ``lhs rel_op rhs``.
+
+    ``lhs`` and ``rhs`` are the sides as stored, which can differ from the
+    order they were given in: the sides of an ``Equality`` or ``Unequality``
+    are kept in a canonical order, and there are no ``GreaterThan`` and
+    ``StrictGreaterThan`` classes, so ``Ge(a, b)`` is ``LessThan(b, a)`` and
+    ``Gt(a, b)`` is ``StrictLessThan(b, a)``, whose ``lhs`` is ``b``.
+    """
 
     @property
     def is_Relational(self):
         return True
 
+    @property
+    def lhs(self):
+        return self.args[0]
+
+    @property
+    def rhs(self):
+        return self.args[1]
+
     def __bool__(self):
         raise TypeError("cannot determine truth value of Relational")
 
@@ -1724,6 +1748,8 @@ Rel = Relational
 
 class Equality(Relational):
 
+    rel_op = "=="
+
     def __new__(cls, *args):
         return eq(*args)
 
@@ -1741,16 +1767,14 @@ class Equality(Relational):
     def is_Equality(self):
         return True
 
-    @property
-    def func(self):
-        return self.__class__
-
 
 Eq = Equality
 
 
 class Unequality(Relational):
 
+    rel_op = "!="
+
     def __new__(cls, *args):
         return ne(*args)
 
@@ -1764,16 +1788,14 @@ class Unequality(Relational):
         s = self.args_as_sage()
         return sage.ne(*s)
 
-    @property
-    def func(self):
-        return self.__class__
-
 
 Ne = Unequality
 
 
 class LessThan(Relational):
 
+    rel_op = "<="
+
     def __new__(cls, *args):
         return le(*args)
 
@@ -1793,6 +1815,8 @@ Le = LessThan
 
 class StrictLessThan(Relational):
 
+    rel_op = "<"
+
     def __new__(cls, *args):
         return lt(*args)
 
@@ -1859,6 +1883,39 @@ cdef class Number(Expr):
     def imag(self):
         return S.Zero
 
+    def round(self, n=None):
+        """Round to the nearest Integer, ties to even, or, if ``n`` is given,
+        to ``n`` decimal places, keeping the kind of number (an exact Rational
+        stays exact). Like SymPy's ``Number.round``.
+
+        A complex number is rounded in its real and imaginary parts; infinities
+        and nan are returned unchanged.
+        """
+        if isinstance(self, (Infinity, NegativeInfinity, ComplexInfinity, NaN)):
+            return self
+        if self.is_complex:
+            return self.real_part().round(n) + I * self.imaginary_part().round(n)
+        if n is not None:
+            scale = Integer(10) ** int(n)
+            result = (self * scale).round() / scale
+            if isinstance(self, Float):
+                bits = int(self.get_prec()) if isinstance(self, RealMPFR) else 53
+                return evalf(result, bits, real=True)
+            return result
+        floored = floor(self)
+        excess = self - floored
+        half = Rational(1, 2)
+        if excess > half:
+            return floored + 1
+        if excess < half:
+            return floored
+        # A tie (compared by value: a Float 0.5 is not equal to 1/2 as an
+        # expression) goes to the even neighbour.
+        return floored + floored % 2
+
+    def __round__(self, ndigits=None):
+        return self.round(ndigits)
+
 
 class Rational(Number):
 
@@ -2545,9 +2602,9 @@ class Function(Expr):
     def is_Function(self):
         return True
 
-    def func(self, *values):
-        import sys
-        return getattr(sys.modules[__name__], self.__class__.__name__.lower())(*values)
+    @property
+    def func(self):
+        return self.__class__
 
 
 class OneArgFunction(Function):
@@ -2973,18 +3030,32 @@ class FunctionSymbol(Function):
         name = self.get_name()
         return sage.function(name)(*self.args_as_sage())
 
-    def func(self, *values):
-        name = self.get_name()
-        return function_symbol(name, *values)
+    @property
+    def func(self):
+        # Like SymPy's, the undefined function itself: f(x).func is f.
+        return UndefFunction(self.get_name())
 
 
 class UndefFunction(object):
     def __init__(self, name):
         self.name = name
 
+    @property
+    def __name__(self):
+        return self.name
+
     def __call__(self, *values):
         return function_symbol(self.name, *values)
 
+    def __eq__(self, other):
+        return isinstance(other, UndefFunction) and self.name == other.name
+
+    def __hash__(self):
+        return hash((UndefFunction, self.name))
+
+    def __repr__(self):
+        return self.name
+
 
 cdef rcp_const_basic pynumber_to_symengine(PyObject* o1):
     cdef Basic X = sympify(<object>o1)
@@ -3079,9 +3150,42 @@ class PyFunction(FunctionSymbol):
         import sage.all as sage
         return sage.SR(self.pyobject())
 
+    @property
+    def func(self):
+        """Rebuilds the wrapped object with new arguments through its own
+        ``func``, so the result is again a PyFunction of the same class."""
+        pyfunc = getattr(self.pyobject(), "func", None)
+        if pyfunc is None:
+            return UndefFunction(self.get_name())
+        return _PyFunctionFunc(pyfunc)
+
     def pyobject(Basic self):
         return <object>deref(symengine.rcp_static_cast_PyFunction(self.thisptr)).get_py_object()
 
+class _PyFunctionFunc(object):
+    """``func`` of a PyFunction: its wrapped (SymPy) class, applied to
+    SymEngine arguments and converted back."""
+
+    def __init__(self, pyfunc):
+        self.pyfunc = pyfunc
+
+    @property
+    def __name__(self):
+        return self.pyfunc.__name__
+
+    def __call__(self, *values):
+        return sympify(self.pyfunc(*[_sympify(v)._sympy_() for v in values]))
+
+    def __eq__(self, other):
+        return isinstance(other, _PyFunctionFunc) and self.pyfunc == other.pyfunc
+
+    def __hash__(self):
+        return hash((_PyFunctionFunc, self.pyfunc))
+
+    def __repr__(self):
+        return repr(self.pyfunc)
+
+
 cdef class PyFunctionClass(object):
 
     def __cinit__(self, function, PyModule module not None):
@@ -3231,9 +3335,24 @@ class Subs(Expr):
 
 class Piecewise(Function):
 
-    def __new__(self, *args):
+    def __new__(cls, *args):
+        # Also accept the flattened form of ``args`` (expr1, cond1, expr2,
+        # cond2, ...), so that ``p.func(*p.args) == p``.
+        if args and not isinstance(args[0], collections.abc.Iterable):
+            if len(args) % 2:
+                raise TypeError("Piecewise takes (expr, cond) pairs or an even "
+                                "number of arguments alternating expr and cond")
+            args = zip(args[0::2], args[1::2])
         return piecewise(*args)
 
+    @property
+    def func(self):
+        return self.__class__
+
+    @property
+    def is_Piecewise(self):
+        return True
+
     def _sympy_(self):
         import sympy
         a = self.args
diff --git a/symengine/tests/test_functions.py b/symengine/tests/test_functions.py
index 7c4573df..906387ec 100644
--- a/symengine/tests/test_functions.py
+++ b/symengine/tests/test_functions.py
@@ -54,6 +54,22 @@ def test_f():
     assert f == g
 
 
+def test_func():
+    x = Symbol("x")
+    y = Symbol("y")
+    for e in (sin(x), log(x), log(x, 2), floor(x), ceiling(x), LambertW(x),
+              KroneckerDelta(x, y), zeta(x, y), exp(x), sign(x)):
+        assert e.func is type(e)
+        assert e.func(*e.args) == e
+    # An undefined function's func is the function itself, as in SymPy
+    f = function_symbol("f", x, y)
+    assert f.func(*f.args) == f
+    assert f.func == function_symbol("f", y).func
+    assert f.func != function_symbol("g", x).func
+    assert f.func.__name__ == "f"
+    assert f.func(y) == function_symbol("f", y)
+
+
 def test_derivative():
     x = Symbol("x")
     y = Symbol("y")
diff --git a/symengine/tests/test_logic.py b/symengine/tests/test_logic.py
index 3d0c3391..06ce0915 100644
--- a/symengine/tests/test_logic.py
+++ b/symengine/tests/test_logic.py
@@ -33,6 +33,25 @@ def test_relationals():
     assert ne.func(*ne.args) == ne
 
 
+def test_relational_attributes():
+    for rel, rel_op in ((Eq(x, y), "=="), (Ne(x, y), "!="), (Le(x, y), "<="),
+                        (Lt(x, y), "<")):
+        assert rel.rel_op == rel_op
+        assert (rel.lhs, rel.rhs) == rel.args
+        assert rel.func(*rel.args) == rel
+    assert (Le(x, 1).lhs, Le(x, 1).rhs) == (x, 1)
+    # Ge and Gt are stored as LessThan and StrictLessThan with the sides swapped
+    assert (Ge(x, 1).rel_op, Ge(x, 1).lhs, Ge(x, 1).rhs) == ("<=", 1, x)
+    assert (Gt(x, 1).rel_op, Gt(x, 1).lhs, Gt(x, 1).rhs) == ("<", 1, x)
+    assert Ge(x, 1).func(*Ge(x, 1).args) == Ge(x, 1)
+
+
+def test_boolean_func():
+    for b in (And(x < 1, y < 1), Or(x < 1, y < 1), Not(Contains(x, Interval(0, 1))),
+              Xor(x < 1, y < 1), Contains(x, Interval(0, 1)), true, false):
+        assert b.func(*b.args) == b
+
+
 def test_rich_cmp():
     assert (x < y) == Lt(x, y)
     assert (x <= y) == Le(x, y)
@@ -117,6 +136,18 @@ def test_Piecewise():
     assert p.diff(x) == q
 
 
+def test_Piecewise_func():
+    p = Piecewise((x, x < 1), (y, True))
+    assert p.is_Piecewise
+    assert not x.is_Piecewise
+    assert p.func is Piecewise
+    assert p.args == (x, x < 1, y, true)
+    assert p.func(*p.args) == p
+    assert Piecewise(x, x < 1, y, True) == p
+    assert Piecewise(1, x < 1, 0, True) == Piecewise((1, x < 1), (0, True))
+    raises(TypeError, lambda: Piecewise(x, x < 1, y))
+
+
 def test_Contains():
     assert Contains(x, FiniteSet(0)) != false
     assert Contains(x, Interval(1, 1)) != false
diff --git a/symengine/tests/test_number.py b/symengine/tests/test_number.py
index f61189a4..fa2f5891 100644
--- a/symengine/tests/test_number.py
+++ b/symengine/tests/test_number.py
@@ -1,6 +1,8 @@
+from fractions import Fraction
+
 from symengine.test_utilities import raises
 
-from symengine import Integer, I, S, Symbol, pi, Rational
+from symengine import Integer, I, S, Symbol, pi, Rational, RealDouble, oo, nan
 from symengine.lib.symengine_wrapper import (perfect_power, is_square, integer_nthroot)
 
 
@@ -185,3 +187,22 @@ def test_is_nonnegative():
 def test_is_real():
     assert Rational(1, 2).is_real
     assert Symbol('x').is_real is None
+
+
+def test_round():
+    # To the nearest Integer, ties to even, like round() and SymPy
+    for v in (0.5, 1.5, 2.5, -0.5, -1.5, -2.5, 2.4999, 2.6, -2.6, 7.0):
+        assert RealDouble(v).round() == round(v)
+        assert RealDouble(v).round().is_Integer
+    for p, q in ((5, 2), (7, 2), (-5, 2), (-7, 2), (1, 3), (-2, 3), (4, 1)):
+        assert Rational(p, q).round() == round(Fraction(p, q))
+        assert Rational(p, q).round().is_Integer
+    assert round(Rational(5, 2)) == 2
+    # To n decimal places, keeping the kind of number
+    assert Rational(7, 3).round(2) == Rational(233, 100)
+    assert RealDouble(2.345).round(1) == RealDouble(2.3)
+    assert RealDouble(2.345).round(1).is_Float
+    assert Integer(1250).round(-2) == 1200
+    assert (Rational(5, 2) + I * Rational(7, 2)).round() == 2 + 4 * I
+    assert oo.round() == oo
+    assert nan.round() == nan
diff --git a/symengine/tests/test_sympy_conv.py b/symengine/tests/test_sympy_conv.py
index 598a5ff4..2b513433 100644
--- a/symengine/tests/test_sympy_conv.py
+++ b/symengine/tests/test_sympy_conv.py
@@ -863,3 +863,25 @@ def test_evalf_symbolic_pyfunction():
     f = sympify(sympy.besselj(x, 1))
     assert isinstance(f, PyFunction)
     assert f.n() == f
+
+
+@unittest.skipIf(not have_sympy, "SymPy not installed")
+def test_pyfunction_func():
+    class price(sympy.Function):
+        @staticmethod
+        def _imp_(a):
+            return 3 * a
+
+    x = Symbol("x")
+    y = Symbol("y")
+    p = sympify(price(sympy.Symbol("x")))
+    assert isinstance(p, PyFunction)
+    assert p.func.__name__ == "price"
+    assert p.func(*p.args) == p
+    # Rebuilding with new arguments keeps the wrapped SymPy class
+    for q in (p.func(y), p.subs(x, y),
+              sympify(price(sympy.floor(sympy.Symbol("y")))).replace(
+                  lambda e: isinstance(e, floor), lambda e: e.args[0])):
+        assert isinstance(q, PyFunction)
+        assert q == sympify(price(sympy.Symbol("y")))
+        assert type(q.pyobject())._imp_(2) == 6