test_that("Univariate Teissier functions work correctly", {
  xi <- 2.0
  z <- 0.5
  
  # PDF
  dens <- dteissier(z, xi)
  expect_true(is.numeric(dens) && dens > 0)
  expect_equal(dteissier(0, xi), 0)
  expect_equal(dteissier(-1, xi), 0)
  expect_equal(dteissier(z, xi, log = TRUE), log(dens))
  
  # CDF and Survival
  p_val <- pteissier(z, xi)
  s_val <- steissier(z, xi)
  expect_equal(p_val + s_val, 1)
  expect_equal(pteissier(0, xi), 0)
  expect_equal(steissier(0, xi), 1)
  expect_equal(pteissier(z, xi, lower.tail = FALSE), s_val)
  
  # Quantile function roundtrip
  probs <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99)
  q_vals <- qteissier(probs, xi)
  expect_equal(pteissier(q_vals, xi), probs, tolerance = 1e-10)
  expect_equal(qteissier(0, xi), 0)
  expect_equal(qteissier(1, xi), Inf)
  
  # Random generation
  set.seed(42)
  n <- 500
  samp <- rteissier(n, xi)
  expect_equal(length(samp), n)
  expect_true(all(samp > 0))
  
  # Missing values and edge cases
  expect_true(is.na(dteissier(NA, xi)))
  expect_true(is.na(pteissier(NA, xi)))
  expect_true(is.na(qteissier(NA, xi)))
  expect_true(is.na(steissier(NA, xi)))
  
  expect_warning(res_bad <- dteissier(0.5, -1))
  expect_true(is.nan(res_bad))
})