From 7500d7130461aa781527ba30a63726cd9cea9224 Mon Sep 17 00:00:00 2001
From: Jan Huwald <jh@sotun.de>
Date: Mon, 3 Jun 2013 15:52:21 +0200
Subject: count garden of eden states


diff --git a/cacount.cpp b/cacount.cpp
index 26b2069..295a40a 100644
--- a/cacount.cpp
+++ b/cacount.cpp
@@ -65,7 +65,11 @@ void init(Trans &t) {
     thread::hardware_concurrency());
 }
 
-void findCycle(Trans &t, Trans &c) {
+void findCycle(Trans &t, Trans &c, bitset<maxState> &reachable) {
+  // compute reachability
+  iterState([&](int s) {
+      reachable[t[s]] = 1;
+    });
   // forward to t=maxState; now every state is in a cycle
   iterTrans(logState, [&](int s) {
       t[s] = t[t[s]]; }, (char*) "fwd time: ");
@@ -78,15 +82,20 @@ void findCycle(Trans &t, Trans &c) {
       t[s] = t[t[s]]; }, (char*) "cycles: ");
 }
 
-void cycleStat(Trans &c) {
+void cycleStat(Trans &c, bitset<maxState> &reachable) {
   struct Stat {
-    State basin, len, minState;
-    Stat() : basin(0), len(0), minState(maxState) {}
+    State basin, len, eden, minState;
+    Stat() : basin(0), len(0), eden(0), minState(maxState) {}
   };
   map<State, Stat> cyc;
 
-  // how big is the basin of attraction?
-  iterState([&](int s) { cyc[c[s]].basin++; });
+  // How big is the basin of attraction?
+  // How many garden of eden states does it contain?
+  iterState([&](int s) {
+      cyc[c[s]].basin++;
+      if (!reachable[s])
+	cyc[c[s]].eden++;
+    });
 
   // how long is the cycle, what is the actual minimal state
   for (auto i : cyc) {
@@ -119,6 +128,7 @@ void cycleStat(Trans &c) {
 	 << i.second.size() << "\t"
 	 << s.len << "\t"
 	 << s.basin << "\t"
+	 << s.eden  << "\t"
 	 << i.first << endl;
   }
 }
@@ -155,9 +165,10 @@ int main(int argc, char **argv) {
   if (!strcmp(argv[2], "cycle")) {
     Trans *t = mmalloc<Trans>(),
       *c = mmalloc<Trans>();
+    bitset<maxState> *r = new bitset<maxState>(0);
     init(*t);
-    findCycle(*t, *c);
-    cycleStat(*c);
+    findCycle(*t, *c, *r);
+    cycleStat(*c, *r);
   }
   return 0;
 }
-- 
cgit v0.10.1