diff --git a/src/main/java/org/joml/AxisAngle4d.java b/src/main/java/org/joml/AxisAngle4d.java
index 0d3e221c..96f56d04 100644
--- a/src/main/java/org/joml/AxisAngle4d.java
+++ b/src/main/java/org/joml/AxisAngle4d.java
@@ -76,7 +76,7 @@ public AxisAngle4d(AxisAngle4d a) {
         x = a.x;
         y = a.y;
         z = a.z;
-        angle = (a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI);
+        angle = (a.angle < 0.0 ? Math.PI_TIMES_2 + a.angle % Math.PI_TIMES_2 : a.angle) % Math.PI_TIMES_2;
     }
 
     /**
@@ -89,7 +89,7 @@ public AxisAngle4d(AxisAngle4f a) {
         x = a.x;
         y = a.y;
         z = a.z;
-        angle = (a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI);
+        angle = (a.angle < 0.0 ? Math.PI_TIMES_2 + a.angle % Math.PI_TIMES_2 : a.angle) % Math.PI_TIMES_2;
     }
 
     /**
@@ -158,7 +158,7 @@ public AxisAngle4d(double angle, double x, double y, double z) {
         this.x = x;
         this.y = y;
         this.z = z;
-        this.angle = (angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI);
+        this.angle = (angle < 0.0 ? Math.PI_TIMES_2 + angle % Math.PI_TIMES_2 : angle) % Math.PI_TIMES_2;
     }
 
     /**
@@ -192,7 +192,7 @@ public AxisAngle4d set(AxisAngle4d a) {
         x = a.x;
         y = a.y;
         z = a.z;
-        angle = (a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI);
+        angle = (a.angle < 0.0 ? Math.PI_TIMES_2 + a.angle % Math.PI_TIMES_2 : a.angle) % Math.PI_TIMES_2;
         return this;
     }
 
@@ -207,7 +207,7 @@ public AxisAngle4d set(AxisAngle4f a) {
         x = a.x;
         y = a.y;
         z = a.z;
-        angle = (a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI);
+        angle = (a.angle < 0.0 ? Math.PI_TIMES_2 + a.angle % Math.PI_TIMES_2 : a.angle) % Math.PI_TIMES_2;
         return this;
     }
 
@@ -228,7 +228,7 @@ public AxisAngle4d set(double angle, double x, double y, double z) {
         this.x = x;
         this.y = y;
         this.z = z;
-        this.angle = (angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI);
+        this.angle = (angle < 0.0 ? Math.PI_TIMES_2 + angle % Math.PI_TIMES_2 : angle) % Math.PI_TIMES_2;
         return this;
     }
 
@@ -345,24 +345,27 @@ public AxisAngle4d set(Matrix3fc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                double invX = 1.0 / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                double invY = 1.0 / y;
+                x = xy * invY;
+                z = yz * invY;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                double invZ = 1.0 / z;
+                x = xz * invZ;
+                y = yz * invZ;
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -405,24 +408,27 @@ public AxisAngle4d set(Matrix3dc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                double invX = 1.0 / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                double invY = 1.0 / y;
+                x = xy * invY;
+                z = yz * invY;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                double invZ = 1.0 / z;
+                x = xz * invZ;
+                y = yz * invZ;
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -465,24 +471,27 @@ public AxisAngle4d set(Matrix4fc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                double invX = 1.0 / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                double invY = 1.0 / y;
+                x = xy * invY;
+                z = yz * invY;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                double invZ = 1.0 / z;
+                x = xz * invZ;
+                y = yz * invZ;
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -525,24 +534,27 @@ public AxisAngle4d set(Matrix4x3fc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                double invX = 1.0 / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                double invY = 1.0 / y;
+                x = xy * invY;
+                z = yz * invY;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                double invZ = 1.0 / z;
+                x = xz * invZ;
+                y = yz * invZ;
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -585,24 +597,27 @@ public AxisAngle4d set(Matrix4dc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                double invX = 1.0 / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                double invY = 1.0 / y;
+                x = xy * invY;
+                z = yz * invY;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                double invZ = 1.0 / z;
+                x = xz * invZ;
+                y = yz * invZ;
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -744,7 +759,7 @@ public AxisAngle4d normalize() {
      */
     public AxisAngle4d rotate(double ang) {
         angle += ang;
-        angle = (angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI);
+        angle = (angle < 0.0 ? Math.PI_TIMES_2 + angle % Math.PI_TIMES_2 : angle) % Math.PI_TIMES_2;
         return this;
     }
 
@@ -868,7 +883,7 @@ public int hashCode() {
         final int prime = 31;
         int result = 1;
         long temp;
-        temp = Double.doubleToLongBits((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        temp = Double.doubleToLongBits((angle < 0.0 ? Math.PI_TIMES_2 + angle % Math.PI_TIMES_2 : angle) % Math.PI_TIMES_2);
         result = prime * result + (int) (temp ^ (temp >>> 32));
         temp = Double.doubleToLongBits(x);
         result = prime * result + (int) (temp ^ (temp >>> 32));
@@ -887,8 +902,8 @@ public boolean equals(Object obj) {
         if (getClass() != obj.getClass())
             return false;
         AxisAngle4d other = (AxisAngle4d) obj;
-        if (Double.doubleToLongBits((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI)) != 
-                Double.doubleToLongBits((other.angle < 0.0 ? Math.PI + Math.PI + other.angle % (Math.PI + Math.PI) : other.angle) % (Math.PI + Math.PI)))
+        if (Double.doubleToLongBits((angle < 0.0 ? Math.PI_TIMES_2 + angle % Math.PI_TIMES_2 : angle) % Math.PI_TIMES_2) !=
+                Double.doubleToLongBits((other.angle < 0.0 ? Math.PI_TIMES_2 + other.angle % Math.PI_TIMES_2 : other.angle) % Math.PI_TIMES_2))
             return false;
         if (Double.doubleToLongBits(x) != Double.doubleToLongBits(other.x))
             return false;
diff --git a/src/main/java/org/joml/AxisAngle4f.java b/src/main/java/org/joml/AxisAngle4f.java
index 5a48ff0e..715587cb 100644
--- a/src/main/java/org/joml/AxisAngle4f.java
+++ b/src/main/java/org/joml/AxisAngle4f.java
@@ -76,7 +76,7 @@ public AxisAngle4f(AxisAngle4f a) {
         x = a.x;
         y = a.y;
         z = a.z;
-        angle = (float) ((a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI));
+        angle = (a.angle < 0.0 ? Math.PI_TIMES_2_f + a.angle % Math.PI_TIMES_2_f : a.angle) % Math.PI_TIMES_2_f;
     }
 
     /**
@@ -120,7 +120,7 @@ public AxisAngle4f(float angle, float x, float y, float z) {
         this.x = x;
         this.y = y;
         this.z = z;
-        this.angle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        this.angle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % Math.PI_TIMES_2_f : angle) % Math.PI_TIMES_2_f;
     }
 
     /**
@@ -145,7 +145,7 @@ public AxisAngle4f set(AxisAngle4f a) {
         y = a.y;
         z = a.z;
         angle = a.angle;
-        angle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        angle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % Math.PI_TIMES_2_f : angle) % Math.PI_TIMES_2_f;
         return this;
     }
 
@@ -161,7 +161,7 @@ public AxisAngle4f set(AxisAngle4d a) {
         y = (float) a.y;
         z = (float) a.z;
         angle = (float) a.angle;
-        angle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        angle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % Math.PI_TIMES_2_f : angle) % Math.PI_TIMES_2_f;
         return this;
     }
 
@@ -182,7 +182,7 @@ public AxisAngle4f set(float angle, float x, float y, float z) {
         this.x = x;
         this.y = y;
         this.z = z;
-        this.angle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        this.angle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % Math.PI_TIMES_2_f : angle) % Math.PI_TIMES_2_f;
         return this;
     }
 
@@ -286,24 +286,27 @@ public AxisAngle4f set(Matrix3fc m) {
             float yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                float invX = 1.0f / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                float invZ = 1.0f / z;
+                x = xy * invZ;
+                z = yz * invZ;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                float invY = 1.0f / y;
+                x = xz * invY;
+                y = yz * invY;
             }
             return this;
         }
-        float s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        float s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -346,24 +349,27 @@ public AxisAngle4f set(Matrix3dc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = (float) Math.sqrt(xx);
-                y = (float) (xy / x);
-                z = (float) (xz / x);
+                float invX = 1.0f / x;
+                y = (float) (xy * invX);
+                z = (float) (xz * invX);
             } else if (yy > zz) {
                 y = (float) Math.sqrt(yy);
-                x = (float) (xy / y);
-                z = (float) (yz / y);
+                float invY = 1.0f / y;
+                x = (float) (xy * invY);
+                z = (float) (yz * invY);
             } else {
                 z = (float) Math.sqrt(zz);
-                x = (float) (xz / z);
-                y = (float) (yz / z);
+                float invZ = 1.0f / z;
+                x = (float) (xz * invZ);
+                y = (float) (yz * invZ);
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = (float) Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (float) ((nm12 - nm21) / s);
-        y = (float) ((nm20 - nm02) / s);
-        z = (float) ((nm01 - nm10) / s);
+        x = (float) ((nm12 - nm21) * s);
+        y = (float) ((nm20 - nm02) * s);
+        z = (float) ((nm01 - nm10) * s);
         return this;
     }
 
@@ -406,24 +412,27 @@ public AxisAngle4f set(Matrix4fc m) {
             float yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                float invX = 1.0f / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                float invZ = 1.0f / z;
+                x = xy * invZ;
+                z = yz * invZ;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                float invY = 1.0f / y;
+                x = xz * invY;
+                y = yz * invY;
             }
             return this;
         }
-        float s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        float s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -466,24 +475,27 @@ public AxisAngle4f set(Matrix4x3fc m) {
             float yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = Math.sqrt(xx);
-                y = xy / x;
-                z = xz / x;
+                float invX = 1.0f / x;
+                y = xy * invX;
+                z = xz * invX;
             } else if (yy > zz) {
                 y = Math.sqrt(yy);
-                x = xy / y;
-                z = yz / y;
+                float invZ = 1.0f / z;
+                x = xy * invZ;
+                z = yz * invZ;
             } else {
                 z = Math.sqrt(zz);
-                x = xz / z;
-                y = yz / z;
+                float invY = 1.0f / y;
+                x = xz * invY;
+                y = yz * invY;
             }
             return this;
         }
-        float s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        float s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (nm12 - nm21) / s;
-        y = (nm20 - nm02) / s;
-        z = (nm01 - nm10) / s;
+        x = (nm12 - nm21) * s;
+        y = (nm20 - nm02) * s;
+        z = (nm01 - nm10) * s;
         return this;
     }
 
@@ -526,24 +538,27 @@ public AxisAngle4f set(Matrix4dc m) {
             double yz = (nm21 + nm12) / 4;
             if ((xx > yy) && (xx > zz)) {
                 x = (float) Math.sqrt(xx);
-                y = (float) (xy / x);
-                z = (float) (xz / x);
+                float invX = 1.0f / x;
+                y = (float) (xy * invX);
+                z = (float) (xz * invX);
             } else if (yy > zz) {
                 y = (float) Math.sqrt(yy);
-                x = (float) (xy / y);
-                z = (float) (yz / y);
+                float invY = 1.0f / y;
+                x = (float) (xy * invY);
+                z = (float) (yz * invY);
             } else {
                 z = (float) Math.sqrt(zz);
-                x = (float) (xz / z);
-                y = (float) (yz / z);
+                float invZ = 1.0f / z;
+                x = (float) (xz * invZ);
+                y = (float) (yz * invZ);
             }
             return this;
         }
-        double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
+        double s = Math.invsqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
         angle = (float) Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
-        x = (float) ((nm12 - nm21) / s);
-        y = (float) ((nm20 - nm02) / s);
-        z = (float) ((nm01 - nm10) / s);
+        x = (float) ((nm12 - nm21) * s);
+        y = (float) ((nm20 - nm02) * s);
+        z = (float) ((nm01 - nm10) * s);
         return this;
     }
 
@@ -685,7 +700,7 @@ public AxisAngle4f normalize() {
      */
     public AxisAngle4f rotate(float ang) {
         angle += ang;
-        angle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        angle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % (Math.PI_TIMES_2_f) : angle) % (Math.PI_TIMES_2_f);
         return this;
     }
 
@@ -777,7 +792,7 @@ public String toString(NumberFormat formatter) {
     public int hashCode() {
         final int prime = 31;
         int result = 1;
-        float nangle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
+        float nangle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % Math.PI_TIMES_2_f : angle) % Math.PI_TIMES_2_f;
         result = prime * result + Float.floatToIntBits(nangle);
         result = prime * result + Float.floatToIntBits(x);
         result = prime * result + Float.floatToIntBits(y);
@@ -793,8 +808,8 @@ public boolean equals(Object obj) {
         if (getClass() != obj.getClass())
             return false;
         AxisAngle4f other = (AxisAngle4f) obj;
-        float nangle = (float) ((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
-        float nangleOther = (float) ((other.angle < 0.0 ? Math.PI + Math.PI + other.angle % (Math.PI + Math.PI) : other.angle) % (Math.PI + Math.PI));
+        float nangle = (angle < 0.0 ? Math.PI_TIMES_2_f + angle % Math.PI_TIMES_2_f : angle) % Math.PI_TIMES_2_f;
+        float nangleOther = (other.angle < 0.0 ? Math.PI_TIMES_2_f + other.angle % Math.PI_TIMES_2_f : other.angle) % Math.PI_TIMES_2_f;
         if (Float.floatToIntBits(nangle) != Float.floatToIntBits(nangleOther))
             return false;
         if (Float.floatToIntBits(x) != Float.floatToIntBits(other.x))
diff --git a/src/main/java/org/joml/Intersectiond.java b/src/main/java/org/joml/Intersectiond.java
index ead23d19..7ffec528 100644
--- a/src/main/java/org/joml/Intersectiond.java
+++ b/src/main/java/org/joml/Intersectiond.java
@@ -1058,7 +1058,7 @@ public static double intersectRayPlane(Vector3dc origin, Vector3dc dir, Vector3d
     public static double intersectRayPlane(double originX, double originY, double originZ, double dirX, double dirY, double dirZ,
             double a, double b, double c, double d, double epsilon) {
         double denom = a * dirX + b * dirY + c * dirZ;
-        if (denom < 0.0) {
+        if (denom < epsilon) {
             double t = -(a * originX + b * originY + c * originZ + d) / denom;
             if (t >= 0.0)
                 return t;
@@ -1257,6 +1257,7 @@ public static double findClosestPointsLineSegments(
         double d2x = b1X - b0X, d2y = b1Y - b0Y, d2z = b1Z - b0Z;
         double rX = a0X - b0X, rY = a0Y - b0Y, rZ = a0Z - b0Z;
         double a = d1x * d1x + d1y * d1y + d1z * d1z;
+        double invA = 1.0 / a;
         double e = d2x * d2x + d2y * d2y + d2z * d2z;
         double f = d2x * rX + d2y * rY + d2z * rZ;
         double EPSILON = 1E-8;
@@ -1277,7 +1278,7 @@ public static double findClosestPointsLineSegments(
             if (e <= EPSILON) {
                 // Second segment degenerates into a point
                 t = 0.0;
-                s = Math.min(Math.max(-c / a, 0.0), 1.0);
+                s = Math.min(Math.max(-c * invA, 0.0), 1.0);
             } else {
                 // The general nondegenerate case starts here
                 double b = d1x * d2x + d1y * d2y + d1z * d2z;
@@ -1296,10 +1297,10 @@ public static double findClosestPointsLineSegments(
                 // and clamp s to [0, 1]
                 if (t < 0.0) {
                     t = 0.0;
-                    s = Math.min(Math.max(-c / a, 0.0), 1.0);
+                    s = Math.min(Math.max(-c * invA, 0.0), 1.0);
                 } else if (t > 1.0) {
                     t = 1.0;
-                    s = Math.min(Math.max((b - c) / a, 0.0), 1.0);
+                    s = Math.min(Math.max((b - c) * invA, 0.0), 1.0);
                 }
             }
         }
@@ -1612,25 +1613,30 @@ public static Vector3d findClosestPointOnRectangle(
         double qX = aX, qY = aY, qZ = aZ;
         double dist = dX * abX + dY * abY + dZ * abZ;
         double maxdist = abX * abX + abY * abY + abZ * abZ;
+        double invMaxdist = 1.0 / maxdist;
+        double distTimesInvMaxDist;
         if (dist >= maxdist) {
             qX += abX;
             qY += abY;
             qZ += abZ;
         } else if (dist > 0.0) {
-            qX += (dist / maxdist) * abX;
-            qY += (dist / maxdist) * abY;
-            qZ += (dist / maxdist) * abZ;
+            distTimesInvMaxDist = dist * invMaxdist;
+            qX += distTimesInvMaxDist * abX;
+            qY += distTimesInvMaxDist * abY;
+            qZ += distTimesInvMaxDist * abZ;
         }
         dist = dX * acX + dY * acY + dZ * acZ;
         maxdist = acX * acX + acY * acY + acZ * acZ;
+        invMaxdist = 1.0 / maxdist;
         if (dist >= maxdist) {
             qX += acX;
             qY += acY;
             qZ += acZ;
         } else if (dist > 0.0) {
-            qX += (dist / maxdist) * acX;
-            qY += (dist / maxdist) * acY;
-            qZ += (dist / maxdist) * acZ;
+            distTimesInvMaxDist = dist * invMaxdist;
+            qX += distTimesInvMaxDist * acX;
+            qY += distTimesInvMaxDist * acY;
+            qZ += distTimesInvMaxDist * acZ;
         }
         res.x = qX;
         res.y = qY;
@@ -1710,12 +1716,13 @@ public static int intersectSweptSphereTriangle(
         double invLen = Math.invsqrt(a * a + b * b + c * c);
         double signedDist = (a * centerX + b * centerY + c * centerZ + d) * invLen;
         double dot = (a * velX + b * velY + c * velZ) * invLen;
+        double invDot = 1.0 / dot;
         if (dot < epsilon && dot > -epsilon)
             return 0;
-        double pt0 = (radius - signedDist) / dot;
+        double pt0 = (radius - signedDist) * invDot;
         if (pt0 > maxT)
             return 0;
-        double pt1 = (-radius - signedDist) / dot;
+        double pt1 = (-radius - signedDist) * invDot;
         double p0X = centerX - radius * a * invLen + velX * pt0;
         double p0Y = centerY - radius * b * invLen + velY * pt0;
         double p0Z = centerZ - radius * c * invLen + velZ * pt0;
@@ -1736,7 +1743,8 @@ public static int intersectSweptSphereTriangle(
         double centerV0Y = centerY - v0Y;
         double centerV0Z = centerZ - v0Z;
         double B0 = 2.0 * (velX * centerV0X + velY * centerV0Y + velZ * centerV0Z);
-        double C0 = centerV0X * centerV0X + centerV0Y * centerV0Y + centerV0Z * centerV0Z - radius2;
+        double baseTo0Len = centerV0X * centerV0X + centerV0Y * centerV0Y + centerV0Z * centerV0Z;
+        double C0 = baseTo0Len - radius2;
         double root0 = computeLowestRoot(A, B0, C0, t0);
         if (root0 < t0) {
             pointAndTime.x = v0X;
@@ -1777,10 +1785,9 @@ public static int intersectSweptSphereTriangle(
             t0 = root2;
             isect = POINT_ON_TRIANGLE_VERTEX_2;
         }
-        double velLen = velX * velX + velY * velY + velZ * velZ;
+        double velLen = A;
         // test against edge10
         double len10 = v10X * v10X + v10Y * v10Y + v10Z * v10Z;
-        double baseTo0Len = centerV0X * centerV0X + centerV0Y * centerV0Y + centerV0Z * centerV0Z;
         double v10Vel = (v10X * velX + v10Y * velY + v10Z * velZ);
         double A10 = len10 * -velLen + v10Vel * v10Vel;
         double v10BaseTo0 = v10X * -centerV0X + v10Y * -centerV0Y + v10Z * -centerV0Z;
@@ -1833,7 +1840,6 @@ public static int intersectSweptSphereTriangle(
             pointAndTime.y = v1Y + f21 * v21Y;
             pointAndTime.z = v1Z + f21 * v21Z;
             pointAndTime.w = root21;
-            t0 = root21;
             isect = POINT_ON_TRIANGLE_EDGE_12;
         }
         return isect;
@@ -1859,8 +1865,9 @@ private static double computeLowestRoot(double a, double b, double c, double max
         if (determinant < 0.0)
             return Double.POSITIVE_INFINITY;
         double sqrtD = Math.sqrt(determinant);
-        double r1 = (-b - sqrtD) / (2.0 * a);
-        double r2 = (-b + sqrtD) / (2.0 * a);
+        double invA2 = 1.0 / (2.0 * a);
+        double r1 = (-b - sqrtD) * invA2;
+        double r2 = (-b + sqrtD) * invA2;
         if (r1 > r2) {
             double temp = r2;
             r2 = r1;
@@ -2578,8 +2585,7 @@ public static boolean testRayTriangleFront(double originX, double originY, doubl
         double v = (dirX * qvecX + dirY * qvecY + dirZ * qvecZ);
         if (v < 0.0 || u + v > det)
             return false;
-        double invDet = 1.0 / det;
-        double t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
+        double t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) / det;
         return t >= epsilon;
     }
 
@@ -2796,9 +2802,7 @@ public static double intersectRayTriangleFront(double originX, double originY, d
         double v = dirX * qvecX + dirY * qvecY + dirZ * qvecZ;
         if (v < 0.0 || u + v > det)
             return -1.0;
-        double invDet = 1.0 / det;
-        double t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
-        return t;
+        return (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) / det;
     }
 
     /**
@@ -2908,8 +2912,7 @@ public static double intersectRayTriangle(double originX, double originY, double
         double v = (dirX * qvecX + dirY * qvecY + dirZ * qvecZ) * invDet;
         if (v < 0.0 || u + v > 1.0)
             return -1.0;
-        double t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
-        return t;
+        return (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
     }
 
     /**
diff --git a/src/main/java/org/joml/Intersectionf.java b/src/main/java/org/joml/Intersectionf.java
index e28778f0..6b72f6c6 100644
--- a/src/main/java/org/joml/Intersectionf.java
+++ b/src/main/java/org/joml/Intersectionf.java
@@ -1058,7 +1058,7 @@ public static float intersectRayPlane(Vector3fc origin, Vector3fc dir, Vector3fc
     public static float intersectRayPlane(float originX, float originY, float originZ, float dirX, float dirY, float dirZ,
             float a, float b, float c, float d, float epsilon) {
         float denom = a * dirX + b * dirY + c * dirZ;
-        if (denom < 0.0f) {
+        if (denom < epsilon) {
             float t = -(a * originX + b * originY + c * originZ + d) / denom;
             if (t >= 0.0f)
                 return t;
@@ -1257,6 +1257,7 @@ public static float findClosestPointsLineSegments(
         float d2x = b1X - b0X, d2y = b1Y - b0Y, d2z = b1Z - b0Z;
         float rX = a0X - b0X, rY = a0Y - b0Y, rZ = a0Z - b0Z;
         float a = d1x * d1x + d1y * d1y + d1z * d1z;
+        float invA = 1.0f / a;
         float e = d2x * d2x + d2y * d2y + d2z * d2z;
         float f = d2x * rX + d2y * rY + d2z * rZ;
         float EPSILON = 1E-5f;
@@ -1277,7 +1278,7 @@ public static float findClosestPointsLineSegments(
             if (e <= EPSILON) {
                 // Second segment degenerates into a point
                 t = 0.0f;
-                s = Math.min(Math.max(-c / a, 0.0f), 1.0f);
+                s = Math.min(Math.max(-c * invA, 0.0f), 1.0f);
             } else {
                 // The general nondegenerate case starts here
                 float b = d1x * d2x + d1y * d2y + d1z * d2z;
@@ -1296,10 +1297,10 @@ public static float findClosestPointsLineSegments(
                 // and clamp s to [0, 1]
                 if (t < 0.0) {
                     t = 0.0f;
-                    s = Math.min(Math.max(-c / a, 0.0f), 1.0f);
+                    s = Math.min(Math.max(-c * invA, 0.0f), 1.0f);
                 } else if (t > 1.0) {
                     t = 1.0f;
-                    s = Math.min(Math.max((b - c) / a, 0.0f), 1.0f);
+                    s = Math.min(Math.max((b - c) * invA, 0.0f), 1.0f);
                 }
             }
         }
@@ -1612,25 +1613,30 @@ public static Vector3f findClosestPointOnRectangle(
         float qX = aX, qY = aY, qZ = aZ;
         float dist = dX * abX + dY * abY + dZ * abZ;
         float maxdist = abX * abX + abY * abY + abZ * abZ;
+        float invMaxdist = 1.0f / maxdist;
+        float distTimesInvMaxDist;
         if (dist >= maxdist) {
             qX += abX;
             qY += abY;
             qZ += abZ;
         } else if (dist > 0.0f) {
-            qX += (dist / maxdist) * abX;
-            qY += (dist / maxdist) * abY;
-            qZ += (dist / maxdist) * abZ;
+            distTimesInvMaxDist = dist * invMaxdist;
+            qX += distTimesInvMaxDist * abX;
+            qY += distTimesInvMaxDist * abY;
+            qZ += distTimesInvMaxDist * abZ;
         }
         dist = dX * acX + dY * acY + dZ * acZ;
         maxdist = acX * acX + acY * acY + acZ * acZ;
+        invMaxdist = 1.0f / maxdist;
         if (dist >= maxdist) {
             qX += acX;
             qY += acY;
             qZ += acZ;
         } else if (dist > 0.0f) {
-            qX += (dist / maxdist) * acX;
-            qY += (dist / maxdist) * acY;
-            qZ += (dist / maxdist) * acZ;
+            distTimesInvMaxDist = dist * invMaxdist;
+            qX += distTimesInvMaxDist * acX;
+            qY += distTimesInvMaxDist * acY;
+            qZ += distTimesInvMaxDist * acZ;
         }
         res.x = qX;
         res.y = qY;
@@ -1710,12 +1716,13 @@ public static int intersectSweptSphereTriangle(
         float invLen = Math.invsqrt(a * a + b * b + c * c);
         float signedDist = (a * centerX + b * centerY + c * centerZ + d) * invLen;
         float dot = (a * velX + b * velY + c * velZ) * invLen;
+        float invDot = 1.0f / dot;
         if (dot < epsilon && dot > -epsilon)
             return 0;
-        float pt0 = (radius - signedDist) / dot;
+        float pt0 = (radius - signedDist) * invDot;
         if (pt0 > maxT)
             return 0;
-        float pt1 = (-radius - signedDist) / dot;
+        float pt1 = (-radius - signedDist) * invDot;
         float p0X = centerX - radius * a * invLen + velX * pt0;
         float p0Y = centerY - radius * b * invLen + velY * pt0;
         float p0Z = centerZ - radius * c * invLen + velZ * pt0;
@@ -1736,7 +1743,8 @@ public static int intersectSweptSphereTriangle(
         float centerV0Y = centerY - v0Y;
         float centerV0Z = centerZ - v0Z;
         float B0 = 2.0f * (velX * centerV0X + velY * centerV0Y + velZ * centerV0Z);
-        float C0 = centerV0X * centerV0X + centerV0Y * centerV0Y + centerV0Z * centerV0Z - radius2;
+        float baseTo0Len = centerV0X * centerV0X + centerV0Y * centerV0Y + centerV0Z * centerV0Z;
+        float C0 = baseTo0Len - radius2;
         float root0 = computeLowestRoot(A, B0, C0, t0);
         if (root0 < t0) {
             pointAndTime.x = v0X;
@@ -1777,10 +1785,9 @@ public static int intersectSweptSphereTriangle(
             t0 = root2;
             isect = POINT_ON_TRIANGLE_VERTEX_2;
         }
-        float velLen = velX * velX + velY * velY + velZ * velZ;
+        float velLen = A;
         // test against edge10
         float len10 = v10X * v10X + v10Y * v10Y + v10Z * v10Z;
-        float baseTo0Len = centerV0X * centerV0X + centerV0Y * centerV0Y + centerV0Z * centerV0Z;
         float v10Vel = (v10X * velX + v10Y * velY + v10Z * velZ);
         float A10 = len10 * -velLen + v10Vel * v10Vel;
         float v10BaseTo0 = v10X * -centerV0X + v10Y * -centerV0Y + v10Z * -centerV0Z;
@@ -1833,7 +1840,6 @@ public static int intersectSweptSphereTriangle(
             pointAndTime.y = v1Y + f21 * v21Y;
             pointAndTime.z = v1Z + f21 * v21Z;
             pointAndTime.w = root21;
-            t0 = root21;
             isect = POINT_ON_TRIANGLE_EDGE_12;
         }
         return isect;
@@ -1859,8 +1865,9 @@ private static float computeLowestRoot(float a, float b, float c, float maxR) {
         if (determinant < 0.0f)
             return Float.POSITIVE_INFINITY;
         float sqrtD = (float) Math.sqrt(determinant);
-        float r1 = (-b - sqrtD) / (2.0f * a);
-        float r2 = (-b + sqrtD) / (2.0f * a);
+        float invA2 = 1.0f / (2.0f * a);
+        float r1 = (-b - sqrtD) * invA2;
+        float r2 = (-b + sqrtD) * invA2;
         if (r1 > r2) {
             float temp = r2;
             r2 = r1;
@@ -2578,8 +2585,7 @@ public static boolean testRayTriangleFront(float originX, float originY, float o
         float v = (dirX * qvecX + dirY * qvecY + dirZ * qvecZ);
         if (v < 0.0f || u + v > det)
             return false;
-        float invDet = 1.0f / det;
-        float t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
+        float t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) / det;
         return t >= epsilon;
     }
 
@@ -2796,9 +2802,7 @@ public static float intersectRayTriangleFront(float originX, float originY, floa
         float v = dirX * qvecX + dirY * qvecY + dirZ * qvecZ;
         if (v < 0.0f || u + v > det)
             return -1.0f;
-        float invDet = 1.0f / det;
-        float t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
-        return t;
+        return (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) / det;
     }
 
     /**
@@ -2908,8 +2912,7 @@ public static float intersectRayTriangle(float originX, float originY, float ori
         float v = (dirX * qvecX + dirY * qvecY + dirZ * qvecZ) * invDet;
         if (v < 0.0f || u + v > 1.0f)
             return -1.0f;
-        float t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
-        return t;
+        return (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
     }
 
     /**
diff --git a/src/main/java/org/joml/Matrix4d.java b/src/main/java/org/joml/Matrix4d.java
index 8fc1a282..2bf094c1 100644
--- a/src/main/java/org/joml/Matrix4d.java
+++ b/src/main/java/org/joml/Matrix4d.java
@@ -14274,7 +14274,8 @@ public Matrix4d setFrustumLH(double left, double right, double bottom, double to
     public Matrix4d setFromIntrinsic(double alphaX, double alphaY, double gamma, double u0, double v0, int imgWidth, int imgHeight, double near, double far) {
         double l00 = 2.0 / imgWidth;
         double l11 = 2.0 / imgHeight;
-        double l22 = 2.0 / (near - far);
+        double invNearFar = 1.0 / (near - far);
+        double l22 = 2.0 * invNearFar;
         this.m00 = l00 * alphaX;
         this.m01 = 0.0;
         this.m02 = 0.0;
@@ -14285,7 +14286,7 @@ public Matrix4d setFromIntrinsic(double alphaX, double alphaY, double gamma, dou
         this.m13 = 0.0;
         this.m20 = l00 * u0 - 1.0;
         this.m21 = l11 * v0 - 1.0;
-        this.m22 = l22 * -(near + far) + (far + near) / (near - far);
+        this.m22 = l22 * -(near + far) + (far + near) * invNearFar;
         this.m23 = -1.0;
         this.m30 = 0.0;
         this.m31 = 0.0;
diff --git a/src/main/java/org/joml/Matrix4f.java b/src/main/java/org/joml/Matrix4f.java
index 116acf3b..e96a546b 100644
--- a/src/main/java/org/joml/Matrix4f.java
+++ b/src/main/java/org/joml/Matrix4f.java
@@ -11288,7 +11288,8 @@ public Matrix4f setFrustumLH(float left, float right, float bottom, float top, f
     public Matrix4f setFromIntrinsic(float alphaX, float alphaY, float gamma, float u0, float v0, int imgWidth, int imgHeight, float near, float far) {
         float l00 = 2.0f / imgWidth;
         float l11 = 2.0f / imgHeight;
-        float l22 = 2.0f / (near - far);
+        float invNearFar = 1.0f / (near - far);
+        float l22 = 2.0f * invNearFar;
         return this
         ._m00(l00 * alphaX)
         ._m01(0.0f)
@@ -11300,7 +11301,7 @@ public Matrix4f setFromIntrinsic(float alphaX, float alphaY, float gamma, float
         ._m13(0.0f)
         ._m20(l00 * u0 - 1.0f)
         ._m21(l11 * v0 - 1.0f)
-        ._m22(l22 * -(near + far) + (far + near) / (near - far))
+        ._m22(l22 * -(near + far) + (far + near) * invNearFar)
         ._m23(-1.0f)
         ._m30(0.0f)
         ._m31(0.0f)
diff --git a/src/main/java/org/joml/Quaterniond.java b/src/main/java/org/joml/Quaterniond.java
index 370a3039..406623d7 100644
--- a/src/main/java/org/joml/Quaterniond.java
+++ b/src/main/java/org/joml/Quaterniond.java
@@ -737,10 +737,11 @@ public Quaterniond fromAxisAngleRad(Vector3dc axis, double angle) {
     public Quaterniond fromAxisAngleRad(double axisX, double axisY, double axisZ, double angle) {
         double hangle = angle / 2.0;
         double sinAngle = Math.sin(hangle);
-        double vLength = Math.sqrt(axisX * axisX + axisY * axisY + axisZ * axisZ);
-        x = axisX / vLength * sinAngle;
-        y = axisY / vLength * sinAngle;
-        z = axisZ / vLength * sinAngle;
+        double invVLength = 1.0 / Math.sqrt(axisX * axisX + axisY * axisY + axisZ * axisZ);
+        double invVLengthTimesSinAngle = invVLength * sinAngle;
+        x = axisX * invVLengthTimesSinAngle;
+        y = axisY * invVLengthTimesSinAngle;
+        z = axisZ * invVLengthTimesSinAngle;
         w = Math.cosFromSin(sinAngle, hangle);
         return this;
     }
diff --git a/src/main/java/org/joml/Quaternionf.java b/src/main/java/org/joml/Quaternionf.java
index ad0f38b6..233fab6e 100644
--- a/src/main/java/org/joml/Quaternionf.java
+++ b/src/main/java/org/joml/Quaternionf.java
@@ -946,10 +946,11 @@ public Quaternionf fromAxisAngleRad(Vector3fc axis, float angle) {
     public Quaternionf fromAxisAngleRad(float axisX, float axisY, float axisZ, float angle) {
         float hangle = angle / 2.0f;
         float sinAngle = Math.sin(hangle);
-        float vLength = Math.sqrt(axisX * axisX + axisY * axisY + axisZ * axisZ);
-        x = axisX / vLength * sinAngle;
-        y = axisY / vLength * sinAngle;
-        z = axisZ / vLength * sinAngle;
+        float invVLength = 1.0f / Math.sqrt(axisX * axisX + axisY * axisY + axisZ * axisZ);
+        float invVLengthTimesSinAngle = invVLength * sinAngle;
+        x = axisX * invVLengthTimesSinAngle;
+        y = axisY * invVLengthTimesSinAngle;
+        z = axisZ * invVLengthTimesSinAngle;
         w = Math.cosFromSin(sinAngle, hangle);
         return this;
     }
diff --git a/src/main/java/org/joml/sampling/Convolution.java b/src/main/java/org/joml/sampling/Convolution.java
index 72fa72a3..78cd200e 100644
--- a/src/main/java/org/joml/sampling/Convolution.java
+++ b/src/main/java/org/joml/sampling/Convolution.java
@@ -65,15 +65,17 @@ public static void gaussianKernel(int rows, int cols, float sigma, FloatBuffer d
         }
         float sum = 0.0f;
         int pos = dest.position();
+        float a = (float) (1.0 / (2.0 * sigma * sigma));
         for (int i = 0, y = -(rows - 1) / 2; y <= (rows - 1) / 2; y++) {
             for (int x = -(cols - 1) / 2; x <= (cols - 1) / 2; x++, i++) {
-                float k = (float) Math.exp(-(y * y + x * x) / (2.0 * sigma * sigma));
+                float k = (float) Math.exp(-(y * y + x * x) * a);
                 dest.put(pos + i, k);
                 sum += k;
             }
         }
+        sum = 1.0f / sum;
         for (int i = 0; i < rows * cols; i++) {
-            dest.put(pos + i, dest.get(pos + i) / sum);
+            dest.put(pos + i, dest.get(pos + i) * sum);
         }
     }
 //#endif
@@ -105,15 +107,17 @@ public static void gaussianKernel(int rows, int cols, float sigma, float[] dest)
             throw new IllegalArgumentException("dest must have a size of at least " + (rows * cols));
         }
         float sum = 0.0f;
+        float a = (float) (1.0 / (2.0 * sigma * sigma));
         for (int i = 0, y = -(rows - 1) / 2; y <= (rows - 1) / 2; y++) {
             for (int x = -(cols - 1) / 2; x <= (cols - 1) / 2; x++, i++) {
-                float k = (float) Math.exp(-(y * y + x * x) / (2.0 * sigma * sigma));
+                float k = (float) Math.exp(-(y * y + x * x) * a);
                 dest[i] = k;
                 sum += k;
             }
         }
+        sum = 1.0f / sum;
         for (int i = 0; i < rows * cols; i++) {
-            dest[i] = dest[i] / sum;
+            dest[i] = dest[i] * sum;
         }
     }
 
diff --git a/src/main/java/org/joml/sampling/PoissonSampling.java b/src/main/java/org/joml/sampling/PoissonSampling.java
index b762d234..6b4159f7 100644
--- a/src/main/java/org/joml/sampling/PoissonSampling.java
+++ b/src/main/java/org/joml/sampling/PoissonSampling.java
@@ -127,8 +127,9 @@ private void compute(int k, Callback2d callback) {
         }
 
         private boolean searchNeighbors(float px, float py) {
-            int row = (int) ((py + diskRadius) / cellSize);
-            int col = (int) ((px + diskRadius) / cellSize);
+            float invCellSize = 1.0f / cellSize;
+            int row = (int) ((py + diskRadius) * invCellSize);
+            int col = (int) ((px + diskRadius) * invCellSize);
             if (grid[row * numCells + col] != null)
                 return true;
             int minX = Math.max(0, col - 1);
@@ -151,8 +152,9 @@ private boolean searchNeighbors(float px, float py) {
         }
 
         private void insert(Vector2f p) {
-            int row = (int) ((p.y + diskRadius) / cellSize);
-            int col = (int) ((p.x + diskRadius) / cellSize);
+            float invCellSize = 1.0f / cellSize;
+            int row = (int) ((p.y + diskRadius) * invCellSize);
+            int col = (int) ((p.x + diskRadius) * invCellSize);
             grid[row * numCells + col] = p;
         }
 
diff --git a/src/main/java/org/joml/sampling/SpiralSampling.java b/src/main/java/org/joml/sampling/SpiralSampling.java
index 80b42b2b..c157ff1c 100644
--- a/src/main/java/org/joml/sampling/SpiralSampling.java
+++ b/src/main/java/org/joml/sampling/SpiralSampling.java
@@ -59,9 +59,11 @@ public SpiralSampling(long seed) {
      *            will be called for each sample generated
      */
     public void createEquiAngle(float radius, int numRotations, int numSamples, Callback2d callback) {
+        float invNumSamples = 1.0f / numSamples;
+        float iNSMinusOne = 1.0f / (numSamples - 1);
         for (int sample = 0; sample < numSamples; sample++) {
-            float angle = 2.0f * (float) Math.PI * (sample * numRotations) / numSamples;
-            float r = radius * sample / (numSamples - 1);
+            float angle = 2.0f * (float) Math.PI * (sample * numRotations) * invNumSamples;
+            float r = radius * sample * iNSMinusOne;
             float x = (float) Math.sin_roquen_9(angle + 0.5f * (float) Math.PI) * r;
             float y = (float) Math.sin_roquen_9(angle) * r;
             callback.onNewSample(x, y);
@@ -89,9 +91,12 @@ public void createEquiAngle(float radius, int numRotations, int numSamples, Call
      */
     public void createEquiAngle(float radius, int numRotations, int numSamples, float jitter, Callback2d callback) {
         float spacing = radius / numRotations;
+        float spacingTimesJitter = spacing * jitter;
+        float invNumSamples = 1.0f / numSamples;
+        float iNSMinusOne = 1.0f / (numSamples - 1);
         for (int sample = 0; sample < numSamples; sample++) {
-            float angle = 2.0f * (float) Math.PI * (sample * numRotations) / numSamples;
-            float r = radius * sample / (numSamples - 1) + (rnd.nextFloat() * 2.0f - 1.0f) * spacing * jitter;
+            float angle = 2.0f * (float) Math.PI * (sample * numRotations) * invNumSamples;
+            float r = radius * sample * iNSMinusOne + (rnd.nextFloat() * 2.0f - 1.0f) * spacingTimesJitter;
             float x = (float) Math.sin_roquen_9(angle + 0.5f * (float) Math.PI) * r;
             float y = (float) Math.sin_roquen_9(angle) * r;
             callback.onNewSample(x, y);
diff --git a/src/main/java/org/joml/sampling/StratifiedSampling.java b/src/main/java/org/joml/sampling/StratifiedSampling.java
index 1556513d..6c938ced 100644
--- a/src/main/java/org/joml/sampling/StratifiedSampling.java
+++ b/src/main/java/org/joml/sampling/StratifiedSampling.java
@@ -56,10 +56,11 @@ public StratifiedSampling(long seed) {
      *            will be called for each generated sample position
      */
     public void generateRandom(int n, Callback2d callback) {
+        float invN = 1.0f / n;
         for (int y = 0; y < n; y++) {
             for (int x = 0; x < n; x++) {
-                float sampleX = (rnd.nextFloat() / n + (float) x / n) * 2.0f - 1.0f;
-                float sampleY = (rnd.nextFloat() / n + (float) y / n) * 2.0f - 1.0f;
+                float sampleX = (rnd.nextFloat() * invN + (float) x * invN) * 2.0f - 1.0f;
+                float sampleY = (rnd.nextFloat() * invN + (float) y * invN) * 2.0f - 1.0f;
                 callback.onNewSample(sampleX, sampleY);
             }
         }
@@ -81,10 +82,11 @@ public void generateRandom(int n, Callback2d callback) {
     public void generateCentered(int n, float centering, Callback2d callback) {
         float start = centering * 0.5f;
         float end = 1.0f - centering;
+        float invN = 1.0f / n;
         for (int y = 0; y < n; y++) {
             for (int x = 0; x < n; x++) {
-                float sampleX = ((start + rnd.nextFloat() * end) / n + (float) x / n) * 2.0f - 1.0f;
-                float sampleY = ((start + rnd.nextFloat() * end) / n + (float) y / n) * 2.0f - 1.0f;
+                float sampleX = ((start + rnd.nextFloat() * end) * invN + (float) x * invN) * 2.0f - 1.0f;
+                float sampleY = ((start + rnd.nextFloat() * end) * invN + (float) y * invN) * 2.0f - 1.0f;
                 callback.onNewSample(sampleX, sampleY);
             }
         }