test_that("discrete divergence handles identities and zero support", {
    p <- matrix(c(0.1, 0.2, 0.3, 0.4), 2)
    expect_equal(biv_kld_discrete(p, p), 0)
    q <- matrix(rev(as.numeric(p)), nrow = 2)
    expect_true(biv_kld_discrete(p, q) > 0)
    expect_equal(biv_kld_discrete(c(1, 0), c(0, 1)), Inf)
    expect_equal(biv_kld_discrete(c(2, 2), c(1, 3), normalize = TRUE),
                 0.5 * log(4 / 3))
})

test_that("normal divergence matches manuscript value", {
    s1 <- matrix(c(1, 1, 1, 2), 2)
    s2 <- matrix(c(3, 3, 3, 6), 2)
    value <- biv_kld_normal(c(-2, 2), s1, c(-2, 2), s2)
    expect_equal(value, 0.4319456, tolerance = 1e-6)
    expect_equal(biv_kld_normal(c(0, 0), diag(2), c(0, 0), diag(2)), 0)
})

test_that("paper model formulas reproduce stated values", {
    expect_equal(biv_kld_pareto2(1, 2), 0.4013877, tolerance = 1e-6)
    expect_equal(biv_kld_pareto2(2, 2), 0)
    expect_equal(biv_kld_independent_weibull(c(1, 2), c(1, 2)), 0)
})

test_that("invalid exact inputs are rejected", {
    expect_error(biv_kld_discrete(c(1, 1), c(1, 1)), "sum to one")
    expect_error(biv_kld_normal(c(0, 0), matrix(1, 2, 2),
                                c(0, 0), diag(2)), "positive-definite")
    expect_error(biv_kld_pareto2(-1, 2), "positive")
})