From 61e27754d91185729a5570d2e9147ee6f85a4faf Mon Sep 17 00:00:00 2001
From: Tyge Løvset <tylo@norceresearch.no>
Date: Mon, 29 Mar 2021 09:20:03 +0200
Subject: Fixed bugs in cbits. (bits_count() and is_disjoint()). Renamed
 is_disjoined() to disjoined() and cbits_subset() to cbits_subset_of().
 Impoved prime.c sieveOfEratosthenes() now 2.5x faster.

---
 examples/prime.c | 47 ++++++++++++++++++++++++++---------------------
 stc/cbits.h      | 25 ++++++++++++++-----------
 2 files changed, 40 insertions(+), 32 deletions(-)

diff --git a/examples/prime.c b/examples/prime.c
index 5c5d3969..f11f1627 100644
--- a/examples/prime.c
+++ b/examples/prime.c
@@ -1,39 +1,44 @@
 #include <stdio.h>
+#include <math.h>
+#include <time.h>
 #include <stc/cbits.h>
 
-static inline cbits sieveOfEratosthenes(size_t n)
+cbits sieveOfEratosthenes(size_t n)
 {
-    cbits primes = cbits_with_size(n + 1, true);
-    cbits_reset(&primes, 0);
-    cbits_reset(&primes, 1);
-
-    c_forrange (i, size_t, 2, n+1) {
-        // If primes[i] is not changed, then it is a prime
-        if (cbits_test(primes, i) && i*i <= n) {
-            c_forrange (j, size_t, i*i, n+1, i) {
-                cbits_reset(&primes, j);
+    cbits bits = cbits_with_pattern(n, 0xAAAAAAAAAAAAAAAA);
+    size_t q = sqrt(n);
+    cbits_reset(&bits, 1);
+    cbits_set(&bits, 2);
+
+    for (size_t j, i = 3; i <= q; i += 2) {
+        for (j = i; j < n; j += 2) {
+            if (cbits_test(bits, j)) {
+                i = j;
+                break;
             }
         }
+        for (j = i*i; j < n; j += i*2)
+            cbits_reset(&bits, j);
     }
-    return primes;
-} 
+    return bits;
+}
 
 
 int main(void)
 {
-    int n = 100000000;
+    size_t n = 100000000;
     printf("computing prime numbers up to %u\n", n);
-    
-    cbits primes = sieveOfEratosthenes(n);
-    puts("done");
-    
+
+    clock_t t1 = clock();
+    cbits primes = sieveOfEratosthenes(n + 1);
     size_t np = cbits_count(primes);
-    printf("number of primes: %zu\n", np);
+    clock_t t2 = clock();
+
+    printf("number of primes: %zu, time: %f\n", np, (t2 - t1)/(float)CLOCKS_PER_SEC);
 
-    printf("2 ");
-    c_forrange (i, int, 3, 1001, 2) {
+    for (size_t i=0; i<=1000; ++i)
        if (cbits_test(primes, i)) printf("%d ", i);
-    }
     puts("");
+
     cbits_del(&primes);
 }
\ No newline at end of file
diff --git a/stc/cbits.h b/stc/cbits.h
index d5618bba..c1ff05ac 100644
--- a/stc/cbits.h
+++ b/stc/cbits.h
@@ -64,8 +64,8 @@ STC_API cbits_t     cbits_clone(cbits_t other);
 STC_API void        cbits_resize(cbits_t* self, size_t size, bool value);
 STC_API cbits_t*    cbits_assign(cbits_t* self, cbits_t other);
 STC_API size_t      cbits_count(cbits_t set);
-STC_API bool        cbits_is_subset(cbits_t set, cbits_t other);
-STC_API bool        cbits_is_disjoint(cbits_t set, cbits_t other);
+STC_API bool        cbits_subset_of(cbits_t set, cbits_t other);
+STC_API bool        cbits_disjoint(cbits_t set, cbits_t other);
 
 STC_INLINE cbits_t  cbits_init() { cbits_t set = {NULL, 0}; return set; }
 STC_INLINE void     cbits_clear(cbits_t* self) { self->size = 0; }
@@ -180,6 +180,11 @@ STC_DEF cbits_t cbits_with_size(size_t size, bool value) {
     cbits_set_all(&set, value);
     return set;
 }
+STC_DEF cbits_t cbits_with_pattern(size_t size, uint64_t pattern) {
+    cbits_t set = {(uint64_t *) c_malloc(((size + 63) >> 6) * 8), size};
+    cbits_set_all64(&set, pattern);
+    return set;
+}
 STC_DEF cbits_t cbits_from_str(const char* str) {
     const char* p = str; while (*p) ++p;
     cbits_t set = cbits_with_size(p - str, false);
@@ -197,26 +202,24 @@ STC_DEF cbits_t cbits_clone(cbits_t other) {
     return set;
 }
 STC_DEF size_t cbits_count(cbits_t s) {
-    size_t count = 0, n = ((s.size + 63) >> 6) - 1;
-    if (s.size > 0) {
-        for (size_t i=0; i<n; ++i) count += cpopcount64(s._arr[i]);
-        count += cpopcount64(s._arr[n] & ((1ull << (s.size & 63)) - 1));
-    }
+    size_t count = 0, n = s.size >> 6;
+    for (size_t i = 0; i < n; ++i) count += cpopcount64(s._arr[i]);
+    count += cpopcount64(s._arr[n] & ((1ull << (s.size & 63)) - 1));
     return count;
 }
 
 #define _cbits_SETOP(OPR, x) \
     assert(s.size == other.size); \
     if (s.size == 0) return false; /* ? */ \
-    size_t n = ((s.size + 63) >> 6) - 1; \
-    for (size_t i=0; i<n; ++i) \
+    size_t n = s.size >> 6; \
+    for (size_t i = 0; i < n; ++i) \
         if ((s._arr[i] OPR other._arr[i]) != x) \
             return false; \
     uint64_t i = n, m = (1ull << (s.size & 63)) - 1; \
     return ((s._arr[i] & m) OPR (other._arr[i] & m)) == (x & m)
 
-STC_DEF bool cbits_is_subset(cbits_t s, cbits_t other) { _cbits_SETOP(|, s._arr[i]); }
-STC_DEF bool cbits_is_disjoint(cbits_t s, cbits_t other) { _cbits_SETOP(^, ~0ull); }
+STC_DEF bool cbits_subset_of(cbits_t s, cbits_t other) { _cbits_SETOP(|, s._arr[i]); }
+STC_DEF bool cbits_disjoint(cbits_t s, cbits_t other) { _cbits_SETOP(&, 0); }
 
 #endif
 #endif
-- 
cgit v1.2.3