add pseudocode

diff --git a/vector_ops.mdwn b/vector_ops.mdwn
index e1982bf6092baba4951d10bbcda5db4002fbced9..11235f24e5fca2d54ba7e2d4678603969226b9ab 100644 (file)
--- a/vector_ops.mdwn
+++ b/vector_ops.mdwn
@@ -20,7 +20,7 @@ All the operands must be vectors of 3 components of a floating-point type.
 
 ## Vector dot product
 
-Computes the dot product of two vectors. Internal accuracy must be greater than the 
+Computes the dot product of two vectors. Internal accuracy must be greater than the
 input vectors and the result.
 
 ## Vector length
@@ -39,48 +39,63 @@ Known as **fmix** in GLSL.
 
 <https://en.m.wikipedia.org/wiki/Linear_interpolation>
 
-## Vector SLERP
+Pseudocode:
 
-<https://en.m.wikipedia.org/wiki/Slerp>
+    // Imprecise method, which does not guarantee v = v1 when t = 1,
+    // due to floating-point arithmetic error.
+    // This form may be used when the hardware has a native fused 
+    // multiply-add instruction.
+    float lerp(float v0, float v1, float t) {
+      return v0 + t * (v1 - v0);
+    }
 
-"""
-Quaternion slerp(Quaternion v0, Quaternion v1, double t) {
-    // Only unit quaternions are valid rotations.
-    // Normalize to avoid undefined behavior.
-    v0.normalize();
-    v1.normalize();
-
-    // Compute the cosine of the angle between the two vectors.
-    double dot = dot_product(v0, v1);
-
-    // If the dot product is negative, slerp won't take
-    // the shorter path. Note that v1 and -v1 are equivalent when
-    // the negation is applied to all four components. Fix by 
-    // reversing one quaternion.
-    if (dot < 0.0f) {
-        v1 = -v1;
-        dot = -dot;
-    }  
-
-    const double DOT_THRESHOLD = 0.9995;
-    if (dot > DOT_THRESHOLD) {
-        // If the inputs are too close for comfort, linearly interpolate
-        // and normalize the result.
-
-        Quaternion result = v0 + t*(v1 - v0);
-        result.normalize();
-        return result;
+    // Precise method, which guarantees v = v1 when t = 1.
+    float lerp(float v0, float v1, float t) {
+      return (1 - t) * v0 + t * v1;
     }
 
-    // Since dot is in range [0, DOT_THRESHOLD], acos is safe
-    double theta_0 = acos(dot);        // theta_0 = angle between input vectors
-    double theta = theta_0*t;          // theta = angle between v0 and result
-    double sin_theta = sin(theta);     // compute this value only once
-    double sin_theta_0 = sin(theta_0); // compute this value only once
+## Vector SLERP
 
-    double s0 = cos(theta) - dot * sin_theta / sin_theta_0;  // == sin(theta_0 - theta) / sin(theta_0)
-    double s1 = sin_theta / sin_theta_0;
+<https://en.m.wikipedia.org/wiki/Slerp>
 
-    return (s0 * v0) + (s1 * v1);
-}
-"""
+Pseudocode:
+
+    Quaternion slerp(Quaternion v0, Quaternion v1, double t) {
+        // Only unit quaternions are valid rotations.
+        // Normalize to avoid undefined behavior.
+        v0.normalize();
+        v1.normalize();
+
+        // Compute the cosine of the angle between the two vectors.
+        double dot = dot_product(v0, v1);
+
+        // If the dot product is negative, slerp won't take
+        // the shorter path. Note that v1 and -v1 are equivalent when
+        // the negation is applied to all four components. Fix by
+        // reversing one quaternion.
+        if (dot < 0.0f) {
+            v1 = -v1;
+            dot = -dot;
+        }
+
+        const double DOT_THRESHOLD = 0.9995;
+        if (dot > DOT_THRESHOLD) {
+            // If the inputs are too close for comfort, linearly interpolate
+            // and normalize the result.
+
+            Quaternion result = v0 + t*(v1 - v0);
+            result.normalize();
+            return result;
+        }
+
+        // Since dot is in range [0, DOT_THRESHOLD], acos is safe
+        double theta_0 = acos(dot);        // theta_0 = angle between input vectors
+        double theta = theta_0*t;          // theta = angle between v0 and result
+        double sin_theta = sin(theta);     // compute this value only once
+        double sin_theta_0 = sin(theta_0); // compute this value only once
+
+        double s0 = cos(theta) - dot * sin_theta / sin_theta_0;  // == sin(theta_0 - theta) / sin(theta_0)
+        double s1 = sin_theta / sin_theta_0;
+
+        return (s0 * v0) + (s1 * v1);
+    }