source("boot.R")
# ── Laurent polynomials in Z[t, t^{-1}] (numeric coefficients) ──

#' @export
laurent <- function(terms = NULL) {
  if (is.null(terms)) return(structure(integer(0), class = "laurent"))
  terms <- terms[terms != 0]
  if (length(terms) == 0) return(structure(integer(0), class = "laurent"))
  structure(terms, class = "laurent")
}

#' @export
laurent_add <- function(a, b) {
  if (length(a) == 0) return(b)
  if (length(b) == 0) return(a)
  e <- union(names(a), names(b))
  r <- numeric(length(e)); names(r) <- e
  for (x in names(a)) r[[x]] <- r[[x]] + a[[x]]
  for (x in names(b)) r[[x]] <- r[[x]] + b[[x]]
  r <- r[r != 0]
  if (length(r) == 0) return(laurent())
  laurent(r)
}

#' @export
laurent_sub <- function(a, b = NULL) {
  if (is.null(a) && is.null(b)) return(laurent())
  if (is.null(b)) { r <- -unclass(a); class(r) <- "laurent"; return(r) }
  if (length(a) == 0) { r <- -unclass(b); class(r) <- "laurent"; return(r) }
  if (length(b) == 0) return(a)
  e <- union(names(a), names(b))
  r <- numeric(length(e)); names(r) <- e
  for (x in names(a)) r[[x]] <- r[[x]] + a[[x]]
  for (x in names(b)) r[[x]] <- r[[x]] - b[[x]]
  r <- r[r != 0]
  if (length(r) == 0) return(laurent())
  laurent(r)
}

#' @export
laurent_mul <- function(a, b) {
  if (length(a) == 0 || length(b) == 0) return(laurent())
  r <- numeric(0)
  for (x in names(a)) for (y in names(b)) {
    e <- as.integer(x) + as.integer(y)
    en <- as.character(e)
    r[[en]] <- (if (en %in% names(r)) r[[en]] else 0) + a[[x]] * b[[y]]
  }
  r <- r[r != 0]
  if (length(r) == 0) return(laurent())
  laurent(r)
}

#' @export
laurent_equal <- function(a, b) {
  if (length(a) != length(b)) return(FALSE)
  if (length(a) == 0) return(TRUE)
  a <- a[order(as.integer(names(a)))]
  b <- b[order(as.integer(names(b)))]
  identical(a, b)
}

#' @export
t <- laurent(c("1" = 1))

#' @export
neg_t <- laurent(c("1" = -1))

#' @export
one <- laurent(c("0" = 1))

#' @export
zero <- laurent(c())

#' @export
one_minus_t <- laurent_add(one, laurent_sub(NULL, t))

# ── Braid words ───────────────────────────────────────────────

#' @export
braid_word <- function(word = integer(0)) structure(word, class = "braid_word")

#' @export
braid_concat <- function(a, b) braid_word(c(unclass(a), unclass(b)))

#' @export
braid_inverse <- function(x) braid_word(rev(-unclass(x)))

#' @export
braid_reduce <- function(w) {
  w <- unclass(w)
  stack <- integer(0)
  for (s in w) {
    if (length(stack) > 0 && stack[length(stack)] == -s) {
      stack <- stack[-length(stack)]
    } else {
      stack <- c(stack, s)
    }
  }
  braid_word(stack)
}

#' @export
sigma <- function(i, positive = TRUE) {
  if (i < 1 || i > 3) stop("B4 has generators 1,2,3 only")
  braid_word(if (positive) i else -i)
}

# ── Burau representation ──────────────────────────────────────

#' @export
is_identity_matrix <- function(m) {
  for (r in 1:4) for (c in 1:4) {
    e <- if (r == c) one else zero
    if (!laurent_equal(m[[r, c]], e)) return(FALSE)
  }
  TRUE
}

burau_sigma <- function(i) {
  m <- matrix(rep(list(zero), 16), nrow = 4, ncol = 4)
  for (k in 1:4) m[[k, k]] <- one
  if (i == 1) {
    m[[1, 1]] <- one_minus_t; m[[1, 2]] <- t; m[[2, 1]] <- one; m[[2, 2]] <- zero
  } else if (i == 2) {
    m[[2, 2]] <- one_minus_t; m[[2, 3]] <- t; m[[3, 2]] <- one; m[[3, 3]] <- zero
  } else if (i == 3) {
    m[[3, 3]] <- one_minus_t; m[[3, 4]] <- t; m[[4, 3]] <- one; m[[4, 4]] <- zero
  }
  m
}

burau_inverse <- function(i) {
  m <- matrix(rep(list(zero), 16), nrow = 4, ncol = 4)
  for (k in 1:4) m[[k, k]] <- one
  ti <- laurent(c("-1" = 1))
  one_minus_ti <- laurent_add(one, laurent_sub(NULL, ti))
  if (i == 1) {
    m[[1, 1]] <- zero; m[[1, 2]] <- one; m[[2, 1]] <- ti; m[[2, 2]] <- one_minus_ti
  } else if (i == 2) {
    m[[2, 2]] <- zero; m[[2, 3]] <- one; m[[3, 2]] <- ti; m[[3, 3]] <- one_minus_ti
  } else if (i == 3) {
    m[[3, 3]] <- zero; m[[3, 4]] <- one; m[[4, 3]] <- ti; m[[4, 4]] <- one_minus_ti
  }
  m
}

#' @export
mat_mul <- function(A, B) {
  C <- matrix(rep(list(zero), 16), nrow = 4, ncol = 4)
  for (i in 1:4) for (j in 1:4) {
    s <- zero
    for (k in 1:4) {
      if (length(A[[i, k]]) > 0 && length(B[[k, j]]) > 0) {
        s <- laurent_add(s, laurent_mul(A[[i, k]], B[[k, j]]))
      }
    }
    C[[i, j]] <- s
  }
  C
}

#' @export
burau_matrix <- function(w) {
  w <- unclass(w)
  r <- matrix(rep(list(zero), 16), nrow = 4, ncol = 4)
  for (k in 1:4) r[[k, k]] <- one
  for (s in w) {
    g <- abs(s)
    r <- mat_mul(r, if (s > 0) burau_sigma(g) else burau_inverse(g))
  }
  r
}

#' @export
burau_matrix_stepwise <- function(w) {
  w <- unclass(w)
  r <- matrix(rep(list(zero), 16), nrow = 4, ncol = 4)
  for (k in 1:4) r[[k, k]] <- one
  monomials <- list()
  for (s in w) {
    g <- abs(s)
    before <- r[[1, 3]]
    r <- mat_mul(r, if (s > 0) burau_sigma(g) else burau_inverse(g))
    after <- r[[1, 3]]
    for (e in unique(c(names(before), names(after)))) {
      cb <- if (e %in% names(before)) before[[e]] else 0
      ca <- if (e %in% names(after)) after[[e]] else 0
      d <- ca - cb
      if (d == 0) next
      d_int <- abs(d)
      sgn <- if (d > 0) 1L else -1L
      for (n in seq_len(d_int)) {
        monomials <- c(monomials, list(list(epsilon = sgn, exponent = as.integer(e))))
      }
    }
  }
  list(matrix = r, monomials = monomials)
}

# ── Display ───────────────────────────────────────────────────

#' @export
laurent_str <- function(x) {
  if (length(x) == 0) return("0")
  nms <- names(x)
  if (is.null(nms) || any(is.na(nms))) return("0")
  exps <- as.integer(nms)
  parts <- character(length(x))
  for (i in seq_along(x)) {
    co <- x[i]
    if (is.na(co) || co == 0) next
    e <- exps[i]
    if (is.na(e)) next
    if (e == 0) {
      parts[i] <- as.character(co)
    } else if (e == 1) {
      if (co == 1) parts[i] <- "t"
      else if (co == -1) parts[i] <- "-t"
      else parts[i] <- paste0(co, "t")
    } else {
      if (co == 1) parts[i] <- paste0("t^", e)
      else if (co == -1) parts[i] <- paste0("-t^", e)
      else parts[i] <- paste0(co, "t^", e)
    }
  }
  parts <- parts[parts != "" & !is.na(parts)]
  if (length(parts) == 0) return("0")
  result <- paste(parts, collapse = " + ")
  gsub("\\+ -", "- ", result)
}

#' @export
braid_str <- function(x) {
  if (length(x) == 0) return("1")
  s <- "\u03C3"
  parts <- sapply(x, function(v) {
    g <- abs(v)
    i <- if (v < 0) "\u207B\u00B9" else ""
    paste0(s, g, i)
  })
  paste(parts, collapse = " \u00B7 ")
}