in_order(theta::vector[K])::int = begin
    if theta[1] == negative_infinity()
        reject("first element must be finite")
    end
    for k in 2:K
        if theta[k] <= theta[k - 1]
            reject("bounds and quantiles are not in the right order")
        end
    end
    1
end

clayton_copula_lpdf(u::real, v::real, theta::real)::real = begin
    if theta <= 0.0
        reject("clayton_copula: theta must > 0")
    end
    log1p(theta) - (theta + 1) * (log(u) + log(v)) -
        (1 + 2 * theta) / theta * log(u^(-theta) + v^(-theta) - 1)
end

frank_copula_lpdf(u::vector[N], v::vector[N], theta::real)::real = begin
    log_lik = 0.0
    if theta < 0
        log_lik = N * log(theta * (1 - exp(-theta))) - theta * sum(append_col(u, v)) -
            sum(log((-expm1(-theta) .- expm1(-theta .* u) .* expm1(-theta .* v)) .^ 2))
    else
        log_lik = N * (log(theta) + log1m_exp(-theta)) - theta * sum(append_col(u, v)) -
            2 * sum(log(-expm1(-theta) .- expm1(-theta .* u) .* expm1(-theta .* v)))
    end
    log_lik
end

gumbel_copula_lpdf(u::real, v::real, theta::real)::real = begin
    neg_log_u = -log(u)
    log_neg_log_u = log(neg_log_u)
    neg_log_v = -log(v)
    log_neg_log_v = log(neg_log_v)
    log_temp = log_sum_exp(theta * log_neg_log_u, theta * log_neg_log_v)
    theta_m1 = theta - 1
    if theta < 1
        reject("theta must be >= 1")
    end
    out = 0.0
    if is_inf(theta)
        out = u == v ? 0.0 : negative_infinity()
    else
        out = theta_m1 * log_neg_log_u + theta_m1 * log_neg_log_v + neg_log_u + neg_log_v -
            exp(log_temp / theta) + log_sum_exp(2 * theta_m1 / -theta * log_temp,
                log(theta_m1) + (1 - 2 * theta) / theta * log_temp)
    end
    out
end

normal_copula_lpdf(u::real, v::real, rho::real)::real = begin
    rho_sq = square(rho)
    (0.5 * rho * (-2.0 * u * v + square(u) * rho + square(v) * rho)) /
        (-1.0 + rho_sq) - 0.5 * log1m(rho_sq)
end

normal_copula_vector_lpdf(u::vector[N], v::vector[N], rho::real)::real = begin
    rho_sq = square(rho)
    a1 = 0.5 * rho
    a2 = rho_sq - 1
    a3 = 0.5 * log1m(rho_sq)
    x = -2 * dot_product(u, v) + rho * (dot_self(u) + dot_self(v))
    a1 * x / a2 - N * a3
end

bivariate_normal_copula_cdf(u::real, v::real, rho::real)::real = begin
    a = 1 / sqrt(1 - square(rho))
    avg_uv = (u + v) / 2
    pu = inv_Phi(u)
    pv = inv_Phi(v)
    alpha_u = a * (pv / pu - rho)
    alpha_v = a * (pu / pv - rho)
    d = 0.0
    if v == 0.5 || (u == 0.5 && v != u)
        x = u == 0.5 ? v : u
        if rho < 0
            return 0.5 * bivariate_normal_copula_cdf(x, x, 1 - 2 * rho^2)
        else
            return u - 0.5 * bivariate_normal_copula_cdf(x, x, 1 - 2 * rho^2)
        end
    end
    if v == u
        return u - 2 * owens_t(inv_Phi(u), sqrt((1 - rho) / (1 + rho)))
    end
    if (u < 0.5 && v >= 0.5) || (u >= 0.5 && v < 0.5)
        d = 0.5
    end
    avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d
end

multi_normal_cholesky_copula_lpdf(U::matrix[N,J], L::matrix[J,J])::real = begin
    Gammainv = chol2inv(L)
    -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U))
end

gpareto_lpdf(y::vector[N], ymin::real, k::real, sigma::real)::real = begin
    inv_k = inv(k)
    if k < 0 && max(y .- ymin) / sigma > -inv_k
        reject("k < 0 and max(y - ymin) / sigma > -1 / k; found k, sigma = ", k, sigma)
    end
    if sigma <= 0
        reject("sigma <= 0; found sigma = ", sigma)
    end
    abs(k) > 1e-15 ? -(1 + inv_k) * sum(log1p((y .- ymin) .* (k / sigma))) -
        N * log(sigma) : -sum(y .- ymin) / sigma - N * log(sigma)
end

gpareto_cdf(y::vector[N], ymin::real, k::real, sigma::real)::real = begin
    inv_k = inv(k)
    if k < 0 && max(y .- ymin) / sigma > -inv_k
        reject("k < 0 and max(y - ymin) / sigma > -1 / k; found k, sigma = ", k, sigma)
    end
    if sigma <= 0
        reject("sigma <= 0; found sigma = ", sigma)
    end
    abs(k) > 1e-15 ? exp(sum(log1m_exp((-inv_k) .* log1p((y .- ymin) .* (k / sigma))))) :
        exp(sum(log1m_exp(-(y .- ymin) ./ sigma)))
end

gpareto_lcdf(y::vector[N], ymin::real, k::real, sigma::real)::real = begin
    inv_k = inv(k)
    if k < 0 && max(y .- ymin) / sigma > -inv_k
        reject("k < 0 and max(y - ymin) / sigma > -1 / k; found k, sigma = ", k, sigma)
    end
    if sigma <= 0
        reject("sigma <= 0; found sigma = ", sigma)
    end
    abs(k) > 1e-15 ? sum(log1m_exp((-inv_k) .* log1p((y .- ymin) .* (k / sigma)))) :
        sum(log1m_exp(-(y .- ymin) ./ sigma))
end

gpareto_lccdf(y::vector[N], ymin::real, k::real, sigma::real)::real = begin
    inv_k = inv(k)
    if k < 0 && max(y .- ymin) / sigma > -inv_k
        reject("k < 0 and max(y - ymin) / sigma > -1 / k; found k, sigma = ", k, sigma)
    end
    if sigma <= 0
        reject("sigma <= 0; found sigma = ", sigma)
    end
    abs(k) > 1e-15 ? (-inv_k) * sum(log1p((y .- ymin) .* (k / sigma))) :
        -sum(y .- ymin) / sigma
end

gpareto_rng(ymin::real, k::real, sigma::real)::real = begin
    if sigma <= 0
        reject("sigma <= 0; found sigma = ", sigma)
    end
    u = uniform_rng(0.0, 1.0)
    abs(k) > 1e-15 ? ymin + (u^(-k) - 1) * sigma / k : ymin - sigma * log(u)
end

student_t_lcdf_stan(x::real, df::real)::real = begin
    if df <= 0
        reject("df must be > 0. Found df = ", df)
    end
    lval = 0.0
    if is_inf(x)
        lval = not_a_number()
    elseif is_inf(df)
        lval = normal_lcdf(x, 0.0, 1.0)
    else
        nx = 1 + (x / df) * x
        lval = nx > 1e100 ? -0.5 * df * (2 * log(abs(x)) - log(df)) -
            lbeta(0.5 * df, 0.5) - log(0.5 * df) :
            (df > x * x ? beta_lccdf(x * x / (df + x * x), 0.5, df / 2) :
                beta_lcdf(1 / nx, df / 2, 0.5))
        lval = x <= 0 ? lval - log2() : log1m(0.5 * exp(lval))
    end
    lval
end

student_t_lccdf_stan(x::real, df::real)::real = begin
    if df <= 0
        reject("df must be > 0. Found df = ", df)
    end
    lval = 0.0
    if is_inf(x)
        lval = not_a_number()
    elseif is_inf(df)
        lval = normal_lccdf(x, 0.0, 1.0)
    else
        nx = 1 + (x / df) * x
        lval = nx > 1e100 ? -0.5 * df * (2 * log(abs(x)) - log(df)) -
            lbeta(0.5 * df, 0.5) - log(0.5 * df) :
            (df > x * x ? beta_lccdf(x * x / (df + x * x), 0.5, df / 2) :
                beta_lcdf(1 / nx, df / 2, 0.5))
        lval = x < 0 ? log1m(0.5 * exp(lval)) : lval - log2()
    end
    lval
end

unit_johnson_lpdf(x::vector[N], mu::real, sigma::real)::real = begin
    N * log(sigma) - sum(log(x) + log1m(x) + 0.5 * log1p(square(logit(x)))) +
        std_normal_lpdf(mu + sigma * asinh(logit(x)))
end

unit_johnson_cdf(x::real, mu::real, sigma::real)::real = std_normal_cdf(mu + sigma * asinh(logit(x)))

unit_johnson_lcdf(x::real, mu::real, sigma::real)::real = std_normal_lcdf(mu + sigma * asinh(logit(x)))

unit_johnson_lccdf(x::real, mu::real, sigma::real)::real = std_normal_lccdf(mu + sigma * asinh(logit(x)))

unit_johnson_rng(mu::real, sigma::real)::real = inv_logit(sinh((inv_Phi(uniform_rng(0.0, 1.0)) - mu) / sigma))

chol_kronecker_prod(A::matrix[m,m], B::matrix[p,p])::matrix[m*p,m*p] = begin
    C = rep_matrix(0.0, m * p, m * p)
    for i in 1:m
        for j in 1:i
            if abs(A[i, j]) > 1e-12
                row_start = (i - 1) * p + 1
                row_end = i * p
                col_start = (j - 1) * p + 1
                col_end = j * p
                C[row_start:row_end, col_start:col_end] = A[i, j] * B
            end
        end
    end
    C
end

angle2chol(angle_mat::matrix[N,N])::matrix[N,N] = begin
    inv_mat = identity_matrix(N)
    inv_mat[:, 1] = cos(angle_mat[:, 1])
    for i in 2:N
        inv_mat[i, i] = prod(sin(angle_mat[i, 1:(i - 1)]))
        if i > 2
            for j in 2:(i - 1)
                inv_mat[i, j] = cos(angle_mat[i, j]) * prod(sin(angle_mat[i, 1:(j - 1)]))
            end
        end
    end
    inv_mat
end

angle_vec2angle_mat(angle::vector[N], K::int)::matrix[K,K] = begin
    mat = add_diag(identity_matrix(K), rep_vector(-1.0, K))
    count = K - 1
    pos = 1
    for k in 1:(K - 1)
        mat[(k + 1):K, k] = segment(angle, pos, count)
        pos += count
        count -= 1
    end
    mat
end

lower_elements(mat::matrix[N,N])::vector[(N*(N-1))÷2] = begin
    out = rep_vector(0.0, (N * (N - 1)) ÷ 2)
    count = N
    pos = 1
    for k in 1:(N - 1)
        count -= 1
        out[pos:(pos + count - 1)] = sub_col(mat, k + 1, k, count)
        pos += count
    end
    out
end

lognormal_qf(p::real, mu::real, sigma::real)::real = begin
    if is_nan(p) || p < 0 || p > 1
        reject("lognormal_icdf: p must be between 0 and 1; found p = ", p)
    end
    if is_nan(mu) || is_inf(mu)
        reject("lognormal_icdf: mu must be finite; found mu = ", mu)
    end
    if is_nan(sigma) || is_inf(sigma) || sigma <= 0
        reject("lognormal_icdf: sigma must be finite and > 0; found sigma = ", sigma)
    end
    exp(mu + sigma * inv_Phi(p))
end

unit_johnson_qf(p::real, mu::real, sigma::real)::real = begin
    if is_nan(p) || p < 0 || p > 1
        reject("unit_johnson_icdf: p must be between 0 and 1; found p = ", p)
    end
    if is_nan(mu) || is_inf(mu)
        reject("unit_johnson_icdf: mu must be finite; found mu = ", mu)
    end
    if is_nan(sigma) || is_inf(sigma) || sigma <= 0
        reject("unit_johnson_icdf: sigma must be finite and > 0; found sigma = ", sigma)
    end
    inv_logit(sinh((inv_Phi(p) - mu) / sigma))
end

inv_erf(x::real)::real = begin
    if is_nan(x) || is_inf(x) || x < 0 || x > 1
        reject("inv_erf: x must be finite and between 0 and 1; found x = ", x)
    end
    w = -log1m(square(x))
    s = sqrt(w)
    p = 0.0
    if w < 5.0
        w -= 2.5
        p = 2.81022636e-08
        p = fma(p, w, 3.43273939e-07)
        p = fma(p, w, -3.5233877e-06)
        p = fma(p, w, -4.39150654e-06)
        p = fma(p, w, 0.00021858087)
        p = fma(p, w, -0.00125372503)
        p = fma(p, w, -0.00417768164)
        p = fma(p, w, 0.246640727)
        p = fma(p, w, 1.50140941)
    else
        w = s - 3.0
        p = -0.000200214257
        p = fma(p, w, 0.000100950558)
        p = fma(p, w, 0.00134934322)
        p = fma(p, w, -0.00367342844)
        p = fma(p, w, 0.00573950773)
        p = fma(p, w, -0.0076224613)
        p = fma(p, w, 0.00943887047)
        p = fma(p, w, 1.00167406)
        p = fma(p, w, 2.83297682)
    end
    p * x
end

inv_erf_nsqrt(x::real)::real = begin
    if is_nan(x) || is_inf(x) || x < 0 || x > 1
        reject("inv_erf_nsqrt: x must be finite and between 0 and 1; found x = ", x)
    end
    t = log(fma(x, -x, 1.0))
    p = 0.0
    if abs(t) > 6.125
        p = 3.03697567e-10
        p = fma(p, t, 2.93243101e-8)
        p = fma(p, t, 1.22150334e-6)
        p = fma(p, t, 2.84108955e-5)
        p = fma(p, t, 3.93552968e-4)
        p = fma(p, t, 3.02698812e-3)
        p = fma(p, t, 4.83185798e-3)
        p = fma(p, t, -2.64646143e-1)
        p = fma(p, t, 8.40016484e-1)
    else
        p = 5.43877832e-9
        p = fma(p, t, 1.43285448e-7)
        p = fma(p, t, 1.22774793e-6)
        p = fma(p, t, 1.12963626e-7)
        p = fma(p, t, -5.61530760e-5)
        p = fma(p, t, -1.47697632e-4)
        p = fma(p, t, 2.31468678e-3)
        p = fma(p, t, 1.15392581e-2)
        p = fma(p, t, -2.32015476e-1)
        p = fma(p, t, 8.86226892e-1)
    end
    p * x
end