diff --git a/pyquil/quilatom.py b/pyquil/quilatom.py
index 7cd7e7c53..06bf1afe7 100644
--- a/pyquil/quilatom.py
+++ b/pyquil/quilatom.py
@@ -482,13 +482,16 @@ def _convert_to_py_expressions(
 
 class Expression(object):
     """
-    Expression involving some unbound parameters. Parameters in Quil are represented as a label
-    like '%x' for the parameter named 'x'. An example expression therefore may be '%x*(%y/4)'.
+    Expression that may involve unbound parameters or the Quil constant "pi".
+    Parameters in Quil are represented as a label like '%x' for the parameter named 'x'.
+    An example expression therefore may be '%x*(%y/4)'.
 
     Expressions may also have function calls, supported functions in Quil are sin, cos, sqrt, exp,
     and cis.
 
-    This class overrides all the Python operators that are supported by Quil.
+    This class overrides all the Python operators that are supported by Quil. Casts to integer,
+    floats, and complex numbers are also supported. These are all implemented by simplifying the
+    expression and returning the result, if the expression contains no unbound parameters.
     """
 
     def __str__(self) -> str:
@@ -537,6 +540,37 @@ def __neg__(self) -> "Mul":
     def _substitute(self, d: Any) -> ExpressionDesignator:
         return self
 
+    def _simplify(self) -> quil_rs_expr.Expression:
+        rs_expression = _convert_to_rs_expression(self)
+        rs_expression.simplify()
+        return rs_expression
+
+    def __int__(self) -> int:
+        return int(float(self))
+
+    def __float__(self) -> float:
+        expr = self._simplify()
+        if expr.is_number():
+            n = expr.to_number()
+            if n.imag == 0:
+                return float(n.real)
+            else:
+                raise TypeError(f"Expression {self} simplifies to a complex number with a non-zero imaginary part, can't cast to a float.")
+        if expr.is_pi():
+            return np.pi
+        raise ValueError(f"Expression {self} contains unbound parameters, can't cast to a float.")
+
+    def __complex__(self) -> complex:
+        expr = self._simplify()
+        if expr.is_number():
+            return expr.to_number()
+        if expr.is_pi():
+            return complex(np.pi)
+        raise ValueError(f"Expression {self} contains unbound parameters, can't cast to a float.")
+
+    def __array__(self, dtype: Optional[type] = None):
+        return np.asarray(float(self), dtype=dtype)
+
 
 ParameterSubstitutionsMapDesignator = Mapping[Union["Parameter", "MemoryReference"], ExpressionValueDesignator]
 
@@ -568,6 +602,12 @@ def substitute_array(a: Union[Sequence[Expression], np.ndarray], d: ParameterSub
     return np.array([substitute(v, d) for v in a.flat]).reshape(a.shape)
 
 
+class Pi(Expression):
+    """
+    Represents the Quil constant "pi".
+    """
+
+
 class Parameter(QuilAtom, Expression):
     """
     Parameters in Quil are represented as a label like '%x' for the parameter named 'x'.
@@ -800,6 +840,8 @@ def _expression_to_string(expression: ExpressionDesignator) -> str:
         return expression.name + "(" + _expression_to_string(expression.expression) + ")"
     elif isinstance(expression, Parameter):
         return str(expression)
+    elif isinstance(expression, Pi):
+        return "pi"
     else:
         return format_parameter(expression)
 
diff --git a/test/unit/test_expressions.py b/test/unit/test_expressions.py
new file mode 100644
index 000000000..75b2960e9
--- /dev/null
+++ b/test/unit/test_expressions.py
@@ -0,0 +1,24 @@
+import numpy as np
+
+from pyquil.gates import RX
+from pyquil.quil import Program
+from pyquil.quilatom import Pi
+
+
+class TestPi:
+    def test_print(self):
+        assert str(Pi()) == "pi"
+        expr = Pi() / 2
+        print(type(expr))
+        assert str(expr) == "pi/2"
+
+    def test_program(self):
+        program = Program(RX(Pi(), 0), RX(Pi() / 2, 1))
+        assert program.out() == "RX(pi) 0\nRX(pi/2) 1\n"
+        assert program[0] == RX(Pi(), 0)
+        assert program[1] == RX(Pi() / 2, 1)
+
+    def test_numpy(self):
+        pi_class_matrix = np.asmatrix([[Pi(), Pi() / 2], [Pi() / 3, Pi() / 4]])
+        pi_float_matrix = np.asmatrix([[np.pi, np.pi / 2], [np.pi / 3, np.pi / 4]])
+        assert np.allclose(pi_class_matrix, pi_float_matrix)