diff --git a/qutip/numpy_backend.py b/qutip/core/numpy_backend.py
similarity index 99%
rename from qutip/numpy_backend.py
rename to qutip/core/numpy_backend.py
index 5139de4bf2..88e2c162d2 100644
--- a/qutip/numpy_backend.py
+++ b/qutip/core/numpy_backend.py
@@ -1,5 +1,6 @@
 from .settings import settings
 
+
 class NumpyBackend:
     @property
     def backend(self):
@@ -9,5 +10,6 @@ def __getattr__(self, name):
         backend = object.__getattribute__(self, 'backend')
         return getattr(backend, name)
 
+
 # Initialize the numpy backend
 np = NumpyBackend()
diff --git a/qutip/core/options.py b/qutip/core/options.py
index d3bc1ec7ca..c9f156729d 100644
--- a/qutip/core/options.py
+++ b/qutip/core/options.py
@@ -128,6 +128,7 @@ class CoreOptions(QutipOptions):
         # Expect, trace, etc. will return real for hermitian matrices.
         # Hermiticity checks can be slow, stop jitting, etc.
         "auto_real_casting": True,
+        # Default backend is numpy
         "numpy_backend": np
     }
     _settings_name = "core"
diff --git a/qutip/entropy.py b/qutip/entropy.py
index bbb35f87de..3b96e0f1d0 100644
--- a/qutip/entropy.py
+++ b/qutip/entropy.py
@@ -2,7 +2,7 @@
            'concurrence', 'entropy_conditional', 'entangling_power',
            'entropy_relative']
 
-from .numpy_backend import np
+from .core.numpy_backend import np
 from .partial_transpose import partial_transpose
 from . import (ptrace, tensor, sigmay, ket2dm,
                expand_operator)
@@ -39,9 +39,9 @@ def entropy_vn(rho, base=np.e, sparse=False):
     vals = rho.eigenenergies(sparse=sparse)
     nzvals = np.clip(vals, 1e-17, None)
     if base == 2:
-        logvals = np.lib.scimath.log2(nzvals)
+        logvals = np.log2(nzvals)
     elif base == np.e:
-        logvals = np.lib.scimath.log(nzvals)
+        logvals = np.log(nzvals)
     else:
         raise ValueError("Base must be 2 or e.")
     return np.real(-sum(nzvals * logvals))
@@ -144,7 +144,7 @@ def negativity(rho, subsys, method='tracenorm', logarithmic=False):
 
 # Return the negativity value (or its logarithm if specified)
     if logarithmic:
-        return np.lib.scimath.log2(2 * N + 1)
+        return np.log2(2 * N + 1)
     else:
         return N
 
@@ -252,9 +252,9 @@ def entropy_relative(rho, sigma, base=np.e, sparse=False, tol=1e-12):
     if rho.dims != sigma.dims:
         raise ValueError("Inputs must have the same shape and dims.")
     if base == 2:
-        log_base = np.lib.scimath.log2
+        log_base = np.log2
     elif base == np.e:
-        log_base = np.lib.scimath.log
+        log_base = np.log
     else:
         raise ValueError("Base must be 2 or e.")
     # S(rho || sigma) = sum_i(p_i log p_i) - sum_ij(p_i P_ij log q_i)