diff --git a/UsageExamples.ipynb b/Demo.ipynb
similarity index 86%
rename from UsageExamples.ipynb
rename to Demo.ipynb
index 0f17115..f2318e9 100644
--- a/UsageExamples.ipynb
+++ b/Demo.ipynb
@@ -18,12 +18,20 @@
    "cell_type": "code",
    "execution_count": 1,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0000000.000 W Numba unavailable. Falling back to pure Python.\n"
+     ]
+    }
+   ],
    "source": [
     "from kyupy import bench\n",
     "\n",
-    "# parse a file\n",
-    "b01 = bench.parse('tests/b01.bench')\n",
+    "# load a file\n",
+    "b01 = bench.load('tests/b01.bench')\n",
     "\n",
     "# ... or specify the circuit as string \n",
     "mycircuit = bench.parse('input(a,b) output(o1,o2,o3) x=buf(a) o1=not(x) o2=buf(x) o3=buf(x)')"
@@ -44,7 +52,7 @@
     {
      "data": {
       "text/plain": [
-       "<Circuit 'tests/b01' with 92 nodes, 130 lines, 4 ports>"
+       "<Circuit 'tests/b01.bench' with 92 nodes, 130 lines, 4 ports>"
       ]
      },
      "execution_count": 2,
@@ -373,7 +381,7 @@
    "source": [
     "from kyupy import verilog\n",
     "\n",
-    "b14 = verilog.parse('tests/b14.v.gz')\n",
+    "b14 = verilog.load('tests/b14.v.gz')\n",
     "b14"
    ]
   },
@@ -456,7 +464,7 @@
    "source": [
     "from kyupy import verilog\n",
     "\n",
-    "b14 = verilog.parse('tests/b14.v.gz')\n",
+    "b14 = verilog.load('tests/b14.v.gz')\n",
     "b14"
    ]
   },
@@ -567,11 +575,11 @@
    "outputs": [],
    "source": [
     "from kyupy import verilog, stil\n",
+    "from kyupy.logic import MVArray, BPArray\n",
     "from kyupy.logic_sim import LogicSim\n",
-    "from kyupy.packed_vectors import PackedVectors\n",
     "\n",
-    "b14 = verilog.parse('tests/b14.v.gz')\n",
-    "s = stil.parse('tests/b14.stuck.stil.gz')\n",
+    "b14 = verilog.load('tests/b14.v.gz')\n",
+    "s = stil.load('tests/b14.stuck.stil.gz')\n",
     "stuck_tests = s.tests(b14)\n",
     "stuck_responses = s.responses(b14)"
    ]
@@ -580,7 +588,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Tests and responses are instances of `PackedVectors`. Its length is the number of test vectors stored (`nvectors`), its `width` is the number of values in a vector, and its `vdim` is the number of bits used for storing one value. By default, the stil parser returns 4-valued test vectors (`vdim=2`)."
+    "Tests and responses are instances of `MVArray`. Its `length` is the number of test vectors stored, its `width` is the number of values in a vector. By default, the stil parser returns 8-valued test vectors (`m=8`)."
    ]
   },
   {
@@ -591,7 +599,7 @@
     {
      "data": {
       "text/plain": [
-       "<PackedVectors nvectors=1081, width=306, vdim=2>"
+       "<MVArray length=1081 width=306 m=8 nbytes=330786>"
       ]
      },
      "execution_count": 19,
@@ -607,7 +615,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The data is stored in a bit-parallel fashion. This internal storage (an `ndarray` of `uint8`) is accessible via `bits`. The first axis is the width, the second axis is `vdim`, the last axis goes along the test set. This last axis is about `nvectors / 8` in length. "
+    "The internal storage (an `ndarray` of `uint8`) is accessible via `data`. The first axis is the width, and the last axis goes along the test set."
    ]
   },
   {
@@ -618,7 +626,7 @@
     {
      "data": {
       "text/plain": [
-       "(306, 2, 136)"
+       "(306, 1081)"
       ]
      },
      "execution_count": 20,
@@ -627,14 +635,14 @@
     }
    ],
    "source": [
-    "stuck_tests.bits.shape"
+    "stuck_tests.data.shape"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The subscript accessor returns a string representation of the given test vector number. Possible values are '0', '1', '-', and 'X'."
+    "The subscript accessor returns a string representation of the given test vector number. Possible values are '0', '1', '-', 'X', 'R', 'F', 'P', and 'N'."
    ]
   },
   {
@@ -645,7 +653,7 @@
     {
      "data": {
       "text/plain": [
-       "'-0--------------------11011111011001100111010101011101----------------------------------00-10111011010110011101110010111010111011101100010000110101111111011010101001010101010101010101001010110101001010101010101010110100000111111111111111011010100100101010010010101101010101001010100111010001010010000011100'"
+       "'P0--------------------11011111011001100111010101011101----------------------------------00-10111011010110011101110010111010111011101100010000110101111111011010101001010101010101010101001010110101001010101010101010110100000111111111111111011010100100101010010010101101010101001010100111010001010010000011100'"
       ]
      },
      "execution_count": 21,
@@ -681,25 +689,80 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The order of values in the vectors correspond to the circuit's interface followed by the scan flip-flops as they appear in `b14.cells`. The test data can be used directly in the simulators as they use the same ordering convention. The following code performs a 4-valued logic simulation and stores the results in a new instance of `PackedVectors`."
+    "The order of values in the vectors correspond to the circuit's interface followed by the scan flip-flops as they appear in `b14.cells`.\n",
+    "The test data can be used directly in the simulators as they use the same ordering convention.\n",
+    "\n",
+    "The logic simulator uses bit-parallel storage of logic values, but our loaded test data uses one `uint8` per logic value.\n",
+    "To convert the storage layout, we instanciate a `BPArray` for the input stimuli.\n",
+    "The storage layout is more compact, but individual values cannot be easily accessed anymore."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 23,
    "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<BPArray length=1081 width=306 m=8 bytes=124848>"
+      ]
+     },
+     "execution_count": 23,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "stuck_tests_bp = BPArray(stuck_tests)\n",
+    "stuck_tests_bp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(306, 3, 136)"
+      ]
+     },
+     "execution_count": 24,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "stuck_tests_bp.data.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The following code performs a 8-valued logic simulation and stores the results in a new instance of `BPArray`.\n",
+    "The packed array is unpacked into an `MVArray` for value access."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
    "outputs": [],
    "source": [
-    "responses = PackedVectors(len(stuck_tests), stuck_tests.width, 2)\n",
-    "simulator = LogicSim(b14, len(responses), 2)\n",
-    "simulator.assign(stuck_tests)\n",
+    "responses_bp = BPArray((stuck_tests_bp.width, len(stuck_tests_bp)))\n",
+    "simulator = LogicSim(b14, sims=len(stuck_tests_bp))\n",
+    "simulator.assign(stuck_tests_bp)\n",
     "simulator.propagate()\n",
-    "simulator.capture(responses)"
+    "simulator.capture(responses_bp)\n",
+    "responses = MVArray(responses_bp)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 26,
    "metadata": {},
    "outputs": [
     {
@@ -708,7 +771,7 @@
        "'--10000010010100010111--------------------------------0101010010101010110101001001010100--011111110011011111000111010101010111011101100010000110101111111011010101001010101010101010101001010110101001010101010101010110100000111111111111111011010100100101010010010101101010101001010101000111111111111111011101'"
       ]
      },
-     "execution_count": 24,
+     "execution_count": 26,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -726,7 +789,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 27,
    "metadata": {},
    "outputs": [
     {
@@ -752,39 +815,32 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Transition faults require test vector pairs for testing. These pairs are generated by `tests8v`, assuming a launch-on-capture scheme (two functional clock cycles after scan-in)."
+    "Transition faults require test vector pairs for testing. These pairs are generated by `tests_loc`, assuming a launch-on-capture scheme (two functional clock cycles after scan-in)."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [],
    "source": [
-    "s = stil.parse('tests/b14.transition.stil.gz')\n",
-    "trans_tests = s.tests8v(b14)\n",
+    "s = stil.load('tests/b14.transition.stil.gz')\n",
+    "trans_tests = s.tests_loc(b14)\n",
     "trans_responses = s.responses(b14)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The returned test data is now 8-valued (`vdim=3`)"
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 29,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "<PackedVectors nvectors=1392, width=306, vdim=3>"
+       "<MVArray length=1392 width=306 m=8 nbytes=425952>"
       ]
      },
-     "execution_count": 27,
+     "execution_count": 29,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -802,16 +858,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "'-0--------------------RRRRRRFRRRRRRRRRRRFFRFRRRRRRRRRR----------------------------------00-00000001110100011111011010000000000000000011001001100101111110101110110001000100010100110111111101101000000111110011100010111000111R1111111111111111111111110001100100000110100000111010101110RFF00F000F0F00F00000FF01F'"
+       "'00--------------------RRRRRRFRRRRRRRRRRRFFRFRRRRRRRRRR----------------------------------00-00000001110100011111011010000000000000000011001001100101111110101110110001000100010100110111111101101000000111110011100010111000111R1111111111111111111111110001100100000110100000111010101110RFF00F000F0F00F00000FF01F'"
       ]
      },
-     "execution_count": 28,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -829,20 +885,22 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 31,
    "metadata": {},
    "outputs": [],
    "source": [
-    "responses = PackedVectors(len(trans_tests), trans_tests.width, 3)\n",
-    "simulator = LogicSim(b14, len(responses), 3)\n",
-    "simulator.assign(trans_tests)\n",
+    "trans_tests_bp = BPArray(trans_tests)\n",
+    "responses_bp = BPArray((trans_tests_bp.width, len(trans_tests_bp)))\n",
+    "simulator = LogicSim(b14, sims=len(trans_tests_bp))\n",
+    "simulator.assign(trans_tests_bp)\n",
     "simulator.propagate()\n",
-    "simulator.capture(responses)"
+    "simulator.capture(responses_bp)\n",
+    "responses = MVArray(responses_bp)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 32,
    "metadata": {},
    "outputs": [
     {
@@ -851,7 +909,7 @@
        "'--F00000F00F0F000F00FF--------------------------------01110101011100000101100000100110R0--0RRRRRRRNNNRNRPRNNNNNRFFRFRRRRRRR000000000011001001100101111110101110110001000100010100110111111101101000000111110011100010111000NNNNNNNNNNNNNNNNNNNNNNNNNNNNP0011001000001101000001110101011101RRRRRRRRRRRRRRRRRRRRP01R'"
       ]
      },
-     "execution_count": 30,
+     "execution_count": 32,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -869,7 +927,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 33,
    "metadata": {},
    "outputs": [
     {
@@ -907,14 +965,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 34,
    "metadata": {},
    "outputs": [],
    "source": [
     "from kyupy import sdf\n",
     "from kyupy.saed import pin_index\n",
     "\n",
-    "df = sdf.parse('tests/b14.sdf.gz')\n",
+    "df = sdf.load('tests/b14.sdf.gz')\n",
     "lt = df.annotation(b14, pin_index, dataset=0, interconnect=False)"
    ]
   },
@@ -927,7 +985,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 35,
    "metadata": {},
    "outputs": [
     {
@@ -936,7 +994,7 @@
        "(46891, 2, 2)"
       ]
      },
-     "execution_count": 33,
+     "execution_count": 35,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -954,7 +1012,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 36,
    "metadata": {},
    "outputs": [
     {
@@ -963,7 +1021,7 @@
        "119676"
       ]
      },
-     "execution_count": 34,
+     "execution_count": 36,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -991,11 +1049,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": 37,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from kyupy.wave_sim_cuda import WaveSimCuda, TMAX\n",
+    "from kyupy.wave_sim import WaveSimCuda, TMAX\n",
     "import numpy as np\n",
     "\n",
     "wsim = WaveSimCuda(b14, lt, sims=32, wavecaps=16)"
@@ -1010,7 +1068,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 38,
    "metadata": {},
    "outputs": [
     {
@@ -1043,23 +1101,23 @@
    "metadata": {},
    "source": [
     "This is a typical simulation loop where the number of patterns is larger than the number of simulators available.\n",
-    "We simulate `trans_tests`.\n",
-    "The timing simulator accepts 4-valued and 8-valued `PackedVectors`, but it will return response (capture) data in a different format."
+    "We simulate `trans_tests_bp`.\n",
+    "The timing simulator accepts 8-valued `BPArray`s, but it will return response (capture) data in a different format."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 39,
    "metadata": {},
    "outputs": [],
    "source": [
-    "nvectors = 128  # len(trans_tests)  # Feel free to simulate all tests if CUDA is set up correctly.\n",
+    "sims = 128  # len(trans_tests_bp)  # Feel free to simulate all tests if CUDA is set up correctly.\n",
     "\n",
-    "cdata = np.zeros((len(wsim.interface), nvectors, 7))  # space to store all capture data\n",
+    "cdata = np.zeros((len(wsim.interface), sims, 7))  # space to store all capture data\n",
     "\n",
-    "for offset in range(0, nvectors, wsim.sims):\n",
-    "    wsim.assign(trans_tests, offset=offset)\n",
-    "    wsim.propagate(sims=nvectors-offset)\n",
+    "for offset in range(0, sims, wsim.sims):\n",
+    "    wsim.assign(trans_tests_bp, offset=offset)\n",
+    "    wsim.propagate(sims=sims-offset)\n",
     "    wsim.capture(time=2.5, cdata=cdata, offset=offset)  # capture at time 2.5"
    ]
   },
@@ -1079,7 +1137,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 40,
    "metadata": {},
    "outputs": [
     {
@@ -1088,7 +1146,7 @@
        "(306, 128, 7)"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 40,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1106,7 +1164,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 41,
    "metadata": {},
    "outputs": [
     {
@@ -1139,7 +1197,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 40,
+   "execution_count": 42,
    "metadata": {},
    "outputs": [
     {
@@ -1148,7 +1206,7 @@
        "2.0610005855560303"
       ]
      },
-     "execution_count": 40,
+     "execution_count": 42,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1166,7 +1224,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 41,
+   "execution_count": 43,
    "metadata": {},
    "outputs": [
     {
@@ -1175,7 +1233,7 @@
        "0.0"
       ]
      },
-     "execution_count": 41,
+     "execution_count": 43,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1193,7 +1251,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 42,
+   "execution_count": 44,
    "metadata": {},
    "outputs": [
     {
@@ -1202,7 +1260,7 @@
        "0.0"
       ]
      },
-     "execution_count": 42,
+     "execution_count": 44,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1222,39 +1280,26 @@
     "If there is an error related to `nvvm`, you probably need to set up some environment variables:\n",
     "```\n",
     "%env LD_LIBRARY_PATH=/usr/local/cuda/lib64\n",
-    "%env NUMBAPRO_NVVM=/usr/local/cuda/nvvm/lib64/libnvvm.so\n",
-    "%env NUMBAPRO_LIBDEVICE=/usr/local/cuda/nvvm/libdevice\n",
+    "%env CUDA_HOME=/usr/local/cuda\n",
     "```\n",
     "If problems persist, refer to documentations for numba and cuda. "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 45,
    "metadata": {},
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Found 1 CUDA devices\n",
-      "id 0              b'TITAN V'                              [SUPPORTED]\n",
-      "                      compute capability: 7.0\n",
-      "                           pci device id: 0\n",
-      "                              pci bus id: 2\n",
-      "Summary:\n",
-      "\t1/1 devices are supported\n"
+     "ename": "ModuleNotFoundError",
+     "evalue": "No module named 'numba'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-45-ffb50e0958fc>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcuda\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mcuda\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdetect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'numba'"
      ]
-    },
-    {
-     "data": {
-      "text/plain": [
-       "True"
-      ]
-     },
-     "execution_count": 43,
-     "metadata": {},
-     "output_type": "execute_result"
     }
    ],
    "source": [
@@ -1287,7 +1332,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.9.1"
   }
  },
  "nbformat": 4,
diff --git a/README.rst b/README.rst
index 1291483..739e7b8 100644
--- a/README.rst
+++ b/README.rst
@@ -21,7 +21,7 @@ Although optional, `numba`_ should be installed for best performance.
 GPU/CUDA support in numba may `require some additional setup <https://numba.pydata.org/numba-doc/latest/cuda/index.html>`_.
 If numba is not available, KyuPy will automatically fall back to slow, pure Python execution.
 
-The Jupyter Notebook `UsageExamples.ipynb <https://github.com/s-holst/kyupy/blob/main/UsageExamples.ipynb>`_ contains some useful examples to get familiar with the API.
+The Jupyter Notebook `Demo.ipynb <https://github.com/s-holst/kyupy/blob/main/Demo.ipynb>`_ contains some useful examples to get familiar with the API.
 
 To work with the latest pre-release source code, clone the `KyuPy GitHub repository <https://github.com/s-holst/kyupy>`_.
 Run ``pip3 install --user -e .`` within your local checkout to make the package available in your Python environment.
diff --git a/setup.py b/setup.py
index fb7aade..9a0bb1b 100644
--- a/setup.py
+++ b/setup.py
@@ -15,7 +15,7 @@ setup(
     author_email='mail@s-holst.de',
     python_requires='>=3.6',
     install_requires=[
-        'numpy>=1.15.0',
+        'numpy>=1.17.0',
         'lark-parser>=0.8.0'
     ],
     extras_requires={
diff --git a/src/kyupy/logic.py b/src/kyupy/logic.py
index f12b293..97446a5 100644
--- a/src/kyupy/logic.py
+++ b/src/kyupy/logic.py
@@ -181,6 +181,21 @@ def mv_xor(x1, x2, out=None):
     return out
 
 
+def mv_transition(init, final, out=None):
+    m = mv_getm(init, final)
+    init, final = mv_cast(init, final, m=m)
+    init = init.data
+    final = final.data
+    out = out or MVArray(np.broadcast(init, final).shape, m=8)
+    out.data[...] = (init & 0b010) | (final & 0b001)
+    out.data[...] |= ((out.data << 1) ^ (out.data << 2)) & 0b100
+    unknown = (init == UNKNOWN) | (init == UNASSIGNED) | (final == UNKNOWN) | (final == UNASSIGNED)
+    unassigned = (init == UNASSIGNED) & (final == UNASSIGNED)
+    np.putmask(out.data, unknown, UNKNOWN)
+    np.putmask(out.data, unassigned, UNASSIGNED)
+    return out
+
+
 class MVArray:
     """An n-dimensional array of m-valued logic values.
 
@@ -188,10 +203,10 @@ class MVArray:
     interpreting 2-valued, 4-valued, and 8-valued logic values.
     Each logic value is stored as an uint8, value manipulations are cheaper than in BPArray.
 
-    An MVArray always has 2 or more axes:
+    An MVArray always has 2 axes:
 
-    * Second-last axis is PI/PO/FF position, the length of this axis is called "width".
-    * Last axis is vector/pattern, the length of this axis is called "length".
+    * Axis 0 is PI/PO/FF position, the length of this axis is called "width".
+    * Axis 1 is vector/pattern, the length of this axis is called "length".
 
     """
 
@@ -203,18 +218,22 @@ class MVArray:
         if isinstance(a, MVArray):
             self.data = a.data.copy()
             self.m = m or a.m
+        elif hasattr(a, 'data'):  # assume it is a BPArray. Can't use isinstance() because BPArray isn't declared yet.
+            self.data = np.zeros((a.width, a.length), dtype=np.uint8)
+            self.m = m or a.m
+            for i in range(a.data.shape[-2]):
+                self.data[...] <<= 1
+                self.data[...] |= np.unpackbits(a.data[..., -i-1, :], axis=1)[:, :a.length]
+            if a.data.shape[-2] == 1:
+                self.data *= 3
         elif isinstance(a, int):
             self.data = np.full((a, 1), UNASSIGNED, dtype=np.uint8)
         elif isinstance(a, tuple):
             self.data = np.full(a, UNASSIGNED, dtype=np.uint8)
         else:
-            if isinstance(a, str):
-                a = [a]
+            if isinstance(a, str): a = [a]
             self.data = np.asarray(interpret(a), dtype=np.uint8)
-            if self.data.ndim == 1:
-                self.data = self.data[:, np.newaxis]
-            else:
-                self.data = np.moveaxis(self.data, -2, -1)
+            self.data = self.data[:, np.newaxis] if self.data.ndim == 1 else np.moveaxis(self.data, -2, -1)
 
         # Cast data to m-valued logic.
         if self.m == 2:
@@ -231,8 +250,166 @@ class MVArray:
         return f'<MVArray length={self.length} width={self.width} m={self.m} nbytes={self.data.nbytes}>'
 
     def __str__(self):
+        return str([self[idx] for idx in range(self.length)])
+
+    def __getitem__(self, vector_idx):
         chars = ["0", "X", "-", "1", "P", "R", "F", "N"]
-        return str([''.join(chars[v] for v in self.data[:, idx]) for idx in range(self.length)])
+        return ''.join(chars[v] for v in self.data[:, vector_idx])
+
+    def __len__(self):
+        return self.length
+
+
+def bp_buf(out, inp):
+    md = out.shape[-2]
+    assert md == inp.shape[-2]
+    if md > 1:
+        unknown = inp[..., 0, :] ^ inp[..., 1, :]
+        if md > 2: unknown &= ~inp[..., 2, :]
+        out[..., 0, :] = inp[..., 0, :] | unknown
+        out[..., 1, :] = inp[..., 1, :] & ~unknown
+        if md > 2: out[..., 2, :] = inp[..., 2, :] & ~unknown
+    else:
+        out[..., 0, :] = inp[..., 0, :]
+
+
+def bp_not(out, inp):
+    md = out.shape[-2]
+    assert md == inp.shape[-2]
+    if md > 1:
+        unknown = inp[..., 0, :] ^ inp[..., 1, :]
+        if md > 2: unknown &= ~inp[..., 2, :]
+        out[..., 0, :] = ~inp[..., 0, :] | unknown
+        out[..., 1, :] = ~inp[..., 1, :] & ~unknown
+        if md > 2: out[..., 2, :] = inp[..., 2, :] & ~unknown
+    else:
+        out[..., 0, :] = ~inp[..., 0, :]
+
+
+def bp_or(out, *ins):
+    # 4-valued:
+    #    0 1 - X
+    # 0  0 1 X X
+    # 1  1 1 1 1
+    # -  X 1 X X
+    # X  X 1 X X
+
+    # 8-valued:           o[0]:            o[1]:            o[2]:
+    #    0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N
+    # 0  0 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  0 0 0 0 1 1 1 1
+    # 1  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0
+    # -  X 1 X X X X X X  1 1 1 1 1 1 1 1  1 0 1 1 1 1 1 1  0 0 0 0 0 0 0 0
+    # X  X 1 X X X X X X  1 1 1 1 1 1 1 1  1 0 1 1 1 1 1 1  0 0 0 0 0 0 0 0
+    # R  R 1 X X R N R R  0 1 1 1 0 1 0 0  0 0 1 1 0 0 0 0  1 0 0 0 1 1 1 1
+    # F  F 1 X X N F F F  1 1 1 1 1 1 1 1  1 0 1 1 0 1 1 1  1 0 0 0 1 1 1 1
+    # P  P 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  1 0 0 0 1 1 1 1
+    # N  N 1 X X R F N N  1 1 1 1 0 1 1 1  0 0 1 1 0 1 0 0  1 0 0 0 1 1 1 1
+    md = out.shape[-2]
+    for inp in ins: assert md == inp.shape[-2]
+    out[...] = 0
+    if md == 1:
+        for inp in ins: out[..., 0, :] |= inp[..., 0, :]
+    elif md == 2:
+        any_unknown = ins[0][..., 0, :] ^ ins[0][..., 1, :]
+        for inp in ins[1:]: any_unknown |= inp[..., 0, :] ^ inp[..., 1, :]
+        any_one = ins[0][..., 0, :] & ins[0][..., 1, :]
+        for inp in ins[1:]: any_one |= inp[..., 0, :] & inp[..., 1, :]
+        for inp in ins:
+            out[..., 0, :] |= inp[..., 0, :] | any_unknown
+            out[..., 1, :] |= inp[..., 1, :] & (~any_unknown | any_one)
+    else:
+        any_unknown = (ins[0][..., 0, :] ^ ins[0][..., 1, :]) & ~ins[0][..., 2, :]
+        for inp in ins[1:]: any_unknown |= (inp[..., 0, :] ^ inp[..., 1, :]) & ~inp[..., 2, :]
+        any_one = ins[0][..., 0, :] & ins[0][..., 1, :] & ~ins[0][..., 2, :]
+        for inp in ins[1:]: any_one |= inp[..., 0, :] & inp[..., 1, :] & ~inp[..., 2, :]
+        for inp in ins:
+            out[..., 0, :] |= inp[..., 0, :] | any_unknown
+            out[..., 1, :] |= inp[..., 1, :] & (~any_unknown | any_one)
+            out[..., 2, :] |= inp[..., 2, :] & (~any_unknown | any_one) & ~any_one
+
+
+def bp_and(out, *ins):
+    # 4-valued:
+    #    0 1 - X
+    # 0  0 0 0 0
+    # 1  0 1 X X
+    # -  0 X X X
+    # X  0 X X X
+
+    # 8-valued:
+    #    0 1 - X R F P N
+    # 0  0 0 0 0 0 0 0 0
+    # 1  0 1 X X R F P N
+    # -  0 X X X X X X X
+    # X  0 X X X X X X X
+    # R  0 R X X R R P R
+    # F  0 F X X R F P F
+    # P  0 P X X P P P P
+    # N  0 N X X R F P N
+    md = out.shape[-2]
+    for inp in ins: assert md == inp.shape[-2]
+    out[...] = 0xff
+    if md == 1:
+        for inp in ins: out[..., 0, :] &= inp[..., 0, :]
+    elif md == 2:
+        any_unknown = ins[0][..., 0, :] ^ ins[0][..., 1, :]
+        for inp in ins[1:]: any_unknown |= inp[..., 0, :] ^ inp[..., 1, :]
+        any_zero = ~ins[0][..., 0, :] & ~ins[0][..., 1, :]
+        for inp in ins[1:]: any_zero |= ~inp[..., 0, :] & ~inp[..., 1, :]
+        for inp in ins:
+            out[..., 0, :] &= inp[..., 0, :] | (any_unknown & ~any_zero)
+            out[..., 1, :] &= inp[..., 1, :] & ~any_unknown
+    else:
+        any_unknown = (ins[0][..., 0, :] ^ ins[0][..., 1, :]) & ~ins[0][..., 2, :]
+        for inp in ins[1:]: any_unknown |= (inp[..., 0, :] ^ inp[..., 1, :]) & ~inp[..., 2, :]
+        any_zero = ~ins[0][..., 0, :] & ~ins[0][..., 1, :] & ~ins[0][..., 2, :]
+        for inp in ins[1:]: any_zero |= ~inp[..., 0, :] & ~inp[..., 1, :] & ~inp[..., 2, :]
+        out[..., 2, :] = 0
+        for inp in ins:
+            out[..., 0, :] &= inp[..., 0, :] | (any_unknown & ~any_zero)
+            out[..., 1, :] &= inp[..., 1, :] & ~any_unknown
+            out[..., 2, :] |= inp[..., 2, :] & (~any_unknown | any_zero) & ~any_zero
+
+
+def bp_xor(out, *ins):
+    # 4-valued:
+    #    0 1 - X
+    # 0  0 1 X X
+    # 1  1 0 X X
+    # -  X X X X
+    # X  X X X X
+
+    # 8-valued:           o[0]:            o[1]:            o[2]:
+    #    0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N
+    # 0  0 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  0 0 0 0 1 1 1 1
+    # 1  1 0 X X F R N P  1 0 1 1 1 0 1 0  0 1 1 1 1 0 0 1  0 0 0 0 1 1 1 1
+    # -  X X X X X X X X  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
+    # X  X X X X X X X X  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
+    # R  R F X X P N R F  0 1 1 1 0 1 0 1  0 1 1 1 1 0 0 1  1 1 0 0 1 1 1 1
+    # F  F R X X N P F R  1 0 1 1 1 0 1 0  1 0 1 1 0 1 1 0  1 1 0 0 1 1 1 1
+    # P  P N X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  1 1 0 0 1 1 1 1
+    # N  N P X X F R N P  1 0 1 1 1 0 1 0  0 1 1 1 1 0 0 1  1 1 0 0 1 1 1 1
+    md = out.shape[-2]
+    for inp in ins: assert md == inp.shape[-2]
+    out[...] = 0
+    if md == 1:
+        for inp in ins: out[..., 0, :] ^= inp[..., 0, :]
+    elif md == 2:
+        any_unknown = ins[0][..., 0, :] ^ ins[0][..., 1, :]
+        for inp in ins[1:]: any_unknown |= inp[..., 0, :] ^ inp[..., 1, :]
+        for inp in ins: out[...] ^= inp
+        out[..., 0, :] |= any_unknown
+        out[..., 1, :] &= ~any_unknown
+    else:
+        any_unknown = (ins[0][..., 0, :] ^ ins[0][..., 1, :]) & ~ins[0][..., 2, :]
+        for inp in ins[1:]: any_unknown |= (inp[..., 0, :] ^ inp[..., 1, :]) & ~inp[..., 2, :]
+        for inp in ins:
+            out[..., 0, :] ^= inp[..., 0, :]
+            out[..., 1, :] ^= inp[..., 1, :]
+            out[..., 2, :] |= inp[..., 2, :]
+        out[..., 0, :] |= any_unknown
+        out[..., 1, :] &= ~any_unknown
+        out[..., 2, :] &= ~any_unknown
 
 
 class BPArray:
@@ -240,9 +417,9 @@ class BPArray:
 
     The primary use of this format is in aiding efficient bit-parallel logic simulation.
     The secondary benefit over MVArray is its memory efficiency.
-    Direct value manipulations are more expensive than with MVArray.
-    It is advised to first construct a MVArray, pack it into a BPArray for simulation and unpack the results
-    back into a MVArray for value access.
+    Accessing individual values is more expensive than with :py:class:`MVArray`.
+    It is advised to first construct a MVArray, pack it into a :py:class:`BPArray` for simulation and unpack the results
+    back into a :py:class:`MVArray` for value access.
 
     The values along the last axis (vectors/patterns) are packed into uint8 words.
     The second-last axis has length ceil(log2(m)) for storing all bits.
@@ -252,6 +429,7 @@ class BPArray:
     def __init__(self, a, m=None):
         if not isinstance(a, MVArray) and not isinstance(a, BPArray):
             a = MVArray(a, m)
+            self.m = a.m
         if isinstance(a, MVArray):
             if m is not None and m != a.m:
                 a = MVArray(a, m)  # cast data
@@ -264,8 +442,12 @@ class BPArray:
                 self.data[..., i, :] = np.packbits((a.data >> i) & 1, axis=-1)
         else:  # we have a BPArray
             self.data = a.data.copy()  # TODO: support conversion to different m
+            self.m = a.m
         self.length = a.length
         self.width = a.width
 
     def __repr__(self):
         return f'<BPArray length={self.length} width={self.width} m={self.m} bytes={self.data.nbytes}>'
+
+    def __len__(self):
+        return self.length
diff --git a/src/kyupy/logic_sim.py b/src/kyupy/logic_sim.py
index e595bc1..cddde47 100644
--- a/src/kyupy/logic_sim.py
+++ b/src/kyupy/logic_sim.py
@@ -1,22 +1,25 @@
+import math
+
 import numpy as np
 
-from . import packed_vectors
+from . import logic
 
 
 class LogicSim:
     """A bit-parallel naïve combinational simulator for 2-, 4-, or 8-valued logic.
     """
-    def __init__(self, circuit, nvectors=1, vdim=1):
+    def __init__(self, circuit, sims=1, m=8):
+        assert m in [2, 4, 8]
+        self.m = m
+        mdim = math.ceil(math.log2(m))
         self.circuit = circuit
-        self.nvectors = nvectors
-        nbytes = (nvectors - 1) // 8 + 1
+        self.sims = sims
+        nbytes = (sims - 1) // 8 + 1
         self.interface = list(circuit.interface) + [n for n in circuit.nodes if 'dff' in n.kind.lower()]
-        self.state = np.zeros((len(circuit.lines), vdim, nbytes), dtype='uint8')
+        self.state = np.zeros((len(circuit.lines), mdim, nbytes), dtype='uint8')
         self.state_epoch = np.zeros(len(circuit.nodes), dtype='int8') - 1
-        self.tmp = np.zeros((5, vdim, nbytes), dtype='uint8')
-        self.zero = np.zeros((vdim, nbytes), dtype='uint8')
-        if vdim > 1:
-            self.zero[1] = 255
+        self.tmp = np.zeros((5, mdim, nbytes), dtype='uint8')
+        self.zero = np.zeros((mdim, nbytes), dtype='uint8')
         self.epoch = 0
 
         self.fork_vd1 = self.fork_vdx
@@ -46,7 +49,7 @@ class LogicSim:
         self.nbuff_vd3 = self.fork_vd3
         self.xor2_vd3 = self.xor_vd3
         
-        known_fct = [(f[:-4], getattr(self, f)) for f in dir(self) if f.endswith(f'_vd{vdim}')]
+        known_fct = [(f[:-4], getattr(self, f)) for f in dir(self) if f.endswith(f'_vd{mdim}')]
         self.node_fct = []
         for n in circuit.nodes:
             t = n.kind.lower().replace('__fork__', 'fork')
@@ -60,9 +63,9 @@ class LogicSim:
 
     def assign(self, stimuli):
         """Assign stimuli to the primary inputs and state-elements (flip-flops)."""
-        if isinstance(stimuli, packed_vectors.PackedVectors):
-            stimuli = stimuli.bits
-        for (stim, node) in zip(stimuli, self.interface):
+        if hasattr(stimuli, 'data'):
+            stimuli = stimuli.data
+        for stim, node in zip(stimuli, self.interface):
             if len(node.outs) == 0: continue
             outputs = [self.state[line.index] if line else self.tmp[3] for line in node.outs]
             self.node_fct[node.index]([stim], outputs)
@@ -80,11 +83,12 @@ class LogicSim:
 
     def capture(self, responses):
         """Capture the current values at the primary outputs and in the state-elements (flip-flops)."""
-        if isinstance(responses, packed_vectors.PackedVectors):
-            responses = responses.bits
-        for (resp, node) in zip(responses, self.interface):
+        if hasattr(responses, 'data'):
+            responses = responses.data
+        for resp, node in zip(responses, self.interface):
             if len(node.ins) == 0: continue
             resp[...] = self.state[node.ins[0].index]
+        # print(responses)
 
     def propagate(self):
         """Propagate the input values towards the outputs (Perform all logic operations in topological order)."""
@@ -98,8 +102,7 @@ class LogicSim:
                 self.state_epoch[line.reader.index] = self.epoch
         self.epoch = (self.epoch + 1) % 128
 
-    @staticmethod
-    def fork_vdx(inputs, outputs):
+    def fork_vdx(self, inputs, outputs):
         for o in outputs: o[...] = inputs[0]
     
     def const0_vdx(self, _, outputs):
@@ -107,40 +110,34 @@ class LogicSim:
 
     # 2-valued simulation
 
-    @staticmethod
-    def not_vd1(inputs, outputs):
+    def not_vd1(self, inputs, outputs):
         outputs[0][0] = ~inputs[0][0]
 
     def const1_vd1(self, _, outputs):
         for o in outputs: o[...] = self.zero
         self.not_vd1(outputs, outputs)
 
-    @staticmethod
-    def and_vd1(inputs, outputs):
+    def and_vd1(self, inputs, outputs):
         o = outputs[0]
         o[0] = inputs[0][0]
         for i in inputs[1:]: o[0] &= i[0]
 
-    @staticmethod
-    def or_vd1(inputs, outputs):
+    def or_vd1(self, inputs, outputs):
         o = outputs[0]
         o[0] = inputs[0][0]
         for i in inputs[1:]: o[0] |= i[0]
 
-    @staticmethod
-    def xor_vd1(inputs, outputs):
+    def xor_vd1(self, inputs, outputs):
         o = outputs[0]
         o[0] = inputs[0][0]
         for i in inputs[1:]: o[0] ^= i[0]
 
-    @staticmethod
-    def sdff_vd1(inputs, outputs):
+    def sdff_vd1(self, inputs, outputs):
         outputs[0][0] = inputs[0][0]
         if len(outputs) > 1:
             outputs[1][0] = ~inputs[0][0]
 
-    @staticmethod
-    def dff_vd1(inputs, outputs):
+    def dff_vd1(self, inputs, outputs):
         outputs[0][0] = inputs[0][0]
         if len(outputs) > 1:
             outputs[1][0] = ~inputs[0][0]
@@ -158,93 +155,26 @@ class LogicSim:
         self.not_vd1(outputs, outputs)
 
     # 4-valued simulation
-    # sym [0] [1] (value, care)
-    #  0   0   1
-    #  1   1   1
-    #  -   0   0
-    #  X   1   0
-
-    @staticmethod
-    def not_vd2(inputs, outputs):
-        # 4-valued not:
-        # i: 0 1 - X
-        # o: 1 0 X X
-        # o0 1 0 1 1
-        # o1 1 1 0 0
-
-        outputs[0][0] = ~inputs[0][0] | ~inputs[0][1]  # value = 0 or DC
-        outputs[0][1] = inputs[0][1]  # care = C
+
+    def not_vd2(self, inputs, outputs):
+        logic.bp_not(outputs[0], inputs[0])
 
     def and_vd2(self, inputs, outputs):
-        # 4-valued:    o[0]:     o[1]:
-        #    0 1 - X   0 1 - X   0 1 - X
-        # 0  0 0 0 0   0 0 0 0   1 1 1 1
-        # 1  0 1 X X   0 1 1 1   1 1 0 0
-        # -  0 X X X   0 1 1 1   1 0 0 0
-        # X  0 X X X   0 1 1 1   1 0 0 0
-
-        i = inputs[0]
-        any0 = self.tmp[0]
-        anyd = self.tmp[1]
-        any0[0] = ~i[0] & i[1]
-        anyd[0] = ~i[1]
-        for i in inputs[1:]:
-            any0[0] |= ~i[0] & i[1]
-            anyd[0] |= ~i[1]
-        o = outputs[0]
-        o[0] = ~any0[0]  # value = no0
-        o[1] = any0[0] | ~anyd[0]  # care = any0 or noDC
+        logic.bp_and(outputs[0], *inputs)
 
     def or_vd2(self, inputs, outputs):
-        # 4-valued:    o[0]:     o[1]:
-        #    0 1 - X   0 1 - X   0 1 - X
-        # 0  0 1 X X   0 1 1 1   1 1 0 0
-        # 1  1 1 1 1   1 1 1 1   1 1 1 1
-        # -  X 1 X X   1 1 1 1   0 1 0 0
-        # X  X 1 X X   1 1 1 1   0 1 0 0
-
-        i = inputs[0]
-        any1 = self.tmp[0]
-        anyd = self.tmp[1]
-        any1[0] = i[0] & i[1]
-        anyd[0] = ~i[1]
-        for i in inputs[1:]:
-            any1[0] |= i[0] & i[1]
-            anyd[0] |= ~i[1]
-        o = outputs[0]
-        o[0] = any1[0] | anyd[0]  # value = any1 or anyDC
-        o[1] = any1[0] | ~anyd[0]  # care = any1 or noDC
+        logic.bp_or(outputs[0], *inputs)
 
     def xor_vd2(self, inputs, outputs):
-        # 4-valued:    o[0]:     o[1]:
-        #    0 1 - X   0 1 - X   0 1 - X
-        # 0  0 1 X X   0 1 1 1   1 1 0 0
-        # 1  1 0 X X   1 0 1 1   1 1 0 0
-        # -  X X X X   1 1 1 1   0 0 0 0
-        # X  X X X X   1 1 1 1   0 0 0 0
-
-        i = inputs[0]
-        odd1 = self.tmp[0]
-        anyd = self.tmp[1]
-        odd1[0] = i[0] & i[1]
-        anyd[0] = ~i[1]
-        for i in inputs[1:]:
-            odd1[0] ^= i[0] & i[1]
-            anyd[0] |= ~i[1]
-        o = outputs[0]
-        o[0] = odd1[0] | anyd[0]  # value = odd1 or anyDC
-        o[1] = ~anyd[0]  # care = noDC
+        logic.bp_xor(outputs[0], *inputs)
 
     def sdff_vd2(self, inputs, outputs):
         self.dff_vd2(inputs, outputs)
         if len(outputs) > 1:
-            outputs[1][0] = ~inputs[0][0] | ~inputs[0][1]  # value = 0 or DC
-            outputs[1][1] = inputs[0][1]  # care = C
+            logic.bp_not(outputs[1], inputs[0])
 
-    @staticmethod
-    def dff_vd2(inputs, outputs):
-        outputs[0][0] = inputs[0][0] | ~inputs[0][1]  # value = 1 or DC
-        outputs[0][1] = inputs[0][1]  # care = C
+    def dff_vd2(self, inputs, outputs):
+        logic.bp_buf(outputs[0], inputs[0])
 
     def nand_vd2(self, inputs, outputs):
         self.and_vd2(inputs, outputs)
@@ -263,149 +193,26 @@ class LogicSim:
         self.not_vd2(outputs, outputs)
 
     # 8-valued simulation
-    # sym [0] [1] [2] (initial value, ~final value, toggles present?)
-    #  0   0   1   0
-    #  1   1   0   0
-    #  -   0   0   0
-    #  X   1   1   0
-    #  R   0   0   1  _/"
-    #  F   1   1   1  "\_
-    #  P   0   1   1  _/\_
-    #  N   1   0   1  "\/"
 
     def not_vd3(self, inputs, outputs):
-        # 8-valued not:
-        # i: 0 1 - X R F P N
-        # i0 0 1 0 1 0 1 0 1
-        # i1 1 0 0 1 0 1 1 0
-        # i2 0 0 0 0 1 1 1 1
-        # o: 1 0 X X F R N P
-        # o0 1 0 1 1 1 0 1 0
-        # o1 0 1 1 1 1 0 0 1
-        # o2 0 0 0 0 1 1 1 1
-        i = inputs[0]
-        dc = self.tmp[0]
-        dc[0] = ~(i[0] ^ i[1]) & ~i[2]
-        dc = self.tmp[0]
-        outputs[0][0] = ~i[0] | dc[0]  # init.v = ~i0 or DC
-        outputs[0][1] = ~i[1] | dc[0]  # init.v = ~i1 or DC
-        outputs[0][2] = i[2]  # toggles = i2
+        logic.bp_not(outputs[0], inputs[0])
 
     def and_vd3(self, inputs, outputs):
-        # 8-valued:           o[0]:            o[1]:            o[2]:
-        #    0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N
-        # 0  0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # 1  0 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  0 0 0 0 1 1 1 1
-        # -  0 X X X X X X X  0 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # X  0 X X X X X X X  0 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # R  0 R X X R R P R  0 0 1 1 0 0 0 0  1 0 1 1 0 0 1 0  0 1 0 0 1 1 1 1
-        # F  0 F X X R F P F  0 1 1 1 0 1 0 1  1 1 1 1 0 1 1 1  0 1 0 0 1 1 1 1
-        # P  0 P X X P P P P  0 0 1 1 0 0 0 0  1 1 1 1 1 1 1 1  0 1 0 0 1 1 1 1
-        # N  0 N X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  0 1 0 0 1 1 1 1
-        i = inputs[0]
-        anyi0 = self.tmp[0]
-        anyf0 = self.tmp[1]
-        anyd = self.tmp[2]
-        any0 = self.tmp[3]
-        any_t = self.tmp[4]
-        anyd[0] = ~(i[0] ^ i[1]) & ~i[2]
-        anyi0[0] = ~i[0] & ~anyd[0]
-        anyf0[0] = i[1] & ~anyd[0]
-        any_t[0] = i[2]
-        any0[0] = anyi0[0] & anyf0[0] & ~i[2]
-        for i in inputs[1:]:
-            dc = ~(i[0] ^ i[1]) & ~i[2]
-            anyd[0] |= dc
-            anyi0[0] |= ~i[0] & ~dc
-            anyf0[0] |= i[1] & ~dc
-            any_t[0] |= i[2]
-            any0[0] |= ~i[0] & ~dc & i[1] & ~i[2]
-        o = outputs[0]
-        o[0] = (~anyi0[0] | anyd[0]) & ~any0[0]  # initial = no_i0 or DC
-        o[1] = anyf0[0] | anyd[0]  # ~final = ~no_f0 or DC
-        o[2] = any_t[0] & ~(anyd[0] | any0[0])  # toggle = anyT and noDC and no0
+        logic.bp_and(outputs[0], *inputs)
 
     def or_vd3(self, inputs, outputs):
-        # 8-valued:           o[0]:            o[1]:            o[2]:
-        #    0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N
-        # 0  0 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  0 0 0 0 1 1 1 1
-        # 1  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0
-        # -  X 1 X X X X X X  1 1 1 1 1 1 1 1  1 0 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # X  X 1 X X X X X X  1 1 1 1 1 1 1 1  1 0 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # R  R 1 X X R N R R  0 1 1 1 0 1 0 0  0 0 1 1 0 0 0 0  1 0 0 0 1 1 1 1
-        # F  F 1 X X N F F F  1 1 1 1 1 1 1 1  1 0 1 1 0 1 1 1  1 0 0 0 1 1 1 1
-        # P  P 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  1 0 0 0 1 1 1 1
-        # N  N 1 X X R F N N  1 1 1 1 0 1 1 1  0 0 1 1 0 1 0 0  1 0 0 0 1 1 1 1
-        i = inputs[0]
-        anyi1 = self.tmp[0]
-        anyf1 = self.tmp[1]
-        anyd = self.tmp[2]
-        any1 = self.tmp[3]
-        any_t = self.tmp[4]
-        anyd[0] = ~(i[0] ^ i[1]) & ~i[2]
-        anyi1[0] = i[0] & ~anyd[0]
-        anyf1[0] = ~i[1] & ~anyd[0]
-        any_t[0] = i[2]
-        any1[0] = (anyi1[0] & anyf1[0]) & ~i[2]
-        for i in inputs[1:]:
-            dc = ~(i[0] ^ i[1]) & ~i[2]
-            anyd[0] |= dc
-            anyi1[0] |= i[0] & ~dc
-            anyf1[0] |= ~i[1] & ~dc
-            any_t[0] |= i[2]
-            any1[0] |= i[0] & ~dc & ~i[1] & ~i[2]
-        o = outputs[0]
-        o[0] = anyi1[0] | anyd[0]  # initial = i1 or DC
-        o[1] = (~anyf1[0] | anyd[0]) & ~any1[0]  # ~final = f1 or DC
-        o[2] = any_t[0] & ~(anyd[0] | any1[0])  # toggle = anyT and no(DC or 1)
+        logic.bp_or(outputs[0], *inputs)
 
     def xor_vd3(self, inputs, outputs):
-        # 8-valued:           o[0]:            o[1]:            o[2]:
-        #    0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N  0 1 - X R F P N
-        # 0  0 1 X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  0 0 0 0 1 1 1 1
-        # 1  1 0 X X F R N P  1 0 1 1 1 0 1 0  0 1 1 1 1 0 0 1  0 0 0 0 1 1 1 1
-        # -  X X X X X X X X  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # X  X X X X X X X X  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  0 0 0 0 0 0 0 0
-        # R  R F X X P N R F  0 1 1 1 0 1 0 1  0 1 1 1 1 0 0 1  1 1 0 0 1 1 1 1
-        # F  F R X X N P F R  1 0 1 1 1 0 1 0  1 0 1 1 0 1 1 0  1 1 0 0 1 1 1 1
-        # P  P N X X R F P N  0 1 1 1 0 1 0 1  1 0 1 1 0 1 1 0  1 1 0 0 1 1 1 1
-        # N  N P X X F R N P  1 0 1 1 1 0 1 0  0 1 1 1 1 0 0 1  1 1 0 0 1 1 1 1
-        i = inputs[0]
-        odd0 = self.tmp[0]
-        odd1 = self.tmp[1]
-        anyd = self.tmp[2]
-        anyt = self.tmp[3]
-        odd0[0] = i[0]
-        odd1[0] = i[1]
-        anyd[0] = ~(i[0] ^ i[1]) & ~i[2]
-        anyt[0] = i[2]
-        for i in inputs[1:]:
-            odd0[0] ^= i[0]
-            odd1[0] ^= i[1]
-            anyd[0] |= ~(i[0] ^ i[1]) & ~i[2]
-            anyt[0] |= i[2]
-        o = outputs[0]
-        o[0] = odd0[0] | anyd[0]
-        o[1] = ~odd1[0] | anyd[0]
-        o[2] = anyt[0] & ~anyd[0]
-        
+        logic.bp_xor(outputs[0], *inputs)
+
     def sdff_vd3(self, inputs, outputs):
         self.dff_vd3(inputs, outputs)
         if len(outputs) > 1:
-            i = inputs[0]
-            dc = self.tmp[0]
-            dc[0] = ~(i[0] ^ i[1]) & ~i[2]
-            outputs[1][0] = ~i[0] | dc[0]  # value = 1 or DC
-            outputs[1][1] = ~i[1] | dc[0]  # value = 1 or DC
-            outputs[1][2] = i[2]  # toggle = T
+            logic.bp_not(outputs[1], inputs[0])
 
     def dff_vd3(self, inputs, outputs):
-        i = inputs[0]
-        dc = self.tmp[0]
-        dc[0] = ~(i[0] ^ i[1]) & ~i[2]
-        outputs[0][0] = i[0] | dc[0]  # value = 1 or DC
-        outputs[0][1] = i[1] | dc[0]  # value = 1 or DC
-        outputs[0][2] = i[2]  # toggle = T
+        logic.bp_buf(outputs[0], inputs[0])
 
     def nand_vd3(self, inputs, outputs):
         self.and_vd3(inputs, outputs)
diff --git a/src/kyupy/packed_vectors.py b/src/kyupy/packed_vectors.py
deleted file mode 100644
index 3581e09..0000000
--- a/src/kyupy/packed_vectors.py
+++ /dev/null
@@ -1,300 +0,0 @@
-import numpy as np
-from . import popcount
-from .logic import bit_in
-
-
-class PackedVectors:
-    def __init__(self, nvectors=8, width=1, vdim=1, from_cache=None):
-        if from_cache is not None:
-            self.bits = np.array(from_cache)
-            self.width, self.vdim, nbytes = self.bits.shape
-        else:
-            self.bits = np.zeros((width, vdim, (nvectors - 1) // 8 + 1), dtype='uint8')
-            self.vdim = vdim
-            self.width = width
-        self.nvectors = nvectors
-        m1 = np.array([2 ** x for x in range(7, -1, -1)], dtype='uint8')
-        m0 = ~m1
-        self.mask = np.rollaxis(np.vstack((m0, m1)), 1)
-
-    @classmethod
-    def from_pair(cls, init, final):
-        assert init.nvectors == final.nvectors
-        assert len(init.bits) == len(final.bits)
-        init_v = init.bits[:, 0]
-        if init.vdim == 3:
-            init_c = (init.bits[:, 0] ^ init.bits[:, 1]) | init.bits[:, 2]
-        elif init.vdim == 2:
-            init_c = init.bits[:, 1]
-        else:
-            init_c = ~np.zeros_like(init.bits[:, 0])
-        final_v = final.bits[:, 0]
-        if final.vdim == 3:
-            final_c = (final.bits[:, 0] ^ final.bits[:, 1]) | final.bits[:, 2]
-            final_v = ~final.bits[:, 1]
-        elif final.vdim == 2:
-            final_c = final.bits[:, 1]
-        else:
-            final_c = ~np.zeros_like(final.bits[:, 0])
-        c = init_c & final_c
-        a0 = init_v & c
-        a1 = ~final_v & c
-        a2 = (init_v ^ final_v) & c
-        p = PackedVectors(init.nvectors, len(init.bits), 3)
-        p.bits[:, 0] = a0
-        p.bits[:, 1] = a1
-        p.bits[:, 2] = a2
-        return p
-        
-    def transition_vectors(self):
-        a = PackedVectors(self.nvectors-1, self.width, 3)
-        for pos in range(self.width):
-            for vidx in range(self.nvectors-1):
-                tr = self.get_value(vidx, pos) + self.get_value(vidx+1, pos)
-                if tr == '00':
-                    a.set_value(vidx, pos, '0')
-                elif tr == '11':
-                    a.set_value(vidx, pos, '1')
-                elif tr == '01':
-                    a.set_value(vidx, pos, 'R')
-                elif tr == '10':
-                    a.set_value(vidx, pos, 'F')
-                elif tr == '--':
-                    a.set_value(vidx, pos, '-')
-                else:
-                    a.set_value(vidx, pos, 'X')
-        return a
-        
-    def __add__(self, other):
-        a = PackedVectors(self.nvectors + other.nvectors, self.width, max(self.vdim, other.vdim))
-        # a.bits[:self.bits.shape[0], 0] = self.bits[:, 0]
-        # if self.vdim == 2:
-        #    a.bits[:self.bits.shape[0], 1] = self.care_bits
-        # elif self.vdim == 3:
-        #    a.bits[:self.bits.shape[0], 1] = ~self.value_bits
-        #    a.bits[:self.bits.shape[0], 2] = self.toggle_bits
-        for i in range(self.nvectors):
-            a[i] = self[i]
-        for i in range(len(other)):
-            a[self.nvectors+i] = other[i]
-        return a
-
-    def __len__(self):
-        return self.nvectors
-    
-    def randomize(self, one_probability=0.5):
-        for data in self.bits:
-            data[0] = np.packbits((np.random.rand(self.nvectors) < one_probability).astype(int))
-            if self.vdim == 2:
-                data[1] = 255
-            elif self.vdim == 3:
-                data[1] = ~np.packbits((np.random.rand(self.nvectors) < one_probability).astype(int))
-                data[2] = data[0] ^ ~data[1]
-            
-    def copy(self, selection_mask=None):
-        if selection_mask is not None:
-            cpy = PackedVectors(popcount(selection_mask), len(self.bits), self.vdim)
-            cur = 0
-            for vidx in range(self.nvectors):
-                if bit_in(selection_mask, vidx):
-                    cpy[cur] = self[vidx]
-                    cur += 1
-        else:
-            cpy = PackedVectors(self.nvectors, len(self.bits), self.vdim)
-            np.copyto(cpy.bits, self.bits)
-        return cpy
-
-    @property
-    def care_bits(self):
-        if self.vdim == 1:
-            return self.bits[:, 0] | 255
-        elif self.vdim == 2:
-            return self.bits[:, 1]
-        elif self.vdim == 3:
-            return (self.bits[:, 0] ^ self.bits[:, 1]) | self.bits[:, 2]
-
-    @property
-    def initial_bits(self):
-        return self.bits[:, 0]
-
-    @property
-    def value_bits(self):
-        if self.vdim == 3:
-            return ~self.bits[:, 1]
-        else:
-            return self.bits[:, 0]
-
-    @property
-    def toggle_bits(self):
-        if self.vdim == 3:
-            return self.bits[:, 2]
-        else:
-            return self.bits[:, 0] & 0
-
-    def get_value(self, vector, position):
-        if vector >= self.nvectors:
-            raise IndexError(f'vector out of range: {vector} >= {self.nvectors}')
-        a = self.bits[position, :, vector // 8]
-        m = self.mask[vector % 8]
-        if self.vdim == 1:
-            return '1' if a[0] & m[1] else '0'
-        elif self.vdim == 2:
-            if a[0] & m[1]:
-                return '1' if a[1] & m[1] else 'X'
-            else:
-                return '0' if a[1] & m[1] else '-'
-        elif self.vdim == 3:
-            if a[2] & m[1]:
-                if a[0] & m[1]:
-                    return 'F' if a[1] & m[1] else 'N'
-                else:
-                    return 'P' if a[1] & m[1] else 'R'
-            else:
-                if a[0] & m[1]:
-                    return 'X' if a[1] & m[1] else '1'
-                else:
-                    return '0' if a[1] & m[1] else '-'                
-
-    def get_values_for_position(self, position):
-        return ''.join(self.get_value(x, position) for x in range(self.nvectors))
-
-    def set_value(self, vector, position, v):
-        if vector >= self.nvectors:
-            raise IndexError(f'vector out of range: {vector} >= {self.nvectors}')
-        a = self.bits[position, :, vector // 8]
-        m = self.mask[vector % 8]
-        if self.vdim == 1:
-            self._set_value_vd1(a, m, v)
-        elif self.vdim == 2:
-            self._set_value_vd2(a, m, v)
-        elif self.vdim == 3:
-            self._set_value_vd3(a, m, v)
-    
-    def set_values(self, vector, v, mapping=None, inversions=None):
-        if vector >= self.nvectors:
-            raise IndexError(f'vector out of range: {vector} >= {self.nvectors}')
-        if not mapping:
-            mapping = [y for y in range(len(v))]
-        if inversions is None:
-            inversions = [False] * len(v)
-        for i, c in enumerate(v):
-            if inversions[i]:
-                if c == '1':
-                    c = '0'
-                elif c == '0':
-                    c = '1'
-                elif c == 'H':
-                    c = 'L'
-                elif c == 'L':
-                    c = 'H'
-                elif c == 'R':
-                    c = 'F'
-                elif c == 'F':
-                    c = 'R'
-            self.set_value(vector, mapping[i], c)
-    
-    def set_values_for_position(self, position, values):
-        for i, v in enumerate(values):
-            self.set_value(i, position, v)
-            
-    def __setitem__(self, vector, value):
-        for i, c in enumerate(value):
-            self.set_value(vector, i, c)
-
-    def __getitem__(self, vector):
-        if isinstance(vector, slice):
-            first = self.get_values_for_position(0)[vector]
-            ret = PackedVectors(len(first), self.width, self.vdim)
-            ret.set_values_for_position(0, first)
-            for pos in range(1, self.width):
-                ret.set_values_for_position(pos, self.get_values_for_position(pos)[vector])
-            return ret
-        return ''.join(self.get_value(vector, pos) for pos in range(len(self.bits)))
-
-    @staticmethod
-    def _set_value_vd1(a, m, v):
-        if v in [True, 1, '1', 'H', 'h']:
-            a[0] |= m[1]
-        else:
-            a[0] &= m[0]
-    
-    @staticmethod
-    def _set_value_vd2(a, m, v):
-        if v in [True, 1, '1', 'H', 'h']:
-            a[0] |= m[1]
-            a[1] |= m[1]
-        elif v in [False, 0, '0', 'L', 'l']:
-            a[0] &= m[0]
-            a[1] |= m[1]
-        elif v in ['X', 'x']:
-            a[0] |= m[1]
-            a[1] &= m[0]
-        else:
-            a[0] &= m[0]
-            a[1] &= m[0]
-
-    #   i fb act
-    # a 0 1 2
-    # - 0 0 0  None, '-'
-    # 0 0 1 0  False, 0, '0', 'l', 'L'
-    # 1 1 0 0  True, 1, '1', 'h', 'H'
-    # X 1 1 0  'x', 'X'
-    # / 0 0 1  '/', 'r', 'R'
-    # ^ 0 1 1  '^', 'p', 'P'
-    # v 1 0 1  'v', 'n', 'N'
-    # \ 1 1 1  '\', 'f', 'F'
-    @staticmethod
-    def _set_value_vd3(a, m, v):
-        if v in [False, 0, '0', 'L', 'l']:
-            a[0] &= m[0]
-            a[1] |= m[1]
-            a[2] &= m[0]
-        elif v in [True, 1, '1', 'H', 'h']:
-            a[0] |= m[1]
-            a[1] &= m[0]
-            a[2] &= m[0]
-        elif v in ['X', 'x']:
-            a[0] |= m[1]
-            a[1] |= m[1]
-            a[2] &= m[0]
-        elif v in ['/', 'r', 'R']:
-            a[0] &= m[0]
-            a[1] &= m[0]
-            a[2] |= m[1]
-        elif v in ['^', 'p', 'P']:
-            a[0] &= m[0]
-            a[1] |= m[1]
-            a[2] |= m[1]
-        elif v in ['v', 'n', 'N']:
-            a[0] |= m[1]
-            a[1] &= m[0]
-            a[2] |= m[1]
-        elif v in ['\\', 'f', 'F']:
-            a[0] |= m[1]
-            a[1] |= m[1]
-            a[2] |= m[1]
-        else:
-            a[0] &= m[0]
-            a[1] &= m[0]
-            a[2] &= m[0]
-                                    
-    def __repr__(self):
-        return f'<PackedVectors nvectors={self.nvectors}, width={self.width}, vdim={self.vdim}>'
-
-    def __str__(self):
-        lst = []
-        for p in range(self.nvectors):
-            lst.append(''.join(self.get_value(p, w) for w in range(len(self.bits))))
-        if len(lst) == 0: return ''
-        if len(lst[0]) > 64:
-            lst = [s[:32] + '...' + s[-32:] for s in lst]
-        if len(lst) <= 16:
-            return '\n'.join(lst)
-        else:
-            return '\n'.join(lst[:8]) + '\n...\n' + '\n'.join(lst[-8:])
-            
-    def diff(self, other, out=None):
-        if out is None:
-            out = np.zeros((self.width, self.bits.shape[-1]), dtype='uint8')
-        out[...] = (self.value_bits ^ other.value_bits) & self.care_bits & other.care_bits
-        return out
diff --git a/src/kyupy/stil.py b/src/kyupy/stil.py
index 6004ea0..5c022ca 100644
--- a/src/kyupy/stil.py
+++ b/src/kyupy/stil.py
@@ -8,13 +8,13 @@ Call :py:func:`StilFile.tests4v`, :py:func:`StilFile.tests8v`, or :py:func:`Stil
 obtain the appropriate vector sets.
 """
 
-from lark import Lark, Transformer
-from collections import namedtuple
 import re
+from collections import namedtuple
 
-from .packed_vectors import PackedVectors
+from lark import Lark, Transformer
+
+from . import readtext, logic
 from .logic_sim import LogicSim
-from . import readtext
 
 
 Call = namedtuple('Call', ['name', 'parameters'])
@@ -86,59 +86,69 @@ class StilFile:
             scan_inversions[chain[-1]] = scan_out_inversion
         return interface, pi_map, po_map, scan_maps, scan_inversions
         
-    def tests4v(self, circuit):
-        """Assembles and returns a scan test pattern set in 4-valued logic for given circuit.
+    def tests(self, circuit):
+        """Assembles and returns a scan test pattern set for given circuit.
 
         This function assumes a static (stuck-at fault) test.
         """
         interface, pi_map, po_map, scan_maps, scan_inversions = self._maps(circuit)
-        tests = PackedVectors(len(self.patterns), len(interface), 2)
+        tests = logic.MVArray((len(interface), len(self.patterns)))
         for i, p in enumerate(self.patterns):
             for si_port in self.si_ports.keys():
-                tests.set_values(i, p.load[si_port], scan_maps[si_port], scan_inversions[si_port])
-            tests.set_values(i, p.launch['_pi'], pi_map)
+                pattern = logic.mv_xor(p.load[si_port], scan_inversions[si_port])
+                tests.data[scan_maps[si_port], i] = pattern.data[:, 0]
+            tests.data[pi_map, i] = logic.MVArray(p.launch['_pi']).data[:, 0]
         return tests
 
-    def tests8v(self, circuit):
-        """Assembles and returns a scan test pattern set in 8-valued logic for given circuit.
+    def tests_loc(self, circuit):
+        """Assembles and returns a LoC scan test pattern set for given circuit.
 
         This function assumes a launch-on-capture (LoC) delay test.
         It performs a logic simulation to obtain the first capture pattern (the one that launches the
         delay test) and assembles the test pattern set from from pairs for initialization- and launch-patterns.
         """
         interface, pi_map, po_map, scan_maps, scan_inversions = self._maps(circuit)
-        init = PackedVectors(len(self.patterns), len(interface), 2)
+        init = logic.MVArray((len(interface), len(self.patterns)), m=4)
+        # init = PackedVectors(len(self.patterns), len(interface), 2)
         for i, p in enumerate(self.patterns):
             # init.set_values(i, '0' * len(interface))
             for si_port in self.si_ports.keys():
-                init.set_values(i, p.load[si_port], scan_maps[si_port], scan_inversions[si_port])
-            init.set_values(i, p.launch['_pi'], pi_map)
-        sim4v = LogicSim(circuit, len(init), 2)
-        sim4v.assign(init)
+                pattern = logic.mv_xor(p.load[si_port], scan_inversions[si_port])
+                init.data[scan_maps[si_port], i] = pattern.data[:, 0]
+            init.data[pi_map, i] = logic.MVArray(p.launch['_pi']).data[:, 0]
+        launch_bp = logic.BPArray(init)
+        sim4v = LogicSim(circuit, len(init), m=4)
+        sim4v.assign(launch_bp)
         sim4v.propagate()
-        launch = init.copy()
-        sim4v.capture(launch)
+        sim4v.capture(launch_bp)
+        launch = logic.MVArray(launch_bp)
         for i, p in enumerate(self.patterns):
             # if there was no launch clock, then init = launch
             if ('P' not in p.launch['_pi']) or ('P' not in p.capture['_pi']):
                 for si_port in self.si_ports.keys():
-                    launch.set_values(i, p.load[si_port], scan_maps[si_port], scan_inversions[si_port])
+                    pattern = logic.mv_xor(p.load[si_port], scan_inversions[si_port])
+                    launch.data[scan_maps[si_port], i] = pattern.data[:, 0]
             if '_pi' in p.capture and 'P' in p.capture['_pi']:
-                launch.set_values(i, p.capture['_pi'], pi_map)
-        
-        return PackedVectors.from_pair(init, launch)
+                launch.data[pi_map, i] = logic.MVArray(p.capture['_pi']).data[:, 0]
+            launch.data[po_map, i] = logic.UNASSIGNED
+
+        return logic.mv_transition(init, launch)
                 
-    def responses4v(self, circuit):
-        """Assembles and returns a scan test response pattern set in 4-valued logic for given circuit."""
+    def responses(self, circuit):
+        """Assembles and returns a scan test response pattern set for given circuit."""
         interface, pi_map, po_map, scan_maps, scan_inversions = self._maps(circuit)
-        resp = PackedVectors(len(self.patterns), len(interface), 2)
+        resp = logic.MVArray((len(interface), len(self.patterns)))
+        # resp = PackedVectors(len(self.patterns), len(interface), 2)
         for i, p in enumerate(self.patterns):
-            if len(p.capture) > 0:
-                resp.set_values(i, p.capture['_po'], po_map)
-            else:
-                resp.set_values(i, p.launch['_po'], po_map)
+            resp.data[po_map, i] = logic.MVArray(p.capture['_po'] if len(p.capture) > 0 else p.launch['_po']).data[:, 0]
+            # if len(p.capture) > 0:
+            #    resp.set_values(i, p.capture['_po'], po_map)
+            # else:
+            #    resp.set_values(i, p.launch['_po'], po_map)
             for so_port in self.so_ports.keys():
-                resp.set_values(i, p.unload[so_port], scan_maps[so_port], scan_inversions[so_port])
+                pattern = logic.mv_xor(p.unload[so_port], scan_inversions[so_port])
+                resp.data[scan_maps[so_port], i] = pattern.data[:, 0]
+                # resp.set_values(i, p.unload[so_port], scan_maps[so_port], scan_inversions[so_port])
         return resp
         
         
diff --git a/src/kyupy/wave_sim.py b/src/kyupy/wave_sim.py
index b421229..2766997 100644
--- a/src/kyupy/wave_sim.py
+++ b/src/kyupy/wave_sim.py
@@ -260,21 +260,26 @@ class WaveSim:
         self.timing[line, 0, polarity] = delay
 
     def assign(self, vectors, time=0.0, offset=0):
-        nvectors = min(vectors.nvectors - offset, self.sims)
+        nvectors = min(len(vectors) - offset, self.sims)
         for i, node in enumerate(self.interface):
             ppi_loc = self.sat[self.ppi_offset + i, 0]
             if ppi_loc < 0: continue
             for p in range(nvectors):
                 vector = p + offset
-                a = vectors.bits[i, :, vector // 8]
+                a = vectors.data[i, :, vector // 8]
                 m = self.mask[vector % 8]
                 toggle = 0
-                if a[0] & m[1]:
-                    self.state[ppi_loc, p] = TMIN
-                    toggle += 1
-                if (len(a) > 2) and (a[2] & m[1]) and ((a[0] & m[1]) == (a[1] & m[1])):
-                    self.state[ppi_loc + toggle, p] = time
-                    toggle += 1
+                if len(a) <= 2:
+                    if a[0] & m[1]:
+                        self.state[ppi_loc, p] = TMIN
+                        toggle += 1
+                else:
+                    if a[1] & m[1]:
+                        self.state[ppi_loc, p] = TMIN
+                        toggle += 1
+                    if (a[2] & m[1]) and ((a[0] & m[1]) != (a[1] & m[1])):
+                        self.state[ppi_loc + toggle, p] = time
+                        toggle += 1
                 self.state[ppi_loc + toggle, p] = TMAX
 
     def propagate(self, sims=None, sd=0.0, seed=1):
@@ -563,12 +568,11 @@ class WaveSimCuda(WaveSim):
     def assign(self, vectors, time=0.0, offset=0):
         assert (offset % 8) == 0
         byte_offset = offset // 8
-        assert byte_offset < vectors.bits.shape[-1]
-        pdim = min(vectors.bits.shape[-1] - byte_offset, self.tdata.shape[-1])
+        assert byte_offset < vectors.data.shape[-1]
+        pdim = min(vectors.data.shape[-1] - byte_offset, self.tdata.shape[-1])
 
-        self.tdata[..., 0:pdim] = vectors.bits[..., byte_offset:pdim + byte_offset]
-        if vectors.vdim == 1:
-            self.tdata[:, 1, 0:pdim] = ~self.tdata[:, 1, 0:pdim]
+        self.tdata[..., 0:pdim] = vectors.data[..., byte_offset:pdim + byte_offset]
+        if vectors.m == 2:
             self.tdata[:, 2, 0:pdim] = 0
         cuda.to_device(self.tdata, to=self.d_tdata)
 
@@ -745,10 +749,10 @@ def assign_kernel(state, sat, ppi_offset, intf_len, tdata, time):
     a2 = tdata[y, 2, vector // 8]
     m = np.uint8(1 << (7 - (vector % 8)))
     toggle = 0
-    if a0 & m:
+    if a1 & m:
         state[line + toggle, x] = TMIN
         toggle += 1
-    if (a2 & m) and ((a0 & m) == (a1 & m)):
+    if (a2 & m) and ((a0 & m) != (a1 & m)):
         state[line + toggle, x] = time
         toggle += 1
     state[line + toggle, x] = TMAX
diff --git a/tests/test_logic.py b/tests/test_logic.py
index 1d6749e..8fb933a 100644
--- a/tests/test_logic.py
+++ b/tests/test_logic.py
@@ -7,10 +7,12 @@ def test_mvarray():
 
     ary = lg.MVArray(4)
     assert ary.length == 1
+    assert len(ary) == 1
     assert ary.width == 4
 
     ary = lg.MVArray((3, 2))
     assert ary.length == 2
+    assert len(ary) == 2
     assert ary.width == 3
 
     # instantiation with single vector
@@ -18,10 +20,14 @@ def test_mvarray():
     ary = lg.MVArray([1, 0, 1])
     assert ary.length == 1
     assert ary.width == 3
+    assert str(ary) == "['101']"
+    assert ary[0] == '101'
 
     ary = lg.MVArray("10X-")
     assert ary.length == 1
     assert ary.width == 4
+    assert str(ary) == "['10X-']"
+    assert ary[0] == '10X-'
 
     ary = lg.MVArray("1")
     assert ary.length == 1
@@ -40,6 +46,8 @@ def test_mvarray():
     ary = lg.MVArray(["000", "001", "110", "---"])
     assert ary.length == 4
     assert ary.width == 3
+    assert str(ary) == "['000', '001', '110', '---']"
+    assert ary[2] == '110'
 
     # casting to 2-valued logic
 
@@ -121,18 +129,86 @@ def test_mv_operations():
     x1_8v = lg.MVArray("00000000XXXXXXXX--------11111111PPPPPPPPRRRRRRRRFFFFFFFFNNNNNNNN", m=8)
     x2_8v = lg.MVArray("0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN", m=8)
 
-    assert str(lg.mv_not(x1_2v)) == "['1100']"
-    assert str(lg.mv_not(x1_4v)) == "['1111XXXXXXXX0000']"
-    assert str(lg.mv_not(x1_8v)) == "['11111111XXXXXXXXXXXXXXXX00000000NNNNNNNNFFFFFFFFRRRRRRRRPPPPPPPP']"
+    assert lg.mv_not(x1_2v)[0] == '1100'
+    assert lg.mv_not(x1_4v)[0] == '1111XXXXXXXX0000'
+    assert lg.mv_not(x1_8v)[0] == '11111111XXXXXXXXXXXXXXXX00000000NNNNNNNNFFFFFFFFRRRRRRRRPPPPPPPP'
 
-    assert str(lg.mv_or(x1_2v, x2_2v)) == "['0111']"
-    assert str(lg.mv_or(x1_4v, x2_4v)) == "['0XX1XXX1XXX11111']"
-    assert str(lg.mv_or(x1_8v, x2_8v)) == "['0XX1PRFNXXX1XXXXXXX1XXXX11111111PXX1PRFNRXX1RRNNFXX1FNFNNXX1NNNN']"
+    assert lg.mv_or(x1_2v, x2_2v)[0] == '0111'
+    assert lg.mv_or(x1_4v, x2_4v)[0] == '0XX1XXX1XXX11111'
+    assert lg.mv_or(x1_8v, x2_8v)[0] == '0XX1PRFNXXX1XXXXXXX1XXXX11111111PXX1PRFNRXX1RRNNFXX1FNFNNXX1NNNN'
 
-    assert str(lg.mv_and(x1_2v, x2_2v)) == "['0001']"
-    assert str(lg.mv_and(x1_4v, x2_4v)) == "['00000XXX0XXX0XX1']"
-    assert str(lg.mv_and(x1_8v, x2_8v)) == "['000000000XXXXXXX0XXXXXXX0XX1PRFN0XXPPPPP0XXRPRPR0XXFPPFF0XXNPRFN']"
+    assert lg.mv_and(x1_2v, x2_2v)[0] == '0001'
+    assert lg.mv_and(x1_4v, x2_4v)[0] == '00000XXX0XXX0XX1'
+    assert lg.mv_and(x1_8v, x2_8v)[0] == '000000000XXXXXXX0XXXXXXX0XX1PRFN0XXPPPPP0XXRPRPR0XXFPPFF0XXNPRFN'
 
-    assert str(lg.mv_xor(x1_2v, x2_2v)) == "['0110']"
-    assert str(lg.mv_xor(x1_4v, x2_4v)) == "['0XX1XXXXXXXX1XX0']"
-    assert str(lg.mv_xor(x1_8v, x2_8v)) == "['0XX1PRFNXXXXXXXXXXXXXXXX1XX0NFRPPXXNPRFNRXXFRPNFFXXRFNPRNXXPNFRP']"
+    assert lg.mv_xor(x1_2v, x2_2v)[0] == '0110'
+    assert lg.mv_xor(x1_4v, x2_4v)[0] == '0XX1XXXXXXXX1XX0'
+    assert lg.mv_xor(x1_8v, x2_8v)[0] == '0XX1PRFNXXXXXXXXXXXXXXXX1XX0NFRPPXXNPRFNRXXFRPNFFXXRFNPRNXXPNFRP'
+
+
+def test_bparray():
+
+    ary = lg.BPArray(4)
+    assert ary.length == 1
+    assert len(ary) == 1
+    assert ary.width == 4
+
+    ary = lg.BPArray((3, 2))
+    assert ary.length == 2
+    assert len(ary) == 2
+    assert ary.width == 3
+
+    assert lg.MVArray(lg.BPArray("01", m=2))[0] == '01'
+    assert lg.MVArray(lg.BPArray("0X-1", m=4))[0] == '0X-1'
+    assert lg.MVArray(lg.BPArray("0X-1PRFN", m=8))[0] == '0X-1PRFN'
+
+    x1_2v = lg.BPArray("0011", m=2)
+    x2_2v = lg.BPArray("0101", m=2)
+    x1_4v = lg.BPArray("0000XXXX----1111", m=4)
+    x2_4v = lg.BPArray("0X-10X-10X-10X-1", m=4)
+    x1_8v = lg.BPArray("00000000XXXXXXXX--------11111111PPPPPPPPRRRRRRRRFFFFFFFFNNNNNNNN", m=8)
+    x2_8v = lg.BPArray("0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN0X-1PRFN", m=8)
+
+    out_2v = lg.BPArray((4, 1), m=2)
+    out_4v = lg.BPArray((16, 1), m=4)
+    out_8v = lg.BPArray((64, 1), m=8)
+
+    lg.bp_buf(out_2v.data, x1_2v.data)
+    lg.bp_buf(out_4v.data, x1_4v.data)
+    lg.bp_buf(out_8v.data, x1_8v.data)
+
+    assert lg.MVArray(out_2v)[0] == '0011'
+    assert lg.MVArray(out_4v)[0] == '0000XXXXXXXX1111'
+    assert lg.MVArray(out_8v)[0] == '00000000XXXXXXXXXXXXXXXX11111111PPPPPPPPRRRRRRRRFFFFFFFFNNNNNNNN'
+
+    lg.bp_not(out_2v.data, x1_2v.data)
+    lg.bp_not(out_4v.data, x1_4v.data)
+    lg.bp_not(out_8v.data, x1_8v.data)
+
+    assert lg.MVArray(out_2v)[0] == '1100'
+    assert lg.MVArray(out_4v)[0] == '1111XXXXXXXX0000'
+    assert lg.MVArray(out_8v)[0] == '11111111XXXXXXXXXXXXXXXX00000000NNNNNNNNFFFFFFFFRRRRRRRRPPPPPPPP'
+
+    lg.bp_or(out_2v.data, x1_2v.data, x2_2v.data)
+    lg.bp_or(out_4v.data, x1_4v.data, x2_4v.data)
+    lg.bp_or(out_8v.data, x1_8v.data, x2_8v.data)
+
+    assert lg.MVArray(out_2v)[0] == '0111'
+    assert lg.MVArray(out_4v)[0] == '0XX1XXX1XXX11111'
+    assert lg.MVArray(out_8v)[0] == '0XX1PRFNXXX1XXXXXXX1XXXX11111111PXX1PRFNRXX1RRNNFXX1FNFNNXX1NNNN'
+
+    lg.bp_and(out_2v.data, x1_2v.data, x2_2v.data)
+    lg.bp_and(out_4v.data, x1_4v.data, x2_4v.data)
+    lg.bp_and(out_8v.data, x1_8v.data, x2_8v.data)
+
+    assert lg.MVArray(out_2v)[0] == '0001'
+    assert lg.MVArray(out_4v)[0] == '00000XXX0XXX0XX1'
+    assert lg.MVArray(out_8v)[0] == '000000000XXXXXXX0XXXXXXX0XX1PRFN0XXPPPPP0XXRPRPR0XXFPPFF0XXNPRFN'
+
+    lg.bp_xor(out_2v.data, x1_2v.data, x2_2v.data)
+    lg.bp_xor(out_4v.data, x1_4v.data, x2_4v.data)
+    lg.bp_xor(out_8v.data, x1_8v.data, x2_8v.data)
+
+    assert lg.MVArray(out_2v)[0] == '0110'
+    assert lg.MVArray(out_4v)[0] == '0XX1XXXXXXXX1XX0'
+    assert lg.MVArray(out_8v)[0] == '0XX1PRFNXXXXXXXXXXXXXXXX1XX0NFRPPXXNPRFNRXXFRPNFFXXRFNPRNXXPNFRP'
diff --git a/tests/test_logic_sim.py b/tests/test_logic_sim.py
index 08c364c..990eec7 100644
--- a/tests/test_logic_sim.py
+++ b/tests/test_logic_sim.py
@@ -1,161 +1,96 @@
 from kyupy.logic_sim import LogicSim
 from kyupy import bench
-from kyupy.packed_vectors import PackedVectors
+from kyupy.logic import MVArray, BPArray
 
 
-def test_vd1():
+def test_2v():
     c = bench.parse('input(x, y) output(a, o, n) a=and(x,y) o=or(x,y) n=not(x)')
-    s = LogicSim(c, 4)
+    s = LogicSim(c, 4, m=2)
     assert len(s.interface) == 5
-    p = PackedVectors(4, len(s.interface))
-    p[0] = '00000'
-    p[1] = '01000'
-    p[2] = '10000'
-    p[3] = '11000'
-    s.assign(p)
+    mva = MVArray(['00000', '01000', '10000', '11000'], m=2)
+    bpa = BPArray(mva)
+    s.assign(bpa)
     s.propagate()
-    s.capture(p)
-    assert p[0] == '00001'
-    assert p[1] == '01011'
-    assert p[2] == '10010'
-    assert p[3] == '11110'
+    s.capture(bpa)
+    mva = MVArray(bpa)
+    assert mva[0] == '00001'
+    assert mva[1] == '01011'
+    assert mva[2] == '10010'
+    assert mva[3] == '11110'
 
 
-def test_vd2():
+def test_4v():
     c = bench.parse('input(x, y) output(a, o, n) a=and(x,y) o=or(x,y) n=not(x)')
-    s = LogicSim(c, 16, 2)
+    s = LogicSim(c, 16, m=4)
     assert len(s.interface) == 5
-    p = PackedVectors(16, len(s.interface), 2)
-    p[0] = '00000'
-    p[1] = '01000'
-    p[2] = '0-000'
-    p[3] = '0X000'
-    p[4] = '10000'
-    p[5] = '11000'
-    p[6] = '1-000'
-    p[7] = '1X000'
-    p[8] = '-0000'
-    p[9] = '-1000'
-    p[10] = '--000'
-    p[11] = '-X000'
-    p[12] = 'X0000'
-    p[13] = 'X1000'
-    p[14] = 'X-000'
-    p[15] = 'XX000'
-    s.assign(p)
+    mva = MVArray(['00000', '01000', '0-000', '0X000',
+                   '10000', '11000', '1-000', '1X000',
+                   '-0000', '-1000', '--000', '-X000',
+                   'X0000', 'X1000', 'X-000', 'XX000'], m=4)
+    bpa = BPArray(mva)
+    s.assign(bpa)
     s.propagate()
-    s.capture(p)
-    assert p[0] == '00001'
-    assert p[1] == '01011'
-    assert p[2] == '0-0X1'
-    assert p[3] == '0X0X1'
-    assert p[4] == '10010'
-    assert p[5] == '11110'
-    assert p[6] == '1-X10'
-    assert p[7] == '1XX10'
-    assert p[8] == '-00XX'
-    assert p[9] == '-1X1X'
-    assert p[10] == '--XXX'
-    assert p[11] == '-XXXX'
-    assert p[12] == 'X00XX'
-    assert p[13] == 'X1X1X'
-    assert p[14] == 'X-XXX'
-    assert p[15] == 'XXXXX'
+    s.capture(bpa)
+    mva = MVArray(bpa)
+    assert mva[0] == '00001'
+    assert mva[1] == '01011'
+    assert mva[2] == '0-0X1'
+    assert mva[3] == '0X0X1'
+    assert mva[4] == '10010'
+    assert mva[5] == '11110'
+    assert mva[6] == '1-X10'
+    assert mva[7] == '1XX10'
+    assert mva[8] == '-00XX'
+    assert mva[9] == '-1X1X'
+    assert mva[10] == '--XXX'
+    assert mva[11] == '-XXXX'
+    assert mva[12] == 'X00XX'
+    assert mva[13] == 'X1X1X'
+    assert mva[14] == 'X-XXX'
+    assert mva[15] == 'XXXXX'
 
     
-def test_vd3():
+def test_8v():
     c = bench.parse('input(x, y) output(a, o, n, xo) a=and(x,y) o=or(x,y) n=not(x) xo=xor(x,y)')
-    s = LogicSim(c, 64, 3)
+    s = LogicSim(c, 64, m=8)
     assert len(s.interface) == 6
-    p = PackedVectors(64, len(s.interface), 3)
-    p[0] = '000010'
-    p[1] = '010111'
-    p[2] = '0-0X1X'
-    p[3] = '0X0X1X'
-    p[4] = '0R0R1R'
-    p[5] = '0F0F1F'
-    p[6] = '0P0P1P'
-    p[7] = '0N0N1N'
-    p[8] = '100101'
-    p[9] = '111100'
-    p[10] = '1-X10X'
-    p[11] = '1XX10X'
-    p[12] = '1RR10F'
-    p[13] = '1FF10R'
-    p[14] = '1PP10N'
-    p[15] = '1NN10P'
-    p[16] = '-00XXX'
-    p[17] = '-1X1XX'
-    p[18] = '--XXXX'
-    p[19] = '-XXXXX'
-    p[20] = '-RXXXX'
-    p[21] = '-FXXXX'
-    p[22] = '-PXXXX'
-    p[23] = '-NXXXX'
-    p[24] = 'X00XXX'
-    p[25] = 'X1X1XX'
-    p[26] = 'X-XXXX'
-    p[27] = 'XXXXXX'
-    p[28] = 'XRXXXX'
-    p[29] = 'XFXXXX'
-    p[30] = 'XPXXXX'
-    p[31] = 'XNXXXX'
-    p[32] = 'R00RFR'
-    p[33] = 'R1R1FF'
-    p[34] = 'R-XXFX'
-    p[35] = 'RXXXFX'
-    p[36] = 'RRRRFP'
-    p[37] = 'RFPNFN'
-    p[38] = 'RPPRFR'
-    p[39] = 'RNRNFF'
-    p[40] = 'F00FRF'
-    p[41] = 'F1F1RR'
-    p[42] = 'F-XXRX'
-    p[43] = 'FXXXRX'
-    p[44] = 'FRPNRN'
-    p[45] = 'FFFFRP'
-    p[46] = 'FPPFRF'
-    p[47] = 'FNFNRR'
-    p[48] = 'P00PNP'
-    p[49] = 'P1P1NN'
-    p[50] = 'P-XXNX'
-    p[51] = 'PXXXNX'
-    p[52] = 'PRPRNR'
-    p[53] = 'PFPFNF'
-    p[54] = 'PPPPNP'
-    p[55] = 'PNPNNN'
-    p[56] = 'N00NPN'
-    p[57] = 'N1N1PP'
-    p[58] = 'N-XXPX'
-    p[59] = 'NXXXPX'
-    p[60] = 'NRRNPF'
-    p[61] = 'NFFNPR'
-    p[62] = 'NPPNPN'
-    p[63] = 'NNNNPP'
-    expect = p.copy()
-    s.assign(p)
+    mva = MVArray(['000010', '010111', '0-0X1X', '0X0X1X', '0R0R1R', '0F0F1F', '0P0P1P', '0N0N1N',
+                   '100101', '111100', '1-X10X', '1XX10X', '1RR10F', '1FF10R', '1PP10N', '1NN10P',
+                   '-00XXX', '-1X1XX', '--XXXX', '-XXXXX', '-RXXXX', '-FXXXX', '-PXXXX', '-NXXXX',
+                   'X00XXX', 'X1X1XX', 'X-XXXX', 'XXXXXX', 'XRXXXX', 'XFXXXX', 'XPXXXX', 'XNXXXX',
+                   'R00RFR', 'R1R1FF', 'R-XXFX', 'RXXXFX', 'RRRRFP', 'RFPNFN', 'RPPRFR', 'RNRNFF',
+                   'F00FRF', 'F1F1RR', 'F-XXRX', 'FXXXRX', 'FRPNRN', 'FFFFRP', 'FPPFRF', 'FNFNRR',
+                   'P00PNP', 'P1P1NN', 'P-XXNX', 'PXXXNX', 'PRPRNR', 'PFPFNF', 'PPPPNP', 'PNPNNN',
+                   'N00NPN', 'N1N1PP', 'N-XXPX', 'NXXXPX', 'NRRNPF', 'NFFNPR', 'NPPNPN', 'NNNNPP'], m=8)
+    bpa = BPArray(mva)
+    s.assign(bpa)
     s.propagate()
-    s.capture(p)
+    resp_bp = BPArray(bpa)
+    s.capture(resp_bp)
+    resp = MVArray(resp_bp)
+
     for i in range(64):
-        assert p[i] == expect[i]
+        assert resp[i] == mva[i]
         
 
 def test_b01(mydir):
     c = bench.load(mydir / 'b01.bench')
 
     # 2-valued
-    s = LogicSim(c, 8)
+    s = LogicSim(c, 8, m=2)
     assert len(s.interface) == 9
-    t = PackedVectors(8, len(s.interface))
-    t.randomize()
-    s.assign(t)
+    mva = MVArray((len(s.interface), 8), m=2)
+    # mva.randomize()
+    bpa = BPArray(mva)
+    s.assign(bpa)
     s.propagate()
-    s.capture(t)
+    s.capture(bpa)
 
     # 8-valued
-    s = LogicSim(c, 8, 3)
-    t = PackedVectors(8, len(s.interface), 3)
-    t.randomize()
-    s.assign(t)
+    s = LogicSim(c, 8, m=8)
+    mva = MVArray((len(s.interface), 8), m=8)
+    # mva.randomize()
+    bpa = BPArray(mva)
+    s.assign(bpa)
     s.propagate()
-    s.capture(t)
+    s.capture(bpa)
diff --git a/tests/test_packed_vectors.py b/tests/test_packed_vectors.py
deleted file mode 100644
index 2f2a4a0..0000000
--- a/tests/test_packed_vectors.py
+++ /dev/null
@@ -1,88 +0,0 @@
-from kyupy.packed_vectors import PackedVectors
-
-
-def test_basic():
-    ba = PackedVectors(8, 1, 1)
-    assert '0\n0\n0\n0\n0\n0\n0\n0' == str(ba)
-    ba.set_value(0, 0, 1)
-    ba.set_value(1, 0, 'H')
-    ba.set_value(2, 0, 'h')
-    ba.set_value(3, 0, True)
-    ba.set_value(4, 0, 0)
-    ba.set_value(5, 0, 'L')
-    ba.set_value(6, 0, 'l')
-    ba.set_value(7, 0, False)
-    assert '1\n1\n1\n1\n0\n0\n0\n0' == str(ba)
-    ba.set_value(1, 0, '0')
-    ba.set_value(5, 0, '1')
-    assert '1\n0\n1\n1\n0\n1\n0\n0' == str(ba)
-    ba = PackedVectors(8, 1, 2)
-    assert '-\n-\n-\n-\n-\n-\n-\n-' == str(ba)
-    ba.set_value(0, 0, 1)
-    ba.set_value(7, 0, 0)
-    ba.set_value(4, 0, 'X')
-    assert '1\n-\n-\n-\nX\n-\n-\n0' == str(ba)
-    ba.set_value(4, 0, '-')
-    assert '1\n-\n-\n-\n-\n-\n-\n0' == str(ba)
-    ba = PackedVectors(8, 2, 2)
-    assert '--\n--\n--\n--\n--\n--\n--\n--' == str(ba)
-    ba.set_value(0, 0, '1')
-    ba.set_value(7, 1, '0')
-    ba.set_values(1, 'XX')
-    assert '1-\nXX\n--\n--\n--\n--\n--\n-0' == str(ba)
-
-
-def test_8v():
-    ba = PackedVectors(1, 8, 3)
-    assert '--------' == str(ba)
-    ba.set_values(0, r'-x01^v\/')
-    assert r'-X01PNFR' == str(ba)
-    ba.set_values(0, '-XLHPNFR')
-    assert r'-X01PNFR' == str(ba)
-    ba.set_values(0, '-xlhpnfr')
-    assert r'-X01PNFR' == str(ba)
-    p1 = PackedVectors(1, 8, 1)
-    p2 = PackedVectors(1, 8, 1)
-    p1.set_values(0, '01010101')
-    p2.set_values(0, '00110011')
-    p = PackedVectors.from_pair(p1, p2)
-    assert r'0FR10FR1' == str(p)
-    p1 = PackedVectors(1, 8, 2)
-    p2 = PackedVectors(1, 8, 2)
-    p1.set_values(0, '0101-X-X')
-    p2.set_values(0, '00110011')
-    p = PackedVectors.from_pair(p1, p2)
-    assert r'0FR1----' == str(p)
-    p1.set_values(0, '0101-X-X')
-    p2.set_values(0, '-X-X--XX')
-    p = PackedVectors.from_pair(p1, p2)
-    assert r'--------' == str(p)
-
-
-def test_slicing():
-    lv = PackedVectors(3, 2, 1)
-    assert '00\n00\n00' == str(lv)
-    lv.set_value(1, 0, '1')
-    lv.set_value(1, 1, '1')
-    assert '00' == lv[0]
-    assert '11' == lv[1]
-    assert 3 == len(lv)
-    lv2 = lv[1:3]
-    assert 2 == len(lv2)
-    assert '11' == lv2[0]
-    assert '00' == lv2[1]
-
-
-def test_copy():
-    lv1 = PackedVectors(8, 1, 1)
-    lv1.set_values_for_position(0, '01010101')
-    lv2 = PackedVectors(8, 1, 1)
-    lv2.set_values_for_position(0, '00100101')
-    diff = lv1.diff(lv2)
-    lv3 = lv1.copy(selection_mask=diff)
-    assert str(lv3) == '1\n0\n1'
-    lv4 = lv1.copy(selection_mask=~diff)
-    assert str(lv4) == '0\n0\n1\n0\n1'
-    lv5 = lv3 + lv4
-    assert str(lv5) == '1\n0\n1\n0\n0\n1\n0\n1'
-
diff --git a/tests/test_wave_sim.py b/tests/test_wave_sim.py
index 7cef055..bea26d3 100644
--- a/tests/test_wave_sim.py
+++ b/tests/test_wave_sim.py
@@ -1,10 +1,10 @@
 import numpy as np
+
 from kyupy.wave_sim import WaveSim, WaveSimCuda, wave_eval, TMIN, TMAX
 from kyupy.logic_sim import LogicSim
-from kyupy import verilog
-from kyupy import sdf
+from kyupy import verilog, sdf, logic
 from kyupy.saed import pin_index
-from kyupy.packed_vectors import PackedVectors
+from kyupy.logic import MVArray, BPArray
 
 
 def test_wave_eval():
@@ -95,24 +95,29 @@ def test_wave_eval():
 
 
 def compare_to_logic_sim(wsim):
-    tests = PackedVectors(wsim.sims, len(wsim.interface), 3)
-    tests.randomize()
-    wsim.assign(tests)
-    wsim.propagate(8)
+    tests = MVArray((len(wsim.interface), wsim.sims))
+    choices = np.asarray([logic.ZERO, logic.ONE, logic.RISE, logic.FALL], dtype=np.uint8)
+    rng = np.random.default_rng(10)
+    tests.data[...] = rng.choice(choices, tests.data.shape)
+    tests_bp = BPArray(tests)
+    wsim.assign(tests_bp)
+    wsim.propagate()
     cdata = wsim.capture()
 
-    resp = tests.copy()
+    resp = MVArray(tests)
 
     for iidx, inode in enumerate(wsim.interface):
         if len(inode.ins) > 0:
             for vidx in range(wsim.sims):
-                resp.set_value(vidx, iidx, 0 if cdata[iidx, vidx, 0] < 0.5 else 1)
+                resp.data[iidx, vidx] = logic.ZERO if cdata[iidx, vidx, 0] < 0.5 else logic.ONE
+                # resp.set_value(vidx, iidx, 0 if cdata[iidx, vidx, 0] < 0.5 else 1)
 
-    lsim = LogicSim(wsim.circuit, len(tests), 3)
-    lsim.assign(tests)
+    lsim = LogicSim(wsim.circuit, len(tests_bp))
+    lsim.assign(tests_bp)
     lsim.propagate()
-    exp = tests.copy()
-    lsim.capture(exp)
+    exp_bp = BPArray(tests_bp)
+    lsim.capture(exp_bp)
+    exp = MVArray(exp_bp)
     
     for i in range(8):
         exp_str = exp[i].replace('R', '1').replace('F', '0').replace('P', '0').replace('N', '1')