diff --git a/parsimony/datasets/simulate/grad.py b/parsimony/datasets/simulate/grad.py
index bbdb6f8..8fa2d25 100644
--- a/parsimony/datasets/simulate/grad.py
+++ b/parsimony/datasets/simulate/grad.py
@@ -10,6 +10,7 @@
 """
 from six import with_metaclass
 import abc
+import math
 
 import numpy as np
 
@@ -333,7 +334,7 @@ def project(self, alpha):
         anorm = ax ** 2 + ay ** 2 + az ** 2
         i = anorm > 1.0
 
-        anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+        anorm_i = math.sqrt(anorm[i])
         ax[i] = np.divide(ax[i], anorm_i)
         ay[i] = np.divide(ay[i], anorm_i)
         az[i] = np.divide(az[i], anorm_i)
@@ -399,7 +400,7 @@ def project(self, a):
             anorm = ax ** 2 + ay ** 2 + az ** 2
             i = anorm > 1.0
 
-            anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+            anorm_i = math.sqrt(anorm[i])
             ax[i] = np.divide(ax[i], anorm_i)
             ay[i] = np.divide(ay[i], anorm_i)
             az[i] = np.divide(az[i], anorm_i)
@@ -429,7 +430,7 @@ def _Nesterov_GroupTV_project(a):
         anorm = ax ** 2 + ay ** 2 + az ** 2
         i = anorm > 1.0
 
-        anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+        anorm_i = math.sqrt(anorm[i])
         ax[i] = np.divide(ax[i], anorm_i)
         ay[i] = np.divide(ay[i], anorm_i)
         az[i] = np.divide(az[i], anorm_i)
@@ -495,7 +496,7 @@ def _Nesterov_TV_project(alpha):
     anorm = ax ** 2 + ay ** 2 + az ** 2
     i = anorm > 1.0
 
-    anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+    anorm_i = math.sqrt(anorm[i])
     ax[i] = np.divide(ax[i], anorm_i)
     ay[i] = np.divide(ay[i], anorm_i)
     az[i] = np.divide(az[i], anorm_i)
diff --git a/parsimony/functions/nesterov/grouptv.py b/parsimony/functions/nesterov/grouptv.py
index 45d8d8e..76a19a0 100644
--- a/parsimony/functions/nesterov/grouptv.py
+++ b/parsimony/functions/nesterov/grouptv.py
@@ -12,6 +12,8 @@
 @email:   lofstedt.tommy@gmail.com
 @license: BSD 3-clause.
 """
+import math
+
 import scipy.sparse as sparse
 import numpy as np
 
@@ -181,7 +183,7 @@ def project(self, a):
             anorm = ax ** 2 + ay ** 2 + az ** 2
             i = anorm > 1.0
 
-            anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+            anorm_i = math.sqrt(anorm[i])
             ax[i] = np.divide(ax[i], anorm_i)
             ay[i] = np.divide(ay[i], anorm_i)
             az[i] = np.divide(az[i], anorm_i)
@@ -225,7 +227,7 @@ def estimate_mu(self, beta):
             ay = A[g + 1].dot(beta_)
             az = A[g + 2].dot(beta_)
 
-            anorm = (ax ** 2 + ay ** 2 + az ** 2) ** 0.5  # Square root is taken here. Faster.
+            anorm = math.sqrt(ax ** 2 + ay ** 2 + az ** 2)
 
             max_norm = max(max_norm, np.max(anorm))  # The overall maximum.
 
diff --git a/parsimony/functions/nesterov/l1tv.py b/parsimony/functions/nesterov/l1tv.py
index 6a99e8c..59e71d1 100644
--- a/parsimony/functions/nesterov/l1tv.py
+++ b/parsimony/functions/nesterov/l1tv.py
@@ -232,7 +232,7 @@ def project(self, a):
             anorm_tv += a[k] ** 2
         i_tv = anorm_tv > 1.0
 
-        anorm_tv_i = anorm_tv[i_tv] ** 0.5  # Square root is taken here. Faster.
+        anorm_tv_i = math.sqrt(anorm_tv[i_tv])
         for k in range(1, len(a)):
             a[k][i_tv] = np.divide(a[k][i_tv], anorm_tv_i)
 
diff --git a/parsimony/functions/nesterov/tv.py b/parsimony/functions/nesterov/tv.py
index f829daf..6ba74ff 100644
--- a/parsimony/functions/nesterov/tv.py
+++ b/parsimony/functions/nesterov/tv.py
@@ -244,7 +244,7 @@ def lambda_max(self):
 #        anorm = ax ** 2 + ay ** 2 + az ** 2
 #        i = anorm > 1.0
 #
-#        anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+#        anorm_i = math.sqrt(anorm[i])
 #        ax[i] = np.divide(ax[i], anorm_i)
 #        ay[i] = np.divide(ay[i], anorm_i)
 #        az[i] = np.divide(az[i], anorm_i)
@@ -261,7 +261,7 @@ def project(self, a):
             anorm += a[k] ** 2
         i = anorm > 1.0
 
-        anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+        anorm_i = math.sqrt(anorm[i])
         for k in range(len(a)):
             a[k][i] = np.divide(a[k][i], anorm_i)
 
diff --git a/parsimony/functions/properties.py b/parsimony/functions/properties.py
index 2458bc7..262a626 100644
--- a/parsimony/functions/properties.py
+++ b/parsimony/functions/properties.py
@@ -15,6 +15,7 @@
 @license: BSD 3-clause.
 """
 import abc
+import math
 from six import with_metaclass
 
 import numpy as np
@@ -1042,7 +1043,7 @@ def Ata(self, a):
 #            anorm = ax ** 2 + ay ** 2 + az ** 2
 #            i = anorm > 1.0
 #
-#            anorm_i = anorm[i] ** 0.5  # Square root is taken here. Faster.
+#            anorm_i = math.sqrt(anorm[i])
 #            ax[i] = np.divide(ax[i], anorm_i)
 #            ay[i] = np.divide(ay[i], anorm_i)
 #            az[i] = np.divide(az[i], anorm_i)