diff --git a/apycula/bitmatrix.py b/apycula/bitmatrix.py
new file mode 100644
index 00000000..2c8332f2
--- /dev/null
+++ b/apycula/bitmatrix.py
@@ -0,0 +1,113 @@
+
+def fliplr(bmp):
+    """
+    Flips the entries in each row in the left/right direction.
+    Returns a new matrix.
+    """
+    return [row[::-1] for row in bmp]
+
+def flipud(bmp):
+    """
+    Reverse the order of elements in each column (up/down).
+    Returns a refence.
+    """
+    return bmp[::-1]
+
+def vstack(bmp_0, bmp_1):
+    """
+    Stack matrices in sequence vertically (row wise).
+    Returns a reference.
+    """
+    return [*bmp_0, *bmp_1]
+
+def hstack(bmp_0, bmp_1):
+    """
+    Stack matrices in sequence horizontally (column wise).
+    Returns a new matrix.
+    """
+    return [bmp[0] + bmp[1] for bmp in zip(bmp_0, bmp_1)]
+
+def shape(bmp):
+    """
+    Return the shape of a matrix.
+    """
+    return [len(bmp), len(bmp[0])]
+
+def ones(rows, cols):
+    """
+    Returns a new matrix of given shape, filled with ones.
+    """
+    return [[1] * cols for i in range(rows)]
+
+def zeros(rows, cols):
+    """
+    Returns a new matrix of given shape, filled with zeros.
+    """
+    return [[0] * cols for i in range(rows)]
+
+def packbits(bmp, axis = None):
+    """
+    Packs the elements of a bitmap into bytes.
+    [1, 1, 0, 0, 0] -> [24]  # [5'b11000]
+    Returns a list of bytes.
+    """
+    byte_list = []
+    if not axis:
+        for bmp_r in bmp:
+            for col in range(shape(bmp)[1] // 8):
+                bcol = col << 3
+                byte_list.append((bmp_r[bcol] << 7) + (bmp_r[bcol + 1] << 6) + (bmp_r[bcol + 2] << 5) +
+                    (bmp_r[bcol + 3] << 4) + (bmp_r[bcol + 4] << 3) + (bmp_r[bcol + 5] << 2) +
+                    (bmp_r[bcol + 6] << 1) + bmp_r[bcol + 7])
+    else:
+        for bmp_r in bmp:
+            byte_list.append([])
+            byte_list_r = byte_list[-1]
+            for col in range(shape(bmp)[1] // 8):
+                bcol = col << 3
+                byte_list_r.append((bmp_r[bcol] << 7) + (bmp_r[bcol + 1] << 6) + (bmp_r[bcol + 2] << 5) +
+                    (bmp_r[bcol + 3] << 4) + (bmp_r[bcol + 4] << 3) + (bmp_r[bcol + 5] << 2) +
+                    (bmp_r[bcol + 6] << 1) + bmp_r[bcol + 7])
+    return byte_list
+
+def xor(bmp_0, bmp_1):
+    """
+    Bitwise XOR
+    Returns a new matrix
+    """
+    return [[ vals[0] ^ vals[1]for vals in zip(row[0], row[1])] for row in zip(bmp_0, bmp_1)]
+
+def histogram(lst, bins):
+    """
+    Compute the histogram of a list.
+    Returns a list of counters.
+    """
+    l_bins = len(bins) - 1
+    r_lst = [0] * l_bins
+    for val in lst:
+        for i in range(l_bins):
+            if val in range(bins[i], bins[i + 1]) or (i == l_bins - 1 and val == bins[-1]):
+                r_lst[i] += 1
+    return r_lst
+
+def any(bmp):
+    """
+    Test whether any matrix element evaluates to True.
+    """
+    for row in bmp:
+        for val in row:
+            if val:
+                return True
+    return False
+
+def nonzero(bmp):
+    """
+    Return the indices of the elements that are non-zero.
+    """
+    res = ([], [])
+    for ri, row in enumerate(bmp):
+        for ci, val in enumerate(row):
+            if val:
+                res[0].append(ri)
+                res[1].append(ci)
+    return res
diff --git a/apycula/bslib.py b/apycula/bslib.py
index af5ddb7f..9eb7ed74 100644
--- a/apycula/bslib.py
+++ b/apycula/bslib.py
@@ -1,6 +1,6 @@
 from math import ceil
-import numpy as np
 import crc
+from apycula import bitmatrix
 
 crc16arc = crc.Configuration(width=16, polynomial=0x8005, reverse_input=True, reverse_output=True)
 
@@ -71,7 +71,7 @@ def read_bitstream(fname):
             bitmap.append(bitarr(line, padding))
             frames = max(0, frames-1)
 
-    return np.fliplr(np.array(bitmap)), hdr, ftr
+    return bitmatrix.fliplr(bitmap), hdr, ftr
 
 def compressLine(line, key8Z, key4Z, key2Z):
     newline = []
@@ -82,26 +82,26 @@ def compressLine(line, key8Z, key4Z, key2Z):
     return newline
 
 def write_bitstream_with_bsram_init(fname, bs, hdr, ftr, compress, bsram_init):
-    new_bs = np.vstack((bs, bsram_init))
+    new_bs = bitmatrix.vstack(bs, bsram_init)
     new_hdr = hdr.copy()
-    frames = int.from_bytes(new_hdr[-1][2:], 'big') + bsram_init.shape[0]
+    frames = int.from_bytes(new_hdr[-1][2:], 'big') + bitmatrix.shape(bsram_init)[0]
     new_hdr[-1][2:] = frames.to_bytes(2, 'big')
     write_bitstream(fname, new_bs, new_hdr, ftr, compress)
 
 def write_bitstream(fname, bs, hdr, ftr, compress):
-    bs = np.fliplr(bs)
+    bs = bitmatrix.fliplr(bs)
     if compress:
-        padlen = (ceil(bs.shape[1] / 64) * 64) - bs.shape[1]
+        padlen = (ceil(bitmatrix.shape(bs)[1] / 64) * 64) - bitmatrix.shape(bs)[1]
     else:
-        padlen = bs.shape[1] % 8
-    pad = np.ones((bs.shape[0], padlen), dtype=np.uint8)
-    bs = np.hstack([pad, bs])
-    assert bs.shape[1] % 8 == 0
-    bs=np.packbits(bs, axis=1)
+        padlen = bitmatrix.shape(bs)[1] % 8
+    pad = bitmatrix.ones(bitmatrix.shape(bs)[0], padlen)
+    bs = bitmatrix.hstack(pad, bs)
+    assert bitmatrix.shape(bs)[1] % 8 == 0
+    bs=bitmatrix.packbits(bs, axis = 1)
 
     if compress:
         # search for smallest values not used in the bitstream
-        lst, _ = np.histogram(bs, bins=[i for i in range(256)])
+        lst = bitmatrix.histogram(bs, bins=[i for i in range(257)]) # 257 iso that the last basket is [255, 256] and not [254, 255]
         [key8Z, key4Z, key2Z] = [i for i,val in enumerate(lst) if val==0][0:3]
 
         # update line 0x51 with keys
@@ -140,7 +140,7 @@ def display(fname, data):
     im = Image.frombytes(
             mode='1',
             size=data.shape[::-1],
-            data=np.packbits(data, axis=1))
+            data=bitmatrix.packbits(data, axis = 1))
     if fname:
         im.save(fname)
     return im
diff --git a/apycula/chipdb.py b/apycula/chipdb.py
index ccb0b3f8..1cc1ec31 100644
--- a/apycula/chipdb.py
+++ b/apycula/chipdb.py
@@ -5,11 +5,11 @@
 import copy
 from functools import reduce
 from collections import namedtuple
-import numpy as np
 from apycula.dat19 import Datfile
 import apycula.fuse_h4x as fuse
 from apycula.wirenames import wirenames, wirenumbers, clknames, clknumbers, hclknames, hclknumbers
 from apycula import pindef
+from apycula import bitmatrix
 
 # the character that marks the I/O attributes that come from the nextpnr
 mode_attr_sep = '&'
@@ -75,7 +75,7 @@ class Device:
     pin_bank: Dict[str, int] = field(default_factory = dict)
     cmd_hdr: List[ByteString] = field(default_factory=list)
     cmd_ftr: List[ByteString] = field(default_factory=list)
-    template: np.ndarray = None
+    template: List[List[int]] = None
     # allowable values of bel attributes
     # {table_name: [(attr_id, attr_value)]}
     logicinfo: Dict[str, List[Tuple[int, int]]] = field(default_factory=dict)
@@ -2093,8 +2093,8 @@ def tile_bitmap(dev, bitmap, empty=False):
         for jdx, td in enumerate(row):
             w = td.width
             h = td.height
-            tile = bitmap[y:y+h,x:x+w]
-            if tile.any() or empty:
+            tile = [row[x:x+w] for row in bitmap[y:y+h]]
+            if bitmatrix.any(tile) or empty:
                 res[(idx, jdx)] = tile
             x+=w
         y+=h
@@ -2102,16 +2102,22 @@ def tile_bitmap(dev, bitmap, empty=False):
     return res
 
 def fuse_bitmap(db, bitmap):
-    res = np.zeros((db.height, db.width), dtype=np.uint8)
+    res = bitmatrix.zeros(db.height, db.width)
     y = 0
     for idx, row in enumerate(db.grid):
         x=0
         for jdx, td in enumerate(row):
             w = td.width
             h = td.height
-            res[y:y+h,x:x+w] = bitmap[(idx, jdx)]
-            x+=w
-        y+=h
+            y0 = y
+            for row in bitmap[(idx, jdx)]:
+                x0 = x
+                for val in row:
+                    res[y0][x0] = val
+                    x0 += 1
+                y0 += 1
+            x += w
+        y += h
 
     return res
 
diff --git a/apycula/clock_fuzzer.py b/apycula/clock_fuzzer.py
index cd03af46..6d60aa23 100644
--- a/apycula/clock_fuzzer.py
+++ b/apycula/clock_fuzzer.py
@@ -9,6 +9,7 @@
 from apycula import fuse_h4x
 from apycula import dat19
 from apycula import gowin_unpack
+from apycula import bitmatrix
 from apycula.wirenames import clknumbers
 
 def dff(mod, cst, row, col, clk=None):
@@ -99,7 +100,7 @@ def quadrants():
 
     res = {}
     for (row, col), (mybs, *_) in zip(idxes, pnr_res):
-        sweep_tiles = fuse_h4x.tile_bitmap(fse, mybs^pnr.bitmap)
+        sweep_tiles = fuse_h4x.tile_bitmap(fse, bitmatrix.xor(mybs, pnr.bitmap))
 
         # find which tap was used
         taps = [r for (r, c, typ), t in sweep_tiles.items() if typ in {13, 14, 15, 16, 18, 19}]
@@ -148,7 +149,7 @@ def center_muxes(ct, rows, cols):
     base = pnr.bitmap
     for i, (bs_sweep, *_) in enumerate(pnr_res):
         pin = true_pins[i]
-        new = base ^ bs_sweep
+        new = bitmatrix.xor(base, bs_sweep)
         tiles = chipdb.tile_bitmap(db, new)
 
         try:
diff --git a/apycula/fuse_h4x.py b/apycula/fuse_h4x.py
index ff763bc3..75e397e7 100644
--- a/apycula/fuse_h4x.py
+++ b/apycula/fuse_h4x.py
@@ -1,6 +1,6 @@
 import sys
-import numpy as np
 import random
+from apycula import bitmatrix
 
 def rint(f, w):
     val = int.from_bytes(f.read(w), 'little', signed=True)
@@ -74,7 +74,7 @@ def readOneFile(f, fuselength):
 def render_tile(d, ttyp):
     w = d[ttyp]['width']
     h = d[ttyp]['height']
-    tile = np.zeros((h, w), np.uint8)#+(255-ttyp)
+    tile = bitmatrix.zeros(h, w)#+(255-ttyp)
     for start, table in [(2, 'shortval'), (2, 'wire'), (16, 'longval'),
                          (1, 'longfuse'), (0, 'const')]:
         if table in d[ttyp]:
@@ -104,17 +104,24 @@ def render_bitmap(d):
     tiles = d['header']['grid'][61]
     width = sum([d[i]['width'] for i in tiles[0]])
     height = sum([d[i[0]]['height'] for i in tiles])
-    bitmap = np.zeros((height, width), np.uint8)
+    bitmap = bitmatrix.zeros(height, width)
     y = 0
     for row in tiles:
         x=0
         for typ in row:
-            #if typ==12: pdb.set_trace()
             td = d[typ]
             w = td['width']
             h = td['height']
             #bitmap[y:y+h,x:x+w] += render_tile(d, typ)
-            bitmap[y:y+h,x:x+w] = typ
+            #bitmap[y:y+h,x:x+w] = typ
+            rtile = render_tile(d, typ)
+            y0 = y
+            for row in rtile:
+                x0 = x
+                for val in row:
+                    bitmap[y0][x0] += val
+                    x0 += 1
+                y0 += 1
             x+=w
         y+=h
 
@@ -122,6 +129,8 @@ def render_bitmap(d):
 
 def display(fname, data):
     from PIL import Image
+    import numpy as np
+    data = np.array(data, dtype = np.uint8)
     im = Image.frombytes(
             mode='P',
             size=data.shape[::-1],
@@ -148,12 +157,11 @@ def tile_bitmap(d, bitmap, empty=False):
     for idx, row in enumerate(tiles):
         x=0
         for jdx, typ in enumerate(row):
-            #if typ==87: pdb.set_trace()
             td = d[typ]
             w = td['width']
             h = td['height']
-            tile = bitmap[y:y+h,x:x+w]
-            if tile.any() or empty:
+            tile = [row[x:x+w] for row in bitmap[y:y+h]]
+            if bitmatrix.any(tile) or empty:
                 res[(idx, jdx, typ)] = tile
             x+=w
         y+=h
@@ -164,7 +172,7 @@ def fuse_bitmap(d, bitmap):
     tiles = d['header']['grid'][61]
     width = sum([d[i]['width'] for i in tiles[0]])
     height = sum([d[i[0]]['height'] for i in tiles])
-    res = np.zeros((height, width), dtype=np.uint8)
+    res = bitmatrix.zeros(height, width)
     y = 0
     for idx, row in enumerate(tiles):
         x=0
@@ -172,7 +180,13 @@ def fuse_bitmap(d, bitmap):
             td = d[typ]
             w = td['width']
             h = td['height']
-            res[y:y+h,x:x+w] = bitmap[(idx, jdx, typ)]
+            y0 = y
+            for row in bitmap[(idx, jdx, typ)]:
+                x0 = x
+                for val in row:
+                    res[y0][x0] = val
+                    x0 += 1
+                y0 += 1
             x+=w
         y+=h
 
@@ -207,7 +221,7 @@ def scan_fuses(d, ttyp, tile):
     w = d[ttyp]['width']
     h = d[ttyp]['height']
     fuses = []
-    rows, cols = np.where(tile==1)
+    rows, cols = bitmatrix.nonzero(tile)
     for row, col in zip(rows, cols):
         # ripe for optimization
         for fnum, fuse in enumerate(d['header']['fuse'][1]):
diff --git a/apycula/gowin_pack.py b/apycula/gowin_pack.py
index eb6d721c..ccfb68c7 100644
--- a/apycula/gowin_pack.py
+++ b/apycula/gowin_pack.py
@@ -5,7 +5,6 @@
 import gzip
 import itertools
 import math
-import numpy as np
 import json
 import argparse
 import importlib.resources
@@ -16,6 +15,7 @@
 from apycula.chipdb import add_attr_val, get_shortval_fuses, get_longval_fuses, get_bank_fuses
 from apycula import attrids
 from apycula import bslib
+from apycula import bitmatrix
 from apycula.wirenames import wirenames, wirenumbers
 
 device = ""
@@ -110,9 +110,9 @@ def store_bsram_init_val(db, row, col, typ, parms, attrs):
     if not has_bsram_init:
         has_bsram_init = True
         # 256 * bsram rows * chip bit width
-        bsram_init_map = np.zeros((256 * len(db.simplio_rows), db.template.shape[1]), dtype=np.uint8)
+        bsram_init_map = bitmatrix.zeros(256 * len(db.simplio_rows), bitmatrix.shape(db.template)[1])
     # 3 BSRAM cells have width 3 * 60
-    loc_map = np.zeros((256, 3 * 60), dtype = np.int8)
+    loc_map = bitmatrix.zeros(256, 3 * 60)
     #print("mapping")
     if not subtype.strip():
         width = 256
@@ -162,8 +162,13 @@ def get_bits(init_data):
     x = 0
     for jdx in range(col):
         x += db.grid[0][jdx].width
-    loc_map = np.flipud(loc_map)
-    bsram_init_map[y:y + height, x:x + 3 * 60] = loc_map
+    loc_map = bitmatrix.flipud(loc_map)
+    for row in loc_map:
+        x0 = x
+        for val in row:
+            bsram_init_map[y][x0] = val
+            x0 += 1
+        y += 1
 
 _bsram_cell_types = {'DP', 'SDP', 'SP', 'ROM'}
 def get_bels(data):
@@ -1245,13 +1250,11 @@ def header_footer(db, bs, compress):
     Currently limited to checksum with
     CRC_check and security_bit_enable set
     """
-    bs = np.fliplr(bs)
-    bs=np.packbits(bs)
+    bs = bitmatrix.fliplr(bs)
+    bs = bitmatrix.packbits(bs)
     # configuration data checksum is computed on all
     # data in 16bit format
-    bb = np.array(bs)
-
-    res = int(bb[0::2].sum() * pow(2,8) + bb[1::2].sum())
+    res = int(sum(bs[0::2]) * pow(2,8) + sum(bs[1::2]))
     checksum = res & 0xffff
 
     if compress:
diff --git a/apycula/gowin_unpack.py b/apycula/gowin_unpack.py
index 9099b3eb..b0a2f714 100644
--- a/apycula/gowin_unpack.py
+++ b/apycula/gowin_unpack.py
@@ -2,7 +2,6 @@
 import os
 import re
 import random
-import numpy as np
 from itertools import chain, count
 import pickle
 import gzip
diff --git a/apycula/tiled_fuzzer.py b/apycula/tiled_fuzzer.py
index dc655b27..87ec00ab 100644
--- a/apycula/tiled_fuzzer.py
+++ b/apycula/tiled_fuzzer.py
@@ -10,7 +10,6 @@
 from random import shuffle, seed
 from warnings import warn
 from math import factorial
-import numpy as np
 from multiprocessing.dummy import Pool
 import pickle
 from shutil import copytree
diff --git a/setup.py b/setup.py
index 5e6dd93f..31b51e6d 100644
--- a/setup.py
+++ b/setup.py
@@ -18,7 +18,6 @@
         "Operating System :: OS Independent",
     ],
     install_requires=[
-        'numpy',
         'crc',
     ],
     python_requires='>=3.6',