diff --git a/apps/web/src/routes/api/misskey/invite-code/+server.ts b/apps/web/src/routes/api/misskey/invite-code/+server.ts
index 79c6509..1e5d47a 100644
--- a/apps/web/src/routes/api/misskey/invite-code/+server.ts
+++ b/apps/web/src/routes/api/misskey/invite-code/+server.ts
@@ -1,4 +1,4 @@
-import { sampleMany } from "@cannorin/utils/array";
+import { sampleMany } from "@cannorin/utils";
 import { type RequestEvent, text } from "@sveltejs/kit";
 import { RateLimiter } from "sveltekit-rate-limiter/server";
 
diff --git a/apps/web/src/routes/kripke/lib/store.ts b/apps/web/src/routes/kripke/lib/store.ts
index f29d6ce..5b91b3c 100644
--- a/apps/web/src/routes/kripke/lib/store.ts
+++ b/apps/web/src/routes/kripke/lib/store.ts
@@ -1,4 +1,4 @@
-import { MultiSet } from "@cannorin/utils/multiset";
+import { MultiSet } from "@cannorin/utils";
 import { persisted } from "svelte-persisted-store";
 import type { Move } from "../components/game.svelte";
 import { date } from "./system";
diff --git a/apps/web/src/routes/kripke/random/[seed]/+page.svelte b/apps/web/src/routes/kripke/random/[seed]/+page.svelte
index f3774e0..f89861e 100644
--- a/apps/web/src/routes/kripke/random/[seed]/+page.svelte
+++ b/apps/web/src/routes/kripke/random/[seed]/+page.svelte
@@ -1,8 +1,7 @@
 <script lang="ts">
 import { Button } from "$lib/components/ui/button";
 import * as Dialog from "$lib/components/ui/dialog";
-import { type Formula, isomorphic } from "@cannorin/kripke";
-import { validWorlds } from "@cannorin/kripke/sat";
+import { type Formula, isomorphic, validWorlds } from "@cannorin/kripke";
 import LuRotateCw from "lucide-svelte/icons/rotate-cw";
 import LuX from "lucide-svelte/icons/x";
 
diff --git a/packages/kripke/package.json b/packages/kripke/package.json
index 4061447..9e096d0 100644
--- a/packages/kripke/package.json
+++ b/packages/kripke/package.json
@@ -3,11 +3,11 @@
   "version": "0.0.1",
   "private": true,
   "exports": {
-    ".": "./index.ts",
-    "./syntax": "./syntax.ts",
-    "./semantics": "./semantics.ts",
-    "./parser": "./parser.ts",
-    "./sat": "./sat.ts"
+    ".": "./src/index.ts",
+    "./syntax": "./src/syntax.ts",
+    "./semantics": "./src/semantics.ts",
+    "./parser": "./src/parser.ts",
+    "./sat": "./src/sat.ts"
   },
   "scripts": {
     "build": "tsc",
diff --git a/packages/kripke/index.ts b/packages/kripke/src/index.ts
similarity index 77%
rename from packages/kripke/index.ts
rename to packages/kripke/src/index.ts
index b2ccf9a..c81309a 100644
--- a/packages/kripke/index.ts
+++ b/packages/kripke/src/index.ts
@@ -1,3 +1,4 @@
 export * from "./syntax";
 export * from "./semantics";
 export * from "./parser";
+export * from "./sat";
diff --git a/packages/kripke/parser.ts b/packages/kripke/src/parser.ts
similarity index 100%
rename from packages/kripke/parser.ts
rename to packages/kripke/src/parser.ts
diff --git a/packages/kripke/sat.ts b/packages/kripke/src/sat.ts
similarity index 64%
rename from packages/kripke/sat.ts
rename to packages/kripke/src/sat.ts
index b176e34..919b3ab 100644
--- a/packages/kripke/sat.ts
+++ b/packages/kripke/src/sat.ts
@@ -1,19 +1,37 @@
-import type {} from "@cannorin/utils/headless";
-import { type Frame, type World, left, right, worlds } from "./semantics";
+import { BitSet, maximal } from "@cannorin/utils";
+import {
+  type Frame,
+  type Model,
+  type World,
+  left,
+  right,
+  worlds,
+} from "./semantics";
 import {
   type Formula,
   type NNFFormula,
-  type PropVar,
   not,
+  propVars,
   simplify,
   tagNNF,
   toNNF,
 } from "./syntax";
-import { maximal } from "./utils";
 
-export type Constraints =
-  | true
-  | ReadonlySet<`${World}${PropVar}` | `!${World}${PropVar}`>;
+const allValuations = worlds.flatMap((w) =>
+  propVars.map((p) => `${w}${p}` as const),
+);
+
+const valuationSet = new BitSet(allValuations);
+
+export type Constraints = {
+  positive: number;
+  negative: number;
+};
+
+const empty: Constraints = {
+  positive: valuationSet.empty,
+  negative: valuationSet.empty,
+};
 
 export type Context = {
   notFml: NNFFormula<number>;
@@ -103,30 +121,25 @@ export function buildContext(frame: Frame, fml: Formula): Context {
   return { notFml, next, memo: new Map(), constant };
 }
 
-const subsumes = (c1: Constraints, c2: Constraints): boolean => {
-  if (c1 === true) return c2 === true;
-  if (c2 === true) return true;
-  for (const v of c2.values()) {
-    if (!c1.has(v)) return false;
-  }
-  return true;
+const isWeaker = (c1: Constraints, c2: Constraints): boolean => {
+  return (
+    valuationSet.isSuperset(c1.positive, c2.positive) &&
+    valuationSet.isSuperset(c1.negative, c2.negative)
+  );
 };
 
 const union = (c1: Constraints, c2: Constraints): Constraints => {
-  if (c1 === true) return c2;
-  if (c2 === true) return c1;
-  return new Set([...c1.values(), ...c2.values()]);
+  return {
+    positive: valuationSet.union(c1.positive, c2.positive),
+    negative: valuationSet.union(c1.negative, c2.negative),
+  };
 };
 
 const consistent = (c: Constraints) => {
-  if (c === true) return true;
-  for (const v of c.values()) {
-    if (!v.startsWith("!") && c.has(`!${v}`)) return false;
-  }
-  return true;
+  return valuationSet.isDisjoint(c.positive, c.negative);
 };
 
-export function sat(
+function sat(
   frame: Frame,
   fml: NNFFormula<number>,
   world: World,
@@ -134,7 +147,7 @@ export function sat(
 ): Constraints[] {
   const key = `${world}:${fml.tag}` as const;
   const c = ctx.constant.get(key);
-  if (c === true) return [true];
+  if (c === true) return [empty];
   if (c === false) return [];
 
   let result = ctx.memo.get(key);
@@ -142,16 +155,26 @@ export function sat(
 
   switch (fml.type) {
     case "top":
-      result = [true];
+      result = [empty];
       break;
     case "bot":
       result = [];
       break;
     case "propvar":
-      result = [new Set([`${world}${fml.name}` as const])];
+      result = [
+        {
+          positive: valuationSet.create(`${world}${fml.name}`),
+          negative: valuationSet.empty,
+        },
+      ];
       break;
     case "not":
-      result = [new Set([`!${world}${fml.fml.name}` as const])];
+      result = [
+        {
+          positive: valuationSet.empty,
+          negative: valuationSet.create(`${world}${fml.fml.name}`),
+        },
+      ];
       break;
     case "box": {
       result = ctx.next[world]
@@ -162,16 +185,16 @@ export function sat(
               current.flatMap((c1) =>
                 prev.map((c2) => union(c1, c2)).filter(consistent),
               ),
-              subsumes,
+              isWeaker,
             ),
-          [true],
+          [empty],
         );
       break;
     }
     case "diamond": {
       result = maximal(
         ctx.next[world].flatMap((w) => sat(frame, fml.fml, w, ctx)),
-        subsumes,
+        isWeaker,
       );
       break;
     }
@@ -185,9 +208,9 @@ export function sat(
             current.flatMap((c1) =>
               prev.map((c2) => union(c1, c2)).filter(consistent),
             ),
-            subsumes,
+            isWeaker,
           ),
-        [true],
+        [empty],
       );
       break;
     }
@@ -197,7 +220,7 @@ export function sat(
           sat(frame, fml.left, world, ctx),
           sat(frame, fml.right, world, ctx),
         ].flat(),
-        subsumes,
+        isWeaker,
       );
       break;
     }
@@ -207,14 +230,53 @@ export function sat(
   return result;
 }
 
-export function validInWorld(frame: Frame, fml: Formula, world: World) {
-  const ctx = buildContext(frame, fml);
+export function findCountermodelsAt(
+  frame: Frame,
+  fml: Formula,
+  world: World,
+  ctx = buildContext(frame, fml),
+): Model[] {
+  return sat(frame, ctx.notFml, world, ctx).map(({ positive }) => ({
+    ...frame,
+    valuations: valuationSet.decode(positive),
+  }));
+}
+
+export function findCountermodels(
+  frame: Frame,
+  fml: Formula,
+  ctx = buildContext(frame, fml),
+): Model[] {
+  return maximal(
+    worlds.flatMap((world) => sat(frame, ctx.notFml, world, ctx)),
+    isWeaker,
+  ).map(({ positive }) => ({
+    ...frame,
+    valuations: valuationSet.decode(positive),
+  }));
+}
+
+export function validInWorld(
+  frame: Frame,
+  fml: Formula,
+  world: World,
+  ctx = buildContext(frame, fml),
+) {
   return sat(frame, ctx.notFml, world, ctx).length === 0;
 }
 
-export function validWorlds(frame: Frame, fml: Formula) {
-  const ctx = buildContext(frame, fml);
-  return worlds.filter(
-    (world) => sat(frame, ctx.notFml, world, ctx).length === 0,
-  );
+export function validWorlds(
+  frame: Frame,
+  fml: Formula,
+  ctx = buildContext(frame, fml),
+) {
+  return worlds.filter((world) => validInWorld(frame, ctx.notFml, world, ctx));
+}
+
+export function validInFrame(
+  frame: Frame,
+  fml: Formula,
+  ctx = buildContext(frame, fml),
+) {
+  return worlds.every((world) => validInWorld(frame, ctx.notFml, world, ctx));
 }
diff --git a/packages/kripke/semantics.ts b/packages/kripke/src/semantics.ts
similarity index 77%
rename from packages/kripke/semantics.ts
rename to packages/kripke/src/semantics.ts
index f49699d..b5f3029 100644
--- a/packages/kripke/semantics.ts
+++ b/packages/kripke/src/semantics.ts
@@ -1,5 +1,5 @@
-import { type Formula, type PropVar, propVars, vars } from "./syntax";
-import { decode, encode, permutations, power } from "./utils";
+import { BitSet, permutate } from "@cannorin/utils";
+import type { Formula, PropVar } from "./syntax";
 
 export const worlds = ["a", "b", "c", "d"] as const;
 
@@ -15,6 +15,8 @@ export const relation: Relation[] = worlds.flatMap((w) =>
   worlds.map((x) => `${w}${x}` as const),
 );
 
+const relationSet = new BitSet(relation);
+
 export interface Frame {
   relations: Set<Relation>;
 }
@@ -23,11 +25,6 @@ export interface Model extends Frame {
   valuations: Set<`${World}${PropVar}`>;
 }
 
-export function valuation(fml?: Formula): `${World}${PropVar}`[] {
-  const vs = fml ? Array.from(vars(fml)) : propVars;
-  return worlds.flatMap((w) => vs.map((p) => `${w}${p}` as const));
-}
-
 export function satisfy(m: Model, w: World, fml: Formula): boolean {
   switch (fml.type) {
     case "top":
@@ -73,37 +70,15 @@ export function satisfy(m: Model, w: World, fml: Formula): boolean {
 export const validInModel = (m: Model, fml: Formula) =>
   worlds.every((w) => satisfy(m, w, fml));
 
-export function validInFrame(f: Frame, fml: Formula) {
-  for (const valuations of power(valuation(fml))) {
-    if (!validInModel({ ...f, valuations }, fml)) return false;
-  }
-  return true;
-}
-
-export function validWorlds(f: Frame, fml: Formula) {
-  const result: World[] = [];
-  for (const w of worlds) {
-    let valid = true;
-    for (const valuations of power(valuation(fml))) {
-      if (!satisfy({ ...f, valuations }, w, fml)) {
-        valid = false;
-        break;
-      }
-    }
-    if (valid) result.push(w);
-  }
-  return result;
-}
-
 export function getFrame(id: number): Frame {
-  return { relations: decode(relation, id) };
+  return { relations: relationSet.decode(id) };
 }
 
 export function getId(frame: Frame) {
-  return encode(relation, frame.relations);
+  return relationSet.encode(frame.relations);
 }
 
-const worldPermutations = permutations(worlds).map(
+const worldPermutations = permutate(worlds).map(
   (perm) => new Map(worlds.map((k, i) => [k, perm[i] as World])),
 );
 
@@ -131,7 +106,7 @@ export function generateAllFrames() {
   for (let id = 0; id < total; id++) {
     if (map.has(id)) continue;
 
-    const relations = decode(relation, id);
+    const relations = relationSet.decode(id);
 
     const frame = { relations };
     const equivalentIds: number[] = [];
@@ -139,7 +114,7 @@ export function generateAllFrames() {
     let canonicalId = id;
     for (const perm of worldPermutations) {
       const permuted = applyPermutation(frame, perm);
-      const permutedId = encode(relation, permuted.relations);
+      const permutedId = relationSet.encode(permuted.relations);
       equivalentIds.push(permutedId);
       if (canonicalId === null || permutedId < canonicalId) {
         canonicalId = permutedId;
diff --git a/packages/kripke/syntax.ts b/packages/kripke/src/syntax.ts
similarity index 100%
rename from packages/kripke/syntax.ts
rename to packages/kripke/src/syntax.ts
diff --git a/packages/kripke/tests/sat.test.ts b/packages/kripke/tests/sat.test.ts
index 4871996..5b09397 100644
--- a/packages/kripke/tests/sat.test.ts
+++ b/packages/kripke/tests/sat.test.ts
@@ -1,7 +1,14 @@
+import { power } from "@cannorin/utils";
 import { expect, test } from "vitest";
-import { validWorlds } from "../sat";
-import { validWorlds as validWorldsNaive } from "../semantics";
-import { prettyPrint } from "../syntax";
+import { validWorlds } from "../src/sat";
+import { type Frame, type World, satisfy, worlds } from "../src/semantics";
+import {
+  type Formula,
+  type PropVar,
+  prettyPrint,
+  propVars,
+  vars,
+} from "../src/syntax";
 import { randomFormula, randomFrame, testFormulas } from "./utils";
 
 const elapsed = <T>(f: () => T) => {
@@ -11,15 +18,43 @@ const elapsed = <T>(f: () => T) => {
   return { value, elapsed: end - start };
 };
 
+export function validWorldsNaive(
+  f: Frame,
+  fml: Formula,
+  allValuations: `${World}${PropVar}`[][],
+) {
+  const result: World[] = [];
+  for (const w of worlds) {
+    let valid = true;
+    for (const v of allValuations) {
+      if (!satisfy({ ...f, valuations: new Set(v) }, w, fml)) {
+        valid = false;
+        break;
+      }
+    }
+    if (valid) result.push(w);
+  }
+  return result;
+}
+
 for (const fml of testFormulas.concat(
-  [...Array(50)].map(() => randomFormula()),
+  [...Array(100)].map(() => randomFormula()),
 )) {
   test(`SAT works for ${prettyPrint(fml)}`, () => {
     const count = 100;
     let diff = 0;
+    const allValuations = power(
+      worlds.flatMap((w) =>
+        (fml ? Array.from(vars(fml)) : propVars).map(
+          (p) => `${w}${p}` as const,
+        ),
+      ),
+    );
     for (let i = 0; i < count; i++) {
       const frame = randomFrame();
-      const expected = elapsed(() => validWorldsNaive(frame, fml));
+      const expected = elapsed(() =>
+        validWorldsNaive(frame, fml, allValuations),
+      );
       const actual = elapsed(() => validWorlds(frame, fml));
 
       expect(actual.value.sort()).toStrictEqual(expected.value.sort());
diff --git a/packages/kripke/tests/utils.ts b/packages/kripke/tests/utils.ts
index 5566871..7b8ed09 100644
--- a/packages/kripke/tests/utils.ts
+++ b/packages/kripke/tests/utils.ts
@@ -1,6 +1,6 @@
-import { sample } from "@cannorin/utils/array";
-import { parse } from "../parser";
-import { getFrame, nontrivials } from "../semantics";
+import { sample } from "@cannorin/utils";
+import { parse } from "../src/parser";
+import { getFrame, nontrivials } from "../src/semantics";
 import {
   type Formula,
   and,
@@ -13,7 +13,7 @@ import {
   propvar,
   to,
   top,
-} from "../syntax";
+} from "../src/syntax";
 
 const formulaTypes = [
   "top",
@@ -52,9 +52,24 @@ export const testFormulas: Formula[] = [
   parse("LLp -> p"),
   parse("p -> LLp"),
 
+  parse("M(p -> Mp)"),
+  parse("r -> M((p -> Mp) & Mr)"),
+  parse("p -> MMp"),
+  parse("M(p -> Lp)"),
+  parse("L0"),
+  parse("ML0"),
+  parse("L(p -> Mp)"),
+  parse("p -> MMMp"),
+  parse("p -> MMMMp"),
+  parse("M1"),
+  parse("ML(p -> Mp)"),
+  parse("Lp -> LMp"),
+  parse("M(p -> LLp)"),
+
   // pathological
   parse("L(M(1 & p) v LM1) -> (s <-> 1)"),
   parse("L(M((q v ~q) & p) v LM(s v ~s)) -> (s <-> (p v ~p))"),
+  parse("~(Lr & ¬Mp) & M~L(q v q)"),
 ];
 
 export function randomFormula(depth = 5): Formula {
diff --git a/packages/kripke/utils.ts b/packages/kripke/utils.ts
deleted file mode 100644
index 7a2d82c..0000000
--- a/packages/kripke/utils.ts
+++ /dev/null
@@ -1,66 +0,0 @@
-const bitIsSet = (num: number, pos: number) => (num & (1 << pos)) !== 0;
-
-export function encode<T>(all: readonly T[], set: Set<T>): number {
-  let flags = 0;
-  for (let i = 0; i < all.length; i++) {
-    if (set.has(all[i] as T)) {
-      flags |= 1 << i;
-    }
-  }
-  return flags;
-}
-
-export function decode<T>(all: readonly T[], flags: number): Set<T> {
-  const total = 2 ** all.length;
-  if (flags < 0 || flags >= total) throw Error("invalid flags");
-  const decoded = new Set<T>();
-  for (let j = 0; j < all.length; j++) {
-    if (bitIsSet(flags, j)) decoded.add(all[j] as T);
-  }
-  return decoded;
-}
-
-export function* power<T>(xs: readonly T[]) {
-  const total = 2 ** xs.length;
-  for (let i = 0; i < total; i++) {
-    yield decode(xs, i);
-  }
-}
-
-export function permutations<T>(arr: readonly T[]): T[][] {
-  if (arr.length === 0) return [[]];
-  const result: T[][] = [];
-  for (let i = 0; i < arr.length; i++) {
-    const current = arr[i];
-    const remaining = [...arr.slice(0, i), ...arr.slice(i + 1)];
-    for (const perm of permutations(remaining)) {
-      result.push([current as T, ...perm]);
-    }
-  }
-  return result;
-}
-
-export function maximal<T>(
-  xs: readonly T[],
-  preorder: (x: T, y: T) => boolean,
-): T[] {
-  const res: T[] = [];
-  outer: for (const x of xs) {
-    for (let i = 0; i < res.length; ) {
-      const y = res[i] as T;
-      const xLeY = preorder(x, y);
-      const yLeX = preorder(y, x);
-
-      if (xLeY && !yLeX) {
-        continue outer;
-      }
-      if (yLeX && !xLeY) {
-        res.splice(i, 1);
-        continue;
-      }
-      i++;
-    }
-    res.push(x);
-  }
-  return res;
-}
diff --git a/packages/utils/index.ts b/packages/utils/index.ts
deleted file mode 100644
index cb0ff5c..0000000
--- a/packages/utils/index.ts
+++ /dev/null
@@ -1 +0,0 @@
-export {};
diff --git a/packages/utils/package.json b/packages/utils/package.json
index c1caa29..0db5db9 100644
--- a/packages/utils/package.json
+++ b/packages/utils/package.json
@@ -3,11 +3,12 @@
   "version": "0.0.1",
   "private": true,
   "exports": {
-    ".": "./index.ts",
+    ".": "./src/index.ts",
     "./headless": "./src/headless.ts",
     "./array": "./src/array.ts",
     "./vector": "./src/vector.ts",
-    "./multiset": "./src/multiset.ts"
+    "./multiset": "./src/multiset.ts",
+    "./bitset": "./src/bitset.ts"
   },
   "scripts": {
     "build": "tsc",
diff --git a/packages/utils/src/array.ts b/packages/utils/src/array.ts
index 077f532..bdf41e4 100644
--- a/packages/utils/src/array.ts
+++ b/packages/utils/src/array.ts
@@ -1,7 +1,9 @@
+/** Returns a random element of the array  */
 export const sample = <T>(arr: readonly T[]): T =>
   arr[Math.floor(Math.random() * arr.length)] as T;
 
-export function sampleMany<T>(arr: T[], n: number = arr.length): T[] {
+/** Returns random N elements of the array  */
+export function sampleMany<T>(arr: readonly T[], n: number = arr.length): T[] {
   if (n > arr.length) return sampleMany(arr, arr.length);
   const copy = [...arr];
   for (let i = 0; i < n; i++) {
@@ -10,3 +12,47 @@ export function sampleMany<T>(arr: T[], n: number = arr.length): T[] {
   }
   return copy.slice(0, n);
 }
+
+/** Returns a powerset of the array */
+export const power = <T>(arr: readonly T[]) =>
+  arr.reduce((a, v) => a.concat(a.map((r) => r.concat(v))), [[]] as T[][]);
+
+/** Returns all the permutations of the array */
+export function permutate<T>(arr: readonly T[]): T[][] {
+  if (arr.length === 0) return [[]];
+  const result: T[][] = [];
+  for (let i = 0; i < arr.length; i++) {
+    const current = arr[i];
+    const remaining = [...arr.slice(0, i), ...arr.slice(i + 1)];
+    for (const perm of permutate(remaining)) {
+      result.push([current as T, ...perm]);
+    }
+  }
+  return result;
+}
+
+/** Returns the maximal elements of the array w.r.t. the preorder */
+export function maximal<T>(
+  arr: readonly T[],
+  preorder: (x: T, y: T) => boolean,
+): T[] {
+  const res: T[] = [];
+  outer: for (const x of arr) {
+    for (let i = 0; i < res.length; ) {
+      const y = res[i] as T;
+      const xLeY = preorder(x, y);
+      const yLeX = preorder(y, x);
+
+      if (xLeY) {
+        continue outer;
+      }
+      if (yLeX) {
+        res.splice(i, 1);
+        continue;
+      }
+      i++;
+    }
+    res.push(x);
+  }
+  return res;
+}
diff --git a/packages/utils/src/bitset.ts b/packages/utils/src/bitset.ts
new file mode 100644
index 0000000..082f593
--- /dev/null
+++ b/packages/utils/src/bitset.ts
@@ -0,0 +1,78 @@
+/** Generates a wrapper for bitset operations.
+ *  Works as an efficient replacement of `Set<T>` if the universal set is known.
+ */
+export class BitSet<T> {
+  readonly universe: readonly T[];
+  readonly empty: number;
+  readonly all: number;
+
+  constructor(universe: readonly T[]) {
+    this.universe = universe;
+    this.empty = 0;
+    this.all = 2 ** universe.length - 1;
+  }
+
+  /** Create a bitset that represents the given values.
+   * The value that is not in the universe is ignored.
+   */
+  create(...values: readonly T[]) {
+    return this.encode(new Set(values));
+  }
+
+  /** Create a bitset from the given set.
+   * The value that is not in the universe is ignored.
+   */
+  encode(set: Set<T>): number {
+    let x = 0;
+    for (let i = 0; i < this.universe.length; i++) {
+      if (set.has(this.universe[i] as T)) {
+        x |= 1 << i;
+      }
+    }
+    return x;
+  }
+
+  /** Create a set from the given bitset. */
+  decode(x: number): Set<T> {
+    if (x < this.empty || x > this.all) throw Error("invalid bits");
+    const decoded = new Set<T>();
+    for (let j = 0; j < this.universe.length; j++) {
+      if ((x & (1 << j)) !== 0) decoded.add(this.universe[j] as T);
+    }
+    return decoded;
+  }
+
+  /** Returns the union of all the given bitsets. */
+  union(...xs: readonly number[]): number {
+    return xs.reduce((x, y) => x | y, this.empty);
+  }
+
+  /** Returns the intersection of all the given bitsets. */
+  intersection(...xs: readonly number[]): number {
+    return xs.reduce((x, y) => x & y, this.all);
+  }
+
+  /** Returns the bitset that is obtained from `x` by removing all the elements in `y`. */
+  difference(x: number, y: number) {
+    return x & (this.all - y);
+  }
+
+  /** Returns if `x` and `y` are disjoint. */
+  isDisjoint(x: number, y: number): boolean {
+    return (x & y) === 0;
+  }
+
+  /**
+   * Returns if `x` is a subset of `y`.
+   */
+  isSubset(x: number, y: number): boolean {
+    return (x & y) === x;
+  }
+
+  /**
+   * Returns if `x` is a superset of `y`.
+   */
+  isSuperset(x: number, y: number): boolean {
+    return (x & y) === y;
+  }
+}
diff --git a/packages/utils/src/index.ts b/packages/utils/src/index.ts
new file mode 100644
index 0000000..dcc5ab3
--- /dev/null
+++ b/packages/utils/src/index.ts
@@ -0,0 +1,9 @@
+export * from "./array";
+
+export * from "./bitset";
+
+export * from "./headless";
+
+export * from "./multiset";
+
+export * from "./vector";
diff --git a/packages/utils/tests/array.test.ts b/packages/utils/tests/array.test.ts
new file mode 100644
index 0000000..1ef30a2
--- /dev/null
+++ b/packages/utils/tests/array.test.ts
@@ -0,0 +1,215 @@
+import { describe, expect, it } from "vitest";
+import { maximal, permutate, power, sample, sampleMany } from "../src/array";
+
+describe("array utilities", () => {
+  describe("sample", () => {
+    it("should return an element from the array", () => {
+      const arr = [1, 2, 3, 4, 5];
+      const result = sample(arr);
+      expect(arr).toContain(result);
+    });
+
+    it("should work with single element array", () => {
+      const arr = [42];
+      expect(sample(arr)).toBe(42);
+    });
+
+    it("should work with different types", () => {
+      const strArr = ["a", "b", "c"];
+      const result = sample(strArr);
+      expect(strArr).toContain(result);
+    });
+  });
+
+  describe("sampleMany", () => {
+    it("should return n random elements", () => {
+      const arr = [1, 2, 3, 4, 5];
+      const result = sampleMany(arr, 3);
+      expect(result).toHaveLength(3);
+      for (const item of result) {
+        expect(arr).toContain(item);
+      }
+    });
+
+    it("should return all elements when n equals array length", () => {
+      const arr = [1, 2, 3];
+      const result = sampleMany(arr, 3);
+      expect(result).toHaveLength(3);
+      expect(result.sort()).toEqual(arr.sort());
+    });
+
+    it("should return all elements when n is greater than array length", () => {
+      const arr = [1, 2, 3];
+      const result = sampleMany(arr, 5);
+      expect(result).toHaveLength(3);
+      expect(result.sort()).toEqual(arr.sort());
+    });
+
+    it("should return all elements when n is not provided", () => {
+      const arr = [1, 2, 3, 4];
+      const result = sampleMany(arr);
+      expect(result).toHaveLength(4);
+      expect(result.sort()).toEqual(arr.sort());
+    });
+
+    it("should return empty array for empty input", () => {
+      const arr: number[] = [];
+      const result = sampleMany(arr, 3);
+      expect(result).toEqual([]);
+    });
+
+    it("should not return duplicate elements", () => {
+      const arr = [1, 2, 3, 4, 5];
+      const result = sampleMany(arr, 5);
+      const uniqueSet = new Set(result);
+      expect(uniqueSet.size).toBe(result.length);
+    });
+  });
+
+  describe("power", () => {
+    it("should return powerset of array", () => {
+      const arr = [1, 2];
+      const result = power(arr);
+      expect(result).toHaveLength(4);
+      expect(result).toContainEqual([]);
+      expect(result).toContainEqual([1]);
+      expect(result).toContainEqual([2]);
+      expect(result).toContainEqual([1, 2]);
+    });
+
+    it("should return [[]] for empty array", () => {
+      const arr: number[] = [];
+      const result = power(arr);
+      expect(result).toEqual([[]]);
+    });
+
+    it("should handle three elements", () => {
+      const arr = [1, 2, 3];
+      const result = power(arr);
+      expect(result).toHaveLength(8);
+      expect(result).toContainEqual([]);
+      expect(result).toContainEqual([1]);
+      expect(result).toContainEqual([2]);
+      expect(result).toContainEqual([3]);
+      expect(result).toContainEqual([1, 2]);
+      expect(result).toContainEqual([1, 3]);
+      expect(result).toContainEqual([2, 3]);
+      expect(result).toContainEqual([1, 2, 3]);
+    });
+
+    it("should work with strings", () => {
+      const arr = ["a", "b"];
+      const result = power(arr);
+      expect(result).toHaveLength(4);
+      expect(result).toContainEqual([]);
+      expect(result).toContainEqual(["a"]);
+      expect(result).toContainEqual(["b"]);
+      expect(result).toContainEqual(["a", "b"]);
+    });
+  });
+
+  describe("permutate", () => {
+    it("should return all permutations of array", () => {
+      const arr = [1, 2];
+      const result = permutate(arr);
+      expect(result).toHaveLength(2);
+      expect(result).toContainEqual([1, 2]);
+      expect(result).toContainEqual([2, 1]);
+    });
+
+    it("should return [[]] for empty array", () => {
+      const arr: number[] = [];
+      const result = permutate(arr);
+      expect(result).toEqual([[]]);
+    });
+
+    it("should handle single element", () => {
+      const arr = [1];
+      const result = permutate(arr);
+      expect(result).toEqual([[1]]);
+    });
+
+    it("should handle three elements", () => {
+      const arr = [1, 2, 3];
+      const result = permutate(arr);
+      expect(result).toHaveLength(6);
+      expect(result).toContainEqual([1, 2, 3]);
+      expect(result).toContainEqual([1, 3, 2]);
+      expect(result).toContainEqual([2, 1, 3]);
+      expect(result).toContainEqual([2, 3, 1]);
+      expect(result).toContainEqual([3, 1, 2]);
+      expect(result).toContainEqual([3, 2, 1]);
+    });
+
+    it("should work with strings", () => {
+      const arr = ["a", "b", "c"];
+      const result = permutate(arr);
+      expect(result).toHaveLength(6);
+      // Check that all permutations are unique
+      const stringified = result.map((p) => p.join(""));
+      const unique = new Set(stringified);
+      expect(unique.size).toBe(6);
+    });
+  });
+
+  describe("maximal", () => {
+    it("should return maximal elements with simple less-than preorder", () => {
+      const arr = [1, 2, 3, 2, 1];
+      const result = maximal(arr, (x, y) => x <= y);
+      expect(result).toEqual([3]);
+    });
+
+    it("should handle incomparable elements", () => {
+      // Divisibility preorder: x <= y if x divides y
+      const arr = [2, 3, 4, 6, 8, 9];
+      const divides = (x: number, y: number) => y % x === 0;
+      const result = maximal(arr, divides);
+      // 8 and 9 are maximal (nothing in the array is divisible by them except themselves)
+      expect(result.sort()).toEqual([6, 8, 9]);
+    });
+
+    it("should return all elements when all are incomparable", () => {
+      const arr = [1, 2, 3];
+      // Always false preorder means all elements are incomparable
+      const result = maximal(arr, () => false);
+      expect(result.sort()).toEqual([1, 2, 3]);
+    });
+
+    it("should handle empty array", () => {
+      const arr: number[] = [];
+      const result = maximal(arr, (x, y) => x <= y);
+      expect(result).toEqual([]);
+    });
+
+    it("should handle single element", () => {
+      const arr = [42];
+      const result = maximal(arr, (x, y) => x <= y);
+      expect(result).toEqual([42]);
+    });
+
+    it("should handle subset preorder", () => {
+      // Subset preorder on sets represented as arrays
+      const arr = [[1], [2], [1, 2], [2, 3], [1, 2, 3]];
+      const isSubset = (x: number[], y: number[]) =>
+        x.every((elem) => y.includes(elem));
+      const result = maximal(arr, isSubset);
+      expect(result).toHaveLength(1);
+      expect(result[0]).toEqual([1, 2, 3]);
+    });
+
+    it("should remove duplicate maximal elements", () => {
+      const arr = [3, 1, 2, 3, 2];
+      const result = maximal(arr, (x, y) => x <= y);
+      expect(result).toEqual([3]);
+    });
+
+    it("should preserve order of maximal elements", () => {
+      // Divisibility preorder: x <= y if x divides y
+      const arr = [2, 6, 4, 3, 9, 8];
+      const divides = (x: number, y: number) => y % x === 0;
+      const result = maximal(arr, divides);
+      // Check that the maximal elements appear in the same order as in the original array
+      expect(result).toEqual([6, 9, 8]);
+    });
+  });
+});
diff --git a/packages/utils/tests/bitset.test.ts b/packages/utils/tests/bitset.test.ts
new file mode 100644
index 0000000..7ca51bd
--- /dev/null
+++ b/packages/utils/tests/bitset.test.ts
@@ -0,0 +1,371 @@
+import { describe, expect, it } from "vitest";
+import { BitSet } from "../src/bitset";
+
+describe("BitSet", () => {
+  // Test with string universe
+  const stringUniverse = ["a", "b", "c", "d", "e"];
+  const stringBitSet = new BitSet<string>(stringUniverse);
+
+  // Test with number universe
+  const numberUniverse = [1, 2, 3, 4, 5, 6, 7, 8];
+  const numberBitSet = new BitSet<number>(numberUniverse);
+
+  describe("constructor", () => {
+    it("should initialize with correct universe", () => {
+      expect(stringBitSet.universe).toEqual(stringUniverse);
+      expect(numberBitSet.universe).toEqual(numberUniverse);
+    });
+
+    it("should calculate empty as 0", () => {
+      expect(stringBitSet.empty).toBe(0);
+      expect(numberBitSet.empty).toBe(0);
+    });
+
+    it("should calculate all correctly", () => {
+      expect(stringBitSet.all).toBe(0b11111); // 2^5 - 1 = 31
+      expect(numberBitSet.all).toBe(0b11111111); // 2^8 - 1 = 255
+    });
+  });
+
+  describe("create", () => {
+    it("should create empty bitset when no values provided", () => {
+      const bitset = stringBitSet.create();
+      expect(bitset).toBe(0);
+    });
+
+    it("should create bitset from single value", () => {
+      const bitset = stringBitSet.create("a");
+      const reference = new Set(["a"]);
+      expect(stringBitSet.decode(bitset)).toEqual(reference);
+    });
+
+    it("should create bitset from multiple values", () => {
+      const bitset = stringBitSet.create("a", "c", "e");
+      const reference = new Set(["a", "c", "e"]);
+      expect(stringBitSet.decode(bitset)).toEqual(reference);
+    });
+
+    it("should ignore values not in universe", () => {
+      const bitset = stringBitSet.create("a", "z", "b", "x");
+      const reference = new Set(["a", "b"]);
+      expect(stringBitSet.decode(bitset)).toEqual(reference);
+    });
+
+    it("should handle duplicate values", () => {
+      const bitset = stringBitSet.create("a", "a", "b", "b");
+      const reference = new Set(["a", "b"]);
+      expect(stringBitSet.decode(bitset)).toEqual(reference);
+    });
+  });
+
+  describe("encode", () => {
+    it("should encode empty set to 0", () => {
+      const set = new Set<string>();
+      expect(stringBitSet.encode(set)).toBe(0);
+    });
+
+    it("should encode full set correctly", () => {
+      const set = new Set(stringUniverse);
+      expect(stringBitSet.encode(set)).toBe(stringBitSet.all);
+    });
+
+    it("should encode subset correctly", () => {
+      const set = new Set(["b", "d"]);
+      const bitset = stringBitSet.encode(set);
+      expect(stringBitSet.decode(bitset)).toEqual(set);
+    });
+
+    it("should ignore elements not in universe", () => {
+      const set = new Set(["a", "x", "y", "c"]);
+      const bitset = stringBitSet.encode(set);
+      const expected = new Set(["a", "c"]);
+      expect(stringBitSet.decode(bitset)).toEqual(expected);
+    });
+  });
+
+  describe("decode", () => {
+    it("should decode 0 to empty set", () => {
+      expect(stringBitSet.decode(0)).toEqual(new Set());
+    });
+
+    it("should decode all bits to full set", () => {
+      expect(stringBitSet.decode(stringBitSet.all)).toEqual(
+        new Set(stringUniverse),
+      );
+    });
+
+    it("should decode specific bits correctly", () => {
+      // bits: 10101 = positions 0, 2, 4 = 'a', 'c', 'e'
+      const bitset = 0b10101;
+      expect(stringBitSet.decode(bitset)).toEqual(new Set(["a", "c", "e"]));
+    });
+
+    it("should throw error for invalid bits (negative)", () => {
+      expect(() => stringBitSet.decode(-1)).toThrow("invalid bits");
+    });
+
+    it("should throw error for invalid bits (too large)", () => {
+      expect(() => stringBitSet.decode(stringBitSet.all + 1)).toThrow(
+        "invalid bits",
+      );
+    });
+  });
+
+  describe("union", () => {
+    it("should return empty for no arguments", () => {
+      expect(stringBitSet.union()).toBe(0);
+    });
+
+    it("should return same value for single argument", () => {
+      const bitset = stringBitSet.create("a", "b");
+      expect(stringBitSet.union(bitset)).toBe(bitset);
+    });
+
+    it("should compute union correctly", () => {
+      const set1 = new Set(["a", "b"]);
+      const set2 = new Set(["b", "c"]);
+      const set3 = new Set(["d"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+      const bitset3 = stringBitSet.encode(set3);
+
+      const unionBitset = stringBitSet.union(bitset1, bitset2, bitset3);
+      const unionSet = stringBitSet.decode(unionBitset);
+
+      // Reference using standard Set
+      const referenceUnion = new Set([...set1, ...set2, ...set3]);
+      expect(unionSet).toEqual(referenceUnion);
+    });
+
+    it("should handle union with empty sets", () => {
+      const bitset1 = stringBitSet.create("a", "b");
+      const empty = stringBitSet.empty;
+
+      expect(stringBitSet.union(bitset1, empty)).toBe(bitset1);
+      expect(stringBitSet.union(empty, bitset1)).toBe(bitset1);
+    });
+  });
+
+  describe("intersection", () => {
+    it("should return all bits for no arguments", () => {
+      // This is because reduce with no initial value uses first element
+      expect(stringBitSet.intersection()).toBe(stringBitSet.all);
+    });
+
+    it("should return same value for single argument", () => {
+      const bitset = stringBitSet.create("a", "b");
+      expect(stringBitSet.intersection(bitset)).toBe(bitset);
+    });
+
+    it("should compute intersection correctly", () => {
+      const set1 = new Set(["a", "b", "c"]);
+      const set2 = new Set(["b", "c", "d"]);
+      const set3 = new Set(["c", "d", "e"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+      const bitset3 = stringBitSet.encode(set3);
+
+      const intersectionBitset = stringBitSet.intersection(
+        bitset1,
+        bitset2,
+        bitset3,
+      );
+      const intersectionSet = stringBitSet.decode(intersectionBitset);
+
+      // Reference using standard Set
+      const referenceIntersection = new Set(
+        [...set1].filter((x) => set2.has(x) && set3.has(x)),
+      );
+      expect(intersectionSet).toEqual(referenceIntersection);
+    });
+
+    it("should handle intersection with empty set", () => {
+      const bitset1 = stringBitSet.create("a", "b");
+      const empty = stringBitSet.empty;
+
+      expect(stringBitSet.intersection(bitset1, empty)).toBe(empty);
+      expect(stringBitSet.intersection(empty, bitset1)).toBe(empty);
+    });
+
+    it("should handle disjoint sets", () => {
+      const bitset1 = stringBitSet.create("a", "b");
+      const bitset2 = stringBitSet.create("c", "d");
+
+      expect(stringBitSet.intersection(bitset1, bitset2)).toBe(0);
+    });
+  });
+
+  describe("difference", () => {
+    it("should compute difference correctly", () => {
+      const set1 = new Set(["a", "b", "c", "d"]);
+      const set2 = new Set(["b", "d"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      const differenceBitset = stringBitSet.difference(bitset1, bitset2);
+      const differenceSet = stringBitSet.decode(differenceBitset);
+
+      // Reference using standard Set
+      const referenceDifference = new Set(
+        [...set1].filter((x) => !set2.has(x)),
+      );
+      expect(differenceSet).toEqual(referenceDifference);
+    });
+
+    it("should return empty when subtracting from empty", () => {
+      const bitset = stringBitSet.create("a", "b");
+      expect(stringBitSet.difference(stringBitSet.empty, bitset)).toBe(0);
+    });
+
+    it("should return original when subtracting empty", () => {
+      const bitset = stringBitSet.create("a", "b");
+      expect(stringBitSet.difference(bitset, stringBitSet.empty)).toBe(bitset);
+    });
+
+    it("should return empty when subtracting self", () => {
+      const bitset = stringBitSet.create("a", "b", "c");
+      expect(stringBitSet.difference(bitset, bitset)).toBe(0);
+    });
+  });
+
+  describe("isDisjoint", () => {
+    it("should return true for disjoint sets", () => {
+      const set1 = new Set(["a", "b"]);
+      const set2 = new Set(["c", "d"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      expect(stringBitSet.isDisjoint(bitset1, bitset2)).toBe(true);
+    });
+
+    it("should return false for overlapping sets", () => {
+      const set1 = new Set(["a", "b", "c"]);
+      const set2 = new Set(["c", "d", "e"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      expect(stringBitSet.isDisjoint(bitset1, bitset2)).toBe(false);
+    });
+
+    it("should return true when one or both sets are empty", () => {
+      const bitset = stringBitSet.create("a", "b");
+      const empty = stringBitSet.empty;
+
+      expect(stringBitSet.isDisjoint(empty, empty)).toBe(true);
+      expect(stringBitSet.isDisjoint(bitset, empty)).toBe(true);
+      expect(stringBitSet.isDisjoint(empty, bitset)).toBe(true);
+    });
+  });
+
+  describe("isSubset", () => {
+    it("should return true for subset", () => {
+      const set1 = new Set(["a", "b"]);
+      const set2 = new Set(["a", "b", "c", "d"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      expect(stringBitSet.isSubset(bitset1, bitset2)).toBe(true);
+    });
+
+    it("should return false for non-subset", () => {
+      const set1 = new Set(["a", "b", "e"]);
+      const set2 = new Set(["a", "b", "c", "d"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      expect(stringBitSet.isSubset(bitset1, bitset2)).toBe(false);
+    });
+
+    it("should return true for equal sets", () => {
+      const bitset = stringBitSet.create("a", "b", "c");
+      expect(stringBitSet.isSubset(bitset, bitset)).toBe(true);
+    });
+
+    it("should handle empty set as subset", () => {
+      const bitset = stringBitSet.create("a", "b");
+      const empty = stringBitSet.empty;
+
+      expect(stringBitSet.isSubset(empty, bitset)).toBe(true);
+      expect(stringBitSet.isSubset(empty, empty)).toBe(true);
+      expect(stringBitSet.isSubset(bitset, empty)).toBe(false);
+    });
+  });
+
+  describe("isSuperset", () => {
+    it("should return true for superset", () => {
+      const set1 = new Set(["a", "b", "c", "d"]);
+      const set2 = new Set(["a", "b"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      expect(stringBitSet.isSuperset(bitset1, bitset2)).toBe(true);
+    });
+
+    it("should return false for non-superset", () => {
+      const set1 = new Set(["a", "b", "c"]);
+      const set2 = new Set(["a", "b", "d"]);
+
+      const bitset1 = stringBitSet.encode(set1);
+      const bitset2 = stringBitSet.encode(set2);
+
+      expect(stringBitSet.isSuperset(bitset1, bitset2)).toBe(false);
+    });
+
+    it("should return true for equal sets", () => {
+      const bitset = stringBitSet.create("a", "b", "c");
+      expect(stringBitSet.isSuperset(bitset, bitset)).toBe(true);
+    });
+
+    it("should handle empty set", () => {
+      const bitset = stringBitSet.create("a", "b");
+      const empty = stringBitSet.empty;
+
+      expect(stringBitSet.isSuperset(bitset, empty)).toBe(true);
+      expect(stringBitSet.isSuperset(empty, empty)).toBe(true);
+      expect(stringBitSet.isSuperset(empty, bitset)).toBe(false);
+    });
+  });
+
+  describe("edge cases and performance", () => {
+    it("should handle large universe efficiently", () => {
+      const largeUniverse = Array.from({ length: 30 }, (_, i) => i);
+      const largeBitSet = new BitSet(largeUniverse);
+
+      expect(largeBitSet.all).toBe(2 ** 30 - 1);
+
+      const set = new Set([0, 5, 10, 15, 20, 25, 29]);
+      const bitset = largeBitSet.encode(set);
+      expect(largeBitSet.decode(bitset)).toEqual(set);
+    });
+
+    it("should maintain consistency across operations", () => {
+      const a = stringBitSet.create("a", "b");
+      const b = stringBitSet.create("b", "c");
+      const c = stringBitSet.create("c", "d");
+
+      // (A ∪ B) ∩ C = (A ∩ C) ∪ (B ∩ C)
+      const left = stringBitSet.intersection(stringBitSet.union(a, b), c);
+      const right = stringBitSet.union(
+        stringBitSet.intersection(a, c),
+        stringBitSet.intersection(b, c),
+      );
+      expect(left).toBe(right);
+
+      // A - (B ∪ C) = (A - B) ∩ (A - C)
+      const left2 = stringBitSet.difference(a, stringBitSet.union(b, c));
+      const right2 = stringBitSet.intersection(
+        stringBitSet.difference(a, b),
+        stringBitSet.difference(a, c),
+      );
+      expect(left2).toBe(right2);
+    });
+  });
+});