anchor.diff:  ;; [TODO] Look into assets (this should be a unit quaternion
   (defun rotate (vec3 quat)
     "Rotate vec3 by quaternion. Assumes quaternion is normalized."
  -  (let* ((v-mult (* 2.0 (+ (* (x quat) (aref vec3 0))
  -                           (* (y quat) (aref vec3 1))
  -                           (* (z quat) (aref vec3 2)))))
  -         (cross-mult (* 2.0 (w quat)))
  -         (p-mult (- (* cross-mult (w quat)) 1.0)))
  -    (v3:make (+ (* p-mult (aref vec3 0))
  -		(* v-mult (x quat))
  -		(* cross-mult
  -		   (- (* (y quat) (aref vec3 2))
  -		      (* (z quat) (aref vec3 1)))))
  -	     (+ (* p-mult (aref vec3 1))
  -		(* v-mult (y quat))
  -		(* cross-mult
  -		   (- (* (z quat) (aref vec3 0))
  -		      (* (x quat) (aref vec3 2)))))
  -	     (+ (* p-mult (aref vec3 2))
  -		(* v-mult (z quat))
  -		(* cross-mult
  -		   (- (* (x quat) (aref vec3 1))
  -		      (* (y quat) (aref vec3 0))))))))
  +  (let* ((v-mult (cl:* 2.0 (cl:+ (cl:* (x quat) (aref vec3 0))
  +				 (cl:* (y quat) (aref vec3 1))
  +				 (cl:* (z quat) (aref vec3 2)))))
  +         (cross-mult (cl:* 2.0 (w quat)))
  +         (p-mult (cl:- (cl:* cross-mult (w quat)) 1.0)))
  +    (v3:make (cl:+ (cl:* p-mult (aref vec3 0))
  +		   (cl:* v-mult (x quat))
  +		   (cl:* cross-mult
  +			 (cl:- (cl:* (y quat) (aref vec3 2))
  +			       (cl:* (z quat) (aref vec3 1)))))
  +	     (cl:+ (cl:* p-mult (aref vec3 1))
  +		   (cl:* v-mult (y quat))
  +		   (cl:* cross-mult
  +			 (cl:- (cl:* (z quat) (aref vec3 0))
  +			       (cl:* (x quat) (aref vec3 2)))))
  +	     (cl:+ (cl:* p-mult (aref vec3 2))
  +		   (cl:* v-mult (z quat))
  +		   (cl:* cross-mult
  +			 (cl:- (cl:* (x quat) (aref vec3 1))
  +			       (cl:* (y quat) (aref vec3 0))))))))
   
   ;; [TODO] Could be faster (see q+1 area)
   (defun lerp (start-quat end-quat pos)
  
positive.code:         ;; if angle between quaternions is less than 90 degrees
          (if (> cos-angle 0f0)
              ;; if angle is greater than zero
              (if (> (- 1.0 cos-angle) 0f0)
                  (let* ((angle (acos cos-angle))
                         (recip-sin-angle (/ 1.0 (sin angle))))
                    (list (* (sin (* (- 1.0 pos) angle))
                             recip-sin-angle)
                          (* (sin (* pos angle))
                             recip-sin-angle)))
                  ;; angle is close to zero
                  (list (- 1.0 pos) pos))
              ;; we take the shorter route
              ;; if angle is less that 180 degrees
              (if (> (+ 1.0 cos-angle) 0f0)
                  (let* ((angle (acos (- cos-angle)))
                         (recip-sin-angle (/ 1.0 (sin angle))))
                    (list (* (sin (* (- pos 1.0) angle))
                             recip-sin-angle)
                          (* (sin (* pos angle))
                             recip-sin-angle)))
                  ;; angle is close to 180 degrees
                  (list (- pos 1.0) pos))))
      (q+1 (q* start-quat start-mult)
           (q* end-quat end-mult))))
  
  (defun approx-slerp (start-quat end-quat pos)
    (let* ((cos-angle (v4:dot start-quat end-quat))
           (factor (expt (- 1.0 (* 0.7878088 cos-angle)) 2.0))
           (k (* 0.5069269 factor))
           (b (* 2.0 k))
           (c (* -3 k))
           (d (+ 1 k))
           (pos (+ (* pos (+ c (* b pos))) d)))
      ;; if angle is less than 90 degrees
      (if (> cos-angle 0f0)
          ;; use standard interp
          (q+1 (q* end-quat pos)
               (q* start-quat (- 1.0 pos)))
          ;; take shorter path
          (q+1 (q* end-quat pos)
               (q* start-quat (- pos 1.0))))))
  
negative.code: 
  (defun 0p (quat)
    (let ((w (w quat)) (x (x quat)) (y (y quat)) (z (z quat)))
      (cl:= 0f0 (cl:+ (cl:* w w) (cl:* x x) (cl:* y y) (cl:* z z)))))
  
  (defun unitp (quat)
    (let ((w (w quat)) (x (x quat)) (y (y quat)) (z (z quat)))
      (cl:= 0f0 (cl:- 1.0 (cl:* w w) (cl:* x x) (cl:* y y) (cl:* z z)))))
  
  (defun identity ()
    (q! 1.0 0.0 0.0 0.0))
  
  (defun identity-p (quat)
    (and (cl:= 0f0 (cl:- 1.0 (w quat)))
         (cl:= 0f0 (x quat))
         (cl:= 0f0 (y quat))
         (cl:= 0f0 (z quat))))
  