diff --git a/gallery/example3.pdf b/gallery/example3.pdf
index 78f7120..ab2412c 100644
Binary files a/gallery/example3.pdf and b/gallery/example3.pdf differ
diff --git a/gallery/example4.pdf b/gallery/example4.pdf
new file mode 100644
index 0000000..12ac42b
Binary files /dev/null and b/gallery/example4.pdf differ
diff --git a/gallery/example4.typ b/gallery/example4.typ
new file mode 100644
index 0000000..a671a61
--- /dev/null
+++ b/gallery/example4.typ
@@ -0,0 +1,34 @@
+#import "/src/lib.typ": *
+
+#let p = poly.poly(0,1,2,3,4)
+#let q = poly.poly(5,6,7)
+
+#let a = poly.add(p, q)
+#let b = poly.sub(p, q)
+
+$
+  p = #poly.display(p) quad q = #poly.display(q)\
+  p + q = #poly.display(a)\
+  p - q = #poly.display(b)\
+$
+
+#let P = poly.poly(7, -3, -3, 1, 1)
+#let Q = poly.poly(-2, 1, 1)
+
+$
+  P = #poly.display(P) quad Q = #poly.display(Q)\
+  P / Q = #poly.display(P) / #poly.display(Q)
+$
+
+#let div = poly.div(P, Q)
+#align(center, grid(
+  columns: 2,
+  align: left,
+  column-gutter: 2em,
+  row-gutter: 1em,
+  [dividend: #poly.display(div.dividend)],
+  grid.cell(rowspan: 4, (div.display)()),
+  [divisor: #poly.display(div.divisor)],
+  [quotient: #poly.display(div.quotient)],
+  [rest: #poly.display(div.rest)],
+))
diff --git a/src/lib.typ b/src/lib.typ
index f90a91b..5b20562 100644
--- a/src/lib.typ
+++ b/src/lib.typ
@@ -1,4 +1,5 @@
 #import "mat.typ"
+#import "poly.typ"
 #import "vec.typ"
 
 #import "gauss.typ"
\ No newline at end of file
diff --git a/src/poly.typ b/src/poly.typ
new file mode 100644
index 0000000..21c23d3
--- /dev/null
+++ b/src/poly.typ
@@ -0,0 +1,233 @@
+#let is-poly(poly) = {
+  if type(poly) != dictionary {return false}
+  if poly.at("type", default: none) != "polynomial" {return false}
+  return true
+}
+
+#let _check(poly) = {
+  if not is-poly(poly) {
+    panic("Argument is not a polynomial")
+  }
+}
+
+#let poly(..args, trim-zeros: true) = {
+  let coefs = args.pos()
+
+  if trim-zeros {
+    let last-i = coefs.rev().position(c => c != 0)
+    if last-i == none {
+      last-i = -1
+    }
+
+    coefs = coefs.slice(0, calc.min(coefs.len(), coefs.len() - last-i))
+  }
+
+  if coefs.len() == 0 {
+    coefs.push(0)
+  }
+
+  return (
+    type: "polynomial",
+    coefs: coefs,
+    deg: coefs.len() - 1
+  )
+}
+
+#let copy(original) = {
+  _check(original)
+  return poly(..original.coefs)
+}
+
+#let display(poly) = {
+  _check(poly)
+
+  let fractions = (2, 3, 4, 5, 6)
+  let beautify-frac(value) = {
+    if calc.round(value) == value {
+      return str(value)
+    }
+    for den in fractions {
+      for num in range(1, 2 * den) {
+        if value == num / den {
+          return str(num) + "/" + str(den)
+        }
+      }
+    }
+    return str(value)
+  }
+
+  let res = ""
+  for (i, coef) in poly.coefs.rev().enumerate() {
+    if coef == 0 and poly.deg > 0 {
+      continue
+    }
+    if i != 0 {
+      res += " "
+      if coef < 0 { res += "-" }
+      else { res += "+" }
+    } else if coef < 0 {
+      res += "-"
+    }
+    
+    if calc.abs(coef) != 1 or poly.deg - i == 0 {
+      res += str(beautify-frac(calc.abs(coef)))
+    }
+    if poly.deg - i > 0 {
+      res += "x"
+    }
+    if poly.deg - i > 1 {
+      res += "^" + str(poly.deg - i)
+    }
+  }
+  return eval("$" + res + "$")
+}
+
+#let pad(poly1, deg) = {
+  let coefs = poly1.coefs
+  coefs += (0,) * calc.max(0, (deg - poly1.deg))
+  return poly(..coefs, trim-zeros: false)
+}
+
+#let shift(poly1, deg) = {
+  let coefs = (0,) * calc.max(0, (deg - poly1.deg))
+  coefs += poly1.coefs
+  return poly(..coefs)
+}
+
+#let add(poly1, poly2) = {
+  let deg = calc.max(poly1.deg, poly2.deg)
+  let poly1 = pad(poly1, deg)
+  let poly2 = pad(poly2, deg)
+
+  let coefs = poly1.coefs.zip(poly2.coefs).map(p => p.sum())
+
+  return poly(..coefs)
+}
+
+#let sub(poly1, poly2) = {
+  let deg = calc.max(poly1.deg, poly2.deg)
+  let poly1 = pad(poly1, deg)
+  let poly2 = pad(poly2, deg)
+
+  let coefs = poly1.coefs.zip(poly2.coefs).map(p => p.at(0) - p.at(1))
+
+  return poly(..coefs)
+}
+
+#let mul(poly1, f) = {
+  let coefs = poly1.coefs.map(c => c * f)
+  return poly(..coefs)
+}
+
+#let to-cells(poly) = {
+  let cells = ()
+  let first = true
+  for (i, coef) in poly.coefs.rev().enumerate() {
+    if coef == 0 and poly.deg > 0 {
+      cells.push("")
+      continue
+    }
+    let cell = ""
+    if i != 0 {
+      if coef < 0 { cell += "-" }
+      else if not first { cell += "+" }
+    } else if coef < 0 {
+      cell += "-"
+    }
+    
+    if calc.abs(coef) != 1 or poly.deg - i == 0 {
+      cell += str(calc.abs(coef))
+    }
+    if poly.deg - i > 0 {
+      cell += "x"
+    }
+    if poly.deg - i > 1 {
+      cell += "^" + str(poly.deg - i)
+    }
+    cells.push(eval("$" + cell + "$"))
+    first = false
+  }
+  if poly.coefs.position(c => c != 0) == none {
+    cells.last() = "0"
+  }
+  return cells
+}
+
+#let div(
+  poly1,
+  poly2,
+  max-steps: 10
+) = {
+  let P = poly1
+  let Q = poly2
+  let p = P
+  let rows = ()
+  let steps = ()
+  let D = poly()
+  let count = 0
+
+  while p.deg >= Q.deg and p.coefs.position(c => c != 0) != none {
+    let f = p.coefs.at(p.deg) / Q.coefs.at(Q.deg)
+    let sub = shift(mul(Q, -f), p.deg)
+    let rem = add(p, sub)
+    let d = shift(poly(f), p.deg - Q.deg)
+    steps.push((
+      p: p,
+      sub: sub,
+      rem: rem,
+      f: f,
+      d: d
+    ))
+    
+    D = add(D, d)
+    p = rem
+    if count >= max-steps {
+      panic("Exceeded maximum number of steps")
+      break
+    }
+    count += 1
+  }
+
+  let make-table(
+    vpad: 5pt,
+    hpad: 5pt,
+    stroke: black + 1pt,
+    debug-grid: false
+  ) = {
+    let cells = ()
+    cells += to-cells(P)
+    cells.push(table.vline(stroke: stroke))
+    cells.push(display(Q))
+    cells.push(table.hline(start: P.deg + 1, stroke: stroke))
+
+    for (i, step) in steps.enumerate() {
+      cells += to-cells(pad(step.sub, P.deg))
+      if i == 0 {
+        cells.push(display(D))
+      } else {
+        cells.push("")
+      }
+      cells.push(table.hline(start: P.deg - step.sub.deg, end: P.deg - step.sub.deg + Q.deg + 1, stroke: stroke))
+      cells += to-cells(pad(step.rem, P.deg))
+      cells.push("")
+    }
+
+    table(
+      columns: P.deg + 2,
+      stroke: if debug-grid {gray + 1pt} else {none},
+      fill: none,
+      align: (x, _) => if x == P.deg + 1 {left} else {right},
+      inset: (left: hpad, right: hpad, top: vpad, bottom: vpad),
+      ..cells
+    )
+  }
+
+  return (
+    dividend: P,
+    divisor: Q,
+    quotient: D,
+    rest: p,
+    steps: steps,
+    display: make-table
+  )
+}
\ No newline at end of file