From df3750364a7077dd95e4798de05748481af0c657 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B6rg=20Thalheim?= <joerg@higgsboson.tk>
Date: Wed, 18 Dec 2013 10:01:06 +0100
Subject: [PATCH] correct fix-point arithmetic

---
 sqrt.c | 42 ++++++++++++++++++++++++++++--------------
 1 file changed, 28 insertions(+), 14 deletions(-)

diff --git a/sqrt.c b/sqrt.c
index 78d951a..d0abeb0 100644
--- a/sqrt.c
+++ b/sqrt.c
@@ -6,9 +6,12 @@
 
 typedef uint32_t u32;
 
+#define COMMA 12
+#define COMMA_FIX (COMMA >> 1)
+
 void heron_sqrt(u32 *mem) {
   int i = mem[0];
-  while (i > 0) {
+  while (i >= 1) {
     u32 s = mem[i];
     if (s > 1) {
       // x_0 = (s + 1) / 2
@@ -24,23 +27,29 @@ void heron_sqrt(u32 *mem) {
         }
         old_x = x;
       }
-      mem[i] = x;
+      mem[i] = x << COMMA_FIX;
     }
     i--;
   }
 }
 
 int main() {
-  u32 data[] = {0,
+  float data[] = {
+    0,
     1,
     2,
-    4,
+    2.5,
+    2.25,
+    2.125,
+    2.0625,
+    2.03125, // precision test
     3<<2,
     1<<3,
     1<<4,
     1<<5,
     1<<6,
     5<<2,
+    11<<7,
     1<<7,
     1<<8,
     1<<16,
@@ -50,21 +59,26 @@ int main() {
     1<<20, // 32 - 14 = 20
     1<<24, // 32 -  8 = 24
     1<<30,
-    0xffffffff};
-  const int length = sizeof(data)/sizeof(u32);
-
-  u32* mem = malloc(sizeof(data) + sizeof(u32));
-  memcpy(&mem[1], data, sizeof(data));
+    0xffffffff
+  };
+  const int length = sizeof(data)/sizeof(float);
+  int i;
 
+  u32* mem = malloc((length + 1) * sizeof(u32));
   mem[0] = length;
+
+  for (i = 1; i <= length; i++) {
+    mem[i] = (u32) (data[i] * (1 << COMMA));
+  }
+
   heron_sqrt(mem);
 
-  int i = 0;
   for (i = 0; i < length; i++) {
-    printf("s             = %u\n", data[i]);
-
-    printf("sqrt(s)       = %f\n", sqrt(data[i]));
-    printf("heron_sqrt(s) = %d\n", mem[i+1]);
+    if (data[i] < (0xfffffff >> (COMMA))) {
+      printf("s             = %f\n", data[i]);
+      printf("sqrt(s)       = %f\n", sqrt(data[i]));
+      printf("heron_sqrt(s) = %f\n", ((float) mem[i]) / (1 << COMMA));
+    }
     puts("\n");
   }
   return 0;