# Dimension p = 1
Sigma1 <- 0.5
Sigma2 <- 1

kl1_12 <- kldcauchy(Sigma1, Sigma2, eps = 1e-16)
kl1_21 <- kldcauchy(Sigma2, Sigma1, eps = 1e-16)

lambda <- 0.5

test_that("kl works (dim 1)", {
  expect_equal(
    round(as.numeric(kl1_21), 15),
    log( (1 + sqrt(lambda))^2 / (4*sqrt(lambda)) )
  )
  expect_equal(
    round(as.numeric(kl1_21), 15),
    log( (1 + sqrt(lambda))^2 / (4*sqrt(lambda)) )
  )
})

# Dimension p = 2

Sigma1 <- diag(0.5, nrow = 2)
Sigma2 <- diag(1, nrow = 2)

kl2_12 <- kldcauchy(Sigma1, Sigma2, eps = 1e-16)
kl2_21 <- kldcauchy(Sigma2, Sigma1, eps = 1e-16)

lambda <- as.complex(0.5)

test_that("kl works (dim 2)", {
  expect_equal(
    round(as.numeric(kl2_12), 15),
    Re(-log(lambda) + 3/sqrt(1-1/lambda) * log(sqrt(lambda) + sqrt(lambda-1)) - 3)
  )
  expect_equal(
    round(as.numeric(kl2_21), 15),
    Re(log(lambda) + 3/sqrt(1-lambda) * log(sqrt(1/lambda) + sqrt(1/lambda-1)) - 3)
  )
})

# Dimension p = 2; 2nd example
Sigma1 <- matrix(c(0.5, 0, 0, 1), nrow = 2)
Sigma2 <- diag(nrow = 2)

lambda <- 0.5

kl2 <- kldcauchy(Sigma1, Sigma2, eps = 1e-16)

test_that("kl works (dim 2, one of the eigenvalues = 1)", {
  expect_equal(
    round(as.numeric(kl2), 15),
    log(lambda) - 3/2 * 1/sqrt(1-lambda) * log((1 - sqrt(1-lambda))/(1 + sqrt(1-lambda))) - 3
  )
})


#Dimension p = 3
Sigma1 <- diag(0.5, nrow = 3)
Sigma2 <- diag(nrow = 3)

lambda <- 0.5

kl3<- kldcauchy(Sigma1, Sigma2, eps = 1e-16)

test_that("kl works (dim 3)", {
  expect_equal(
    round(as.numeric(kl3), 15),
    -3/2*log(lambda) + 4*log(0.5 + sqrt(lambda)/2) - 2*((1 - sqrt(lambda))/(1 + sqrt(lambda)))
  )
})