From e09e88e69c322159fe3d18018b95b9acac1dadab Mon Sep 17 00:00:00 2001
From: itamar <itamar@itamar.site>
Date: Fri, 20 Mar 2026 18:34:42 +0000
Subject: [PATCH] push

---
 bigints.nim   |  192 +++++++
 fenwick.nim   |  112 ++++
 graphs.nim    | 1484 +++++++++++++++++++++++++++++++++++++++++++++++++
 nt.nim        |  152 +++++
 segmented.nim |   45 ++
 segtree.nim   |   98 ++++
 sieve.nim     |  250 +++++++++
 7 files changed, 2333 insertions(+)
 create mode 100644 bigints.nim
 create mode 100644 fenwick.nim
 create mode 100644 graphs.nim
 create mode 100644 nt.nim
 create mode 100644 segmented.nim
 create mode 100644 segtree.nim
 create mode 100644 sieve.nim

diff --git a/bigints.nim b/bigints.nim
new file mode 100644
index 0000000..c6e95b9
--- /dev/null
+++ b/bigints.nim
@@ -0,0 +1,192 @@
+import strutils
+
+const BASE = 1_000_000_000'u64
+const BASE_DIGITS = 9
+
+type
+  BigInt* = object
+    sign*: int        # -1, 0, 1
+    digits*: seq[uint32]  # little endian limbs
+
+
+proc normalize(a: var BigInt) =
+  while a.digits.len > 0 and a.digits[^1] == 0:
+    a.digits.setLen(a.digits.len - 1)
+
+  if a.digits.len == 0:
+    a.sign = 0
+  elif a.sign == 0:
+    a.sign = 1
+
+
+proc big*(x: int): BigInt =
+  if x == 0:
+    return BigInt(sign: 0)
+
+  var v = abs(x)
+  result.sign = if x < 0: -1 else: 1
+
+  while v > 0:
+    result.digits.add(uint32(v mod int(BASE)))
+    v = v div int(BASE)
+
+
+proc parse*(s: string): BigInt =
+  var str = s
+  result.sign = 1
+
+  if str[0] == '-':
+    result.sign = -1
+    str = str[1..^1]
+
+  for i in countdown(str.len, 1, BASE_DIGITS):
+    let start = max(0, i - BASE_DIGITS)
+    let chunk = parseInt(str[start..<i])
+    result.digits.add(uint32(chunk))
+
+  result.normalize()
+
+
+proc cmp(a, b: BigInt): int =
+  if a.digits.len != b.digits.len:
+    return cmp(a.digits.len, b.digits.len)
+
+  for i in countdown(a.digits.len - 1, 0):
+    if a.digits[i] != b.digits[i]:
+      return cmp(a.digits[i], b.digits[i])
+
+  return 0
+
+
+proc `==`*(a, b: BigInt): bool =
+  a.sign == b.sign and a.digits == b.digits
+
+
+proc add(a, b: BigInt): BigInt =
+  var carry: uint64 = 0
+  let n = max(a.digits.len, b.digits.len)
+
+  result.digits = newSeq[uint32](n)
+
+  for i in 0..<n:
+    var x = carry
+
+    if i < a.digits.len: x += a.digits[i]
+    if i < b.digits.len: x += b.digits[i]
+
+    result.digits[i] = uint32(x mod BASE)
+    carry = x div BASE
+
+  if carry > 0:
+    result.digits.add(uint32(carry))
+
+  result.sign = 1
+  result.normalize()
+
+
+proc sub(a, b: BigInt): BigInt =
+  var borrow: int64 = 0
+  result.digits = newSeq[uint32](a.digits.len)
+
+  for i in 0..<a.digits.len:
+    var x = int64(a.digits[i]) - borrow
+    if i < b.digits.len:
+      x -= int64(b.digits[i])
+
+    if x < 0:
+      x += int64(BASE)
+      borrow = 1
+    else:
+      borrow = 0
+
+    result.digits[i] = uint32(x)
+
+  result.sign = 1
+  result.normalize()
+
+
+proc `+`*(a, b: BigInt): BigInt =
+  if a.sign == 0: return b
+  if b.sign == 0: return a
+
+  if a.sign == b.sign:
+    result = add(a, b)
+    result.sign = a.sign
+  else:
+    let c = cmp(a, b)
+
+    if c == 0:
+      return BigInt(sign: 0)
+
+    if c > 0:
+      result = sub(a, b)
+      result.sign = a.sign
+    else:
+      result = sub(b, a)
+      result.sign = b.sign
+
+
+proc `-`*(a, b: BigInt): BigInt =
+  var nb = b
+  nb.sign *= -1
+  a + nb
+
+
+proc `*`*(a, b: BigInt): BigInt =
+  if a.sign == 0 or b.sign == 0:
+    return BigInt(sign: 0)
+
+  result.digits = newSeq[uint32](a.digits.len + b.digits.len)
+
+  for i in 0..<a.digits.len:
+    var carry: uint64 = 0
+
+    for j in 0..<b.digits.len:
+      let k = i + j
+
+      var cur =
+        uint64(result.digits[k]) +
+        uint64(a.digits[i]) * uint64(b.digits[j]) +
+        carry
+
+      result.digits[k] = uint32(cur mod BASE)
+      carry = cur div BASE
+
+    var pos = i + b.digits.len
+
+    while carry > 0:
+      let cur = uint64(result.digits[pos]) + carry
+      result.digits[pos] = uint32(cur mod BASE)
+      carry = cur div BASE
+      inc pos
+
+  result.sign = a.sign * b.sign
+  result.normalize()
+
+
+proc `$`*(a: BigInt): string =
+  if a.sign == 0:
+    return "0"
+
+  var parts: seq[string]
+
+  for d in a.digits:
+    parts.add($d)
+
+  result = parts[^1]
+
+  for i in countdown(parts.len - 2, 0):
+    result.add(parts[i].align(BASE_DIGITS, '0'))
+
+  if a.sign < 0:
+    result = "-" & result
+
+
+when isMainModule:
+  let a = parse("123456789123456789123456789")
+  let b = parse("987654321987654321987654321")
+
+  echo "a = ", a
+  echo "b = ", b
+  echo "a + b = ", a + b
+  echo "a * b = ", a * b
\ No newline at end of file
diff --git a/fenwick.nim b/fenwick.nim
new file mode 100644
index 0000000..24bd434
--- /dev/null
+++ b/fenwick.nim
@@ -0,0 +1,112 @@
+type Fenwick*[T] = ref object
+  arr*: seq[T]
+
+type FenwickMod*[T] = ref object
+  arr*: seq[T]
+  m*: int64
+
+proc Fenwick*[T](f: var Fenwick[T], len: int): void =
+  new(f)
+  f.arr.newSeq(len)
+
+proc Fenwick*[T](f: var Fenwick[T], len: int, default: T): void =
+  ##ializes a fenwick tree with a constant array, f[i] = default for all i.
+  new(f)
+  f.arr.newSeq(len)
+  for i in 0..<len:
+    f.arr[i] = default
+  for i in 1..len:
+    var j = i + (i and (-i))
+    if j<=len: f.arr[j-1] += f.arr[i-1]
+
+proc Fenwick*[T](f: var Fenwick[T], default: seq[T]): void =
+  ##Converts the seq into a fenwick tree in O(default.len) time.
+  new(f)
+  f.arr = default
+  for i in 1..default.len:
+    var j = i + (i and (-i))
+    if j<=default.len: f.arr[j-1] += f.arr[i-1]
+
+proc Fenwick*[T](f: var FenwickMod[T], len: int, m: int64): void =
+  new(f)
+  f.arr.newSeq(len)
+  f.m = m
+
+proc Fenwick*[T](f: var FenwickMod[T], len: int, default: T): void =
+  ##ializes a fenwick tree with a constant array, f[i] = default for all i.
+  new(f)
+  f.arr.newSeq(len)
+  for i in 0..<len:
+    f.arr[i] = default
+  for i in 1..len:
+    var j = i + (i and (-i))
+    if j<=len: f.arr[j-1] = (f.arr[j-1] + f.arr[i-1]) mod f.m
+    
+proc Fenwick*[T](f: var FenwickMod[T], default: seq[T]): void =
+  ##Converts the seq into a fenwick tree in O(default.len) time.
+  new(f)
+  f.arr = default
+  for i in 1..len:
+    var j = i + (i and (-i))
+    if j<=len: f.arr[j-1] = (f.arr[j-1] + f.arr[i-1]) mod f.m
+
+type SomeFenwick*[T] = Fenwick[T] | FenwickMod[T]
+
+
+proc len*[T](f: SomeFenwick[T]): int = f.arr.len
+
+proc sum*[T](f: Fenwick[T], i: SomeInteger): T =
+  ##Returns f[0] + f[1] + ... + f[i]. Time O(log i).
+  var ii = i+1
+  while ii>0:
+    result += f.arr[ii-1]
+    ii -= (ii and (-ii))
+
+proc sum*[T](f: FenwickMod[T], i: SomeInteger): int64 =
+  ##Returns f[0] + f[1] + ... + f[i]. Time O(log i).
+  var ii = i+1
+  while ii>0:
+    result += f.arr[ii-1]
+    result = result mod f.m
+    ii -= (ii and (-ii))
+
+proc addTo*[T](f: Fenwick[T], i: SomeInteger, x: T): void =
+  var ii = i+1
+  while ii<=f.arr.len:
+    f.arr[ii-1] += x
+    ii += (ii and (-ii))
+
+proc addTo*[T](f: var Fenwick[T], vals: seq[T]): void =
+  ##Adds vals[i] to f[i] for each i. Time O(vals.len).
+  var vals = vals
+  for i in 1..vals.len:
+    var j = i + (i and (-i))
+    if j<=vals.len: vals[j-1] = vals[j-1] + vals[i-1]
+    f.arr[i-1] += vals[i-1]
+
+proc addTo*[T](f: FenwickMod[T], i: SomeInteger, x: T): void =
+  var ii = i+1
+  while ii<=f.arr.len:
+    f.arr[ii-1] += x
+    f.arr[ii-1] = (f.arr[ii-1] mod f.m).T
+    ii += (ii and (-ii))
+
+proc addTo*[T](f: var FenwickMod[T], vals: seq[T]): void =
+  ##Adds vals[i] to f[i] for each i. Time O(vals.len).
+  var vals = vals
+  for i in 1..vals.len:
+    var j = i + (i and (-i))
+    if j<=vals.len: vals[j-1] = (vals[j-1] + vals[i-1]) mod f.m
+    f.arr[i-1] = (f.arr[i-1] + vals[i-1]) mod f.m
+
+proc `[]`*[T](f: Fenwick[T], i: SomeInteger): T =
+  if i==0: return f.sum(0)
+  return f.sum(i) - f.sum(i-1)
+
+proc `[]`*[T](f: FenwickMod[T], i: SomeInteger): T =
+  ##Accesses a single element of the base array. O(log i)
+  if i==0: return f.sum(0)
+  return (f.sum(i) - f.sum(i-1) + f.m) mod f.m
+
+proc `[]=`*[T](f: SomeFenwick[T], i: SomeInteger, x: T): void =
+  f.addTo(i, x-f[i])
\ No newline at end of file
diff --git a/graphs.nim b/graphs.nim
new file mode 100644
index 0000000..7bf5153
--- /dev/null
+++ b/graphs.nim
@@ -0,0 +1,1484 @@
+import linear,iops, sets, tables, algorithm, heapqueue, random
+export linear
+
+type GraphRepr = enum
+  AdjMatrix
+  Sparse
+  Dense
+
+type Graph* = ref object
+  #representation of vertices
+  V*: int ##number of vertices.
+  capacity: int #number of possible vertices.
+  E*: int ##number of edges.
+  directed*: bool
+  representation: GraphRepr
+  #AdjMatrix
+  adjMatrix: matrix[int]
+  #Sparse / Dense
+  adjList: seq[HashSet[int]]
+  inverseAdjList: seq[HashSet[int]] #only used for digraph
+
+func contains*(g: Graph, v, w: int): bool =
+  ##Test if v ==> w is an edge in the graph g.
+  if v >= g.V or w >= g.V: return false
+  case g.representation:
+  of AdjMatrix:
+    return g.adjMatrix[v, w] > 0
+  of Sparse:
+    return w in g.adjList[v]
+  of Dense:
+    return not (w in g.adjList[v])
+
+func outDegree*(g: Graph, v: int): int =
+  if v >= g.V: return 0
+  case g.representation:
+  of AdjMatrix:
+    for w in 0..<g.V:
+      result += g.adjMatrix[v, w]
+    return result
+  of Sparse:
+    return g.adjList[v].len
+  of Dense:
+    return g.V - g.adjList[v].len
+
+func inDegree*(g: Graph, v: int): int =
+  if v >= g.V: return 0
+  if not g.directed: return g.outDegree(v) #easier
+  case g.representation:
+  of AdjMatrix:
+    for w in 0..<g.V:
+      result += g.adjMatrix[w, v]
+    return result
+  of Sparse:
+    return g.inverseAdjList[v].len
+  of Dense:
+    return g.V - g.inverseAdjList[v].len
+
+func degree*(g: Graph, v: int): int =
+  if v >= g.V: return 0
+  if g.directed: return g.inDegree(v) + g.outDegree(v)
+  else: return g.outDegree(v)
+
+func edgeDensity*(g: Graph): float =
+  var E = g.E.float
+  var V = g.V.float
+  if g.directed:
+    return E / (V * (V-1))
+  else:
+    return 2*E / (V * (V-1))
+
+func adjacencyMatrix*(g: Graph): matrix[int] =
+  ##Returns the V x V adjacency matrix for g.
+  ##This may be different from the internal adjacency matrix if g.capacity > g.V.
+  result.initMatrix(g.V, g.V)
+  case g.representation
+  of AdjMatrix:
+    #trim it down
+    for v in 0..<g.V:
+      for w in 0..<g.V:
+        if v != w:
+          result[v, w] = g.adjMatrix[v, w]
+  of Sparse:
+    for v in 0..<g.V:
+      for w in g.adjList[v]:
+        if v != w:
+          result[v, w] += 1
+  of Dense:
+    for v in 0..<g.V:
+      for w in 0..<g.V:
+        if v != w:
+          result[v, w] = 1
+      for w in g.adjList[v]:
+        result[v, w] -= 1
+
+func initGraph*(V: int = 0): Graph =
+  ##This returns a graph with V vertices and 0 edges.
+  new(result)
+  result.V = V
+  result.capacity = V
+  result.E = 0
+  result.directed = false
+  #empty graph - thus Sparse.
+  result.representation = Sparse
+  result.adjList = newSeq[HashSet[int]](V)
+  for i in 0..<V:
+    result.adjList[i] = initHashSet[int]()
+
+func initDenseGraph*(V: int = 0): Graph =
+  new(result)
+  result.V = V
+  result.capacity = V
+  result.E = (V*(V-1)) div 2
+  result.directed = false
+  #the densest dense graph is the complete graph.
+  result.representation = Dense
+  result.adjList = newSeq[HashSet[int]](V)
+  for i in 0..<V:
+    result.adjList[i] = initHashSet[int]()
+
+func initDigraph*(V: int = 0): Graph =
+  new(result)
+  result.V = V
+  result.capacity = V
+  result.E = 0
+  result.directed = true
+  #empty graph - thus Sparse.
+  result.representation = Sparse
+  result.adjList = newSeq[HashSet[int]](V)
+  result.inverseAdjList = newSeq[HashSet[int]](V)
+  for i in 0..<V:
+    result.adjList[i] = initHashSet[int]()
+    result.inverseAdjList[i] = initHashSet[int]()
+
+func initDenseDigraph*(V: int = 0): Graph =
+  new(result)
+  result.V = V
+  result.capacity = V
+  result.E = (V*(V-1)) div 2
+  result.directed = true
+  #the densest dense graph is the complete graph.
+  result.representation = Dense
+  result.adjList = newSeq[HashSet[int]](V)
+  result.inverseAdjList = newSeq[HashSet[int]](V)
+  for i in 0..<V:
+    result.adjList[i] = initHashSet[int]()
+    result.inverseAdjList[i] = initHashSet[int]()
+
+proc regularize(g: var Graph) =
+  #actually implement this later lmao
+  return
+  #[
+    AdjMatrix:
+      O(V^2) space
+      O(1) time to check/modify (v, w)
+      O(1) time to get adjacency matrix
+      O(V) time to get deg(v)
+      Converting into:
+        Sparse/Dense:
+          O(V^2) time
+    Sparse/Dense:
+      O(2E) space
+      O(deg(v)) to check/modify (v, w)
+      O(V^2) time to get adjacency matrix
+      O(1) time to get deg(v)
+      Converting into:
+        Dense/Sparse:
+          O(V^2) time
+        AdjMatrix:
+          O(V^2) time
+  ]#
+  # var density = g.edgeDensity()
+  # case g.representation:
+  # of AdjMatrix:
+  #   if density <= 0.05:
+  #     #switch to sparse
+  #     # echo "Switching to sparse (density is ", density, ")"
+  #     g.representation = Sparse
+  #     g.adjList = newSeq[HashSet[int]](g.capacity)
+  #     for i in 0..<g.capacity:
+  #       g.adjList[i] = initHashSet[int]()
+  #     for v in 0..<g.V:
+  #       for w in 0..<g.V:
+  #         for i in 1..g.adjMatrix[v, w]:
+  #           g.adjList[v].incl w
+  #     if g.directed:
+  #       g.inverseAdjList = newSeq[HashSet[int]](g.capacity)
+  #       for i in 0..<g.capacity:
+  #         g.inverseAdjList[i] = initHashSet[int]()
+  #       for v in 0..<g.V:
+  #         for w in 0..<g.V:
+  #           for i in 1..g.adjMatrix[v, w]:
+  #             g.inverseAdjList[w].incl v
+  #     g.adjMatrix = nil
+  #   if density >= 0.95:
+  #     return #TODO TESTING
+  #     #switch to dense
+  #     #automatically reduce capacity as much as possible
+  #     g.representation = Dense
+  #     g.capacity = g.V
+  #     g.adjList = newSeq[HashSet[int]](g.V)
+  #     for i in 0..<g.V:
+  #       g.adjList[i] = initHashSet[int]()
+  #     for v in 0..<g.V:
+  #       for w in 0..<g.V:
+  #         if g.adjMatrix[v, w] == 0:
+  #           g.adjList[v].incl w
+  #     if g.directed:
+  #       g.inverseAdjList = newSeq[HashSet[int]](g.V)
+  #       for i in 0..<g.V:
+  #         g.inverseAdjList[i] = initHashSet[int]()
+  #       for v in 0..<g.V:
+  #         for w in 0..<g.V:
+  #           if g.adjMatrix[v, w] == 0:
+  #             g.inverseAdjList[w].incl v
+  #     g.adjMatrix = nil
+  # of Sparse:
+  #   if density >= 0.10:
+  #     # echo "Switching to adjmatrix (density is ", density, ")"
+  #     #convert to AdjMatrix
+  #     g.adjMatrix = g.adjacencyMatrix()
+  #     g.capacity = g.V #cuts down capacity automatically
+  #     g.adjList = @[]
+  #     if g.directed: g.inverseAdjList = @[]
+  #     g.representation = AdjMatrix
+  # of Dense:
+  #   if density <= 0.90:
+  #     # echo "Switching to adjmatrix (density is ", density, ")"
+  #     #convert to AdjMatrix
+  #     g.adjMatrix = g.adjacencyMatrix()
+  #     g.capacity = g.V #cuts down capacity automatically
+  #     g.adjList = @[]
+  #     if g.directed: g.inverseAdjList = @[]
+  #     g.representation = AdjMatrix
+
+proc enlarge(g: var Graph, capNew: int) =
+  var capNew = capNew
+  if capNew == 0: capNew = 1
+  # echo "Graph (", g.representation, ") expanding. Capacity ", g.capacity, " ==> ", capNew
+  case g.representation:
+  of AdjMatrix:
+    var newMatrix: matrix[int]
+    newMatrix.initMatrix(capNew, capNew)
+    for v in 0..<g.capacity:
+      for w in 0..<g.capacity:
+        newMatrix[v, w] = g.adjMatrix[v, w]
+    g.adjMatrix = newMatrix
+  of Sparse:
+    for v in g.capacity..<capNew:
+      g.adjList.add initHashSet[int]()
+      if g.directed:
+        g.inverseAdjList.add initHashSet[int]()
+  of Dense:
+    var all = initHashSet[int]()
+    for i in 0..<capNew: all.incl i
+    for v in g.capacity..<capNew:
+      g.adjList.add all
+      if g.directed:
+        g.inverseAdjList.add all
+    for i in 0..<g.capacity:
+      for v in g.capacity..<capNew:
+        g.adjList[i].incl v
+        if g.directed:
+          g.inverseAdjList[i].incl v
+  g.capacity = capNew
+
+proc connect*(g: var Graph, v, w: int) =
+  if g.contains(v, w) or v==w: return
+
+  #verify bounds
+  while v >= g.capacity or w >= g.capacity:
+    g.enlarge(g.capacity * 2)
+  if v >= g.V or w >= g.V:
+    #increase g.V
+    g.V = max(v, w) + 1
+  
+  inc g.E
+
+  if g.directed:
+    case g.representation:
+    of AdjMatrix:
+      g.adjMatrix[v, w] = 1
+    of Sparse:
+      g.adjList[v].incl w
+      g.inverseAdjList[w].incl v
+    of Dense: #remove (v, w) if it exists, and remove (w, v) from inverseAdjList
+      g.adjList[v].excl w
+      g.inverseAdjList[w].excl v
+  else: #undirected
+    case g.representation:
+    of AdjMatrix:
+      g.adjMatrix[v, w] = 1
+      g.adjMatrix[w, v] = 1
+    of Sparse:
+      g.adjList[v].incl w
+      g.adjList[w].incl v
+    of Dense: #remove (v, w) and (w, v) if it exists
+      g.adjList[v].excl(w)
+      g.adjList[w].excl(v)
+  g.regularize()
+
+proc disconnect*(g: var Graph, v, w: int) =
+  if not g.contains(v, w): return
+
+  #verify bounds
+  while v >= g.capacity or w >= g.capacity:
+    g.enlarge(g.capacity * 2)
+  if v >= g.V or w >= g.V:
+    #increase g.V
+    g.V = max(v, w) + 1
+
+  dec g.E
+
+  if g.directed:
+    case g.representation:
+    of AdjMatrix:
+      g.adjMatrix[v, w] = 0
+    of Sparse:
+      g.adjList[v].excl w
+      g.inverseAdjList[w].excl v
+    of Dense:
+      g.adjList[v].incl w
+      g.inverseAdjList[w].incl v
+  else: #undirected
+    case g.representation:
+    of AdjMatrix:
+      g.adjMatrix[v, w] = 0
+      g.adjMatrix[w, v] = 0
+    of Sparse:
+      g.adjList[v].excl w
+      g.adjList[w].excl v
+    of Dense: #remove (v, w) and (w, v) if it exists
+      g.adjList[v].incl(w)
+      g.adjList[w].incl(v)
+  g.regularize()
+
+template defineGraph*(graphName: untyped, V: untyped, body: untyped): untyped =
+  var `graphName` {.inject.} = initGraph(V)
+  block:
+    proc `==>`(v, w: int) {.used.} =
+      `graphName`.connect(v, w)
+    proc cycle(vs: openArray[int]) {.used.} = 
+      for i in 0..vs.high:
+        vs[i] ==> vs[(i+1) mod vs.len]
+    body
+
+template defineDigraph*(graphName: untyped, V: untyped, body: untyped): untyped =
+  var `graphName` {.inject.} = initDigraph(V)
+  proc `==>`(v, w: int) {.used.} =
+    `graphName`.connect(v, w)
+  proc cycle(vs: openArray[int]) {.used.} =
+    for i in 0..vs.high:
+      vs[i] ==> vs[(i+1) mod vs.len]
+  body
+
+#[
+  =====================
+   ADVANCED OPERATIONS
+  =====================
+]#
+
+iterator edges*(g: Graph): (int, int) =
+  case g.representation:
+  of AdjMatrix:
+    for v in 0..<g.V:
+      var wStart = 0
+      if not g.directed: wStart = v+1
+      for w in wStart..<g.V:
+        if w == v: continue
+        for i in 1..g.adjMatrix[v, w]: yield (v, w)
+  of Sparse: 
+    for v in countdown(g.V - 1, 0):
+      for w in g.adjList[v]: 
+        if w > v or g.directed: yield (v, w)
+  of Dense:
+    var forbidden = newSeq[bool](g.V)
+    for v in 0..<g.V:
+      for w in 0..<g.V: forbidden[w] = false
+      for w in g.adjList[v]: forbidden[w] = true
+      var wStart = 0
+      if not g.directed: wStart = v+1
+      for w in wStart..<g.V: 
+        if not forbidden[w]: yield (v, w)
+
+iterator edges*(g: Graph, v: int): int =
+  ##Iterates over the edges (v, w), where you specify v.
+  case g.representation:
+  of AdjMatrix:
+    for w in 0..<g.V:
+      if w == v: continue
+      for i in 1..g.adjMatrix[v, w]: yield w
+  of Sparse: 
+    for w in g.adjList[v]: yield w
+  of Dense:
+    var forbidden = newSeq[bool](g.V)
+    for w in 0..<g.V: forbidden[w] = false
+    for w in g.adjList[v]: forbidden[w] = true
+    for w in 0..<g.V: 
+      if not forbidden[w] and w != v: yield w
+
+iterator preimage*(g: Graph, v: int): int =
+  ##Iterates over the edges (w, v), where you specify v.
+  ##This is the same as g.edges(v) if g is undirected.
+  if not g.directed: 
+    for w in g.edges(v): yield w
+
+  case g.representation:
+  of AdjMatrix:
+    for w in 0..<g.V:
+      if w == v: continue
+      for i in 1..g.adjMatrix[w, v]: yield w
+  of Sparse: 
+    for w in 0..<g.V:
+      if v in g.adjList[w]: yield w
+  of Dense:
+    var forbidden = newSeq[bool](g.V)
+    for w in 0..<g.V: forbidden[w] = false
+    for w in g.inverseAdjList[v]: forbidden[w] = true
+    for w in 0..<g.V: 
+      if not forbidden[w] and w != v: yield w
+
+
+proc shortestPathLengths*(g: Graph): matrix[int] =
+  ##result[v, w] is the lowest length path from v to w, or -1 if none exists.
+  ##Uses Floyd-Warshall, runs in O(V^3).
+  result.initMatrix(g.V, g.V)
+  for v in 0..<g.V:
+    for w in 0..<g.V:
+      result[v, w] = -1
+  for (v, w) in g.edges:
+    result[v, w] = 1
+  for v in 0..<g.V:
+    result[v, v] = 0
+  for k in 0..<g.V:
+    for i in 0..<g.V:
+      for j in 0..<g.V:
+        if result[i, k] >= 0 and result[k, j] >= 0 and (result[i, j] < 0 or result[i, j] > result[i, k] + result[k, j]):
+          result[i, j] = result[i, k] + result[k, j]
+
+proc shortestPathCosts*[T](g: Graph, weights: matrix[T]): matrix[T] =
+  ##Given the set of edge weights, result[v, w] is the lowest cost path from v to w.
+  ##Uses Floyd-Warshall, runs in O(V^3).
+  result.initMatrix(g.V, g.V)
+  for v in 0..<g.V:
+    for w in 0..<g.V:
+      result[v, w] = Inf
+  for (v, w) in g.edges:
+    result[v, w] = weights[v, w]
+  for v in 0..<g.V:
+    result[v, v] = weights[0, 0] - weights[0, 0]
+  for k in 0..<g.V:
+    for i in 0..<g.V:
+      for j in 0..<g.V:
+        if result[i, j] > result[i, k] + result[k, j]:
+          result[i, j] = result[i, k] + result[k, j]
+
+proc clone*(g: Graph): Graph =
+  new(result)
+  result.representation = g.representation
+  result.directed = g.directed
+  result.V = g.V
+  result.capacity = g.capacity
+  result.E = g.E
+  case g.representation
+  of AdjMatrix:
+    result.adjMatrix = g.adjMatrix.copy
+  of Sparse, Dense:
+    result.adjList = g.adjList
+
+proc countEdges(g: Graph): int =
+  #used internally
+  case g.representation:
+  of AdjMatrix:
+    for u in 0..<g.V:
+      for v in 0..<g.V:
+        if g.adjMatrix[u, v] > 0: result += 1
+  of Sparse:
+    for u in 0..<g.V:
+      result += g.adjList[u].len
+  of Dense:
+    for u in 0..<g.V:
+      result += g.V - g.adjList[u].len
+  if not g.directed: return result shr 1
+  return result
+
+proc contract*(g: Graph, v, w: int, simple: bool = true): Graph =
+  var gc = g.clone
+  if v == w: return gc
+  var (v, w) = (v, w)
+  if v > w: 
+    (v, w) = (w, v)
+  #relabel w as v
+  #relabel g.V - 1 as w
+  #track number of edges lost
+
+  #[
+    any edge v -> w or w -> v will be lost.
+    any edge u -> w will be turned into an edge u -> v.
+    any edge w -> u will be turned into an edge v -> u.
+  ]#
+  case gc.representation
+  of AdjMatrix:
+    for u in 0..<gc.V:
+      if u == v: continue
+      if gc.adjMatrix[w, u] > 0:
+        gc.adjMatrix[v, u] = 1
+      if gc.adjMatrix[u, w] > 0:
+        gc.adjMatrix[u, w] = 0
+        gc.adjMatrix[u, v] = 1
+    for u in 0..<gc.V:
+      gc.adjMatrix[w, u] = gc.adjMatrix[gc.V - 1, u]
+      gc.adjMatrix[gc.V - 1, u] = 0
+  of Sparse:
+    gc.adjList[w].excl v
+    gc.adjList[v].excl w
+    for u in gc.adjList[w]:
+      if u == v: continue
+      #edge w ==> u
+      #change it to an edge v ==> u
+      gc.adjList[v].incl u
+    for u in 0..<gc.V:
+      if u == v: continue
+      if w in gc.adjList[u]:
+        #edge u ==> w
+        #change it to an edge u ==> v
+        gc.adjList[u].excl w
+        gc.adjList[u].incl v
+    gc.adjList[w] = gc.adjList[gc.V - 1]
+    for u in 0..<gc.V:
+      if gc.V - 1 in gc.adjList[u]:
+        gc.adjList[u].excl(g.V - 1)
+        gc.adjList[u].incl w
+    gc.adjList[gc.V - 1] = initHashSet[int]()
+  of Dense:
+    for u in gc.adjList[w]:
+      #edge w ==> u
+      #change it to an edge v ==> u
+      gc.adjList[v].incl u
+    for u in 0..<g.V:
+      if w in gc.adjList[u]:
+        #edge u ==> w
+        #change it to an edge u ==> v
+        gc.adjList[u].excl w
+        gc.adjList[u].incl v
+    gc.adjList[w] = gc.adjList[gc.V - 1]
+    gc.adjList[gc.V - 1] = initHashSet[int]()
+  #TODO for adj and dense
+  dec gc.V
+  gc.E = gc.countEdges
+  return gc
+
+proc `/`*(g: Graph, vw: (int, int)): Graph = g.contract(vw[0], vw[1])
+
+proc `*`*(g, h: Graph): Graph =
+  ##Returns the direct product of the graphs.
+  if not g.directed and not h.directed:
+    defineGraph(output, g.V * h.V):
+      for e1 in g.edges:
+        for e2 in h.edges:
+          e1[0] * h.V + e2[0] ==> e1[1] * h.V + e2[1]
+          e1[1] * h.V + e2[0] ==> e1[0] * h.V + e2[1]
+    return output
+  elif g.directed and h.directed:
+    defineDigraph(output, g.V * h.V):
+      for e1 in g.edges:
+        for e2 in h.edges:
+          e1[0] * h.V + e2[0] ==> e1[1] * h.V + e2[1]
+    return output
+  else:
+    echo "Graph product not yet implemented."
+
+proc directProduct*(g, h: Graph): Graph = g*h
+
+proc cartesianProduct*(g, h: Graph): Graph =
+  ##Returns the cartesian product of the graphs.
+  if not g.directed and not h.directed:
+    defineGraph(output, g.V * h.V):
+      for e1 in g.edges:
+        for v in 0..<h.V:
+          e1[0] * h.V + v ==> e1[1] * h.V + v
+      for e1 in h.edges:
+        for v in 0..<g.V:
+          e1[0] + v * h.V ==> e1[1] + v * h.V
+    return output
+  elif g.directed and h.directed:
+    defineDigraph(output, g.V * h.V):
+      for e1 in g.edges:
+        for v in 0..<h.V:
+          e1[0] * h.V + v ==> e1[1] * h.V + v
+      for e1 in h.edges:
+        for v in 0..<g.V:
+          e1[0] + v * h.V ==> e1[1] + v * h.V
+    return output
+
+proc lexProduct*(g, h: Graph): Graph =
+  ##Returns the lexicographical product of the graphs.
+  if not g.directed and not h.directed:
+    defineGraph(output, g.V * h.V):
+      for e1 in g.edges:
+        for v in 0..<h.V:
+          for w in 0..<h.V:
+            e1[0] * h.V + v ==> e1[1] * h.V + w
+      for e1 in h.edges:
+        for v in 0..<g.V:
+          e1[0] + v * h.V ==> e1[1] + v * h.V
+    return output
+  elif g.directed and h.directed:
+    defineDigraph(output, g.V * h.V):
+      for e1 in g.edges:
+        for v in 0..<h.V:
+          for w in 0..<h.V:
+            e1[0] * h.V + v ==> e1[1] * h.V + w
+      for e1 in h.edges:
+        for v in 0..<g.V:
+          e1[0] + v * h.V ==> e1[1] + v * h.V
+    return output
+  else:
+    echo "Graph product not yet implemented."
+
+func `+`*(g, h: Graph): Graph =
+  ##Returns the disjoint sum of the graphs g and h.
+  if g.directed == h.directed:
+    result = g.clone
+    for e in h.edges:
+      result.connect e[0] + g.V, e[1] + g.V
+  #else: TODO
+
+proc union*(g, h: Graph): Graph =
+  ##Returns the union of the graphs g and h.
+  if g.directed == h.directed:
+    result = g.clone
+    for e in h.edges:
+      result.connect e[0], e[1]
+  else:
+    echo "Graph union not yet implemented."
+
+func complement*(g: Graph): Graph =
+  ##Returns the graph with the same vertices, but with opposite edges.
+  result = g.clone
+  case g.representation:
+  of AdjMatrix:
+    for v in 0..<g.V:
+      for w in 0..<g.V:
+        if w == v: continue
+        result.adjMatrix[v, w] = 1 - result.adjMatrix[v, w]
+  of Sparse:
+    result.representation = Dense
+  of Dense:
+    result.representation = Sparse
+  result.E = ((g.V * (g.V - 1)) shr 1) - g.E
+
+func toDigraph*(g: Graph): Graph =
+  ##Returns the directed version of the graph g.
+  ##If g is directed it will return a clone of g.
+  if g.directed: return g.clone
+
+  case g.representation:
+  of AdjMatrix:
+    result = g.clone
+    result.directed = true
+  of Sparse:
+    result = initDigraph(g.V)
+    for e in g.edges:
+      result.connect e[0], e[1]
+      result.connect e[1], e[0]
+  of Dense:
+    result = initDenseDigraph(g.V)
+    #do this manually because there isn't a great g.nonEdges
+    for v in countdown(g.V - 1, 0):
+      for w in g.adjList[v]: 
+        result.disconnect v, w
+        result.disconnect w, v
+  return result
+
+func toGraph*(g: Graph): Graph = #mostly copied from toDigraph.
+  ##Returns the UNDIRECTED version of the graph g.
+  ##The edge v <-> w will be present iff v ==> w or w ==> v.
+  ##If g is undirected it will return a clone of g.
+  if not g.directed: return g.clone
+  
+  case g.representation:
+  of AdjMatrix:
+    result = g.clone
+    for v in 0..<g.V:
+      for w in v..<g.V:
+        result.adjMatrix[v, w] = max(result.adjMatrix[v, w], result.adjMatrix[w, v])
+    result.directed = false
+  of Sparse:
+    result = initGraph(g.V)
+    for e in g.edges:
+      result.connect e[0], e[1]
+  of Dense:
+    result = initDenseGraph(g.V)
+    #do this manually because there isn't a great g.nonEdges
+    for v in countdown(g.V - 1, 0):
+      for w in g.adjList[v]: 
+        result.disconnect v, w
+  return result
+
+func wolframFormat*(g: Graph): string =
+  ##Converts the graph into a string parseable by Wolfram Language.
+  result = "Graph[{"
+  var comma = false
+  for e in g.edges:
+    if comma: result &= ","
+    else: comma = true
+    if g.directed: result &= $(e[0]) & "->" & $(e[1])
+    else: result &= $(e[0]) & "<->" & $(e[1])
+  return result & "}]"
+
+func lineGraph*(g: Graph): Graph =
+  ##This returns a graph (isomorphic to) the graph whose vertices are the edges of the original graph,
+  ##and two new vertices are connected if the corresponding edges are incident at a vertex.
+  var g = g.toGraph() #forget direction
+  result = initGraph(g.E)
+  #label each edge in a matrix
+  var edgeLabels: matrix[int]
+  block:
+    var i = 0
+    edgeLabels.initMatrix(g.V, g.V)
+    for r in 0..<g.V:
+      for c in 0..<g.V:
+        edgeLabels[r, c] = -1
+    for e in g.edges:
+      edgeLabels[e[0], e[1]] = i
+      edgeLabels[e[1], e[0]] = i
+      inc i
+  #edges are labeled
+  for v in 0..<g.V:
+    var clique = newSeq[int]()
+    for e in g.edges(v): clique.add e
+    for w in 0..clique.high:
+      for u in w+1..clique.high:
+        result.connect edgeLabels[v, clique[w]], edgeLabels[v, clique[u]]
+
+  return result
+
+func del*(g: Graph, v: int): Graph =
+  ##Returns the graph obtained by removing the vertex v.
+  ##The vertex g.V - 1 will be relabeled.
+  #TODO
+
+func inducedSubgraph*(g: Graph, vertices: openArray[int]): Graph =
+  ##Returns the graph obtained by keeping only those vertices.
+  ##All vertices will be relabeled.
+  #TODO
+
+func degreeSequence*(g: Graph): seq[int] =
+  ##Returns the sequence of degrees of the vertices of g.
+  ##The vertex v has degree result[v].
+  result = newSeq[int](g.V)
+  case g.representation:
+  of AdjMatrix:
+    for v in 0..<g.V:
+      for w in 0..<g.V:
+        result[v] += g.adjacencyMatrix[v, w]
+  of Sparse:
+    for v in 0..<g.V:
+      result[v] = g.adjList[v].len
+      if g.directed: result[v] += g.inverseAdjList[v].len
+  of Dense:
+    for v in 0..<g.V:
+      result[v] = g.V - g.adjList[v].len
+
+func degreeMultiset*(g: Graph): seq[int] =
+  ##Returns the sorted sequence of degrees of the vertices of g.
+  ##In this implementation they are listed from smallest to largest.
+  result = g.degreeSequence
+  result.sort
+
+proc connectedComponents*(g: Graph): seq[Graph] =
+  ##Returns a sequence of connected graphs whose disjoint union is isomorphic to g.
+  result = @[]
+  var ids = newSeq[int](g.V)
+  var seen = newSeq[bool](g.V)
+  for v in 0..<g.V:
+    if seen[v]: continue
+
+    var component: Graph
+    if g.directed: component = initDigraph(1)
+    else: component = initGraph(1)
+
+    var chunk = initHashSet[int]()
+    chunk.incl v
+    #ids[v] = 0
+    var u = 0
+    while chunk.len > 0:
+      var w = chunk.pop
+      if seen[w]: continue
+      ids[w] = u
+      inc u
+      for z in g.edges(w):
+        if seen[z]: component.connect ids[w], ids[z]
+        else: 
+          chunk.incl z
+      for z in g.preimage(w):
+        if seen[z]: component.connect ids[z], ids[w]
+        else: chunk.incl z
+      seen[w] = true
+    result.add component
+
+
+
+
+func spanningForest*(g: Graph): seq[Graph] =
+  ##Returns one spanning tree for each connected component of g.
+  #TODO
+
+func grundyValues*(g: Graph): seq[int] =
+  ##Treat the vertices as states of an impartial game.
+  ##Then this returns the grundy values of the states of the game.
+  ##result[v] = 0 if and only if v is a losing state under normal play.
+  ##This assumes that g is an acyclic directed graph.
+  if not g.directed: 
+    return @[]
+
+  result = newSeq[int](g.V)
+  var following = newSeq[seq[int]](g.V)
+  for v in 0..<g.V: following[v] = @[]
+  var toProcess = initHashSet[int]()
+  for v in 0..<g.V:
+    if g.outDegree(v) == 0: toProcess.incl v
+  while toProcess.len > 0:
+    for v in toProcess:
+      if g.outDegree(v) == following[v].len:
+        #we can compute the grundy value of v and add its preimage
+        var u = 0
+        while u in following[v]: inc u
+        result[v] = u
+        for w in g.preimage(v):
+          following[w].add u
+          toProcess.incl w
+        toProcess.excl v
+        break
+
+#[
+  ===============
+   GRAPH QUERIES
+  ===============
+]#
+
+type partial = tuple
+  assignments: seq[int]
+  remaining: HashSet[int]
+
+func findSubgraph*(g, h: Graph): seq[int] =
+  ##Returns an isomorphic copy of g inside of h, or @[] if none exists.
+  ##This is usually hard to do.
+  if g.V > h.V: return @[]
+  if g.E > h.E: return @[]
+  if g.E == 0: return @[0] #after this point we know h.E > 0
+
+  var states = newSeq[partial]() #a state is an assignment of the first _ vertices of g to vertices of h.
+  block:
+    var initial: partial
+    initial.assignments = @[]
+    initial.remaining = initHashSet[int]()
+    for v in 0..<h.V: initial.remaining.incl v
+    states.add initial
+  var gdegs = g.degreeSequence
+  var hdegs = h.degreeSequence
+  while states.len > 0:
+    var state = states.pop
+    
+    for v in state.remaining:
+      if hdegs[v] >= gdegs[state.assignments.len]:
+        #try to use v next
+        var useable = true
+        for i, u in state.assignments.pairs:
+          if g.contains(i, state.assignments.len) and not h.contains(state.assignments[i], v):
+            useable = false
+            break
+        if useable:
+          var stateNew: partial
+          stateNew.assignments = state.assignments
+          stateNew.assignments.add v
+          if stateNew.assignments.len == g.V:
+            return stateNew.assignments
+          stateNew.remaining = state.remaining
+          stateNew.remaining.excl v
+          states.add stateNew
+  return @[]
+
+func findInducedSubgraph*(g, h: Graph): seq[int] =
+  ##TODO
+  return @[]
+
+func findIsomorphism*(g, h: Graph): seq[int] =
+  ##Returns an isomorphism between g and h, or @[] if none exists.
+  ##This is usually hard to do.
+  if g.V != h.V: return @[]
+  if g.E != h.E: return @[]
+  var gdegs = g.degreeMultiset
+  var hdegs = h.degreeMultiset
+  for i in 0..<g.V:
+    if gDegs[i] != hDegs[i]: return @[]
+
+  var states = newSeq[partial]() #a state is an assignment of the first _ vertices of g to vertices of h.
+  block:
+    var initial: partial
+    initial.assignments = @[]
+    initial.remaining = initHashSet[int]()
+    for v in 0..<h.V: initial.remaining.incl v
+    states.add initial
+  gdegs = g.degreeSequence
+  hdegs = h.degreeSequence
+  while states.len > 0:
+    var state = states.pop
+    
+    for v in state.remaining:
+      if hdegs[v] == gdegs[state.assignments.len]:
+        #try to use v next
+        var useable = true
+        for i, u in state.assignments.pairs:
+          if g.contains(i, state.assignments.len) and not h.contains(state.assignments[i], v):
+            useable = false
+            break
+        if useable:
+          var stateNew: partial
+          stateNew.assignments = state.assignments
+          stateNew.assignments.add v
+          if stateNew.assignments.len == g.V:
+            return stateNew.assignments
+          stateNew.remaining = state.remaining
+          stateNew.remaining.excl v
+          states.add stateNew
+  return @[]
+
+func subgraphQ*(g, h: Graph): bool = findSubgraph(g, h).len > 0
+
+func inducedSubgraphQ*(g, h: Graph): bool = findSubgraph(g, h).len > 0 #TODO
+
+func isomorphicQ*(g, h: Graph): bool = findIsomorphism(g, h).len > 0
+
+proc swapVertices*(g: var Graph, v, w: int) =
+  ##Swaps the labels of vertices v and w.
+  if g.directed:
+    echo "Not implemented for directed graphs."
+    quit(-1)
+  if v == w: return
+  case g.representation:
+  of AdjMatrix:
+    #todo
+    return
+  of Sparse:
+    for u in (g.adjList[v] - g.adjList[w]):
+      if u == w: continue
+      g.adjList[u].excl v
+      g.adjList[u].incl w
+    for u in (g.adjList[w] - g.adjList[v]):
+      if u == v: continue
+      g.adjList[u].excl w
+      g.adjList[u].incl v
+    swap g.adjList[v], g.adjList[w]
+
+  of Dense:
+    #todo
+    return
+
+proc relabel*(g: var Graph) =
+  ##Relabels g so that the degree string is increasing.
+  var list = g.degreeSequence
+  var list2 = list.sorted
+  for i in 0..<g.V:
+    if list[i] != list2[i]:
+      for j in i+1..<g.V:
+        if list[j] == list2[i]:
+          swap list[i], list[j]
+          g.swapVertices(i, j)
+
+
+var chromaticTable* = initTable[string, poly[int64]]()
+# var seenGraphs* = newSeq[seq[(Graph, poly[int64])]](12)
+
+proc chromaticPolynomial*(g: Graph): poly[int64] =
+  if g.directed:
+    echo "Trying to get chrom poly of directed graph. Exiting."
+    quit(-1)
+  
+  if g.E == 0:
+    return createPoly(@[1'i64]).shift(g.V)
+  if g.E == 1: 
+    return createPoly(@[-1'i64, 1'i64]).shift(g.V - 1)
+  if g.E == 2 and g.V == 3: 
+    return createPoly(@[0'i64, 1'i64, -2'i64, 1'i64])
+  var singletons = 0
+  # for v in 0..<g.V:
+  #   if g.degree(v) == 0: inc singletons
+  if g.E == ((g.V - singletons) * (g.V - singletons - 1)) div 2:
+    var factor = createPoly(@[1'i64]).shift(1)
+    result = factor
+    for i in 1..<g.V - singletons:
+      factor = factor - 1
+      result = result * factor
+    return result.shift(singletons)
+  # var g = g
+  # if rand(1.0) > 0.01: g.relabel
+  var id = $g.V
+  for e in g.edges:
+    id = id & $e
+  if chromaticTable.hasKey(id): return chromaticTable[id]
+
+  # if g.V < seenGraphs.len:
+  #   for (h, ph) in seenGraphs[g.V]:
+  #     if g.isomorphicQ(h): 
+  #       # chromaticTable[id] = ph
+  #       return ph
+
+  if g.edgeDensity <= 1.0: #change when dense graphs are implemented:
+    for v in countdown(g.V - 1, 0):
+      for w in g.edges(v):
+        #only need the first edge we find
+        var h = g.clone
+        h.disconnect(v, w)
+        var q = g / (v, w)
+        var hp = h.chromaticPolynomial
+        var qp = q.chromaticPolynomial
+        # quit(-1)
+        result = hp - qp
+        chromaticTable[id] = result
+        # if g.V < seenGraphs.len:
+        #   seenGraphs[g.V].add (g, result)
+        return result
+  else:
+    #find first not-edge
+    for v in countdown(g.V - 1, 0):
+      for w in countdown(v - 1, 0):
+        if not g.contains(v, w):
+          #only need the first edge we find
+          var h = g.clone
+          h.connect(v, w)
+          var q = g / (v, w)
+          var hp = h.chromaticPolynomial
+          var qp = q.chromaticPolynomial
+          result = hp + qp
+          chromaticTable[id] = result
+          return result
+  
+  # echo g.V, ", ", g.E
+
+var chromaticTableMod = initTable[string, poly[zmod]]()
+proc chromaticPolynomialMod*(g: Graph, m: int64): poly[zmod] =
+  if g.directed:
+    quit(-1)
+  
+  if g.E == 0:
+    return createPoly(@[(1'i64, m)]).shift(g.V)
+  if g.E == 1:
+    var x = createPoly(@[(1'i64, m)]).shift(1)
+    if g.V == 2: return x * (x - (1'i64, m))
+    return x.shift(g.V - 2) * (x - (1'i64, m))
+  if g.E == (g.V * (g.V - 1)) div 2:
+    var factor = createPoly(@[(1'i64, m)]).shift(1)
+    result = factor
+    for i in 1..<g.V:
+      factor = factor - (1'i64, m)
+      result = result * factor
+    return result
+  var id = $g.V
+  for e in g.edges: id = id & $e
+  if chromaticTableMod.hasKey(id): return chromaticTableMod[id]
+  if g.edgeDensity <= 0.5:
+    for (v, w) in g.edges:
+      #only need the first edge we find
+      var h = g.clone
+      h.disconnect(v, w)
+      var q = g / (v, w)
+      var hp = h.chromaticPolynomialMod(m)
+      var qp = q.chromaticPolynomialMod(m)
+      # quit(-1)
+      result = hp - qp
+      chromaticTableMod[id] = result
+      return result
+  else:
+    #find first not-edge
+    for v in countdown(g.V - 1, 0):
+      # if g.degree(v) == 0: continue
+      for w in countdown(v - 1, 0):
+        if not g.contains(v, w):
+          #only need the first edge we find
+          var h = g.clone
+          h.connect(v, w)
+          var q = g / (v, w)
+          var hp = h.chromaticPolynomialMod(m)
+          var qp = q.chromaticPolynomialMod(m)
+          # quit(-1)
+          result = hp + qp
+          chromaticTableMod[id] = result
+          return result
+
+proc chromaticNumber*(g: Graph): int =
+  ##Returns the chromatic number of g.
+  ##Will leave behind chromatic polynomial memoization data.
+  ##Use sanitizeChromaticTables() if you wish to clear it.
+  var p = g.chromaticPolynomial
+  result = 1
+  while (p <~ result) == 0: inc result
+
+proc sanitizeChromaticTables* =
+  chromaticTable.clear
+  chromaticTableMod.clear
+
+
+func acyclicQ*(g: Graph): bool =
+  ##Returns whether g is acyclic (a tree).
+  if g.E <= 2: return true
+
+  var unmarked = initHashSet[int]()
+  for v in 0..<g.V:
+    unmarked.incl v
+
+  var traversed: matrix[bool]
+  traversed.initMatrix(g.V, g.V)
+
+  while unmarked.len > 0:
+    var batch = initHashSet[int]()
+    batch.incl unmarked.pop
+    while batch.len > 0:
+      var batchNew = initHashSet[int]()
+      for v in batch:
+        for w in g.edges(v):
+          if not (w in unmarked) and not traversed[v, w]: return false
+          traversed[v, w] = true
+          if not g.directed: traversed[w, v] = true
+          if w in unmarked:
+            batchNew.incl w
+            unmarked.excl w
+      batch = batchNew
+  return true
+ 
+func bipartiteQ*(g: Graph): bool =
+  ##Returns whether g is bipartite.
+  ##This is equivalent to g having no odd cycle.
+  var unmarked = initHashSet[int]()
+  for v in 0..<g.V:
+    unmarked.incl v
+  var red = newSeq[bool](g.V)
+  while unmarked.len > 0:
+    var batch = initHashSet[int]()
+    batch.incl unmarked.pop
+    while batch.len > 0:
+      var batchNew = initHashSet[int]()
+      for v in batch:
+        for w in g.edges(v):
+          if w in unmarked: batchNew.incl w
+          elif red[v] == red[w]: return false
+          red[w] = not red[v]
+          unmarked.excl w
+      batch = batchNew
+  return true
+
+func shortestPaths*(g: Graph): matrix[seq[int]] = 
+  ##Returns a sequence of vertices which form a shortest path from v to w, or an empty seq if none exists.
+  ##result[v, v] will tell you the shortest circuit starting and ending at v.
+  result.initMatrix(g.V, g.V)
+  for v in 0..<g.V:
+    for w in 0..<g.V:
+      result[v, w] = newSeq[int]()
+  for (v, w) in g.edges:
+    result[v, w] = @[v, w]
+    if not g.directed:
+      result[w, v] = @[w, v]
+  #if you uncomment this it will not find the shortest circuits v -> ... -> v.
+  # for v in 0..<g.V:
+  #   result[v, v] = @[v]
+  for k in 0..<g.V:
+    for i in 0..<g.V:
+      for j in 0..<g.V:
+        if result[i, k].len > 0 and result[k, j].len > 0 and (result[i, j].len == 0 or result[i, j].len > result[i, k].len + result[k, j].len):
+          result[i, j] = result[i, k][0..^2] & result[k, j]
+
+proc connectedQ*(g: Graph): bool =
+  ##Returns whether g is connected.
+  if g.E <= 1: return true
+  var lengths = g.shortestPaths
+  for v in 0..<g.V:
+    for w in 0..<g.V:
+      if v == w: continue
+      if lengths[v, w].len == 0: 
+        echo v, ", ", w
+        return false
+  return true
+func completeGraph*(n: int): Graph = initDenseGraph(n)
+
+proc completeBipartiteGraph*(n: int, m: int): Graph =
+  defineGraph(output, n+m):
+    for v in 0..<n:
+      for w in n..<n+m:
+        v ==> w
+  return output
+
+func cycleGraph*(n: int): Graph =
+  defineGraph(output, n):
+    for v in 0..<n:
+      v ==> (v + 1) mod n
+  return output
+
+func starGraph*(n: int): Graph = completeBipartiteGraph(1, n)
+
+func cubicGraph*(n: int): Graph =
+  ##Returns the n dimensional cube graph.
+  defineGraph(output, 1 shl n):
+    for v in 0..<(1 shl n):
+      var z = 1
+      for i in 1..n:
+        v ==> (v xor z)
+        z = z shl 1
+  return output
+
+func petersenGraph*: Graph =
+  defineGraph(output, 10):
+    cycle [0, 1, 2, 3, 4]
+    cycle [5, 7, 9, 6, 8]
+    for i in 0..4:
+      i ==> i+5
+  return output
+
+func octahedralGraph*: Graph =
+  defineGraph(output, 6):
+    cycle [0, 1, 2]
+    cycle [3, 4, 5]
+    for i in 0..2:
+      i ==> i + 3
+      i ==> ((i + 1) mod 3) + 3
+  return output
+
+#generated by Mathematica
+func dodecahedralGraph*: Graph =
+  defineGraph(output, 20):
+    0 ==> 13
+    0 ==> 14
+    0 ==> 15
+    1 ==> 4
+    1 ==> 5
+    1 ==> 12
+    2 ==> 6
+    2 ==> 13
+    2 ==> 18
+    3 ==> 7
+    3 ==> 14
+    3 ==> 19
+    4 ==> 10
+    4 ==> 18
+    5 ==> 11
+    5 ==> 19
+    6 ==> 10
+    6 ==> 15
+    7 ==> 11
+    7 ==> 15
+    8 ==> 9
+    8 ==> 13
+    8 ==> 16
+    9 ==> 14
+    9 ==> 17
+    10 ==> 11
+    12 ==> 16
+    12 ==> 17
+    16 ==> 18
+    17 ==> 19
+  return output
+
+func icosahedralGraph*: Graph =
+  defineGraph(output, 12):
+    0 ==> 2
+    0 ==> 4
+    0 ==> 5
+    0 ==> 8
+    0 ==> 9
+    1 ==> 3
+    1 ==> 6
+    1 ==> 7
+    1 ==> 10
+    1 ==> 11
+    2 ==> 6
+    2 ==> 7
+    2 ==> 8
+    2 ==> 9
+    3 ==> 4
+    3 ==> 5
+    3 ==> 10
+    3 ==> 11
+    4 ==> 5
+    4 ==> 8
+    4 ==> 10
+    5 ==> 9
+    5 ==> 11
+    6 ==> 7
+    6 ==> 8
+    6 ==> 10
+    7 ==> 9
+    7 ==> 11
+    8 ==> 10
+    9 ==> 11
+  return output
+
+func moserSpindleGraph*: Graph =
+  defineGraph(output, 7):
+    0 ==> 1
+    0 ==> 4
+    0 ==> 6
+    1 ==> 2
+    1 ==> 5
+    2 ==> 3
+    2 ==> 5
+    3 ==> 4
+    3 ==> 5
+    3 ==> 6
+    4 ==> 6
+  return output
+
+func wheelGraph*(n: int): Graph =
+  result = cycleGraph(n-1)
+  for v in 0..<(n-1):
+    result.connect(v, n-1)
+
+func pathGraph*(n: int): Graph =
+  defineGraph(output, n):
+    for i in 0..<(n-1):
+      i ==> i+1
+  return output
+
+func gridGraph*(r, c: int): Graph =
+  ##Returns the rectangular r x c grid graph.
+  defineGraph(output, r*c):
+    for i in 0..<r:
+      for j in 0..<c:
+        if i < (r-1):
+          i*c + j ==> (i+1)*c + j
+        if j < (c-1):
+          i*c + j ==> i*c + (j+1)
+  return output
+
+func rookGraph*(m, n: int): Graph = cartesianProduct(completeGraph(m), completeGraph(n)) 
+#also edge graph of completeBipartiteGraph(m, n)
+
+func torusGridGraph*(m, n: int): Graph = cartesianProduct(cycleGraph(m), cycleGraph(n))
+
+func turanGraph*(n, r: int): Graph =
+  result = initGraph(1)
+  var q = n div r
+  var s = n mod r
+  var at = 0
+  for u in 1..s:
+    #each group is size (q+1)
+    for v in at..<at + (q+1):
+      for w in 0..<at:
+        result.connect(v, w)
+      for w in at+(q+1)..<n:
+        result.connect(v, w)
+    at += q+1
+  for u in s+1..r:
+    #each group is size q
+    for v in at..<at + q:
+      for w in 0..<at:
+        result.connect(v, w)
+      for w in at+q..<n:
+        result.connect(v, w)
+    at += q
+
+func prismGraph*(n: int): Graph = pathGraph(2).cartesianProduct(cycleGraph(n))
+
+
+proc `$`*(g: Graph): string =
+  ##This attempts to name the graph.
+  ##If it can't come up with a name then it will describe its properties.
+  if g.directed:
+    #do something else
+    #TODO
+    discard 0
+  result = ""
+  var components = g.connectedComponents
+  proc recognizeComponent(h: Graph): string =
+    #recognize a connected graph
+    if h.V == 1: return "Singleton"
+    elif h.V == 2: return "Path[2]"
+    elif h.V == 3 and h.E == 2: return "Path[3]"
+    elif h.V == 3 and h.E == 3: return "Triangle"
+    elif h.V == 4 and h.E == 4 and h.isomorphicQ(cycleGraph(4)): return "Square"
+    elif h.V == 4 and h.E == 6: return "Tetrahedron (K[4])"
+    elif h.V == 5 and h.E == 5 and h.isomorphicQ(cycleGraph(5)): return "Pentagon"
+    elif h.V == 6 and h.E == 6 and h.isomorphicQ(cycleGraph(6)): return "Hexagon"
+    elif h.V == h.E and h.isomorphicQ(cycleGraph(h.V)): return "C[" & $h.V & "]"
+    elif h.V == 8 and h.E == 12 and h.isomorphicQ(cubicGraph(3)): return "Cube"
+    elif h.V == 12 and h.E == 24 and h.isomorphicQ(cubicGraph(3).lineGraph): return "Cuboctahedron"
+    elif h.V == 16 and h.E == 32 and h.isomorphicQ(cubicGraph(4)): return "Tesseract"
+    #test for other cube graphs
+    for l in 5..8:
+      if h.V == (1 shl l) and h.E == l * (1 shl (l-1)) and h.isomorphicQ(cubicGraph(l)): return $h.V & "-Cube"
+    if h.V == 4 and h.E == 3 and h.isomorphicQ(starGraph(3)): return "Claw"
+    elif h.E == h.V - 1 and h.isomorphicQ(starGraph(h.V - 1)): return "Star[" & $(h.V -  1) & "]"
+    elif h.V == 10 and h.E == 15 and h.isomorphicQ(petersenGraph()): return "Petersen"
+    elif h.V == 6 and h.E == 12 and h.isomorphicQ(octahedralGraph()): return "Octahedron"
+    elif h.E == h.V - 1 and h.isomorphicQ(pathGraph(h.V)): return "Path[" & $h.V & "]"
+    elif h.V == 20 and h.E == 30 and h.isomorphicQ(dodecahedralGraph()): return "Dodecahedron"
+    elif h.V == 30 and h.E == 60 and h.isomorphicQ(dodecahedralGraph().lineGraph): return "Icosidodecahedron"
+    elif h.V == 12 and h.E == 30 and h.isomorphicQ(icosahedralGraph()): return "Icosahedron"
+    elif h.V == 30 and h.E == 120 and h.isomorphicQ(icosahedralGraph().lineGraph): return "Icosahedral Line Graph"
+    elif h.E == 2 * (h.V - 1) and h.isomorphicQ(wheelGraph(h.V)): return "Wheel[" & $h.V & "]"
+    elif h.V == 7 and h.E == 11 and h.isomorphicQ(moserSpindleGraph()): return "Moser Spindle"
+    #test for grid graph
+    for r in 1..h.V:
+      if r*r > h.V: break
+      if h.V mod r == 0:
+        var c = h.V div r
+        if h.E == 2*r*c - r - c and h.isomorphicQ(gridGraph(r, c)):
+          return "Grid[" & $r & ", " & $c & "]"
+    #test for complete graph
+    if 2 * h.E == h.V * (h.V - 1): return "K[" & $h.V & "]"
+    #test for complete bipartite graph
+    if h.V * h.V - 4 * h.E >= 0:
+      var u = h.V * h.V - 4 * h.E
+      var urt = isqrt(u)
+      if urt * urt == u and (h.V + urt) mod 2 == 0:
+        var r = (h.V - urt) div 2
+        var c = (h.V + urt) div 2
+        if r >= 1 and c >= 1 and h.isomorphicQ(completeBipartiteGraph(r.int, c.int)):
+          return "CompleteBipartite[" & $r & ", " & $c & "]"
+    #test for prism
+    if h.V * 3 == h.E * 2 and h.isomorphicQ(prismGraph(h.V div 2)): return "Prism[" & $(h.V div 2) & "]"
+    #test for torus grid graph
+    if h.E == 2 * h.V:
+      for p in 1..h.V:
+        if p*p > h.V: break
+        if h.V mod p == 0:
+          var q = h.V div p
+          if h.isomorphicQ(torusGridGraph(p, q)):
+            return $p & " x " & $q & " Torus Grid Graph / " & $p & " - " & $q & " Duoprism"
+    #do not test for turan graph. sorry :[
+    #test for rook graph
+    for r in 1..h.V:
+      if r*r > h.V: break
+      if h.V mod r == 0:
+        var c = h.V div r
+        if h.E == ((r*c*(r+c)) div 2) - r*c and h.isomorphicQ(rookGraph(r, c)):
+          return "Rook[" & $r & ", " & $c & "]"
+    #hasn't been recognized.. just describe it
+    result = "V = " & $h.V & ", E = " & $h.E & ", "
+    if h.bipartiteQ: result &= "Bipartite, "
+    if h.acyclicQ: result &= "Acyclic, "
+    result &= "with degree multiset " & $h.degreeMultiset & "."
+
+  for i, h in components.pairs:
+    var part = ""
+    if components.len > 1: part = "Component " & $i & ": "
+    part &= recognizeComponent(h)
+    if i != components.high: part &= "\n"
+    result &= part
+
+proc characteristicPoly*(g: Graph): poly[int64] =
+  ##Returns the characteristic polynomial of a graph.
+  var m = g.adjacencyMatrix
+  castCopyMat(m, m64, int64)
+  return m64.characteristicPoly
+
+type candidateMIS = tuple
+  chosen: seq[int]
+  remaining: HashSet[int]
+proc `<`(x, y: candidateMIS): bool = x.chosen.len > y.chosen.len
+proc maximalIndependentSet*(g: Graph): seq[int] =
+  var maxSoFar = 0
+  var selections: seq[int]
+  var stk = initHeapQueue[candidateMIS]()
+  block:
+    var initial = initHashSet[int]()
+    for v in 0..<g.V: initial.incl v
+    stk.push (@[], initial)
+  while stk.len > 0:
+    var current = stk.pop
+    # echo current
+    if maxSoFar >= current.chosen.len + current.remaining.len: continue
+    if maxSoFar < current.chosen.len:
+      maxSoFar = current.chosen.len
+      selections = current.chosen
+    if current.remaining.len == 0: continue
+    var remainingNew = current.remaining
+    block:
+      var v = remainingNew.pop
+      stk.push (current.chosen, remainingNew)
+      for w in g.edges(v):
+        remainingNew.excl w
+      var chosenNew = current.chosen
+      chosenNew.add v
+      stk.push (chosenNew, remainingNew)
+  return selections
+
diff --git a/nt.nim b/nt.nim
new file mode 100644
index 0000000..98b379d
--- /dev/null
+++ b/nt.nim
@@ -0,0 +1,152 @@
+import std/[sequtils]
+
+const MOD* = 1000000007
+type ll* = int64
+
+proc reduce*(x: ll): ll =
+  var r = x mod MOD
+  if r < 0: r += MOD
+  r
+
+proc bexp*(b: ll, p: ll): ll =
+  var
+    curr = reduce(b)
+    power = p
+    res: ll = 1
+  while power > 0:
+    if (power and 1) == 1:
+      res = (res * curr) mod MOD
+    curr = (curr * curr) mod MOD
+    power = power shr 1
+  res
+
+proc inv*(x: ll): ll =
+  bexp(x, MOD - 2)
+
+type
+  Combinatorics*[N: static[int]] = object
+    factorialCache: array[N+1, ll]
+    inverseFactorialCache: array[N+1, ll]
+
+proc initCombinatorics*[N: static[int]](): Combinatorics[N] =
+  var c: Combinatorics[N]
+  c.factorialCache[0] = 1
+
+  for i in 1..N:
+    c.factorialCache[i] = (c.factorialCache[i-1] * i) mod MOD
+
+  c.inverseFactorialCache[N] = inv(c.factorialCache[N])
+
+  for i in countdown(N-1, 0):
+    c.inverseFactorialCache[i] =
+      ((i+1).ll * c.inverseFactorialCache[i+1]) mod MOD
+
+  c
+
+proc fact*[N](c: Combinatorics[N], x: int): ll =
+  c.factorialCache[x]
+
+proc ifact*[N](c: Combinatorics[N], x: int): ll =
+  c.inverseFactorialCache[x]
+
+proc permute*[N](c: Combinatorics[N], n, k: int): ll =
+  if k > n: return 0
+  (c.fact(n) * c.ifact(n-k)) mod MOD
+
+proc choose*[N](c: Combinatorics[N], n, k: int): ll =
+  if k > n: return 0
+  (((c.fact(n) * c.ifact(k)) mod MOD) * c.ifact(n-k)) mod MOD
+
+proc gcd*(a, b: ll): ll =
+  if b == 0: a else: gcd(b, a mod b)
+
+proc lcm*(a, b: ll): ll =
+  (a div gcd(a,b)) * b
+
+type
+  EGCDResult* = object
+    x*: ll
+    y*: ll
+    d*: ll
+
+proc extendedGcd*(a, b: ll): EGCDResult =
+  if b == 0:
+    return EGCDResult(x: 1, y: 0, d: a)
+
+  let next = extendedGcd(b, a mod b)
+
+  EGCDResult(
+    x: next.y,
+    y: next.x - (a div b) * next.y,
+    d: next.d
+  )
+
+type
+  FactorSieve*[N: static[int]] = object
+    sieve: array[N+1, int]
+
+proc initFactorSieve*[N: static[int]](): FactorSieve[N] =
+  var s: FactorSieve[N]
+  s.sieve[1] = 1
+
+  for i in 2..N:
+    if s.sieve[i] == 0:
+      for j in countup(i, N, i):
+        s.sieve[j] = i
+
+  s
+
+proc factor*[N](s: FactorSieve[N], x0: int): seq[(int,int)] =
+  var x = x0
+  var res: seq[(int,int)] = @[]
+
+  while x != 1:
+    let curr = s.sieve[x]
+    var cnt = 0
+
+    while s.sieve[x] == curr:
+      inc cnt
+      x = x div curr
+
+    res.add((curr,cnt))
+
+  res
+
+type
+  BooleanSieve*[N: static[int]] = object
+    isComposite: array[N+1, bool]
+
+proc initBooleanSieve*[N: static[int]](): BooleanSieve[N] =
+  var s: BooleanSieve[N]
+  s.isComposite[1] = true
+
+  for i in 2..N:
+    if not s.isComposite[i]:
+      for j in countup(2*i, N, i):
+        s.isComposite[j] = true
+
+  s
+
+proc isPrime*[N](s: BooleanSieve[N], x: int): bool =
+  not s.isComposite[x]
+
+# ---------------- MAIN ----------------
+
+proc main() =
+  let comb = initCombinatorics[100000]()
+  echo "10 choose 3 = ", comb.choose(10,3)
+
+  echo "gcd(48,18) = ", gcd(48,18)
+  echo "lcm(48,18) = ", lcm(48,18)
+
+  let eg = extendedGcd(48,18)
+  echo "extended gcd: x=", eg.x, " y=", eg.y, " d=", eg.d
+
+  let sieve = initFactorSieve[100000]()
+  echo "Factors of 84: ", sieve.factor(84)
+
+  let bs = initBooleanSieve[100000]()
+  echo "Is 97 prime? ", bs.isPrime(97)
+
+when isMainModule:
+  main()
\ No newline at end of file
diff --git a/segmented.nim b/segmented.nim
new file mode 100644
index 0000000..e56ccd7
--- /dev/null
+++ b/segmented.nim
@@ -0,0 +1,45 @@
+import math, sequtils
+proc simple(limit: int): seq[int] =
+  var isPrime = newSeq[bool](limit + 1)
+  for i in 0..limit:
+    isPrime[i] = true
+  isPrime[0] = false
+  isPrime[1] = false
+
+  for i in 2..int(sqrt(float(limit))):
+    if isPrime[i]:
+      var j = i * i
+      while j <= limit:
+        isPrime[j] = false
+        j += i
+
+  result = @[]
+  for i, prime in isPrime:
+    if prime:
+      result.add(i)
+
+proc segmented(L, R: int): seq[int] =
+  let limit = int(sqrt(float(R)))
+  let primes = simple(limit)
+  var isPrime = newSeq[bool](R - L + 1)
+  for i in 0..<isPrime.len:
+    isPrime[i] = true
+
+  for p in primes:
+    var start = (L div p) * p
+    if start < L: start += p
+    if start == p: start += p  
+
+    var j = start
+    while j <= R:
+      isPrime[j - L] = false
+      j += p
+
+  result = @[]
+  for i in 0..<isPrime.len:
+    if isPrime[i] and (i + L >= 2):
+      result.add(i + L)
+
+let L = 10
+let R = 500000000
+echo segmented(L, R)
diff --git a/segtree.nim b/segtree.nim
new file mode 100644
index 0000000..8c59613
--- /dev/null
+++ b/segtree.nim
@@ -0,0 +1,98 @@
+type SegTree*[T, U] = ref object
+  len*: int
+  arr*: seq[T]
+  shift*: seq[U]
+
+proc SegTree*[T, U](f: var SegTree[T, U], len: int): void =
+  new(f)
+  f.len = len
+  f.arr.newSeq(4*len)
+  f.shift.newSeq(4*len)
+
+proc build[T, U](f: var SegTree[T, U], v, tl, tr: int, a: seq[T]): void =
+    if tl==tr:
+      f.arr[v] = a[tl]
+    else:
+      var tm: int = (tl + tr) div 2
+      f.build(2*v, tl, tm, a)
+      f.build(2*v+1, tm+1, tr, a)
+      f.arr[v] = f.arr[2*v] + f.arr[2*v+1]    
+
+proc build*[T, U](f: var SegTree[T, U], a: seq[T]): void =
+  f.build(1, 0, f.len - 1, a)
+  
+proc modifyRange[T, U](tl, tr: int, sum: T, x: U, additive: bool): T =
+  if additive: return sum @ ((tr-tl+1) * x)
+  else: return sum @ x
+
+proc sum[T, U, R](f: var SegTree[T, U], v, tl, tr, l, r: int, phi: proc (x: T): R, additive: bool): R =
+  if l>r:
+    return
+  if f.shift[v] - f.shift[v] != f.shift[v]: #nonzero
+    f.arr[v] = modifyRange(tl, tr, f.arr[v], f.shift[v], additive)
+    if tl!=tr:
+      f.shift[2*v] = f.shift[2*v] + f.shift[v]
+      f.shift[2*v+1] = f.shift[2*v+1] + f.shift[v]
+    f.shift[v] = f.shift[v] - f.shift[v]
+  if l<=tl and tr<=r: return f.arr[v].phi
+  elif tl!=tr:
+    var tm: int = (tl + tr) div 2
+    return f.sum(2*v, tl, tm, l, min(r, tm), phi, additive) + f.sum(2*v+1, tm+1, tr, max(l, tm+1), r, phi, additive)
+
+proc sum*[T, U, R](f: var SegTree[T, U], i, j: int, phi: proc (x: T): R, additive: bool): R =
+  return f.sum(1, 0, f.len - 1, i, j, phi, additive)
+
+proc showDetails*[T, U](f: var SegTree[T, U], v: int = 1, tl: int = 0, tr: int = f.len - 1, prefix: string = ""): void =
+  echo prefix, "(", tl, ", ", tr, "): ", f.arr[v], ", ", f.shift[v]
+  if tl != tr:
+    var tm: int = (tl + tr) div 2
+    f.showDetails(2*v, tl, tm, prefix & "  ") 
+    f.showDetails(2*v+1, tm+1, tr, prefix & "  ")
+
+
+proc applyTo[T, U](f: var SegTree[T, U], v, tl, tr, l, r: int, x: U, additive: bool): void =
+  if f.shift[v] - f.shift[v] != f.shift[v]: #nonzero
+    f.arr[v] = modifyRange(tl, tr, f.arr[v], f.shift[v], additive)
+    if tl != tr:
+      f.shift[2*v] = f.shift[2*v] + f.shift[v]
+      f.shift[2*v+1] = f.shift[2*v+1] + f.shift[v]
+    f.shift[v] = f.shift[v] - f.shift[v]
+  if tl>tr or tl>r or tr<l: return
+  if l<=tl and tr<=r:
+    f.arr[v] = modifyRange(tl, tr, f.arr[v], x, additive)
+    if tl != tr:
+      f.shift[2*v] = f.shift[2*v] + x
+      f.shift[2*v+1] = f.shift[2*v+1] + x
+  elif tl!=tr:
+    var tm: int = (tl + tr) div 2
+    f.applyTo(2*v, tl, tm, l, min(r, tm), x, additive)
+    f.applyTo(2*v+1, tm+1, tr, max(l, tm+1), r, x, additive)
+    f.arr[v] = f.arr[2*v] + f.arr[2*v+1]
+
+proc addTo*[T, U](f: var SegTree[T, U], i, j: int, x: U): void =
+
+  f.applyTo(1, 0, f.len - 1, i, j, x, true)
+
+proc multTo*[T, U](f: var SegTree[T, U], i, j: int, x: U): void =
+
+  f.applyTo(1, 0, f.len - 1, i, j, x, false)
+
+proc toBaseArray*[T, U, R](f: var SegTree[T, U], phi: proc (x: T): R, additive: bool): seq[R] =
+  result.newSeq(f.len)
+  for i in 0..<f.len:
+    result[i] = f.sum(i, i, phi, additive)
+
+when isMainModule:
+  proc `@`(x, y: int): int = x + y
+  var seg: SegTree[int, int]
+  seg.SegTree(10)
+  var a: seq[int] =  @[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+  seg.build(a)
+  seg.showDetails
+  proc phi(x: int): int = x 
+  seg.addTo(1, 8, 4)
+  seg.showDetails
+  seg.addTo(3, 9, 4)
+  seg.showDetails
+  for i in 0..9:
+    echo i, ": ", seg.sum(i, i, phi, true)
\ No newline at end of file
diff --git a/sieve.nim b/sieve.nim
new file mode 100644
index 0000000..213a6f3
--- /dev/null
+++ b/sieve.nim
@@ -0,0 +1,250 @@
+import math
+
+type
+  Factor* = object
+    prime*: uint64
+    power*: uint64
+
+  Factors* = object
+    numPrimes*: uint64
+    factors*: seq[Factor]
+
+  U64Array* = object
+    elems*: seq[uint64]
+
+# =========================
+# Utilities
+# =========================
+
+proc factorAppend*(f: var Factors, p, k: uint64) =
+  f.factors.add Factor(prime: p, power: k)
+  inc f.numPrimes
+
+# =========================
+# Math helpers
+# =========================
+
+proc u64Pow*(a, e: uint64): uint64 =
+  var base = a
+  var exp = e
+  var r: uint64 = 1
+  while exp > 0:
+    if (exp and 1) == 1:
+      r *= base
+    base *= base
+    exp = exp shr 1
+  r
+
+proc u64NthRoot*(x, n: uint64): uint64 =
+  uint64(pow(float64(x), 1.0 / float64(n)))
+
+proc u64Binom*(n, k: uint64): uint64 =
+  if k > n: return 0
+  var r: uint64 = 1
+  for i in 1..k:
+    r = r * (n - k + i) div i
+  r
+
+proc i64Gcd*(a, b: int64): int64 =
+  var x = a
+  var y = b
+  while y != 0:
+    let t = x mod y
+    x = y
+    y = t
+  x
+
+proc i64Lcm*(a, b: int64): int64 =
+  (a div i64Gcd(a, b)) * b
+
+proc i64Egcd*(a, b: int64, x, y: var int64): int64 =
+  if b == 0:
+    x = 1
+    y = 0
+    return a
+
+  var x1, y1: int64
+  let g = i64Egcd(b, a mod b, x1, y1)
+  x = y1
+  y = x1 - y1 * (a div b)
+  g
+
+# =========================
+# Prime estimate helpers
+# =========================
+
+proc maxPrimeDivs*(max: uint64): uint64 =
+  let primorials = [
+    6'u64,30,210,2310,30030,
+    510510,9699690,223092870
+  ]
+
+  for i, p in primorials:
+    if max < p:
+      return uint64(i+1)
+
+  15
+
+proc maxPrimesLe*(max: uint64): uint64 =
+  uint64(1.25506 * float64(max) / ln(float64(max)))
+
+# =========================
+# Omega sieve
+# =========================
+
+proc sieveOmega*(max: uint64): seq[uint8] =
+  result = newSeq[uint8](max+1)
+
+  for n in 2..max:
+    if result[n] != 0: continue
+    var m = n
+    while m <= max:
+      inc result[m]
+      m += n
+
+# =========================
+# Largest prime factor sieve
+# =========================
+
+proc sieveLargestFactors*(max: uint64): seq[uint64] =
+  result = newSeq[uint64](max+1)
+
+  for p in 2..max:
+    if result[p] != 0: continue
+    var m = p
+    while m <= max:
+      result[m] = p
+      m += p
+
+# =========================
+# Smallest prime factor sieve
+# =========================
+
+proc sieveSmallestFactors*(max: uint64): seq[uint32] =
+  result = newSeq[uint32](max+1)
+
+  let r = uint64(sqrt(float64(max)))
+
+  for p in 2..r:
+    if result[p] != 0: continue
+
+    var m = p*p
+    while m <= max:
+      if result[m] == 0:
+        result[m] = uint32(p)
+      m += p
+
+  for i in 2..max:
+    if result[i] == 0:
+      result[i] = uint32(i)
+
+# =========================
+# Factorization
+# =========================
+
+proc fillFactorsFromSmallest*(res: var Factors, n: uint64, spf: seq[uint32]) =
+  var x = n
+  res.numPrimes = 0
+  res.factors.setLen(0)
+
+  while x > 1:
+    let p = uint64(spf[x])
+    var k: uint64 = 0
+
+    while spf[x] == uint32(p):
+      x = x div p
+      inc k
+
+    factorAppend(res, p, k)
+# =========================
+# Totient sieve
+# =========================
+
+proc sievePhi*(max: uint64): seq[uint64] =
+  result = newSeq[uint64](max+1)
+
+  for i in 1..max:
+    result[i] = i
+
+  for p in 2..max:
+    if result[p] != p: continue
+
+    var m = p
+    while m <= max:
+      result[m] -= result[m] div p
+      m += p
+
+# =========================
+# Sigma_0 (divisor count)
+# =========================
+
+proc sigma0*(max: uint64): seq[uint64] =
+  result = newSeq[uint64](max+1)
+
+  for i in 1..max:
+    var j = i
+    while j <= max:
+      inc result[j]
+      j += i
+
+# =========================
+# Mobius sieve
+# =========================
+
+proc mobius*(max: uint64): seq[int8] =
+  result = newSeq[int8](max+1)
+
+  for i in 1..max:
+    result[i] = 1
+
+  for p in 2..max:
+    if result[p] == 1:
+      var m = p
+      while m <= max:
+        result[m] = -result[m]
+        m += p
+
+      let p2 = p*p
+      m = p2
+      while m <= max:
+        result[m] = 0
+        m += p2
+
+# =========================
+# Prime sieve
+# =========================
+
+proc sievePrimes*(max: uint64): seq[uint64] =
+  var isComp = newSeq[bool](max+1)
+
+  for p in 2..max:
+    if isComp[p]: continue
+
+    result.add p
+
+    if p*p <= max:
+      var m = p*p
+      while m <= max:
+        isComp[m] = true
+        m += p
+
+# =========================
+# Main
+# =========================
+
+when isMainModule:
+
+  let max = 10_000_000'u64
+
+  let primes = sievePrimes(max)
+
+  echo "Primes up to ", max, ":"
+  for p in primes:
+    stdout.write($p & " ")
+  echo ""
+
+  let phi = sievePhi(max)
+  echo "phi(36) = ", phi[36]
+
+  let mu = mobius(max)
+  echo "mu(30) = ", mu[30]