diff --git a/classical/chebyshev.py b/classical/chebyshev.py
index d98bf1da75b3fb79738602b33b6add57e57588ad..c0f89fd0ed693d2ae59aa15d1d0db303800234dd 100644
--- a/classical/chebyshev.py
+++ b/classical/chebyshev.py
@@ -6,15 +6,18 @@ from scipy.special import jn, iv
 from propagators import PolynomialPropagator
 from liouville import LiouvillePowersQP
 
+
 def cheb_coef(bf, nterms, alpha):
     """ chebyshev expansion coefficients for the function exp(z*t) """
     coef = bf(range(nterms), np.full(nterms, alpha))
     coef[1:] = coef[1:]*2.
     return coef
 
-def besf(fpprec, signed, prec=None):
+
+def besself(fpprec, signed, prec=None):
     """ select the proper bessel function according to requested precision
-        and selected axis: if signed if True [-1, +1] otherwise [-i, +i] """
+        and selected axis: if signed is True [-1, +1] otherwise [-i, +i] """
+    assert fpprec in ['fixed', 'multi']
     if fpprec == 'fixed':
         bessel = iv if signed else jn
     elif fpprec == 'multi':
@@ -23,16 +26,14 @@ def besf(fpprec, signed, prec=None):
         def bessel(orders, args):
             bfun = [bess(o, a).evalf(*fargs) for o, a in zip(orders, args)]
             return np.array(bfun)
-    else:
-        bessel = None
     return bessel
 
+
 def chebyshev_scalar(z, t, delta=None, lambda_min=None, nterms=None,
                      signed=False, fpprec='fixed', prec=None):
     """ Chebyshev polynomial expansion of the scalar function exp(z*t) """
-
     assert isinstance(z, np.ndarray)
-    delta = max(abs(z)) if delta is None else delta 
+    delta = max(abs(z)) if delta is None else delta
     alpha = delta*t/2.
     lambda_min = -delta/2. if lambda_min is None else lambda_min
     nterms = int(np.ceil(abs(alpha)))+1 if nterms is None else nterms
@@ -41,10 +42,9 @@ def chebyshev_scalar(z, t, delta=None, lambda_min=None, nterms=None,
     shift = -(2.*lambda_min/delta+1.)
     sign = -1 if signed else 1
     z = z*scale+shift
-    exp = np.exp if fpprec == 'fixed' else lambda x: sp.exp(x).evalf(prec)
+    exp = np.exp if fpprec == 'fixed' else sp.exp
     scaling = exp((lambda_min+delta/2.)*t)
-    bessel_func = besf(fpprec, signed, prec=prec)
-    ch_coef = cheb_coef(bessel_func, nterms, alpha)*scaling
+    ch_coef = cheb_coef(besself(fpprec, signed, prec), nterms, alpha)*scaling
     ch_poly = np.empty((nterms, len(z)), dtype=z.dtype)
     ch_poly[0, :] = 1.
     ch_poly[1, :] = z
@@ -68,7 +68,7 @@ class Chebyshev(PolynomialPropagator):
         self.lambda_min = -self.delta/2. if lambda_min is None else lambda_min
         self.scale = 2./self.delta
         self.shift = -(2.*self.lambda_min/self.delta+1.)
-        bf = besf(fpprec, signed, prec=prec)
+        bf = besself(fpprec, signed, prec=prec)
         sign = -1 if signed else 1
         exp = np.exp if fpprec == 'fixed' else lambda x: sp.exp(x).evalf(prec)
 
diff --git a/classical/tests.py b/classical/tests.py
index 4f3d7c24cb17a09f5db6190d45e66dc5dd75c01d..f6c79d75b4007adcee58015f25139f7c5de97e59 100644
--- a/classical/tests.py
+++ b/classical/tests.py
@@ -470,10 +470,13 @@ class ChebyshevAssessmentTest(unittest.TestCase):
 
         sample_length = 10
         lambdas = (np.random.rand((sample_length))-0.5)*6*np.random.rand()
-        for sign, fact in product((False, True), (1., 1.j)):
-            apr = chebyshev_scalar(lambdas*fact, 1.1, nterms=16, signed=sign)
-            ref = np.exp(lambdas*1.1*fact)
-            self.assertTrue(np.allclose(apr, ref))
+        pars = [(False, True), (1., 1.j), ('fixed', 'multi')]
+        for sign, fact, fpprec in product(*pars):
+            values = lambdas*fact
+            apr = chebyshev_scalar(values, 1.1, nterms=24, signed=sign,
+                                   fpprec=fpprec)
+            ref = np.exp(values*1.1)
+            self.assertTrue(np.allclose(np.array(apr, dtype=ref.dtype), ref))
 
     def chebyshev_symplectic_1p_test(self):
         """ with one-particle system: anharmonic oscillator """