From ad3232936c06880e6465ddbd9f1e9e75c314d9e7 Mon Sep 17 00:00:00 2001
From: Thibaut Horel <thibaut.horel@gmail.com>
Date: Sun, 14 Sep 2014 09:37:42 -0400
Subject: [ps1] add C implementation

---
 ps1/lsb.c | 72 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 72 insertions(+)
 create mode 100644 ps1/lsb.c

diff --git a/ps1/lsb.c b/ps1/lsb.c
new file mode 100644
index 0000000..c0fac2c
--- /dev/null
+++ b/ps1/lsb.c
@@ -0,0 +1,72 @@
+#include <stdio.h>
+
+const unsigned int w = 16;
+const unsigned int l = 4; // sqrt(w)
+const unsigned int b = 0x8888; // 1000 1000 1000 1000
+const unsigned int a = 0x1111; // 0001 0001 0001 0001
+const unsigned int mm = 0xf;   // 0000 0000 0000 1111
+const unsigned int m = 0x2490; // 0010 0100 1001 0000
+const unsigned int c = 0xf731; // 1111 0111 0011 0001
+
+char* itoa(unsigned int i, char* str, unsigned int base) {
+    const char chars[] = "0123456789abcdef";
+    char* p = str; char b;
+    while (i) {
+        *str++ = chars[i % base];
+        i /= base;
+    }
+    *str = '\0';
+    while (p < str) { // reverse str through pairwise swapping
+        b = *p;
+        *p++ = *--str;
+        *str = b; 
+    }
+    return str;
+}
+
+void print_n(unsigned int i) {
+    char buffer[2*w];
+    itoa(i, buffer, 2);
+    printf("%s\n", buffer);
+}
+
+/* assumes that x is 0 everywhere except possibly the leading bit of each block
+ * and computes the number of leading bits equal to one (4 op) */
+unsigned int count_ones(unsigned int x) { return ((x * a) >> (w-1)) & mm; }
+
+/* assumes that x is 0 everywhere except possibly the leading bit of each block
+ * and compress the leading bits of each block into the least significant
+ * block keeping their order (3 op) */
+unsigned int perfect_sketch(unsigned int x) { return ((x * m) >> w) & mm; }
+
+/* assumes that x consists of a single block and repeat it l times (1 op) */
+unsigned int repeat(unsigned int x) { return x * a; }
+
+/* the leading bit of each block indicates whether or not the corresponding
+ * block in x is empty or not (7 op) */
+unsigned int non_empty_clusters(unsigned int x) { return (x & b) | (b & ~( b - (x & ~b))); }
+
+/* the leading bit of each block indicates whether or not the corresponding
+ * block in x is empty or not (6 op) */
+unsigned int empty_clusters(unsigned int x) { return (x ^ b) & (b & ( b - (x & ~b))); }
+
+/* computes the lsb of a single block (10 op) */
+unsigned int lsb_short(unsigned int x) {
+    return count_ones(empty_clusters(repeat(x) & c));
+}
+
+/* computes the lsb of x (34 op) */
+unsigned int lsb(unsigned int x) {
+    int z = lsb_short(perfect_sketch(non_empty_clusters(x)));
+    if (z == l) {
+        return -1;
+    } else {
+        return z*l + lsb_short((x >> z * l) & mm);
+    }
+}
+
+int main() {
+    unsigned int i = 0x8a00;
+    print_n(i);
+    printf("%d\n", lsb(i));
+}
-- 
cgit v1.2.3-70-g09d2