diff --git a/Project.toml b/Project.toml
index 4552d47..060b6e5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,20 +7,28 @@ version = "0.5.1"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 nauty_jll = "55c6dc9b-343a-50ca-8ff2-b71adb3733d5"
 
 [compat]
+Aqua = "0.8"
 Graphs = "1.9"
+Interfaces = "0.3"
+JET = "0.4, 0.9, 0.10"
+LinearAlgebra = "1"
+Pkg = "1"
+Random = "1"
 SHA = "0.7, 1"
-SparseArrays = "1.6"
+Test = "1"
 julia = "1.6"
 nauty_jll = "2.8.9"
 
 [extras]
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Interfaces = "85a1e053-f937-4924-92a5-1367d23b7b87"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "BenchmarkTools"]
+test = ["Aqua", "JET", "Test", "Random", "Pkg", "Interfaces"]
diff --git a/src/NautyGraphs.jl b/src/NautyGraphs.jl
index 1657cd9..e84714d 100644
--- a/src/NautyGraphs.jl
+++ b/src/NautyGraphs.jl
@@ -1,18 +1,17 @@
 module NautyGraphs
 
-using Graphs, SparseArrays, LinearAlgebra, SHA
+using Graphs, LinearAlgebra, SHA
+using Graphs.SimpleGraphs: SimpleEdgeIter
 import nauty_jll
 
-const libnauty = nauty_jll.libnautyTL
-const WORDSIZE = 64
-const WordType = Culong
-#const WordType = WORDSIZE == 32 ? Cuint : WORDSIZE == 64 ? Culong : error("only wordsize 32 or 64 supported") 
 const Cbool = Cint
 const HashType = UInt
 
-include("densenautygraphs.jl")
+abstract type AbstractNautyGraph{T} <: AbstractGraph{T} end
+
 include("utils.jl")
-include("bitutils.jl")
+include("graphset.jl")
+include("densenautygraph.jl")
 include("nauty.jl")
 
 const NautyGraph = DenseNautyGraph{false}
@@ -20,8 +19,9 @@ const NautyDiGraph = DenseNautyGraph{true}
 
 function __init__()
     # global default options to nauty carry a pointer reference that needs to be initialized at runtime
-    libnauty_dispatch = cglobal((:dispatch_graph, libnauty), Cvoid)
-    DEFAULT_OPTIONS.dispatch = libnauty_dispatch
+    DEFAULTOPTIONS16.dispatch = cglobal((:dispatch_graph, libnauty(UInt16)), Cvoid)
+    DEFAULTOPTIONS32.dispatch = cglobal((:dispatch_graph, libnauty(UInt32)), Cvoid)
+    DEFAULTOPTIONS64.dispatch = cglobal((:dispatch_graph, libnauty(UInt64)), Cvoid)
     return
 end
 
@@ -29,6 +29,7 @@ export
     AbstractNautyGraph,
     NautyGraph,
     NautyDiGraph,
+    DenseNautyGraph,
     AutomorphismGroup,
     labels,
     nauty,
diff --git a/src/bitutils.jl b/src/bitutils.jl
deleted file mode 100644
index ef69c83..0000000
--- a/src/bitutils.jl
+++ /dev/null
@@ -1,135 +0,0 @@
-function push_bits(word::WordType, fill::WordType, n::Integer, k::Integer)
-    # Pushes n rightmost bits k steps to the left, overwriting step by step and replacing with values from fill
-    #TODO: OPTIMIZE THIS
-    wmax = one(WordType) << (WORDSIZE - 1)
-    move_mask = reduce(|, one(WordType) << i for i in 0:n-1; init=zero(WordType))
-    fill_mask = reduce(|, wmax >> i for i in 0:k-1; init=zero(WordType))
-    write_mask = ~reduce(|, one(WordType) << i for i in 0:n+k-1; init=zero(WordType))
-
-    move_chunk = (word & move_mask) << k
-    copy_chunk = (fill & fill_mask) >> (WORDSIZE - k)
-
-    return (word & write_mask) | move_chunk | copy_chunk
-end
-function rightshift_set!(set::Vector{WordType}, offset::Integer)
-    # @assert offset > 0
-
-    overflow = zero(WordType)
-    for i in eachindex(set)
-        w = set[i]
-        set[i] = (w >> offset) | overflow
-        overflow = w << (WORDSIZE - offset)
-    end
-    if overflow > 0
-        push!(set, overflow)
-    end
-end
-function transfer_set!(target::Vector{WordType}, set::Vector{WordType}, offset::Integer, m_target::Integer=1, m_set::Integer=1)
-    # @assert offset >= 0
-
-    overflow = zero(WordType)
-    dn = offset ÷ WORDSIZE
-    dk = offset % WORDSIZE
-    for i in eachindex(set)
-        target_idx = dn + 1 + (i - 1) % m_set + m_target * ((i - 1) ÷ m_set)
-        w = set[i]
-        target[target_idx] |= (w >> dk)
-        overflow = w << (WORDSIZE - dk)
-        if overflow > 0
-            target[target_idx+1] |= overflow
-        end
-    end
-    return
-end
-
-
-@inline function has_bit(word::WordType, i::Integer)
-    mask = one(WordType) << (WORDSIZE - i)
-    return (mask & word) != zero(WordType)
-end
-
-@inline function count_bits(word::WordType)
-    return sum(has_bit(word, i) for i in 1:WORDSIZE)
-end
-
-function word_to_bitvec(word::WordType)
-    bitvec = BitArray(undef, WORDSIZE)
-    word_to_bitvec!(word, bitvec)
-    return bitvec
-end
-function word_to_bitvec!(word::WordType, bitvec::BitVector)
-    @inbounds for i in 1:WORDSIZE
-        bitvec[i] = has_bit(word, i)
-    end
-    return
-end
-# From https://discourse.julialang.org/t/parse-an-array-of-bits-bitarray-to-an-integer/42361/23
-function bitvec_to_word(b::AbstractVector)
-    w = WordType(0)
-    v = WordType(1)
-
-    @inbounds for j in length(b):-1:1
-        w += v * convert(WordType, b[j])
-        v <<= 1
-    end
-
-    return w
-end
-
-function word_to_idxs(word::WordType, max_nidxs::Integer=WORDSIZE)
-    idxs = zeros(Int, max_nidxs)
-    k = word_to_idxs!(word, idxs)
-    resize!(idxs, k)
-    return idxs
-end
-function word_to_idxs!(word::WordType, idxs::AbstractVector{<:Integer})
-    k = 1
-    @inbounds for i in 1:WORDSIZE
-        if has_bit(word, i)
-            idxs[k] = i
-            k += 1
-        end
-    end
-    return k - 1
-end
-
-function set_to_idxs(set::AbstractVector{WordType}, shift::Bool=false, max_nidxs::Integer=WORDSIZE)
-    n = length(set)
-    idxs = zeros(Int, n * max_nidxs)
-
-    nidx = set_to_idxs!(set, idxs, shift)
-    resize!(idxs, nidx)
-    return idxs
-end
-function set_to_idxs!(set::AbstractVector{WordType}, idxs::AbstractVector{<:Integer}, shift::Bool=false)
-    n = length(set)
-
-    nidx = 0
-    for i in 1:n
-        k = word_to_idxs!(set[i], @view idxs[nidx+1:end])
-        if shift
-            @views idxs[nidx+1:end] .+= (i - 1) * WORDSIZE
-        end
-        nidx += k
-    end
-    return nidx
-end
-
-function _concatbytes(bytes::AbstractVector{<:UInt8})
-    # @assert length(bytes) == sizeof(HashType)
-    w = HashType(0)
-    for b in bytes
-        w |= b
-        w <<= 8
-    end
-    return w
-end
-
-function _to_matrixidx(idx::Integer, m::Integer)
-    return 1 + (idx - 1) ÷ m, mod1(idx, m)
-end
-
-function _to_vecidx(i::Integer, j::Integer, m::Integer)
-    # @assert j <= m
-    return (i - 1) * m + j
-end
diff --git a/src/densenautygraph.jl b/src/densenautygraph.jl
new file mode 100644
index 0000000..8c1ff1d
--- /dev/null
+++ b/src/densenautygraph.jl
@@ -0,0 +1,276 @@
+"""
+    DenseNautyGraph{D,W}
+
+Memory-efficient graph format compatible with nauty. Can be directed (`D = true`) or undirected (`D = false`).
+This graph format stores the adjacency matrix in bit vector form. `W` is the underlying
+unsigned integer type that holds the individual bits (defaults to `UInt`).
+"""
+mutable struct DenseNautyGraph{D,W<:Unsigned} <: AbstractNautyGraph{Int}
+    graphset::Graphset{W}
+    labels::Vector{Int}
+    ne::Int
+    hashval::Union{HashType,Nothing}
+end
+function DenseNautyGraph{D}(graphset::Graphset{W}; vertex_labels=nothing) where {D,W}
+    if !isnothing(vertex_labels) && graphset.n != length(vertex_labels)
+        throw(ArgumentError("The length of `graphset` is not compatible with length of `vertex_labels`. See the nauty user guide for how to correctly construct `graphset`."))
+    end
+    ne = sum(graphset)
+    !D && (ne ÷= 2)
+
+    if isnothing(vertex_labels)
+        vertex_labels = zeros(Int, graphset.n)
+    end
+    return DenseNautyGraph{D,W}(graphset, vertex_labels, ne, nothing)
+end
+
+
+"""
+    DenseNautyGraph{D}(n::Integer; [vertex_labels]) where {D}
+
+Construct a `DenseNautyGraph` on `n` vertices and 0 edges. 
+Can be directed (`D = true`) or undirected (`D = false`).
+Vertex labels can optionally be specified.
+"""
+function DenseNautyGraph{D,W}(n::Integer; vertex_labels=nothing) where {D,W<:Unsigned}
+    graphset = Graphset{W}(n)
+    return DenseNautyGraph{D}(graphset; vertex_labels)
+end
+DenseNautyGraph{D}(n::Integer; vertex_labels=nothing) where {D} = DenseNautyGraph{D,UInt}(n; vertex_labels)
+
+"""
+    DenseNautyGraph{D}(A::AbstractMatrix; [vertex_labels]) where {D}
+
+Construct a `DenseNautyGraph{D}` from the adjacency matrix `A`.
+If `A[i][j] != 0`, an edge `(i, j)` is inserted. `A` must be a square matrix.
+The graph can be directed (`D = true`) or undirected (`D = false`). If `D = false`, `A` must be symmetric.
+Vertex labels can optionally be specified.
+"""
+function DenseNautyGraph{D,W}(A::AbstractMatrix; vertex_labels=nothing) where {D,W<:Unsigned}
+    D || issymmetric(A) || throw(ArgumentError("Adjacency / distance matrices must be symmetric"))
+    graphset = Graphset{W}(A)
+    return DenseNautyGraph{D}(graphset; vertex_labels)
+end
+DenseNautyGraph{D}(A::AbstractMatrix; vertex_labels=nothing) where {D} = DenseNautyGraph{D,UInt}(A; vertex_labels)
+
+function (::Type{G})(g::AbstractGraph) where {G<:AbstractNautyGraph}
+    ng = G(nv(g))
+    for e in edges(g)
+        add_edge!(ng, e)
+        !is_directed(g) && is_directed(ng) && add_edge!(ng, reverse(e))
+    end
+    return ng
+end
+function (::Type{G})(g::AbstractNautyGraph) where {G<:AbstractNautyGraph}
+    h = invoke(G, Tuple{AbstractGraph}, g)
+    @views h.labels .= g.labels
+    return h
+end
+
+"""
+    DenseNautyGraph{D}(edge_list::Vector{<:AbstractEdge}; [vertex_labels]) where {D}
+
+Construct a `DenseNautyGraph` from a vector of edges.
+The number of vertices is the highest that is used in an edge in `edge_list`.
+The graph can be directed (`D = true`) or undirected (`D = false`).
+Vertex labels can optionally be specified.
+"""
+function DenseNautyGraph{D,W}(edge_list::Vector{<:AbstractEdge}; vertex_labels=nothing) where {D,W<:Unsigned}
+    nvg = 0
+    @inbounds for e in edge_list
+        nvg = max(nvg, src(e), dst(e))
+    end
+
+    g = DenseNautyGraph{D,W}(nvg; vertex_labels)
+    for edge in edge_list
+        add_edge!(g, edge)
+    end
+    return g
+end
+DenseNautyGraph{D}(edge_list::Vector{<:AbstractEdge}; vertex_labels=nothing) where {D} = DenseNautyGraph{D,UInt}(edge_list; vertex_labels)
+
+
+Base.copy(g::G) where {G<:DenseNautyGraph} = G(copy(g.graphset), copy(g.labels), g.ne, g.hashval)
+function Base.copy!(dest::G, src::G) where {G<:DenseNautyGraph}
+    copy!(dest.graphset, src.graphset)
+    copy!(dest.labels, src.labels)
+    dest.ne = src.ne
+    dest.hashval = src.hashval
+    return dest
+end
+
+Base.show(io::Core.IO, g::DenseNautyGraph{false}) = print(io, "{$(nv(g)), $(ne(g))} undirected NautyGraph")
+Base.show(io::Core.IO, g::DenseNautyGraph{true}) = print(io, "{$(nv(g)), $(ne(g))} directed NautyGraph")
+
+Base.hash(g::DenseNautyGraph, h::UInt) = hash(g.labels, hash(g.graphset, h))
+Base.:(==)(g::DenseNautyGraph, h::DenseNautyGraph) = (g.graphset == h.graphset) && (g.labels == h.labels)
+
+# BASIC GRAPH API
+labels(g::AbstractNautyGraph) = g.labels
+Graphs.nv(g::DenseNautyGraph) = g.graphset.n
+Graphs.ne(g::DenseNautyGraph) = g.ne
+Graphs.vertices(g::DenseNautyGraph) = Base.OneTo(nv(g))
+Graphs.has_vertex(g::DenseNautyGraph, v) = v ∈ vertices(g)
+function Graphs.has_edge(g::DenseNautyGraph, s::Integer, d::Integer)
+    has_vertex(g, s) && has_vertex(g, s) || return false
+    return g.graphset[s, d]
+end
+function Graphs.outdegree(g::DenseNautyGraph, v::Integer)
+    return sum(adjrow(g, v))
+end
+function Graphs.outneighbors(g::DenseNautyGraph, v::Integer)
+    return findall(adjrow(g, v))
+end
+function adjrow(g::DenseNautyGraph, v::Integer)
+    return @view g.graphset[v, :]
+end
+
+function Graphs.indegree(g::DenseNautyGraph, v::Integer)
+    return sum(adjcol(g, v))
+end
+function Graphs.inneighbors(g::DenseNautyGraph, v::Integer)
+    return findall(adjcol(g, v))
+end
+function adjcol(g::DenseNautyGraph, v::Integer)
+    return @view g.graphset[:, v]
+end
+
+function Graphs.edges(g::DenseNautyGraph)
+    return SimpleEdgeIter(g)
+end
+eltype(::Type{SimpleEdgeIter{<:DenseNautyGraph{true}}}) = Graphs.SimpleGraphEdge{Int}
+eltype(::Type{SimpleEdgeIter{<:DenseNautyGraph{false}}}) = Graphs.SimpleDiGraphEdge{Int}
+function Base.iterate(eit::SimpleEdgeIter{G}, state=(1, 1)) where {G<:DenseNautyGraph}
+    g = eit.g
+    n = nv(g)
+    i0, j0 = state
+
+    jstart = j0
+    for i in i0:n
+        for j in jstart:n
+            g.graphset[i, j] && return Graphs.SimpleEdge(i, j), (i, j+1)
+        end
+        jstart = is_directed(g) ? 1 : i+1
+    end
+    return nothing
+end
+function Base.:(==)(e1::SimpleEdgeIter{<:DenseNautyGraph}, e2::SimpleEdgeIter{<:DenseNautyGraph})
+    g = e1.g
+    h = e2.g
+    ne(g) == ne(h) || return false
+    m = min(nv(g), nv(h))
+
+    g.graphset[1:m, 1:m] == h.graphset[1:m, 1:m] || return false
+    nv(g) == nv(h) && return true
+
+    g.graphset[m+1:end, :] == 0 || return false
+    is_directed(g) || g.graphset[m+1:end, 1:m] == 0 || return false
+
+    h.graphset[m+1:end, :] == 0 || return false
+    is_directed(h) || h.graphset[m+1:end, 1:m] == 0 || return false
+    return true
+end
+function Base.:(==)(e1::SimpleEdgeIter{<:DenseNautyGraph}, e2::SimpleEdgeIter{<:Graphs.SimpleGraphs.AbstractSimpleGraph})
+    g = e1.g
+    h = e2.g
+    ne(g) == ne(h) || return false
+    is_directed(g) == is_directed(h) || return false
+
+    m = min(nv(g), nv(h))
+    for i in 1:m
+        outneighbors(g, i) == Graphs.SimpleGraphs.fadj(h, i) || return false
+        if is_directed(h)
+            inneighbors(g, i) == Graphs.SimpleGraphs.badj(h, i) || return false
+        end
+    end
+    nv(g) == nv(h) && return true
+
+    g.graphset[m+1:end, :] == 0 || return false
+    is_directed(g) || g.graphset[m+1:end, 1:m] == 0 || return false
+
+    for i in (m + 1):nv(h)
+        isempty(Graphs.SimpleGraphs.fadj(h, i)) || return false
+        if is_directed(h)
+            isempty(Graphs.SimpleGraphs.badj(h, i)) || return false
+        end
+    end
+    return true
+end
+Base.:(==)(e1::SimpleEdgeIter{<:Graphs.SimpleGraphs.AbstractSimpleGraph}, e2::SimpleEdgeIter{<:DenseNautyGraph}) = e2 == e1
+
+Graphs.is_directed(::Type{<:DenseNautyGraph{D}}) where {D} = D
+Graphs.edgetype(::AbstractNautyGraph) = Graphs.SimpleGraphs.SimpleEdge{Int}
+Base.eltype(::AbstractNautyGraph{T}) where {T} = T
+Base.zero(::G) where {G<:AbstractNautyGraph} = G(0)
+Base.zero(::Type{G}) where {G<:AbstractNautyGraph} = G(0)
+
+function _induced_subgraph(g::DenseNautyGraph, iter)
+    h, vmap = invoke(Graphs.induced_subgraph, Tuple{AbstractGraph,typeof(iter)}, g, iter)
+    @views h.labels .= g.labels[vmap]
+    return h, vmap
+end
+Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{<:Integer}) = _induced_subgraph(g::DenseNautyGraph, iter)
+Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{Bool}) = _induced_subgraph(g::DenseNautyGraph, iter)
+Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{<:AbstractEdge}) = _induced_subgraph(g::DenseNautyGraph, iter)
+
+# GRAPH MODIFY METHODS
+function Graphs.add_edge!(g::DenseNautyGraph, e::Edge)
+    has_vertex(g, e.src) && has_vertex(g, e.dst) || return false
+    edge_present = g.graphset[e.src, e.dst]
+    edge_present && return false
+
+    g.graphset[e.src, e.dst] = 1
+    is_directed(g) || (g.graphset[e.dst, e.src] = 1)
+    g.ne += 1
+    g.hashval = nothing
+    return true
+end
+Graphs.add_edge!(g::AbstractNautyGraph, i::Integer, j::Integer) = Graphs.add_edge!(g, edgetype(g)(i, j))
+
+function Graphs.rem_edge!(g::DenseNautyGraph, e::Edge)
+    has_vertex(g, e.src) && has_vertex(g, e.dst) || return false
+    edge_present = g.graphset[e.src, e.dst]
+    edge_present || return false
+
+    g.graphset[e.src, e.dst] = 0
+    is_directed(g) || (g.graphset[e.dst, e.src] = 0)
+    g.ne -= 1
+    g.hashval = nothing
+    return true
+end
+Graphs.rem_edge!(g::AbstractNautyGraph, i::Integer, j::Integer) = Graphs.rem_edge!(g, edgetype(g)(i, j))
+
+function Graphs.add_vertices!(g::DenseNautyGraph, n::Integer, vertex_labels=0)
+    vertex_labels isa Number || n != length(vertex_labels) && throw(ArgumentError("Incompatible length: trying to add `n` vertices, but`vertex_labels` has length $(length(vertex_labels))."))
+    ng = nv(g)
+    _add_vertices!(g.graphset, n)
+    resize!(g.labels, ng + n)
+    g.labels[ng+1:end] .= vertex_labels
+    g.hashval = nothing
+    return n
+end
+Graphs.add_vertex!(g::DenseNautyGraph, vertex_label::Integer=0) = Graphs.add_vertices!(g, 1, vertex_label) > 0
+
+function Graphs.rem_vertices!(g::DenseNautyGraph, inds)
+    all(i->has_vertex(g, i), inds) || return false
+
+    _rem_vertices!(g.graphset, inds)
+
+    ne = sum(g.graphset)
+    is_directed(g) || (ne ÷= 2)
+    g.ne = ne
+
+    g.hashval = nothing
+    return true
+end
+Graphs.rem_vertex!(g::DenseNautyGraph, i::Integer) = rem_vertices!(g, (i,))
+
+function Graphs.blockdiag(g::DenseNautyGraph{D1,W}, h::DenseNautyGraph{D2}) where {D1,D2,W}
+    ng, nh = nv(g), nv(h)
+
+    gset = Graphset{wordtype(g.graphset)}(ng+nh)
+    gset[1:ng, 1:ng] .= g.graphset
+    gset[ng+1:end, ng+1:end] .= h.graphset
+    D = D1 || D2
+    return DenseNautyGraph{D,W}(gset; vertex_labels=vcat(g.labels, h.labels))
+end
\ No newline at end of file
diff --git a/src/densenautygraphs.jl b/src/densenautygraphs.jl
deleted file mode 100644
index ceda4c8..0000000
--- a/src/densenautygraphs.jl
+++ /dev/null
@@ -1,347 +0,0 @@
-abstract type AbstractNautyGraph <: AbstractGraph{Cint} end
-# TODO: abstract type AbstractSparseNautyGraph <: AbstractNautyGraph end
-
-mutable struct DenseNautyGraph{D} <: AbstractNautyGraph
-    graphset::Vector{WordType}
-    n_vertices::Cint
-    n_edges::Cint
-    n_words::Cint
-    labels::Vector{Cint}
-    hashval::Union{HashType,Nothing}
-
-    function DenseNautyGraph(graphset, labels, directed::Bool)
-        n_vertices = length(labels)
-        n_words = n_vertices > 0 ? length(graphset) ÷ n_vertices : 1
-
-        if length(graphset) != n_words * n_vertices
-            error("length of `graphset` is not compatible with length of `labels`. See the nauty user guide for how to correctly construct `graphset`.")
-        end
-
-        n_edges = 0
-        for word in graphset
-            n_edges += count_bits(word)
-        end
-        if !directed
-            n_edges ÷= 2
-        end
-        return new{directed}(graphset, n_vertices, n_edges, n_words, labels, nothing)
-    end
-    function DenseNautyGraph{D}(graphset, n_vertices, n_edges, n_words, labels, hashval) where {D}
-        return new{D}(graphset, n_vertices, n_edges, n_words, labels, hashval)
-    end
-end
-
-function DenseNautyGraph{D}(n::Integer, vertex_labels=nothing) where {D}
-    if !isnothing(vertex_labels) && n != length(vertex_labels)
-        error("Incompatible length: `vertex_labels` has length $(length(vertex_labels)) instead of `n=$n`.")
-    end
-    m = ceil(Cint, n / WORDSIZE)
-    graphset = zeros(WordType, Int(n * m))
-    labels = isnothing(vertex_labels) ? zeros(Cint, n) : convert(Vector{Cint}, vertex_labels)
-    return DenseNautyGraph(graphset, labels, D)
-end
-function DenseNautyGraph{D}(adjmx::AbstractMatrix, vertex_labels=nothing) where {D}
-    n, n2 = size(adjmx)
-
-    if !D
-        !issymmetric(adjmx) && error("Cannot create an undirected NautyGraph from a non-symmetric adjacency matrix. Make sure the adjacency matrix is square symmetric!")
-    else
-        n != n2 && error("Cannot create a NautyGraph from a rectangular matrix. Make sure the adjacency matrix is square!")
-    end
-    
-    graphset = _adjmatrix_to_graphset(adjmx)
-    labels = isnothing(vertex_labels) ? zeros(Cint, n) : convert(Vector{Cint}, vertex_labels)
-
-    return DenseNautyGraph(graphset, labels, D)
-end
-
-function (::Type{G})(g::AbstractGraph, vertex_labels=nothing) where {G<:AbstractNautyGraph}
-    if is_directed(g) != is_directed(G)
-        error("Cannot create an undirected NautyGraph from a directed graph (or vice versa). Please make sure the directedness of the graph types is matching.")
-    end
-
-    ng = G(nv(g), vertex_labels)
-
-    for e in edges(g)
-        add_edge!(ng, e)
-    end
-    return ng
-end
-
-Base.copy(g::G) where {G<:DenseNautyGraph} = G(copy(g.graphset), g.n_vertices, g.n_edges, g.n_words, copy(g.labels), g.hashval)
-function Base.copy!(dest::G, src::G) where {G<:DenseNautyGraph}
-    copy!(dest.graphset, src.graphset)
-    copy!(dest.labels, src.labels)
-    dest.n_vertices = src.n_vertices
-    dest.n_edges = src.n_edges
-    dest.n_words = src.n_words
-    dest.hashval = src.hashval
-    return dest
-end
-
-Base.show(io::Core.IO, g::DenseNautyGraph{false}) = print(io, "{$(nv(g)), $(ne(g))} undirected NautyGraph")
-Base.show(io::Core.IO, g::DenseNautyGraph{true}) = print(io, "{$(nv(g)), $(ne(g))} directed NautyGraph")
-
-Base.hash(g::DenseNautyGraph, h::UInt) = hash(g.labels, hash(g.graphset, h))
-Base.:(==)(g::DenseNautyGraph, h::DenseNautyGraph) = (g.graphset == h.graphset) && (g.labels == h.labels)
-
-begin # BASIC GRAPH API
-    labels(g::AbstractNautyGraph) = g.labels
-    Graphs.nv(g::DenseNautyGraph) = g.n_vertices
-    Graphs.ne(g::DenseNautyGraph) = g.n_edges
-    Graphs.vertices(g::DenseNautyGraph) = Base.OneTo(nv(g))
-    Graphs.has_vertex(g::DenseNautyGraph, v) = v ∈ vertices(g)
-    function Graphs.has_edge(g::DenseNautyGraph, s::Integer, d::Integer)
-        i_row = s
-        i_col = 1 + (d - 1) ÷ WORDSIZE
-        k = mod1(d, WORDSIZE)
-        i_vec = _to_vecidx(i_row, i_col, g.n_words)
-        return has_bit(g.graphset[i_vec], k)
-    end
-    function Graphs.outdegree(g::DenseNautyGraph, v::Integer)
-        m = g.n_words
-        i = _to_vecidx(v, 1, m)
-        return sum(count_bits, @view g.graphset[i:i+m-1])
-    end
-    function Graphs.outneighbors(g::DenseNautyGraph, v::Integer)
-        neighs = zeros(Int, nv(g))
-        k = outneighbors!(neighs, g, v)
-        resize!(neighs, k)
-        return neighs
-    end
-    function outneighbors!(neighs::AbstractVector{<:Integer}, g::DenseNautyGraph, v::Integer)
-        m = g.n_words
-        i = _to_vecidx(v, 1, m)
-        return set_to_idxs!((@view g.graphset[i:i+m-1]), neighs, true)
-    end
-    function Graphs.indegree(g::DenseNautyGraph, v::Integer)
-        m = g.n_words
-        i_col = 1 + (v - 1) ÷ WORDSIZE
-        vmod = mod1(v, WORDSIZE)
-
-        idx_convert(i) = _to_vecidx(i, i_col, m)
-        counter(i) = has_bit(g.graphset[idx_convert(i)], vmod)
-        return sum(counter, 1:nv(g))
-    end
-    function Graphs.inneighbors(g::DenseNautyGraph, v::Integer)
-        i_col = 1 + (v - 1) ÷ WORDSIZE
-        idxs = _to_vecidx.(1:nv(g), Ref(i_col), Ref(g.n_words))
-
-        neighs = zeros(Cint, nv(g))
-        k = 1
-        vmod = mod1(v, WORDSIZE)
-        for (i, word) in enumerate(@view g.graphset[idxs])
-            if has_bit(word, vmod)
-                neighs[k] = i
-                k += 1
-            end
-        end
-        resize!(neighs, k - 1)
-        return neighs
-    end
-    function Graphs.edges(g::DenseNautyGraph{true})
-        #TODO: define edgeiterator
-        return _directed_edges(g)
-    end
-    function Graphs.edges(g::DenseNautyGraph{false})
-        edges = _directed_edges(g)
-        filter!(e -> e.src <= e.dst, edges) #TODO: optimize
-        return edges
-    end
-
-    Graphs.is_directed(::Type{<:DenseNautyGraph{D}}) where {D} = D
-    Graphs.edgetype(::AbstractNautyGraph) = Edge{Cint}
-    Base.eltype(::AbstractNautyGraph) = Cint
-    Base.zero(::G) where {G<:AbstractNautyGraph} = G(0)
-    Base.zero(::Type{G}) where {G<:AbstractNautyGraph} = G(0)
-
-    function _induced_subgraph(g::DenseNautyGraph, iter)
-        h, vmap = invoke(Graphs.induced_subgraph, Tuple{AbstractGraph,typeof(iter)}, g, iter)
-        @views h.labels .= g.labels[vmap]
-        return h, vmap
-    end
-
-    Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{<:Integer}) = _induced_subgraph(g::DenseNautyGraph, iter)
-    Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{Bool}) = _induced_subgraph(g::DenseNautyGraph, iter)
-    Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{<:AbstractEdge}) = _induced_subgraph(g::DenseNautyGraph, iter)
-end
-
-begin # GRAPH MODIFY METHODS
-    function Graphs.add_edge!(g::DenseNautyGraph{true}, e::Edge)
-        edge_added = _modify_edge!(g, e, true)
-        if edge_added
-            g.n_edges += 1
-            g.hashval = nothing
-        end
-        return edge_added
-    end
-    function Graphs.add_edge!(g::DenseNautyGraph{false}, e::Edge)
-        fwd_edge_added = _modify_edge!(g, e, true)
-
-        if e.src != e.dst
-            bwd_edge_added = _modify_edge!(g, reverse(e), true)
-            edge_added = fwd_edge_added && bwd_edge_added
-        else
-            edge_added = fwd_edge_added
-        end
-
-        if edge_added
-            g.n_edges += 1
-            g.hashval = nothing
-        end
-        return edge_added
-    end
-    Graphs.add_edge!(g::AbstractNautyGraph, i::Integer, j::Integer) = Graphs.add_edge!(g, Graphs.Edge{Cint}(i, j))
-
-    function Graphs.rem_edge!(g::DenseNautyGraph{true}, e::Edge)
-        edge_removed = _modify_edge!(g, e, false)
-        if edge_removed
-            g.n_edges -= 1
-            g.hashval = nothing
-        end
-        return edge_removed
-    end
-    function Graphs.rem_edge!(g::DenseNautyGraph{false}, e::Edge)
-        fwd_edge_removed = _modify_edge!(g, e, false)
-        bwd_edge_removed = _modify_edge!(g, reverse(e), false)
-        edge_removed = fwd_edge_removed && bwd_edge_removed
-
-        if edge_removed
-            g.n_edges -= 1
-            g.hashval = nothing
-        end
-        return edge_removed
-    end
-    Graphs.rem_edge!(g::AbstractNautyGraph, i::Integer, j::Integer) = Graphs.rem_edge!(g, Graphs.Edge{Cint}(i, j))
-
-    function Graphs.add_vertex!(g::DenseNautyGraph, vertex_label::Union{<:Integer,Nothing}=nothing)
-        if isnothing(vertex_label)
-            labeled = !all(iszero, g.labels)
-            labeled && error("Cannot add an unlabeled vertex to a labeled nautygraph. Use `add_vertex!(g, label)` to add a labeled vertex.")
-        end
-
-        n = nv(g)
-        m = g.n_words
-
-        if n + 1 > m * WORDSIZE
-            m = g.n_words += 1
-            for i in m:m:n*m
-                insert!(g.graphset, i, zero(WordType))
-            end
-        end
-
-        append!(g.graphset, fill(zero(WordType), m))
-
-        if isnothing(vertex_label)
-            push!(g.labels, zero(Cint))
-        else
-            push!(g.labels, convert(Cint, vertex_label))
-        end
-
-        g.hashval = nothing
-        g.n_vertices += 1
-        return true
-    end
-    function Graphs.add_vertices!(g::AbstractNautyGraph, n::Integer, vertex_labels::Union{AbstractVector,Nothing}=nothing)
-        if isnothing(vertex_labels)
-            return sum(_->add_vertex!(g), 1:n)
-        end
-
-        n != length(vertex_labels) && error("Incompatible length: `vertex_labels` has length $(length(vertex_labels)) instead of `n=$n`.")
-        for l in vertex_labels
-            add_vertex!(g, l)
-        end
-        return n
-    end
-
-    function Graphs.rem_vertices!(g::DenseNautyGraph, inds::AbstractVector{<:Integer})
-        n = nv(g)
-        sort!(inds)
-        if last(inds) > n
-            return false
-        end
-        
-        m = g.n_words
-
-        i_vecs = sort!(_to_vecidx.(inds, Ref(1), m))
-        deleteat!(g.graphset, Iterators.flatten(i:i+m-1 for i in i_vecs))
-        # Shift all bits such that the vertices are renamed correctly
-        d = 0
-        for i in inds
-            i -= d
-
-            i_col = 1 + (i - 1) ÷ WORDSIZE
-            k = i - (i_col - 1) * WORDSIZE
-            for j in eachindex(g.graphset)
-                _, j_col = _to_matrixidx(j, m)
-
-                if j_col < i_col
-                    continue
-                end
-                offset = j_col == i_col ? k : 0
-
-                next_word = j_col != m ? g.graphset[j+1] : zero(WordType)
-                g.graphset[j] = push_bits(g.graphset[j], next_word, WORDSIZE - offset, 1)
-            end
-
-            d += 1
-        end
-
-        deleteat!(g.labels, inds)
-
-        g.hashval = nothing
-        g.n_vertices -= length(inds)
-
-        # Remove extra words if they are no longer needed
-        # TODO: make this optional
-        m_new = ceil(Cint, g.n_vertices / WORDSIZE)
-        if m_new < m
-            deleteat!(g.graphset, Iterators.flatten(i+m_new:i+m-1 for i in 1:m:g.n_vertices*m))
-            g.n_words = m_new
-        end
-
-        # Recount edges
-        n_edges = 0
-        for word in g.graphset
-            n_edges += count_bits(word)
-        end
-        if !is_directed(g)
-            n_edges ÷= 2
-        end
-        g.n_edges = n_edges
-        return true
-    end
-    Graphs.rem_vertex!(g::DenseNautyGraph, i::Integer) = rem_vertices!(g, [i])
-end
-
-function Graphs.blockdiag(g::G, h::G) where {G<:DenseNautyGraph}
-    ng, nh = g.n_vertices, h.n_vertices
-    vl = vcat(g.labels, h.labels)
-
-    k = G(Int(ng + nh), vl)
-
-    transfer_set!(k.graphset, g.graphset, 0, k.n_words, g.n_words)
-    transfer_set!(k.graphset, h.graphset, ng * k.n_words * WORDSIZE + ng, k.n_words, h.n_words)
-    k.n_edges = g.n_edges + h.n_edges
-    return k
-end
-
-"""
-    blockdiag!(g::G, h::G) where {G<:DenseNautyGraph}
-
-Compute `blockdiag(g, h)` and store it in `g`, whose size is increased to accomodate it.
-"""
-function blockdiag!(g::G, h::G) where {G<:DenseNautyGraph}
-    @assert g !== h # Make sure g and h don't alias the same memory.
-
-    for i in vertices(h)
-        add_vertex!(g, h.labels[i])
-    end
-    ng, mg = g.n_vertices, g.n_words
-    nh, mh = h.n_vertices, h.n_words
-
-    transfer_set!(g.graphset, h.graphset, (ng - nh) * mg * WORDSIZE + (ng - nh), mg, mh)
-    g.n_edges += h.n_edges
-    g.hashval = nothing
-    return g
-end
diff --git a/src/graphset.jl b/src/graphset.jl
new file mode 100644
index 0000000..1d8162b
--- /dev/null
+++ b/src/graphset.jl
@@ -0,0 +1,165 @@
+"""
+    Graphset{W}
+
+A graphset is a special bit matrix used to represent the adjacency matrix 
+of a nauty graph in dense format. For a graph on `n` vertices, the graphset
+contains `n*m` "words", i.e. unsigned integers that contain the bits of 
+the adjacency matrix, where `m` is the number of words per vertex.
+
+The organization of words is as follows:
+
+   ------------ m words per vertex ----->  
+   | 0x00000000, 0x00000000, 0x00000000, ...  
+ n | 0x00000000, 0x00000000, 0x00000000, ...  
+   | 0x00000000, 0x00000000, 0x00000000, ...  
+   v  
+"""
+mutable struct Graphset{W<:Unsigned} <: AbstractMatrix{Bool}
+    words::Vector{W}
+    n::Int
+    m::Int
+
+    function Graphset{W}(n::Integer, m=cld(n, wordsize(W))) where {W}
+        if n > m * wordsize(W)
+            throw(ArgumentError("Not enough words to hold n=$n vertices. Increase m or use a larger word type."))
+        end
+        words = zeros(W, n*m)
+        return new{W}(words, n, m)
+    end
+end
+function Graphset{W}(A::AbstractMatrix, m=cld(size(A,1), wordsize(W))) where {W}
+    n1, n2 = size(A)
+    n1 == n2 || throw(ArgumentError("Adjacency / distance matrices must be square"))
+    gset = Graphset{W}(n1, m)
+    gset .= A
+    return gset
+end
+Graphset(args...) = Graphset{UInt}(args...)
+
+Base.size(gset::Graphset) = (gset.n, gset.n)
+Base.IndexStyle(::Type{Graphset}) = IndexCartesian()
+Base.similar(gset::Graphset{W}) where {W} = Graphset{W}(gset.n, gset.m)
+
+function Base.copy!(dest::Graphset{W}, src::Graphset{W}) where {W}
+    dest.n == src.n || throw(ArgumentError("graphsets must have the same size for copy!"))
+    if dest.m == src.m
+        copyto!(dest.words, src.words)
+    else
+        dest .== src
+    end
+    return dest
+end
+
+@inline wordtype(::Graphset{W}) where {W} = W
+@inline wordsize(u::Unsigned) = 8 * sizeof(u)
+@inline wordsize(T::Type{<:Unsigned}) = 8 * sizeof(T)
+@inline wordsize(::Graphset{W}) where {W} = wordsize(W)
+@inline logwordsize(::Type{UInt16}) = 4
+@inline logwordsize(::Type{UInt32}) = 5
+@inline logwordsize(::Type{UInt64}) = 6
+@inline logwordsize(::Type{UInt128}) = 7
+@inline logwordsize(::Graphset{W}) where {W} = logwordsize(W)
+
+@inline mod1pow2(p2, x) = 1 + (x - 1) & (p2 - 1)
+@inline divlogpow2(lp2, x) = x >> lp2
+
+@inline function bitaddress(gset::Graphset, i, j)
+    ws = wordsize(gset)
+    lws = logwordsize(gset)
+    wordidx = (i - 1) * gset.m + 1 + divlogpow2(lws, (j - 1))
+    bitidx = mod1pow2(ws, j)
+    return wordidx, bitidx
+end
+
+@inline function getbit(word::W, i::Integer) where {W<:Unsigned}
+    ws = wordsize(W)
+    @boundscheck checkindex(Bool, Base.OneTo(ws), i) || throw(BoundsError(word, i))
+    mask = one(W) << (ws - i)
+    return (mask & word) != zero(W)
+end
+
+@inline function setbit(word::W, x::Bool, i::Integer) where {W<:Unsigned}
+    ws = wordsize(W)
+    @boundscheck checkindex(Bool, Base.OneTo(ws), i) || throw(BoundsError(word, i))
+    mask = one(W) << (ws - i)
+    return ifelse(x, word | mask, word & ~mask)
+end
+
+@inline function Base.getindex(gset::Graphset, inds::Vararg{Int,2})
+    i, j = inds
+    @boundscheck checkbounds(gset, i, j)
+    wordidx, bitidx = bitaddress(gset, i, j)
+    word = gset.words[wordidx]
+    return getbit(word, bitidx)
+end
+
+@inline function Base.setindex!(gset::Graphset, x, inds::Vararg{Int,2})
+    i, j = inds
+    @boundscheck checkbounds(gset, i, j)
+    wordidx, bitidx = bitaddress(gset, i, j)
+    gset.words[wordidx] = setbit(gset.words[wordidx], convert(Bool, x), bitidx)
+    return gset
+end
+
+function _increase_padding!(gset::Graphset{W}, m::Integer=1) where {W}
+    # TODO: optimize this
+    for _ in Base.OneTo(m)
+        gset.m += 1
+        for i in Base.OneTo(gset.n)
+            insert!(gset.words, i * gset.m, zero(W))
+        end
+    end
+    return gset
+end
+
+@inline function partial_leftshift(word::Unsigned, n::Integer, start::Integer, fillword::Unsigned=zero(word))
+    # Starting from the `start`th bit from the left of `word`, shift all bits to the left `n` times,
+    # refill right-most bits with bits from `fillword`. The `n` bits to the left of `start` are
+    # overwritten.
+
+    # Select all to-be-moved bits
+    mask = typemax(word) >> (start - 1)
+
+    unshifted_part = word & (~mask << n)
+    shifted_part = (word & mask) << n
+    filled_part = fillword >> (wordsize(typeof(word)) - n)
+
+    return unshifted_part | shifted_part | filled_part
+end
+
+function _add_vertices!(gset::Graphset{W}, n::Integer) where {W} # TODO think of a better name
+    _increase_padding!(gset, cld(gset.n + n, wordsize(gset)) - gset.m)
+    append!(gset.words, fill(zero(W), n*gset.m))
+    gset.n += n
+    return gset
+end
+_add_vertex!(gset::Graphset) = _add_vertices!(gset, 1)
+
+function _rem_vertices!(gset::Graphset{W}, inds) where {W}
+    nrv = length(inds)
+
+    deleteat!(gset.words, Iterators.flatten(1+(i-1)*gset.m:i*gset.m for i in inds))
+    gset.n -= nrv
+
+    wordblock = reshape(gset.words, gset.m, gset.n)'
+
+    δ = 0
+    lastind = 0
+    for ind in inds
+        ind < lastind && throw(ArgumentError("indices must be unique and sorted"))
+
+        wordidx, bitidx = bitaddress(gset, 1, ind - δ)
+        for i in axes(wordblock, 1)
+            fillword = wordidx == size(wordblock, 2) ? zero(W) : wordblock[i, wordidx+1]
+            wordblock[i, wordidx] = partial_leftshift(wordblock[i, wordidx], 1, bitidx+1, fillword)
+
+            for j in wordidx+1:size(wordblock, 2)
+                fillword = j == size(wordblock, 2) ? zero(W) : wordblock[i, j+1]
+                wordblock[i, j] = partial_leftshift(wordblock[i, j], 1, 1, fillword)
+            end
+        end
+        δ += 1
+    end
+    return gset
+end
+_rem_vertex!(gset::Graphset, i::Integer) = _rem_vertices!(gset, (i,))
\ No newline at end of file
diff --git a/src/nauty.jl b/src/nauty.jl
index f860abc..a8414df 100644
--- a/src/nauty.jl
+++ b/src/nauty.jl
@@ -1,3 +1,8 @@
+libnauty(::Type{UInt16}) = nauty_jll.libnautyTS
+libnauty(::Type{UInt32}) = nauty_jll.libnautyTW
+libnauty(::Type{UInt64}) = nauty_jll.libnautyTL
+libnauty(::DenseNautyGraph{D,W}) where {D,W} = libnauty(W)
+
 mutable struct NautyOptions
     getcanon::Cint # Warning: setting getcanon to false means that nauty will NOT compute the canonical representative, which may lead to unexpected results.
     digraph::Cbool # This needs to be true if the graph is directed or has loops. Disabling this option for undirected graphs with no loops may increase performance.
@@ -25,16 +30,25 @@ mutable struct NautyOptions
     schreier::Cbool
     extra_options::Ptr{Cvoid}
 
-    NautyOptions(; digraph_or_loops=true, ignorelabels=false, groupinfo=false) = new(
-        1, digraph_or_loops, groupinfo, false, ignorelabels, false, 78,
-        C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL,
-        100, 0, 1, 0,
-        cglobal((:dispatch_graph, libnauty), Cvoid),
-        false, C_NULL
+    function NautyOptions(dispatch_pointer::Ptr{Cvoid}; digraph_or_loops, ignorelabels, groupinfo)
+        return new(1, digraph_or_loops, groupinfo, false, ignorelabels, false, 78,
+                C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL,
+                100, 0, 1, 0,
+                dispatch_pointer,
+                false, C_NULL
     )
+    end
+end
+@generated function NautyOptions(::Type{W}; digraph_or_loops=true, ignorelabels=false, groupinfo=false) where {W}
+    return :(NautyOptions(cglobal((:dispatch_graph, $(libnauty(W))), Cvoid); digraph_or_loops, ignorelabels, groupinfo))
 end
 
-const DEFAULT_OPTIONS = NautyOptions(digraph_or_loops=true)
+const DEFAULTOPTIONS16 = NautyOptions(C_NULL; digraph_or_loops=true, ignorelabels=false, groupinfo=false)
+const DEFAULTOPTIONS32 = NautyOptions(C_NULL; digraph_or_loops=true, ignorelabels=false, groupinfo=false)
+const DEFAULTOPTIONS64 = NautyOptions(C_NULL; digraph_or_loops=true, ignorelabels=false, groupinfo=false)
+default_options(::DenseNautyGraph{D,UInt16}) where {D} = DEFAULTOPTIONS16
+default_options(::DenseNautyGraph{D,UInt32}) where {D} = DEFAULTOPTIONS32
+default_options(::DenseNautyGraph{D,UInt64}) where {D} = DEFAULTOPTIONS64
 
 mutable struct NautyStatistics
     grpsize1::Cdouble
@@ -62,30 +76,34 @@ struct AutomorphismGroup
     # generators::Vector{Vector{Cint}} #TODO: not implemented
 end
 
-function _densenauty(g::DenseNautyGraph, options::NautyOptions=DEFAULT_OPTIONS, statistics::NautyStatistics=NautyStatistics())
+function _densenauty(g::DenseNautyGraph{D,W}, options::NautyOptions=default_options(g), statistics::NautyStatistics=NautyStatistics()) where {D,W}
     # TODO: allow the user to pass pre-allocated arrays for lab, ptn, orbits, canong in a safe way.
-    n, m = g.n_vertices, g.n_words
+    n, m = g.graphset.n, g.graphset.m
 
-    lab, ptn = _vertexlabels_to_labptn(g.labels)
+    lab, ptn = vertexlabels2labptn(g.labels)
     orbits = zeros(Cint, n)
-    canong = zero(g.graphset)
+    canong = Graphset{W}(n, m)
 
-    @ccall libnauty.densenauty(
-        g.graphset::Ref{WordType},
+    _ccall_densenauty(g, lab, ptn, orbits, options, statistics, canong)
+
+    canonperm = (lab .+= 1)
+    return canong, canonperm, orbits, statistics
+end
+
+@generated function _ccall_densenauty(g::DenseNautyGraph{D,W}, lab, ptn, orbits, options, statistics, canong) where {D,W}
+    return quote @ccall $(libnauty(W)).densenauty(
+        g.graphset.words::Ref{W},
         lab::Ref{Cint},
         ptn::Ref{Cint},
         orbits::Ref{Cint},
         Ref(options)::Ref{NautyOptions},
         Ref(statistics)::Ref{NautyStatistics},
-        m::Cint,
-        n::Cint,
-        canong::Ref{WordType})::Cvoid
-
-    canonperm = (lab .+= 1)
-    return canong, canonperm, orbits, statistics
+        g.graphset.m::Cint,
+        g.graphset.n::Cint,
+        canong.words::Ref{W})::Cvoid end
 end
 
-function _sethash!(g::DenseNautyGraph, canong, canonperm)
+function _sethash!(g::DenseNautyGraph, canong::Graphset, canonperm)
     # Base.hash skips elements in arrays of length >= 8192
     # Use SHA in these cases
     canong_hash = length(canong) >= 8192 ? hash_sha(canong) : hash(canong)
@@ -95,8 +113,8 @@ function _sethash!(g::DenseNautyGraph, canong, canonperm)
     g.hashval = hashval
     return
 end
-function _canonize!(g::DenseNautyGraph, canong, canonperm)
-    copyto!(g.graphset, canong)
+function _canonize!(g::DenseNautyGraph, canong::Graphset, canonperm)
+    copy!(g.graphset, canong)
     permute!(g.labels, canonperm)
     return
 end
@@ -109,7 +127,7 @@ Compute a graph `g`'s canonical form and automorphism group.
 """
 function nauty(::AbstractNautyGraph, ::NautyOptions; kwargs...) end
 
-function nauty(g::DenseNautyGraph, options::NautyOptions=DEFAULT_OPTIONS; canonize=false, compute_hash=true)
+function nauty(g::DenseNautyGraph, options::NautyOptions=default_options(g); canonize=false, compute_hash=true)
     if is_directed(g) && !isone(options.digraph)
         error("Nauty options need to match the directedness of the input graph. Make sure to instantiate options with `digraph=true` if the input graph is directed.")
     end
diff --git a/src/utils.jl b/src/utils.jl
index ba89904..afb2932 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,87 +1,32 @@
-function _adjmatrix_to_graphset(A::AbstractMatrix{<:Integer})
-    n, _n = size(A)
-
-    m = ceil(Cint, n / WORDSIZE)
-
-    G = zeros(Cint, n, m * WORDSIZE)
-    G[:, 1:n] .= A
-    # Reshape the padded adjacency matrix from 
-    # (a b c)
-    # (d f g)
-    # to
-    # (a b c d f g)'
-    G = reshape(G', Int(WORDSIZE), Int(n * m))'
-
-    graphset = zeros(WordType, Int(n * m))
-    for i in eachindex(graphset)
-        graphset[i] = bitvec_to_word(@view G[i, :])
-    end
-    return graphset
-end
-
-function _directed_edges(g::DenseNautyGraph)
-    n, m = g.n_vertices, g.n_words
-
-    edges = Edge{Cint}[]
-    js = zeros(Cint, WORDSIZE)
-    for set_idx in eachindex(g.graphset)
-        i = 1 + (set_idx - 1) ÷ m
-        j0 = ((set_idx - 1) % m) * WORDSIZE
-
-        k = word_to_idxs!(g.graphset[set_idx], js)
-        js .+= j0
-        for j in 1:k
-            push!(edges, Edge{Cint}(i, js[j]))
-        end
-    end
-
-    return edges
-end
-
-function _vertexlabels_to_labptn(labels::Vector{<:Integer})
+function vertexlabels2labptn(labels::Vector{<:Integer})
     n = length(labels)
     lab = zeros(Cint, n)
     ptn = zeros(Cint, n)
-    return _vertexlabels_to_labptn!(lab, ptn, labels)
+    return vertexlabels2labptn!(lab, ptn, labels)
 end
-function _vertexlabels_to_labptn!(lab::Vector{<:Integer}, ptn::Vector{<:Integer}, labels::Vector{<:Integer})
+function vertexlabels2labptn!(lab::Vector{<:Integer}, ptn::Vector{<:Integer}, labels::Vector{<:Integer})
     lab .= 1:length(labels)
-    sort!(lab, alg=QuickSort, by=k -> labels[k])
+    sort!(lab, alg=QuickSort, by=k->labels[k])
     @views lab .-= 1
 
     for i in 1:length(lab)-1
-        ptn[i] = labels[lab[i+1]+1] == labels[lab[i]+1] ? 1 : 0
+        ptn[i] = ifelse(labels[lab[i+1]+1] == labels[lab[i]+1], 1, 0)
     end
     return lab, ptn
 end
 
-function _modify_edge!(g::AbstractNautyGraph, e::Edge, add::Bool)
-    # Adds or removes the edge e
-    n, m = g.n_vertices, g.n_words
-    i, j = e.src, e.dst
-    if i > n || j > n
-        return false
-    end
-
-    set_idx = 1 + (i - 1) * m + (j - 1) ÷ WORDSIZE
-    # left = 1000000000.... 
-    left = one(WordType) << (WORDSIZE - 1)
-    new_edge = left >> ((j - 1) % WORDSIZE)
-
-    word_old = g.graphset[set_idx]
-    if add
-        g.graphset[set_idx] |= new_edge
-        counter = +1
-    else
-        g.graphset[set_idx] &= ~new_edge
-        counter = -1
-    end
-
-    return g.graphset[set_idx] != word_old
-end
-
 function hash_sha(x)
     io = IOBuffer()
     write(io, x)
-    return _concatbytes(sha256(take!(io))[1:8])
-end
\ No newline at end of file
+    return concatbytes(@view sha256(take!(io))[1:8])
+end
+
+function concatbytes(W::Type{<:Unsigned}, bytes)
+    w = zero(W)
+    for b in bytes
+        w |= b
+        w <<= 8
+    end
+    return w
+end
+concatbytes(bytes) = concatbytes(UInt64, bytes)
\ No newline at end of file
diff --git a/test/aqua.jl b/test/aqua.jl
new file mode 100644
index 0000000..0b54ba3
--- /dev/null
+++ b/test/aqua.jl
@@ -0,0 +1,5 @@
+using Aqua
+
+@testset "aqua" begin 
+   Aqua.test_all(NautyGraphs)
+end
diff --git a/test/densenautygraphs.jl b/test/densenautygraph.jl
similarity index 65%
rename from test/densenautygraphs.jl
rename to test/densenautygraph.jl
index d814d6f..8d3d25d 100644
--- a/test/densenautygraphs.jl
+++ b/test/densenautygraph.jl
@@ -2,34 +2,50 @@ rng = Random.Random.MersenneTwister(0) # Use MersenneTwister for Julia 1.6 compa
 symmetrize_adjmx(A) = (A = convert(typeof(A), (A + A') .> 0); for i in axes(A, 1); A[i, i] = 0; end; A)
 
 @testset "modify" begin
-    nverts = [5, 10, 20, 31, 32, 33, 50, 63, 64, 65, 100, 200, 500, 1000]
-    rvs = [sort(unique(rand(rng, 1:i, 4))) for i in nverts]
-    As = [symmetrize_adjmx(rand(rng, [0, 1], i, i)) for i in nverts]
-
-    gs = [Graph(A) for A in As]
-    ngs = [NautyGraph(A) for A in As]
+    nverts = [1, 2, 3, 4, 5, 10, 20, 31, 32, 33, 50, 63, 64, 
+              65, 100, 122, 123, 124, 125, 126, 200, 500, 1000]
+    As = [rand(rng, [0, 1], i, i) for i in nverts]
+
+    wtypes = [UInt16, UInt32, UInt64]
+
+    gs = []
+    ngs = []
+    for A in As, wt in wtypes
+        Asym = symmetrize_adjmx(A)
+        push!(gs, Graph(Asym))
+        push!(gs, DiGraph(A))
+        push!(ngs, NautyGraph{wt}(Asym))
+        push!(ngs, NautyDiGraph{wt}(A))
+    end
 
-    for (g, ng, rv) in zip(gs, ngs, rvs)
+    for (g, ng) in zip(gs, ngs)
         g, ng = copy(g), copy(ng)
 
         @test adjacency_matrix(g) == adjacency_matrix(ng)
+        @test edges(ng) == edges(g)
+        @test collect(edges(g)) == collect(edges(ng))
+
+        rv = sort(unique(rand(rng, 1:nv(ng), 4)))
 
         rem_vertices!(g, rv, keep_order=true)
         rem_vertices!(ng, rv)
         @test adjacency_matrix(g) == adjacency_matrix(ng)
     end
 
-    for (g, ng, rv) in zip(gs, ngs, rvs)
+    for (g, ng) in zip(gs, ngs)
         g, ng = copy(g), copy(ng)
 
-        edge = collect(edges(g))[end]
+        es = edges(g)
+        if !isempty(es)
+            edge = last(collect(es))
 
-        rem_edge!(g, edge)
-        rem_edge!(ng, edge)
-        @test adjacency_matrix(g) == adjacency_matrix(ng)
+            rem_edge!(g, edge)
+            rem_edge!(ng, edge)
+            @test adjacency_matrix(g) == adjacency_matrix(ng)
+        end
     end
 
-    for (g, ng, rv) in zip(gs, ngs, rvs)
+    for (g, ng) in zip(gs, ngs)
         g, ng = copy(g), copy(ng)
 
         add_vertex!(g)
@@ -39,7 +55,7 @@ symmetrize_adjmx(A) = (A = convert(typeof(A), (A + A') .> 0); for i in axes(A, 1
         @test adjacency_matrix(g) == adjacency_matrix(ng)
     end
 
-    for (g, ng, rv) in zip(gs, ngs, rvs)
+    for (g, ng) in zip(gs, ngs)
         g, ng = copy(g), copy(ng)
 
         add_vertices!(g, 500)
@@ -88,7 +104,8 @@ end
     
     g0 = erdos_renyi(70, 100; rng=rng)
     rand_g = NautyGraph(g0)
-    @test_throws ErrorException (rand_g_dir = NautyDiGraph(g0))
+    rand_g_dir = NautyDiGraph(g0)
+    @test rand_g_dir.graphset == rand_g.graphset
 
     @test nv(rand_g) == 70
     @test ne(rand_g) == 100
@@ -106,10 +123,15 @@ end
         @test indegree(rand_g, vertex) == length(inneighbors(rand_g, vertex))
     end
 
+    es = [Edge(1, 2), Edge(2, 3), Edge(2, 4)]
     g = NautyDiGraph(4)
-    add_edge!(g, 1, 2)
-    add_edge!(g, 2, 3)
-    add_edge!(g, 2, 4)
+    for e in es
+        add_edge!(g, e)
+    end
+
+    k = NautyDiGraph(es)
+    
+    @test g == k
 
     g2 = copy(g)
     g3 = copy(g)
@@ -137,23 +159,30 @@ end
     h = NautyDiGraph(4)
     copy!(h, g)
     @test h.graphset == g.graphset
-    @test h.n_vertices == g.n_vertices
-    @test h.n_edges == g.n_edges
-    @test h.n_words == g.n_words
+    @test h.graphset.n == g.graphset.n
+    @test h.ne == g.ne
+    @test h.graphset.m == g.graphset.m
     @test h.labels == g.labels
     @test h.hashval == g.hashval
 
 
-    g = NautyGraph(5, 1:5)
-    add_edge!(g, 1, 2)
-    add_edge!(g, 1, 3)
-    add_edge!(g, 1, 4)
-    add_edge!(g, 1, 5)
-    add_edge!(g, 2, 5)
+    glab = NautyGraph(5; vertex_labels=1:5)
+    add_edge!(glab, 1, 2)
+    add_edge!(glab, 1, 3)
+    add_edge!(glab, 1, 4)
+    add_edge!(glab, 1, 5)
+    add_edge!(glab, 2, 5)
 
-    gind1 = g[[1, 5, 2]]
+    gind1 = glab[[1, 5, 2]]
     @test gind1.labels == [1, 5, 2]
 
-    gind2 = g[[Edge(1, 2), Edge(1, 4)]]
+    gind2 = glab[[Edge(1, 2), Edge(1, 4)]]
     @test gind2.labels == [1, 2, 4]
+
+    gb = NautyGraph(g0)
+    vg = DiGraph(g)
+
+    bb_ng = blockdiag(gb, g)
+    bb_g = NautyDiGraph(blockdiag(DiGraph(g0), vg))
+    @test bb_ng == bb_g
 end
\ No newline at end of file
diff --git a/test/graphset.jl b/test/graphset.jl
new file mode 100644
index 0000000..d0c7802
--- /dev/null
+++ b/test/graphset.jl
@@ -0,0 +1,39 @@
+rng = Random.Random.MersenneTwister(0) # Use MersenneTwister for Julia 1.6 compat
+
+function test_graphsets(A; mfacts)
+    n, _ = size(A)
+    for mf in mfacts
+        g16 = Graphset{UInt16}(A, mf * cld(n, NautyGraphs.wordsize(UInt16)))
+        @test g16 == A
+        
+        g32 = Graphset{UInt32}(A, mf * cld(n, NautyGraphs.wordsize(UInt32)))
+        @test g32 == A
+
+        g64 = Graphset{UInt64}(A, mf * cld(n, NautyGraphs.wordsize(UInt64)))
+        @test g64 == A
+    end
+    return
+end
+
+@testset "graphset" begin
+    ns = [1, 2, 5, 15, 16, 17, 31, 32, 33, 63, 64, 65, 500]
+    As = [rand(rng, Bool, n, n) for n in ns]
+    test_graphsets.(As; mfacts=1:3)
+
+    gs1 = Graphset{UInt64}(3, 1)
+    @test_throws BoundsError gs1[1, 4]
+
+    gs2 = Graphset{UInt64}(3, 2)
+    @test_throws BoundsError gs2[1, 4]
+
+    A = [1 0 0; 1 1 0; 0 0 1]
+    gs1 .= A
+    gs2 .= A
+
+    @test gs1 == gs2
+
+    NautyGraphs._rem_vertex!(gs1, 3)
+    NautyGraphs._rem_vertex!(gs2, 1)
+
+    @test gs1 != gs2
+end
\ No newline at end of file
diff --git a/test/interface.jl b/test/interface.jl
new file mode 100644
index 0000000..87515d9
--- /dev/null
+++ b/test/interface.jl
@@ -0,0 +1,12 @@
+using Pkg; Pkg.add(url="https://github.com/JuliaGraphs/GraphsInterfaceChecker.jl")
+using GraphsInterfaceChecker, Interfaces
+
+@testset "interface" begin 
+    test_graphs = [NautyGraph(0), 
+                NautyDiGraph(0), 
+                NautyGraph([1 0 1; 0 0 0; 1 0 1]),
+                NautyDiGraph([0 0 1; 1 0 0; 1 1 1])]
+
+    @implements AbstractGraphInterface{(:mutation)} DenseNautyGraph test_graphs
+    @test Interfaces.test(AbstractGraphInterface, DenseNautyGraph)
+end
diff --git a/test/jet.jl b/test/jet.jl
new file mode 100644
index 0000000..1eebbbc
--- /dev/null
+++ b/test/jet.jl
@@ -0,0 +1,31 @@
+using JET
+
+@testset "JET" begin 
+   test_package(NautyGraphs; target_modules=(NautyGraphs,))
+
+   @test_opt target_modules=(NautyGraphs,) NautyGraph(10)
+   @test_opt target_modules=(NautyGraphs,) NautyDiGraph(10)
+
+   A = [1 0 1; 0 1 0; 1 1 1]
+   @test_opt target_modules=(NautyGraphs,) NautyDiGraph(A)
+
+   g = NautyDiGraph(5; vertex_labels=1:5)
+   h = NautyGraph(g)
+   @test_opt target_modules=(NautyGraphs,) copy(g)
+   @test_opt target_modules=(NautyGraphs,) NautyDiGraph(g)
+
+   add_edge!(g, 2, 5)
+
+   @test_opt target_modules=(NautyGraphs,) add_edge!(g, 1, 2)
+   @test_opt target_modules=(NautyGraphs,) add_vertex!(g)
+   @test_opt target_modules=(NautyGraphs,) add_vertex!(g, 5)
+   @test_opt target_modules=(NautyGraphs,) rem_vertex!(g, 3)
+   @test_opt target_modules=(NautyGraphs,) rem_edge!(g, 2, 5)
+   @test_opt target_modules=(NautyGraphs,) outneighbors(g, 1)
+   @test_opt target_modules=(NautyGraphs,) inneighbors(g, 1)
+   @test_opt target_modules=(NautyGraphs,) collect(edges(g))
+   @test_opt target_modules=(NautyGraphs,) blockdiag(g, h)
+
+   @test_opt target_modules=(NautyGraphs,) nauty(g)
+   @test_opt target_modules=(NautyGraphs,) canonize!(g)
+end
diff --git a/test/nauty.jl b/test/nauty.jl
index 9fc5f79..8598984 100644
--- a/test/nauty.jl
+++ b/test/nauty.jl
@@ -1,5 +1,3 @@
-using Base.Threads
-
 @testset "nauty" begin
     verylarge_g = NautyGraph(50)
 
@@ -19,12 +17,35 @@ using Base.Threads
     @test g1 ≃ h1
     @test ghash(g1) == ghash(h1)
 
-    g2 = NautyGraph(4, [0, 0, 1, 1])
+    g1_16 = NautyGraph{UInt16}(g1)
+    @test g1_16 == g1
+    @test g1_16 ≃ g1
+    @test ghash(g1_16) == ghash(g1)
+
+    g1_32 = NautyGraph{UInt32}(g1)
+    @test g1_32 == g1_16
+    @test g1_32 ≃ g1_16
+    @test ghash(g1_32) == ghash(g1_16)
+    @test g1_32 == g1
+    @test g1_32 ≃ g1
+    @test ghash(g1_32) == ghash(g1)
+
+    k1 = copy(g1)
+    rem_edge!(k1, 2, 3)
+    @test !(k1 ≃ h1)
+    @test ghash(k1) != ghash(h1)
+
+    f1 = copy(g1)
+    rem_vertex!(f1, 2)
+    @test !(f1 ≃ h1)
+    @test ghash(f1) != ghash(h1)
+
+    g2 = NautyGraph(4; vertex_labels=[0, 0, 1, 1])
     add_edge!(g2, 1, 2)
     add_edge!(g2, 2, 3)
     add_edge!(g2, 2, 4)
 
-    h2 = NautyGraph(4, [1, 1, 0, 0])
+    h2 = NautyGraph(4; vertex_labels=[1, 1, 0, 0])
     add_edge!(h2, 3, 4)
     add_edge!(h2, 4, 1)
     add_edge!(h2, 4, 2)
@@ -32,8 +53,20 @@ using Base.Threads
     @test g2 ≃ h2
     @test ghash(g2) == ghash(h2)
 
+    g2_16 = NautyGraph{UInt16}(g2)
+    @test g2_16 == g2
+    @test g2_16 ≃ g2
+    @test ghash(g2_16) == ghash(g2)
+
+    g2_32 = NautyGraph{UInt32}(g2)
+    @test g2_32 == g2_16
+    @test g2_32 ≃ g2_16
+    @test ghash(g2_32) == ghash(g2_16)
+    @test g2_32 == g2
+    @test g2_32 ≃ g2
+    @test ghash(g2_32) == ghash(g2)
 
-    k2 = NautyGraph(4, [1, 0, 0, 1])
+    k2 = NautyGraph(4; vertex_labels=[1, 0, 0, 1])
     add_edge!(k2, 3, 4)
     add_edge!(k2, 4, 1)
     add_edge!(k2, 4, 2)
@@ -41,6 +74,10 @@ using Base.Threads
     @test !(g2 ≃ k2)
     @test ghash(g2) != ghash(k2)
 
+    @test !(g2_16 ≃ k2)
+    @test !(g2_32 ≃ k2)
+    @test ghash(g2_16) != ghash(k2)
+    @test ghash(g2_32) != ghash(k2)
 
     g3 = NautyGraph(6)
     add_edge!(g3, 1, 2)
@@ -70,8 +107,8 @@ using Base.Threads
     @test adjacency_matrix(m3) == adjacency_matrix(h3)
 
 
-    g4 = NautyGraph(3, [1, 2, 3])
-    h4 = NautyGraph(3, [1, 2, 3])
+    g4 = NautyGraph(3; vertex_labels=[1, 2, 3])
+    h4 = NautyGraph(3; vertex_labels=[1, 2, 3])
     add_edge!(g4, 1, 2)
     add_edge!(h4, 1, 2)
 
@@ -83,7 +120,7 @@ using Base.Threads
     @test Base.hash(g4) == Base.hash(h4)
 
 
-    g5 = NautyGraph(10, 10:-1:1)
+    g5 = NautyGraph(10; vertex_labels=10:-1:1)
     add_edge!(g5, 1, 2)
     add_edge!(g5, 5, 2)
     add_edge!(g5, 6, 7)
@@ -96,13 +133,16 @@ using Base.Threads
     canonperm5 = canonical_permutation(g5)
     @test canon5.labels == g5.labels[canonperm5]
 
-    thread_gs = fill(copy(g4), 10)
-    vals = []
-    @threads for i in eachindex(thread_gs)
-        push!(vals, nauty(thread_gs[i]))
+    # Just test that multithreading doesnt lead to errors
+    thread_gs = [copy(g4) for i in 1:10]
+    vals = Any[nothing for i in 1:10]
+    @threads for i in eachindex(vals, thread_gs)
+        for j in 1:20
+            vals[i] = nauty(thread_gs[i])
+            sleep(0.01)
+        end
     end
-    @test length(vals) == length(thread_gs)
-
+    @test true
 
     gnoloop = NautyGraph(5)
     add_edge!(gnoloop, 1, 2)
@@ -125,4 +165,9 @@ using Base.Threads
 
     @test_nowarn nauty(gdiloop)
     @test !is_isomorphic(gdinoloop, gdiloop)
+
+    # Test that ghash doesnt error for large graphs
+    glarge = NautyGraph(200)
+    ghash(glarge)
+    @test true
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 0646184..1a29d37 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,8 +1,15 @@
 using NautyGraphs, Graphs
+using NautyGraphs: Graphset
 using Test
 using Random, LinearAlgebra
+using Base.Threads
 
 @testset verbose=true "NautyGraphs" begin
-    include("densenautygraphs.jl")
+    include("densenautygraph.jl")
     include("nauty.jl")
+    include("graphset.jl")
+    include("utils.jl")
+    VERSION >= v"1.9-" && include("interface.jl")
+    include("aqua.jl")
+    VERSION >= v"1.10" && include("jet.jl")
 end
diff --git a/test/utils.jl b/test/utils.jl
new file mode 100644
index 0000000..156e38f
--- /dev/null
+++ b/test/utils.jl
@@ -0,0 +1,7 @@
+
+@testset "utils" begin
+   # Test that sha works without error
+   A = ones(10_000)
+   NautyGraphs.hash_sha(ones(10_000))
+   @test true
+end