diff --git a/src/mynonlinearity.cc b/src/mynonlinearity.cc
index 4efdd603c12ac0d2f3d0f748745dbab141268a77..3e9f9028ec5be6cf9eb85bfb57e173e988b42f36 100644
--- a/src/mynonlinearity.cc
+++ b/src/mynonlinearity.cc
@@ -24,15 +24,22 @@ class MyNonlinearity {
   typedef SmallVector VectorType;
   typedef SmallMatrix MatrixType;
 
+  // directional subdifferential: at u on the line u + t*v
   void directionalSubDiff(VectorType u, VectorType v, Interval<double>& D) {
     if (u == SmallVector(0.0)) {
       D[0] = D[1] = func_.rightDifferential(0);
-    } else if (u * v > 0) {
-      D[1] = (v * u) * func_.rightDifferential(u.two_norm()) / u.two_norm();
-      D[0] = (v * u) * func_.leftDifferential(u.two_norm()) / u.two_norm();
+      return;
+    }
+    double const un = u.two_norm();
+    double const ndotp = (u * v) / (un * v.two_norm());
+    if (ndotp > 0) {
+      // If we had |v| = 1, this would be
+      // <f'_pm(|u|)*u/|u|,v>
+      D[1] = ndotp * func_.rightDifferential(un);
+      D[0] = ndotp * func_.leftDifferential(un);
     } else {
-      D[1] = (v * u) * func_.leftDifferential(u.two_norm()) / u.two_norm();
-      D[0] = (v * u) * func_.rightDifferential(u.two_norm()) / u.two_norm();
+      D[1] = ndotp * func_.leftDifferential(un);
+      D[0] = ndotp * func_.rightDifferential(un);
     }
   }