From ec630d85034ddf239bf420d19183d1fd78434484 Mon Sep 17 00:00:00 2001
From: Jens Nockert <jens@nockert.se>
Date: Sat, 17 Mar 2012 14:28:41 +0100
Subject: [PATCH] Minor refactoring mistake

---
 lib/csmath.js   | 51 +++++++++++++++++++++++++++++++------------------
 src/norm.erb.js |  2 +-
 2 files changed, 33 insertions(+), 20 deletions(-)

diff --git a/lib/csmath.js b/lib/csmath.js
index 7ce3e44..f08f5af 100644
--- a/lib/csmath.js
+++ b/lib/csmath.js
@@ -237,7 +237,7 @@ function(global) {
     CSMath.min = Math.min
     CSMath.max = Math.max
 
-    CSMath.hypot = norm2
+    CSMath.hypot = CSMath.norm2
     /* TODO: The trigonometric functions have low precision */
 
     CSMath.round = Math.round
@@ -249,34 +249,47 @@ function(global) {
     if (CSMath.accurate['sqrt']) {
         CSMath.sqrt = Math.sqrt
     } else {
-        /* Algorithm from 'Divide, Square Root and Remainder Algorithms for the IA64 Architecture */
+        /*
+         * Algorithm from 'Divide, Square Root and Remainder Algorithms for the IA64
+         * Architecture' by Intel
+         *
+         *  - It is proven correct under a few assumptions
+         *      1. That Math.sqrt is at least as accurate as the IA64 approximation
+         *         instruction.
+         *      2. That we're using the IA64 extended floating point format. (Which
+         *         we are not)
+         *  - I haven't checked if this algorithm actually works in double-precision,
+         *    but my intuition says that it should be fine.
+         */
 
-        Math.sqrt = function(a) {
-            var y0 = 1.0 / Math.sqrta(a)
+        var sqrtApproximation = Math.sqrt
 
-            var H0 = 0.5 * y0,
-                S0 = (a * y0)
+        CSMath.sqrt = function(a) {
+                var y0 = 1.0 / sqrtApproximation(a)
 
-                var d0 = -Math.fms(S0, H0, 0.5)
+                var H0 = 0.5 * y0,
+                    S0 = (a * y0)
 
-            var H1 = Math.fma(d0, H0, H0),
-                S1 = Math.fma(d0, S0, S0)
+                    var d0 = -CSMath.fms(S0, H0, 0.5)
 
-                var d1 = -Math.fms(S1, H1, 0.5)
+                var H1 = CSMath.fma(d0, H0, H0),
+                    S1 = CSMath.fma(d0, S0, S0)
 
-            var H2 = Math.fma(d1, H1, H1),
-                S2 = Math.fma(d1, S1, S1)
+                    var d1 = -CSMath.fms(S1, H1, 0.5)
 
-                var d2 = -Math.fms(S2, H2, 0.5),
-                e2 = -Math.fms(S2, S2, a)
+                var H2 = CSMath.fma(d1, H1, H1),
+                    S2 = CSMath.fma(d1, S1, S1)
 
-                var H3 = Math.fma(d2, H2, H2),
-                S3 = Math.fma(e2, H2, S2)
+                    var d2 = -CSMath.fms(S2, H2, 0.5),
+                    e2 = -CSMath.fms(S2, S2, a)
 
-                var e3 = -Math.fms(S3, S3, a)
+                    var H3 = CSMath.fma(d2, H2, H2),
+                    S3 = CSMath.fma(e2, H2, S2)
 
-            return Math.fma(e3, H3, S3)
-        }
+                    var e3 = -CSMath.fms(S3, S3, a)
+
+                return CSMath.fma(e3, H3, S3)
+            }
     }
 
     /* TODO: The root functions have low precision */
diff --git a/src/norm.erb.js b/src/norm.erb.js
index b3b19a6..97495b8 100644
--- a/src/norm.erb.js
+++ b/src/norm.erb.js
@@ -51,4 +51,4 @@ CSMath.hypot2 = function() {
 CSMath.min = Math.min
 CSMath.max = Math.max
 
-CSMath.hypot = norm2
\ No newline at end of file
+CSMath.hypot = CSMath.norm2
\ No newline at end of file