diff --git a/nutils/__init__.py b/nutils/__init__.py
index 4cad96d50..1d57638dc 100644
--- a/nutils/__init__.py
+++ b/nutils/__init__.py
@@ -1,4 +1,4 @@
 'Numerical Utilities for Finite Element Analysis'
 
-__version__ = version = '9a8'
+__version__ = version = '9a9'
 version_name = 'jook-sing'
diff --git a/nutils/evaluable.py b/nutils/evaluable.py
index 97d968840..36c3ba098 100644
--- a/nutils/evaluable.py
+++ b/nutils/evaluable.py
@@ -62,9 +62,9 @@ def simplified(value: 'Evaluable'):
     return value.simplified
 
 
-_type_order = bool, int, float, complex
-asdtype = lambda arg: arg if any(arg is dtype for dtype in _type_order) else {'f': float, 'i': int, 'b': bool, 'c': complex}[numpy.dtype(arg).kind]
-Dtype = typing.Union[_type_order]
+_array_dtypes = bool, int, float, complex
+asdtype = lambda arg: arg if any(arg is dtype for dtype in _array_dtypes) else {'f': float, 'i': int, 'b': bool, 'c': complex}[numpy.dtype(arg).kind]
+Dtype = typing.Union[_array_dtypes]
 
 
 def asarray(arg):
@@ -626,6 +626,9 @@ def dot(a, b, axes):
     Contract ``a`` and ``b`` along ``axes``.
     '''
 
+    a, b = _numpy_align(a, b)
+    if a.dtype == bool or b.dtype == bool:
+        raise ValueError('The boolean dot product is not supported.')
     return multiply(a, b).sum(axes)
 
 
@@ -637,9 +640,6 @@ def conjugate(arg):
         return arg
 
 
-conjugate
-
-
 def real(arg):
     arg = asarray(arg)
     if arg.dtype == complex:
@@ -838,7 +838,7 @@ class Array(Evaluable, metaclass=_ArrayMeta):
 
     def __init__(self, args, shape: typing.Tuple['Array', ...], dtype: Dtype):
         assert isinstance(shape, tuple) and all(_isindex(n) for n in shape), f'shape={shape!r}'
-        assert isinstance(dtype, type) and dtype in _type_order, f'dtype={dtype!r}'
+        assert isinstance(dtype, type) and dtype in _array_dtypes, f'dtype={dtype!r}'
         self.shape = shape
         self.dtype = dtype
         super().__init__(args=args)
@@ -1062,8 +1062,8 @@ class Orthonormal(Array):
     'make a vector orthonormal to a subspace'
 
     def __init__(self, basis: Array, vector: Array):
-        assert isinstance(basis, Array) and basis.ndim >= 2 and basis.dtype != complex, f'basis={basis!r}'
-        assert isinstance(vector, Array) and vector.ndim >= 1 and vector.dtype != complex, f'vector={vector!r}'
+        assert isinstance(basis, Array) and basis.ndim >= 2 and basis.dtype not in (bool, complex), f'basis={basis!r}'
+        assert isinstance(vector, Array) and vector.ndim >= 1 and vector.dtype not in (bool, complex), f'vector={vector!r}'
         assert equalshape(basis.shape[:-1], vector.shape)
         self._basis = basis
         self._vector = vector
@@ -1278,11 +1278,11 @@ def _derivative(self, var, seen):
 
     def _sum(self, i):
         if i == self.ndim - 1:
-            return self.func * self.length
+            return self.func * astype(self.length, self.func.dtype)
         return InsertAxis(sum(self.func, i), self.length)
 
     def _product(self):
-        return self.func**self.length
+        return self.func**astype(self.length, self.func.dtype)
 
     def _power(self, n):
         unaligned1, unaligned2, where = unalign(self, n)
@@ -1570,7 +1570,7 @@ def evalf(arr):
 
     def _derivative(self, var, seen):
         grad = derivative(self.func, var, seen)
-        funcs = Product(insertaxis(self.func, -2, self.func.shape[-1]) + Diagonalize(1 - self.func))  # replace diagonal entries by 1
+        funcs = Product(insertaxis(self.func, -2, self.func.shape[-1]) + Diagonalize(astype(1, self.func.dtype) - self.func))  # replace diagonal entries by 1
         return einsum('Ai,AiB->AB', funcs, grad)
 
     def _take(self, indices, axis):
@@ -1587,7 +1587,7 @@ class Inverse(Array):
     '''
 
     def __init__(self, func: Array):
-        assert isinstance(func, Array) and func.ndim >= 2 and _equals_simplified(func.shape[-1], func.shape[-2]), f'func={func!r}'
+        assert isinstance(func, Array) and func.dtype != bool and func.ndim >= 2 and _equals_simplified(func.shape[-1], func.shape[-2]), f'func={func!r}'
         self.func = func
         super().__init__(args=(func,), shape=func.shape, dtype=complex if func.dtype == complex else float)
 
@@ -1628,7 +1628,7 @@ def _unravel(self, axis, shape):
 class Determinant(Array):
 
     def __init__(self, func: Array):
-        assert isarray(func) and func.ndim >= 2 and _equals_simplified(func.shape[-1], func.shape[-2])
+        assert isarray(func) and func.dtype != bool and func.ndim >= 2 and _equals_simplified(func.shape[-1], func.shape[-2])
         self.func = func
         super().__init__(args=(func,), shape=func.shape[:-2], dtype=complex if func.dtype == complex else float)
 
@@ -1722,7 +1722,7 @@ def _add(self, other):
             return multiply(*common) * add(multiply(*factors), multiply(*other_factors))
         nz = factors or other_factors
         if not nz: # self equals other (up to factor ordering)
-            return self * 2
+            return self * astype(2, self.dtype)
         if len(nz) == 1 and tuple(nz)[0]._const_uniform == -1:
             # Since the subtraction x - y is stored as x + -1 * y, this handles
             # the simplification of x - x to 0. While we could alternatively
@@ -1740,7 +1740,7 @@ def _determinant(self, axis1, axis2):
             unaligned, where = unalign(fi)
             if axis1 not in where and axis2 not in where:
                 det = determinant(multiply(*factors[:i], *factors[i+1:]), (axis1, axis2))
-                scale = align(unaligned**self.shape[axis1], [i-(i > axis1)-(i > axis2) for i in where if i not in (axis1, axis2)], det.shape)
+                scale = align(unaligned**astype(self.shape[axis1], unaligned.dtype), [i-(i > axis1)-(i > axis2) for i in where if i not in (axis1, axis2)], det.shape)
                 return det * scale
 
     def _product(self):
@@ -1906,7 +1906,7 @@ def _loopsum(self, index):
         for f in self._terms:
             (dep if index in f.arguments else indep).append(f)
         if indep:
-            return add(*indep) * index.length + loop_sum(add(*dep), index)
+            return add(*indep) * astype(index.length, self.dtype) + loop_sum(add(*dep), index)
 
     def _multiply(self, other):
         f_other = [f._multiply(other) or other._multiply(f) or (f._inflations or f._diagonals) and f * other for f in self._terms]
@@ -2117,7 +2117,7 @@ def _intbounds_impl(self):
 class Power(Array):
 
     def __init__(self, func: Array, power: Array):
-        assert isinstance(func, Array), f'func={func!r}'
+        assert isinstance(func, Array) and func.dtype != bool, f'func={func!r}'
         assert isinstance(power, Array) and (power.dtype in (float, complex) or power.dtype == int and power._intbounds[0] >= 0), f'power={power!r}'
         assert equalshape(func.shape, power.shape) and func.dtype == power.dtype, 'Power({}, {})'.format(func, power)
         self.func = func
@@ -2156,7 +2156,7 @@ def _derivative(self, var, seen):
         # ln self = power * ln func
         # self` / self = power` * ln func + power * func` / func
         # self` = power` * ln func * self + power * func` * func**(power-1)
-        return einsum('A,A,AB->AB', self.power, power(self.func, self.power - 1), derivative(self.func, var, seen)) \
+        return einsum('A,A,AB->AB', self.power, power(self.func, self.power - astype(1, self.power.dtype)), derivative(self.func, var, seen)) \
             + einsum('A,A,AB->AB', ln(self.func), self, derivative(self.power, var, seen))
 
     def _power(self, n):
@@ -2164,7 +2164,7 @@ def _power(self, n):
             return
         func = self.func
         newpower = multiply(self.power, n)
-        if iszero(self.power % 2) and not iszero(newpower % 2):
+        if iszero(self.power % astype(2, self.power.dtype)) and not iszero(newpower % astype(2, newpower.dtype)):
             func = abs(func)
         return Power(func, newpower)
 
@@ -2184,7 +2184,6 @@ class Pointwise(Array):
     '''
 
     deriv = None
-    complex_deriv = None
     return_type = None
 
     def __init__(self, *args: Array, **params):
@@ -2237,9 +2236,7 @@ def _optimized_for_numpy(self):
             return constant(retval)
 
     def _derivative(self, var, seen):
-        if self.complex_deriv is not None:
-            return util.sum(einsum('A,AB->AB', deriv(*self.args, **self.params), derivative(arg, var, seen)) for arg, deriv in zip(self.args, self.complex_deriv))
-        elif self.dtype == complex or var.dtype == complex:
+        if self.dtype == complex or var.dtype == complex:
             raise NotImplementedError('The complex derivative is not implemented.')
         elif self.deriv is not None:
             return util.sum(einsum('A,AB->AB', deriv(*self.args, **self.params), derivative(arg, var, seen)) for arg, deriv in zip(self.args, self.deriv))
@@ -2256,9 +2253,29 @@ def _unravel(self, axis, shape):
         return self._newargs(*[unravel(arg, axis, shape) for arg in self.args])
 
 
-class Reciprocal(Pointwise):
+class Holomorphic(Pointwise):
+    '''
+    Abstract base class for holomorphic array functions.
+    '''
+
+    @staticmethod
+    def return_type(*dtypes, **params):
+        return_type = dtypes[-1]
+        if not all(dtype == return_type for dtype in dtypes[:-1]):
+            raise ValueError('All arguments must have the same dtype but got {} and {}.'.format(', '.join(map(str, dtypes[:-1])), return_type))
+        if return_type not in (float, complex):
+            raise ValueError(f'{self.__class__.__name__} is not defined for arguments of dtype {return_type}')
+        return return_type
+
+    def _derivative(self, var, seen):
+        if self.deriv is not None:
+            return util.sum(einsum('A,AB->AB', deriv(*self.args, **self.params), derivative(arg, var, seen)) for arg, deriv in zip(self.args, self.deriv))
+        else:
+            return super()._derivative(var, seen)
+
+
+class Reciprocal(Holomorphic):
     evalf = staticmethod(numpy.reciprocal)
-    return_type = lambda T: complex if T == complex else float
 
 
 class Negative(Pointwise):
@@ -2275,12 +2292,21 @@ def _intbounds_impl(self):
 
 class FloorDivide(Pointwise):
     evalf = staticmethod(numpy.floor_divide)
-    return_type = lambda T1, T2: complex if complex in (T1, T2) else float if float in (T1, T2) else int
+    def return_type(dividend, divisor):
+        if dividend != divisor:
+            raise ValueError(f'All arguments must have the same dtype but got {dividend} and {divisor}.')
+        if dividend == bool:
+            raise ValueError(f'The boolean floor division is not supported.')
+        return dividend
 
 
 class Absolute(Pointwise):
     evalf = staticmethod(numpy.absolute)
-    return_type = lambda T: float if T in (float, complex) else int
+
+    def return_type(T):
+        if T == bool:
+            raise ValueError('The boolean absolute value is not implemented.')
+        return float if T == complex else T
 
     def _intbounds_impl(self):
         lower, upper = self.args[0]._intbounds
@@ -2291,11 +2317,10 @@ def _intbounds_impl(self):
             return min(extrema), max(extrema)
 
 
-class Cos(Pointwise):
+class Cos(Holomorphic):
     'Cosine, element-wise.'
     evalf = staticmethod(numpy.cos)
-    complex_deriv = lambda x: -Sin(x),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: -Sin(x),
 
     def _simplified(self):
         arg, = self.args
@@ -2304,11 +2329,10 @@ def _simplified(self):
         return super()._simplified()
 
 
-class Sin(Pointwise):
+class Sin(Holomorphic):
     'Sine, element-wise.'
     evalf = staticmethod(numpy.sin)
-    complex_deriv = Cos,
-    return_type = lambda T: complex if T == complex else float
+    deriv = Cos,
 
     def _simplified(self):
         arg, = self.args
@@ -2317,86 +2341,80 @@ def _simplified(self):
         return super()._simplified()
 
 
-class Tan(Pointwise):
+class Tan(Holomorphic):
     'Tangent, element-wise.'
     evalf = staticmethod(numpy.tan)
-    complex_deriv = lambda x: Cos(x)**-2,
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: Cos(x)**astype(-2, x.dtype),
 
 
-class ArcSin(Pointwise):
+class ArcSin(Holomorphic):
     'Inverse sine, element-wise.'
     evalf = staticmethod(numpy.arcsin)
-    complex_deriv = lambda x: reciprocal(sqrt(1-x**2)),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: reciprocal(sqrt(astype(1, x.dtype)-x**astype(2, x.dtype))),
 
 
-class ArcCos(Pointwise):
+class ArcCos(Holomorphic):
     'Inverse cosine, element-wise.'
     evalf = staticmethod(numpy.arccos)
-    complex_deriv = lambda x: -reciprocal(sqrt(1-x**2)),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: -reciprocal(sqrt(astype(1, x.dtype)-x**astype(2, x.dtype))),
 
 
-class ArcTan(Pointwise):
+class ArcTan(Holomorphic):
     'Inverse tangent, element-wise.'
     evalf = staticmethod(numpy.arctan)
-    complex_deriv = lambda x: reciprocal(1+x**2),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: reciprocal(astype(1, x.dtype)+x**astype(2, x.dtype)),
 
 
-class Sinc(Pointwise):
+class Sinc(Holomorphic):
     evalf = staticmethod(numeric.sinc)
-    complex_deriv = lambda x, n: Sinc(x, n=n+1),
-    return_type = lambda T, n: complex if T == complex else float
+    deriv = lambda x, n: Sinc(x, n=n+1),
 
 
-class CosH(Pointwise):
+class CosH(Holomorphic):
     'Hyperbolic cosine, element-wise.'
     evalf = staticmethod(numpy.cosh)
-    complex_deriv = lambda x: SinH(x),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: SinH(x),
 
 
-class SinH(Pointwise):
+class SinH(Holomorphic):
     'Hyperbolic sine, element-wise.'
     evalf = staticmethod(numpy.sinh)
-    complex_deriv = CosH,
-    return_type = lambda T: complex if T == complex else float
+    deriv = CosH,
 
 
-class TanH(Pointwise):
+class TanH(Holomorphic):
     'Hyperbolic tangent, element-wise.'
     evalf = staticmethod(numpy.tanh)
-    complex_deriv = lambda x: 1 - TanH(x)**2,
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: astype(1, x.dtype) - TanH(x)**astype(2, x.dtype),
 
 
-class ArcTanH(Pointwise):
+class ArcTanH(Holomorphic):
     'Inverse hyperbolic tangent, element-wise.'
     evalf = staticmethod(numpy.arctanh)
-    complex_deriv = lambda x: reciprocal(1-x**2),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: reciprocal(astype(1, x.dtype)-x**astype(2, x.dtype)),
 
 
-class Exp(Pointwise):
+class Exp(Holomorphic):
     evalf = staticmethod(numpy.exp)
-    complex_deriv = lambda x: Exp(x),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: Exp(x),
 
 
-class Log(Pointwise):
+class Log(Holomorphic):
     evalf = staticmethod(numpy.log)
-    complex_deriv = lambda x: reciprocal(x),
-    return_type = lambda T: complex if T == complex else float
+    deriv = lambda x: reciprocal(x),
 
 
 class Mod(Pointwise):
     evalf = staticmethod(numpy.mod)
-    def return_type(T1, T2):
-        if T1 == complex or T2 == complex:
-            raise ValueError('mod is not defined for complex numbers')
-        return float if float in (T1, T2) else int
+
+    def return_type(dividend, divisor):
+        if dividend != divisor:
+            raise ValueError(f'All arguments must have the same dtype but got {dividend} and {divisor}.')
+        if dividend == bool:
+            raise ValueError(f'The boolean floor division is not supported.')
+        if dividend == complex:
+            raise ValueError(f'The complex floor division is not supported.')
+        return dividend
 
     def _intbounds_impl(self):
         dividend, divisor = self.args
@@ -2421,7 +2439,7 @@ def _simplified(self):
 
 class ArcTan2(Pointwise):
     evalf = staticmethod(numpy.arctan2)
-    deriv = lambda x, y: y / (x**2 + y**2), lambda x, y: -x / (x**2 + y**2)
+    deriv = lambda x, y: y / (x**astype(2, x.dtype) + y**astype(2, x.dtype)), lambda x, y: -x / (x**astype(2, x.dtype) + y**astype(2, x.dtype))
     def return_type(T1, T2):
         if T1 == complex or T2 == complex:
             raise ValueError('arctan2 is not defined for complex numbers')
@@ -2431,14 +2449,21 @@ def return_type(T1, T2):
 class Greater(Pointwise):
     evalf = staticmethod(numpy.greater)
     def return_type(T1, T2):
-        if T1 == complex or T2 == complex:
+        if T1 != T2:
+            raise ValueError('Cannot compare different dtypes.')
+        elif T1 == complex:
             raise ValueError('Complex numbers have no total order.')
+        elif T1 == bool:
+            raise ValueError('Use logical operators to compare booleans.')
         return bool
 
 
 class Equal(Pointwise):
     evalf = staticmethod(numpy.equal)
-    return_type = lambda T1, T2: bool
+    def return_type(T1, T2):
+        if T1 != T2:
+            raise ValueError('Cannot compare different dtypes.')
+        return bool
 
     def _simplified(self):
         a1, a2 = self.args
@@ -2456,8 +2481,12 @@ def _simplified(self):
 class Less(Pointwise):
     evalf = staticmethod(numpy.less)
     def return_type(T1, T2):
-        if T1 == complex or T2 == complex:
+        if T1 != T2:
+            raise ValueError('Cannot compare different dtypes.')
+        elif T1 == complex:
             raise ValueError('Complex numbers have no total order.')
+        elif T1 == bool:
+            raise ValueError('Use logical operators to compare booleans.')
         return bool
 
 
@@ -2511,7 +2540,10 @@ def _intbounds_impl(self):
 
 class Conjugate(Pointwise):
     evalf = staticmethod(numpy.conjugate)
-    return_type = lambda T: int if T == bool else T
+    def return_type(T):
+        if T != complex:
+            raise ValueError(f'Conjugate is not defined for arguments of type {T}')
+        return complex
 
     def _simplified(self):
         retval = self.args[0]._conjugate()
@@ -2522,7 +2554,10 @@ def _simplified(self):
 
 class Real(Pointwise):
     evalf = staticmethod(numpy.real)
-    return_type = lambda T: float if T == complex else T
+    def return_type(T):
+        if T != complex:
+            raise ValueError(f'Real is not defined for arguments of type {T}')
+        return float
 
     def _simplified(self):
         retval = self.args[0]._real()
@@ -2533,7 +2568,10 @@ def _simplified(self):
 
 class Imag(Pointwise):
     evalf = staticmethod(numpy.imag)
-    return_type = lambda T: float if T == complex else T
+    def return_type(T):
+        if T != complex:
+            raise ValueError(f'Real is not defined for arguments of type {T}')
+        return float
 
     def _simplified(self):
         retval = self.args[0]._imag()
@@ -2641,12 +2679,14 @@ def _derivative(self, var, seen):
 
 def astype(arg, dtype):
     arg = asarray(arg)
-    i = _type_order.index(arg.dtype)
-    j = _type_order.index(dtype)
-    if i > j:
+    if arg.dtype == bool and dtype != bool:
+        arg = BoolToInt(arg)
+    if arg.dtype == int and dtype != int:
+        arg = IntToFloat(arg)
+    if arg.dtype == float and dtype != float:
+        arg = FloatToComplex(arg)
+    if arg.dtype != dtype:
         raise TypeError('Downcasting is forbidden.')
-    for cast in (BoolToInt, IntToFloat, FloatToComplex)[i:j]:
-        arg = cast(arg)
     return arg
 
 
@@ -3268,7 +3308,7 @@ class Argument(DerivativeTargetBase):
     >>> from nutils import evaluable
     >>> a = evaluable.Argument('x', ())
     >>> b = evaluable.Argument('y', ())
-    >>> f = a**3 + b**2
+    >>> f = a**3. + b**2.
     >>> evaluable.derivative(f, b).simplified == 2.*b
     True
 
@@ -3896,7 +3936,7 @@ def _derivative(self, var, seen):
         d = numpy.zeros((self._degree+1,)*2, dtype=int)
         for i in range(self._degree+1):
             d[i, i+1::2] = 2*i+1
-        dself = einsum('Ai,ij->Aj', self, constant(d))
+        dself = einsum('Ai,ij->Aj', self, astype(d, self.dtype))
         return einsum('Ai,AB->AiB', dself, derivative(self._x, var, seen))
 
     def _simplified(self):
@@ -4292,7 +4332,7 @@ def _simplified(self):
         if iszero(self.func):
             return zeros_like(self)
         elif self.index not in self.func.arguments:
-            return self.func * self.index.length
+            return self.func * astype(self.index.length, self.func.dtype)
         return self.func._loopsum(self.index)
 
     def _takediag(self, axis1, axis2):
@@ -4588,11 +4628,6 @@ def _numpy_align(a, b):
 
     a = asarray(a)
     b = asarray(b)
-    if a.dtype != b.dtype:
-        if _type_order.index(a.dtype) < _type_order.index(b.dtype):
-            a = astype(a, b.dtype)
-        else:
-            b = astype(b, a.dtype)
     if not a.ndim:
         return _inflate_scalar(a, b.shape), b
     if not b.ndim:
@@ -4692,11 +4727,18 @@ def ones_like(arr):
 
 
 def reciprocal(arg):
-    return power(arg, astype(-1, float))
+    arg = asarray(arg)
+    if arg.dtype in (bool, int):
+        raise ValueError('The boolean or integer reciprocal is not supported.')
+    return power(arg, astype(-1, arg.dtype))
 
 
 def negative(arg):
-    return multiply(arg, -1)
+    arg = asarray(arg)
+    if arg.dtype == bool:
+        raise ValueError('The boolean negative is not supported.')
+    else:
+        return multiply(arg, astype(-1, arg.dtype))
 
 
 def sin(x):
@@ -4746,15 +4788,20 @@ def mod(arg1, arg2):
 
 
 def log2(arg):
-    return ln(arg) / constant(numpy.log(2))
+    arg = asarray(arg)
+    return ln(arg) / astype(numpy.log(2), arg.dtype)
 
 
 def log10(arg):
-    return ln(arg) / constant(numpy.log(10))
+    arg = asarray(arg)
+    return ln(arg) / astype(numpy.log(10), arg.dtype)
 
 
 def sqrt(arg):
-    return power(arg, constant(.5))
+    arg = asarray(arg)
+    if arg.dtype in (bool, int):
+        raise ValueError('The boolean or integer square root is not supported.')
+    return power(arg, astype(.5, arg.dtype))
 
 
 def arctan2(arg1, arg2):
@@ -4763,7 +4810,7 @@ def arctan2(arg1, arg2):
 
 def abs(arg):
     if arg.dtype == complex:
-        return sqrt(arg.real**2 + arg.imag**2)
+        return sqrt(arg.real**2. + arg.imag**2.)
     else:
         return arg * sign(arg)
 
@@ -4822,6 +4869,9 @@ def get(arg, iax, item):
 
 
 def determinant(arg, axes=(-2, -1)):
+    arg = asarray(arg)
+    if arg.dtype == bool:
+        raise ValueError('The boolean determinant is not supported.')
     return Determinant(Transpose.to_end(arg, *axes))
 
 
@@ -4840,6 +4890,9 @@ def sqrt_abs_det_gram(arg, axes=(-2, -1)):
 
 
 def inverse(arg, axes=(-2, -1)):
+    arg = asarray(arg)
+    if arg.dtype == bool:
+        raise ValueError('The boolean inverse is not supported.')
     return Transpose.from_end(Inverse(Transpose.to_end(arg, *axes)), *axes)
 
 
diff --git a/nutils/function.py b/nutils/function.py
index 9649d8903..fb666ac3b 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -803,7 +803,7 @@ def _derivative(self, var: evaluable.Array, seen: Dict[evaluable.Array, evaluabl
             fpd_expected_shape = tuple(n.__index__() for n in self.shape[self.points_dim:] + arg.shape[self.points_dim:])
             if fpd.shape != fpd_expected_shape:
                 raise ValueError('`partial_derivative` to argument {} returned an array with shape {} but {} was expected.'.format(iarg, fpd.shape, fpd_expected_shape))
-            epd = evaluable.appendaxes(fpd.lower(self.lower_args), var.shape)
+            epd = evaluable.astype(evaluable.appendaxes(fpd.lower(self.lower_args), var.shape), self.dtype)
             eda = evaluable.derivative(arg, var, seen)
             eda = evaluable.Transpose.from_end(evaluable.appendaxes(eda, self.shape[self.points_dim:]), *range(self.points_dim, self.ndim))
             result += (epd * eda).sum(range(self.ndim, self.ndim + arg.ndim - self.points_dim))
@@ -1063,7 +1063,7 @@ def lower(self, args: LowerArgs) -> evaluable.Array:
         dfunc_dref = evaluable.concatenate([evaluable.derivative(func, ref) for ref in refs], axis=-1)
         dgeom_dref = evaluable.concatenate([evaluable.derivative(geom, ref) for ref in refs], axis=-1)
         dref_dgeom = evaluable.inverse(dgeom_dref)
-        return evaluable.einsum('Ai,Aij->Aj', dfunc_dref, dref_dgeom)
+        return evaluable.einsum('Ai,Aij->Aj', dfunc_dref, evaluable.astype(dref_dgeom, dfunc_dref.dtype))
 
 
 class _SurfaceGradient(Array):
@@ -1089,7 +1089,7 @@ def lower(self, args: LowerArgs) -> evaluable.Array:
         dfunc_dref = evaluable.concatenate([evaluable.derivative(func, ref) for ref in refs], axis=-1)
         dgeom_dref = evaluable.concatenate([evaluable.derivative(geom, ref) for ref in refs], axis=-1)
         dref_dgeom = evaluable.einsum('Ajk,Aik->Aij', dgeom_dref, evaluable.inverse(evaluable.grammium(dgeom_dref)))
-        return evaluable.einsum('Ai,Aij->Aj', dfunc_dref, dref_dgeom)
+        return evaluable.einsum('Ai,Aij->Aj', dfunc_dref, evaluable.astype(dref_dgeom, dfunc_dref.dtype))
 
 
 class _Jacobian(Array):
@@ -1174,7 +1174,7 @@ def lower(self, args: LowerArgs) -> evaluable.Array:
             normal = evaluable.Take(rgrad[..., 0], i) * evaluable.Take(rgrad[..., 1], j) - evaluable.Take(rgrad[..., 1], i) * evaluable.Take(rgrad[..., 0], j)
         else:
             raise NotImplementedError
-        return normal / evaluable.InsertAxis(evaluable.sqrt(evaluable.Sum(normal**2)), normal.shape[-1])
+        return normal / evaluable.InsertAxis(evaluable.sqrt(evaluable.Sum(normal**2.)), normal.shape[-1])
 
 
 class _Concatenate(Array):
@@ -2650,7 +2650,7 @@ def f_dofs_coeffs(self, index: evaluable.Array) -> Tuple[evaluable.Array, evalua
     def lower(self, args: LowerArgs) -> evaluable.Array:
         index = _WithoutPoints(self.index).lower(args)
         coords = self.coords.lower(args)
-        leg = evaluable.Legendre(evaluable.get(coords, coords.ndim-1, evaluable.constant(0)) * 2 - 1, self._degree)
+        leg = evaluable.Legendre(evaluable.get(coords, coords.ndim-1, evaluable.constant(0)) * 2. - 1., self._degree)
         dofs = evaluable.Range(evaluable.constant(self._degree+1)) + index * (self._degree+1)
         return evaluable.Inflate(leg, dofs, evaluable.constant(self.ndofs))
 
diff --git a/nutils/sample.py b/nutils/sample.py
index 85d0b1aaf..ca992e8af 100644
--- a/nutils/sample.py
+++ b/nutils/sample.py
@@ -917,8 +917,8 @@ def __init__(self, integrand: function.Array, sample: Sample) -> None:
 
     def lower(self, args: function.LowerArgs) -> evaluable.Array:
         ielem = evaluable.loop_index('_sample_' + '_'.join(self._sample.spaces), self._sample.nelems)
-        weights = self._sample.get_evaluable_weights(ielem)
-        integrand = self._integrand.lower(args | self._sample.get_lower_args(ielem))
+        weights = evaluable.astype(self._sample.get_evaluable_weights(ielem), self.dtype)
+        integrand = evaluable.astype(self._integrand.lower(args | self._sample.get_lower_args(ielem)), self.dtype)
         elem_integral = evaluable.einsum('B,ABC->AC', weights, integrand, B=weights.ndim, C=self.ndim)
         return evaluable.loop_sum(elem_integral, ielem)
 
diff --git a/nutils/solver.py b/nutils/solver.py
index 9b827e831..a3e874b8c 100644
--- a/nutils/solver.py
+++ b/nutils/solver.py
@@ -644,7 +644,7 @@ def __init__(self, target, residual, inertia, timestep: float, constrain, argume
         self.old_new.append((timetarget+historysuffix, timetarget))
         subs0 = {new: evaluable.Argument(old, tuple(map(evaluable.constant, self.lhs0[new].shape))) for old, new in self.old_new}
         dt = evaluable.Argument(timetarget, ()) - subs0[timetarget]
-        self.residuals = tuple(res * theta + evaluable.replace_arguments(res, subs0) * (1-theta) + (inert - evaluable.replace_arguments(inert, subs0)) / dt for res, inert in zip(residual, inertia))
+        self.residuals = tuple(res * evaluable.astype(theta, res.dtype) + evaluable.replace_arguments(res, subs0) * evaluable.astype(1-theta, res.dtype) + (inert - evaluable.replace_arguments(inert, subs0)) / evaluable.astype(dt, res.dtype) for res, inert in zip(residual, inertia))
         self.jacobians = _derivative(self.residuals, target)
 
     def _step(self, lhs0, dt):
diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py
index 07973a2c2..4478ea750 100644
--- a/tests/test_evaluable.py
+++ b/tests/test_evaluable.py
@@ -25,7 +25,14 @@ def setUp(self):
         self.actual = self.op(*self.args)
         self.desired = self.n_op(*self.arg_values)
         assert numpy.isfinite(self.desired).all(), 'something is wrong with the design of this unit test'
-        self.other = numpy.random.normal(size=self.desired.shape)
+        if self.actual.dtype == bool:
+            self.other = numpy.random.normal(size=self.desired.shape) < 0
+        elif self.actual.dtype == int:
+            self.other = numpy.random.randint(0, 100, size=self.desired.shape)
+        elif self.actual.dtype == float:
+            self.other = numpy.random.normal(size=self.desired.shape)
+        elif self.actual.dtype == complex:
+            self.other = numpy.random.normal(size=self.desired.shape) + 1j * numpy.random.normal(size=self.desired.shape)
         self.pairs = [(i, j) for i in range(self.actual.ndim-1) for j in range(i+1, self.actual.ndim) if self.actual.shape[i] == self.actual.shape[j]]
         _builtin_warnings.simplefilter('ignore', evaluable.ExpensiveEvaluationWarning)
 
@@ -141,6 +148,8 @@ def test_eig(self):
                                             desired=(numpy.expand_dims(A, ax2) * numpy.expand_dims(V, ax2+1).swapaxes(ax1, ax2+1)).sum(ax2+1))
 
     def test_inv(self):
+        if self.actual.dtype == bool:
+            return
         for ax1, ax2 in self.pairs:
             trans = [i for i in range(self.desired.ndim) if i not in (ax1, ax2)] + [ax1, ax2]
             invtrans = list(map(trans.index, range(len(trans))))
@@ -149,6 +158,8 @@ def test_inv(self):
                                            actual=evaluable.inverse(self.actual, axes=(ax1, ax2)))
 
     def test_determinant(self):
+        if self.actual.dtype == bool:
+            return
         for ax1, ax2 in self.pairs:
             self.assertFunctionAlmostEqual(decimal=11,
                                            desired=numpy.linalg.det(self.desired.transpose([i for i in range(self.desired.ndim) if i not in (ax1, ax2)] + [ax1, ax2])),
@@ -245,33 +256,46 @@ def test_sumaxis(self):
                                            actual=self.actual.sum(idim))
 
     def test_add(self):
+        if self.actual.dtype == bool:
+            return
         self.assertFunctionAlmostEqual(decimal=14,
                                        desired=self.desired + self.other,
                                        actual=(self.actual + self.other))
 
     def test_multiply(self):
+        if self.actual.dtype == bool:
+            return
         self.assertFunctionAlmostEqual(decimal=14,
                                        desired=self.desired * self.other,
                                        actual=(self.actual * self.other))
 
     def test_dot(self):
+        if self.actual.dtype == bool:
+            return
         for iax in range(self.actual.ndim):
             self.assertFunctionAlmostEqual(decimal=13,
                                            desired=numeric.contract(self.desired, self.other, axis=iax),
                                            actual=evaluable.dot(self.actual, self.other, axes=iax))
 
     def test_pointwise(self):
+        if self.actual.dtype == bool:
+            return
+        dtype = complex if self.actual.dtype == complex else float
         self.assertFunctionAlmostEqual(decimal=14,
-                                       desired=numpy.sin(self.desired).astype(complex if self.desired.dtype.kind == 'c' else float),  # "astype" necessary for boolean operations (float16->float64)
-                                       actual=evaluable.sin(self.actual))
+                                       desired=numpy.sin(self.desired).astype(dtype),  # "astype" necessary for boolean operations (float16->float64)
+                                       actual=evaluable.sin(evaluable.astype(self.actual, dtype)))
 
     def test_power(self):
+        if self.actual.dtype == bool:
+            return
         self.assertFunctionAlmostEqual(decimal=13,
                                        desired=self.desired**3,
-                                       actual=(self.actual**3))
+                                       actual=(self.actual**self.actual.dtype(3)))
 
     def test_power0(self):
-        power = (numpy.arange(self.desired.size) % 2).reshape(self.desired.shape)
+        if self.actual.dtype == bool:
+            return
+        power = (numpy.arange(self.desired.size) % 2).reshape(self.desired.shape).astype(self.actual.dtype)
         self.assertFunctionAlmostEqual(decimal=13,
                                        desired=self.desired**power,
                                        actual=self.actual**power)
@@ -482,7 +506,6 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2):
 _check('sign', evaluable.sign, numpy.sign, ANY(4, 4), zerograd=True)
 _check('power', evaluable.power, numpy.power, POS(4, 4), ANY(4, 4))
 _check('power-complex', evaluable.power, numpy.power, NZC(4, 4), (2-2j)*ANC(4, 4), hasgrad=False)
-_check('power-complex-int', lambda a: evaluable.power(a, 2), lambda a: numpy.power(a, 2), ANC(4, 4))
 _check('negative', evaluable.negative, numpy.negative, ANY(4, 4))
 _check('negative-complex', evaluable.negative, numpy.negative, ANC(4, 4))
 _check('reciprocal', evaluable.reciprocal, numpy.reciprocal, NZ(4, 4))
@@ -526,7 +549,7 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2):
 _check('multiply-complex', evaluable.multiply, numpy.multiply, ANC(4, 4), 1-1j+ANC(4, 4))
 _check('dot', lambda a, b: evaluable.dot(a, b, axes=1), lambda a, b: (a*b).sum(1), ANY(4, 2, 4), ANY(4, 2, 4))
 _check('divide', evaluable.divide, lambda a, b: a * b**-1, ANY(4, 4), NZ(4, 4))
-_check('divide2', lambda a: evaluable.asarray(a)/2, lambda a: a/2, ANY(4, 1))
+_check('divide2', lambda a: evaluable.asarray(a)/2., lambda a: a/2., ANY(4, 1))
 _check('divide-complex', evaluable.divide, lambda a, b: a * b**-1, ANC(4, 4), NZC(4, 4).T)
 _check('add', evaluable.add, numpy.add, ANY(4, 4), ANY(4, 4))
 _check('add-complex', evaluable.add, numpy.add, ANC(4, 4), numpy.exp(.2j)*ANC(4, 4))
@@ -569,9 +592,9 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2):
 _check('loopsum4', lambda: evaluable.loop_sum(evaluable.Inflate(evaluable.loop_index('index', 3), evaluable.constant(0), evaluable.constant(2)), evaluable.loop_index('index', 3)), lambda: numpy.array([3, 0]))
 _check('loopsum5', lambda: evaluable.loop_sum(evaluable.loop_index('index', 1), evaluable.loop_index('index', 1)), lambda: numpy.array(0))
 _check('loopsum6', lambda: evaluable.loop_sum(evaluable.Guard(evaluable.constant(1) + evaluable.loop_index('index', 4)), evaluable.loop_index('index', 4)) * evaluable.loop_sum(evaluable.loop_index('index', 4), evaluable.loop_index('index', 4)), lambda: numpy.array(60))
-_check('loopconcatenate1', lambda a: evaluable.loop_concatenate(a+evaluable.prependaxes(evaluable.loop_index('index', 3), a.shape), evaluable.loop_index('index', 3)), lambda a: a+numpy.arange(3)[None], ANY(3, 1))
+_check('loopconcatenate1', lambda a: evaluable.loop_concatenate(a+evaluable.prependaxes(evaluable.astype(evaluable.loop_index('index', 3), float), a.shape), evaluable.loop_index('index', 3)), lambda a: a+numpy.arange(3)[None], ANY(3, 1))
 _check('loopconcatenate2', lambda: evaluable.loop_concatenate(evaluable.Elemwise(tuple(types.arraydata(numpy.arange(48).reshape(4, 4, 3)[:, :, a:b]) for a, b in util.pairwise([0, 2, 3])), evaluable.loop_index('index', 2), int), evaluable.loop_index('index', 2)), lambda: numpy.arange(48).reshape(4, 4, 3))
-_check('loopconcatenatecombined', lambda a: evaluable.loop_concatenate_combined([a+evaluable.prependaxes(evaluable.loop_index('index', 3), a.shape)], evaluable.loop_index('index', 3))[0], lambda a: a+numpy.arange(3)[None], ANY(3, 1), hasgrad=False)
+_check('loopconcatenatecombined', lambda a: evaluable.loop_concatenate_combined([a+evaluable.prependaxes(evaluable.astype(evaluable.loop_index('index', 3), float), a.shape)], evaluable.loop_index('index', 3))[0], lambda a: a+numpy.arange(3)[None], ANY(3, 1), hasgrad=False)
 _check('legendre', lambda a: evaluable.Legendre(evaluable.asarray(a), 5), lambda a: numpy.moveaxis(numpy.polynomial.legendre.legval(a, numpy.eye(6)), 0, -1), ANY(3, 4, 3))
 
 _check('polyval_1d_p0', lambda c, x: evaluable.Polyval(c, x), poly.eval_outer, POS(1), ANY(4, 1), ndim=1)