R/RcppExports.R

In sphunif: Uniformity Tests on the Circle, Sphere, and Hypersphere

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' @title Surface area of the intersection of two hyperspherical caps
#'
#' @description Computation of
#' \deqn{A_x(\theta_{ij}) := \frac{1}{\omega_p}
#' \int_{S^{p - 1}} 1_{\{{\bf X}_i'\boldsymbol\gamma \le x,
#' {\bf X}_j'\boldsymbol\gamma \le x\}}\,\mathrm{d}\boldsymbol\gamma,}{
#' A_x(\theta_{ij}) := \frac{1}{\omega_p} \int_{S^{p - 1}}
#' 1_{X_i'\gamma \le x, X_j'\gamma \le x} d\gamma,}
#' where \eqn{\theta_{ij} := \cos^{-1}({\bf X}_i'{\bf X}_j)
#' \in [0, \pi]}{\theta_{ij} := \cos^{-1}(X_i'X_j) \in [0, \pi]},
#' \eqn{x \in [-1, 1]}, and \eqn{\omega_{p}} is the surface area of
#' \eqn{S^{p - 1}}. \eqn{A_x(\theta_{ij})} is the proportion of surface area
#' of \eqn{S^{p - 1}} covered by the intersection of two hyperspherical caps
#' centered at \eqn{{\bf X}_i}{X_i} and \eqn{{\bf X}_j}{X_j} and with
#' common solid angle \eqn{\pi - \cos^{-1}(x)}.
#'
#' @param theta vector with values in \eqn{[0, \pi]}.
#' @param x vector with values in \eqn{[-1, 1]}.
#' @inheritParams r_unif
#' @param N number of points used in the
#' \link[=Gauss_Legen_nodes]{Gauss-Legendre quadrature}. Defaults to
#' \code{160}.
#' @param as_matrix return a matrix with the values of \eqn{A_x(\theta)} on
#' the grid formed by \code{theta} and \code{x}? If \code{FALSE},
#' \eqn{A_x(\theta)} is evaluated on \code{theta} and \code{x} if they equal
#' in size. Defaults to \code{TRUE}.
#' @return A matrix of size \code{c(length(theta), length(x))} containing the
#' evaluation of \eqn{A_x(\theta)} if \code{as_matrix = TRUE}. Otherwise,
#' a vector of size \code{c(length(theta)} if \code{theta} and \code{x} equal
#' in size.
#' @details
#' See García-Portugués et al. (2023) for more details about the
#' \eqn{A_x(\theta)} function.
#' @references
#' García-Portugués, E., Navarro-Esteban, P., Cuesta-Albertos, J. A. (2023)
#' On a projection-based class of uniformity tests on the hypersphere.
#' \emph{Bernoulli}, 29(1):181--204. \doi{10.3150/21-BEJ1454}.
#' @examples
#' # Plot A_x(theta) for several dimensions and x's
#' A_lines <- function(x, th = seq(0, pi, l = 200)) {
#'
#'   plot(th, A_theta_x(theta = th, x = x, p = 2), type = "l",
#'        col = 1, ylim = c(0, 1.25), main = paste("x =", x),
#'        ylab = expression(A[x](theta)),
#'        xlab = expression(theta), axes = FALSE)
#'   axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
#'        labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
#'   axis(2); box()
#'   abline(h = c(0, 1), lty = 2)
#'   lines(th, A_theta_x(theta = th, x = x, p = 3), col = 2)
#'   lines(th, A_theta_x(theta = th, x = x, p = 4), col = 3)
#'   lines(th, A_theta_x(theta = th, x = x, p = 5), col = 4)
#'   legend("top", lwd = 2, legend = paste("p =", 2:5),
#'          col = 1:4, cex = 0.75, horiz = TRUE)
#'
#' }
#' old_par <- par(mfrow = c(2, 3))
#' A_lines(x = -0.75)
#' A_lines(x = -0.25)
#' A_lines(x = 0)
#' A_lines(x = 0.25)
#' A_lines(x = 0.5)
#' A_lines(x = 0.75)
#' par(old_par)
#'
#' # As surface of (theta, x) for several dimensions
#' A_surf <- function(p, x = seq(-1, 1, l = 201), th = seq(0, pi, l = 201)) {
#'
#'   col <- c("white", viridisLite::viridis(20))
#'   breaks <- c(-1, seq(1e-15, 1, l = 21))
#'   A <- A_theta_x(theta = th, x = x, p = p)
#'   image(th, x, A, main = paste("p =", p), col = col, breaks = breaks,
#'         xlab = expression(theta), axes = FALSE)
#'   axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
#'        labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
#'   axis(2); box()
#'   contour(th, x, A, levels = breaks, add = TRUE)
#'
#' }
#' old_par <- par(mfrow = c(2, 2))
#' A_surf(p = 2)
#' A_surf(p = 3)
#' A_surf(p = 4)
#' A_surf(p = 5)
#' par(old_par)
#'
#' # No matrix return
#' th <- seq(0, pi, l = 5)
#' x <- seq(-1, 1, l = 5)
#' diag(A_theta_x(theta = th, x = x, p = 2))
#' A_theta_x(theta = th, x = x, p = 2, as_matrix = FALSE)
#' @export
A_theta_x <- function(theta, x, p, N = 160L, as_matrix = TRUE) {
    .Call('_sphunif_A_theta_x', PACKAGE = 'sphunif', theta, x, p, N, as_matrix)
}

#' @title Shortest angles matrix
#'
#' @description Efficient computation of the shortest angles matrix
#' \eqn{\boldsymbol\Psi}{\Psi}, defined as
#' \deqn{\Psi_{ij}:=\cos^{-1}({\bf X}_i'{\bf X}_j),\quad
#' i,j=1,\ldots,n,}{\Psi_{ij} = \cos^{-1}(X_i'X_j), i, j = 1, \ldots, n,}
#' for a sample \eqn{{\bf X}_1,\ldots,{\bf X}_n\in S^{p-1}:=\{{\bf x}\in
#' R^p:||{\bf x}||=1\}}{X_1, \ldots, X_n \in
#' S^{p - 1} := \{x \in R^p : ||x|| = 1\}}, \eqn{p\ge 2}.
#'
#' For a circular sample \eqn{\Theta_1, \ldots, \Theta_n \in [0, 2\pi)},
#' \eqn{\boldsymbol\Psi}{\Psi} can be expressed as
#' \deqn{\Psi_{ij}=\pi-|\pi-|\Theta_i-\Theta_j||,\quad
#' i,j=1,\ldots,n.}{\Psi_{ij}=\pi-|\pi-|\Theta_i-\Theta_j||, i,j=1,\ldots,n.}
#'
#' @param data an array of size \code{c(n, p, M)} containing the Cartesian
#' coordinates of \code{M} samples of size \code{n} of directions on
#' \eqn{S^{p-1}}. Alternatively if \code{p = 2}, an array of size
#' \code{c(n, 1, M)} containing the angles on \eqn{[0, 2\pi)} of the \code{M}
#' circular samples of size \code{n} on \eqn{S^{1}}. Must not contain
#' \code{NA}'s.
#' @param ind_tri if \code{use_ind_tri = TRUE}, the vector of 0-based indexes
#' provided by \code{upper_tri_ind(n)}, which allows to extract the upper
#' triangular part of the matrix \eqn{\boldsymbol\Psi}{\Psi}. See the examples.
#' @param use_ind_tri use the already computed vector index \code{ind_tri}? If
#' \code{FALSE} (default), \code{ind_tri} is computed internally.
#' @param scalar_prod return the scalar products
#' \eqn{{\bf X}_i'{\bf X}}{X_i'X_j} instead of the shortest angles? Only taken
#' into account for data in \emph{Cartesian} form. Defaults to
#' \code{FALSE}.
#' @param angles_diff return the (unwrapped) angles difference
#' \eqn{\Theta_i-\Theta_j} instead of the shortest angles? Only taken into
#' account for data in \emph{angular} form. Defaults to \code{FALSE}.
#' @param n sample size, used to determine the index vector that gives the
#' upper triangular part of \eqn{\boldsymbol\Psi}{\Psi}.
#' @return
#' \itemize{
#'   \item \code{Psi_mat}: a matrix of size
#'   \code{c(n * (n - 1) / 2, M)} containing, for each column, the vector
#'   half of \eqn{\boldsymbol\Psi}{\Psi} for each of the \code{M} samples.
#'   \item \code{upper_tri_ind}: a matrix of size \code{n * (n - 1) / 2}
#'   containing the 0-based linear indexes for extracting the upper triangular
#'   matrix of a matrix of size \code{c(n, n)}, diagonal excluded, assuming
#'   column-major order.
#' }
#' @section Warning:
#' Be careful on avoiding the next bad usages of \code{Psi_mat}, which will
#' produce spurious results:
#' \itemize{
#'   \item The directions in \code{data} do \emph{not} have unit norm when
#'   Cartesian coordinates are employed.
#'   \item The entries of \code{data} are \emph{not} in \eqn{[0, 2\pi)} when
#'   polar coordinates are employed.
#'   \item \code{ind_tri} is a vector of size \code{n * (n - 1) / 2} that
#'   does \emph{not} contain the indexes produced by \code{upper_tri_ind(n)}.
#' }
#' @examples
#' # Shortest angles
#' n <- 5
#' X <- r_unif_sph(n = n, p = 2, M = 2)
#' Theta <- X_to_Theta(X)
#' dim(Theta) <- c(n, 1, 2)
#' Psi_mat(X)
#' Psi_mat(Theta)
#'
#' # Precompute ind_tri
#' ind_tri <- upper_tri_ind(n)
#' Psi_mat(X, ind_tri = ind_tri, use_ind_tri = TRUE)
#'
#' # Compare with R
#' A <- acos(tcrossprod(X[, , 1]))
#' ind <- upper.tri(A)
#' A[ind]
#'
#' # Reconstruct matrix
#' Psi_vec <- Psi_mat(Theta[, , 1, drop = FALSE])
#' Psi <- matrix(0, nrow = n, ncol = n)
#' Psi[upper.tri(Psi)] <- Psi_vec
#' Psi <- Psi + t(Psi)
#' @name Psi
NULL

#' @rdname Psi
#' @export
Psi_mat <- function(data, ind_tri = 0L, use_ind_tri = FALSE, scalar_prod = FALSE, angles_diff = FALSE) {
    .Call('_sphunif_Psi_mat', PACKAGE = 'sphunif', data, ind_tri, use_ind_tri, scalar_prod, angles_diff)
}

#' @rdname Psi
#' @export
upper_tri_ind <- function(n) {
    .Call('_sphunif_upper_tri_ind', PACKAGE = 'sphunif', n)
}

#' @title Sort the columns of a matrix
#'
#' @description Convenience functions to sort the columns of a matrix in an
#' increasing way.
#'
#' @param A a matrix of arbitrary dimensions.
#' @return A matrix with the same dimensions as \code{A} such that each
#' column is sorted increasingly (for \code{sort_each_col}) or contains the
#' sorting indexes (for \code{sort_index_each_col}).
#' \code{sort_index_each_col}.
#' @keywords internal
sort_each_col <- function(A) {
    .Call('_sphunif_sort_each_col', PACKAGE = 'sphunif', A)
}

#' @rdname sort_each_col
#' @keywords internal
sort_index_each_col <- function(A) {
    .Call('_sphunif_sort_index_each_col', PACKAGE = 'sphunif', A)
}

#' @title Transforming between polar and Cartesian coordinates
#'
#' @description Transformation between a matrix \code{Theta} containing
#' \code{M} circular samples of size \code{n} on \eqn{[0, 2\pi)} and an array
#' \code{X} containing the associated Cartesian coordinates on
#' \eqn{S^1:=\{{\bf x}\in R^2:||{\bf x}||=1\}}{S^1:=\{x\in R^2:||x||=1\}}.
#'
#' @inheritParams cir_stat
#' @param X an \bold{array} of size \code{c(n, 2, M)} containing the Cartesian
#' coordinates of \code{M} samples of size \code{n} of directions on
#' \eqn{S^{1}}. Must not contain \code{NA}'s.
#' @return
#' \itemize{
#'   \item \code{Theta_to_X}: the corresponding \code{X}.
#'   \item \code{X_to_Theta}: the corresponding \code{Theta}.
#' }
#' @examples
#' # Sample
#' Theta <- r_unif_cir(n = 10, M = 2)
#' X <- r_unif_sph(n = 10, p = 2, M = 2)
#'
#' # Check equality
#' sum(abs(X - Theta_to_X(X_to_Theta(X))))
#' sum(abs(Theta - X_to_Theta(Theta_to_X(Theta))))
#' @name cir_coord_conv
NULL

#' @title Low-level utilities for \pkg{sphunif}
#'
#' @description Internal and undocumented low-level utilities for
#' \pkg{sphunif}.
#'
#' @param n_dist a positive integer \eqn{(n - 1) * n / 2} for which \eqn{n}
#' is to be recovered.
#' @param t a vector to evaluate \eqn{t / \sqrt{1 - t^2}}.
#' @name utils
NULL

#' @rdname cir_coord_conv
#' @export
Theta_to_X <- function(Theta) {
    .Call('_sphunif_Theta_to_X', PACKAGE = 'sphunif', Theta)
}

#' @rdname cir_coord_conv
#' @export
X_to_Theta <- function(X) {
    .Call('_sphunif_X_to_Theta', PACKAGE = 'sphunif', X)
}

#' @title Circular gaps
#'
#' @description Computation of the circular gaps of an angular sample
#' \eqn{\Theta_1,\ldots,\Theta_n} on \eqn{[0, 2\pi)}, defined as
#' \deqn{\Theta_{(2)} - \Theta_{(1)},\ldots,\Theta_{(n)} - \Theta_{(n - 1)},
#' 2\pi - \Theta_{(n)} - \Theta_{(1)},}
#' where
#' \deqn{0 \le \Theta_{(1)} \le \Theta_{(2)} \le \ldots \le
#' \Theta_{(n)} \le 2\pi.}
#'
#' @inheritParams cir_stat
#' @param sorted are the columns of \code{Theta} sorted increasingly? If
#' \code{TRUE}, performance is improved. If \code{FALSE} (default), each
#' column of \code{Theta} is sorted internally.
#' @return A matrix of size \code{c(n, M)} containing the \code{n} circular
#' gaps for each of the \code{M} circular samples.
#' @section Warning:
#' Be careful on avoiding the next bad usages of \code{cir_gaps}, which will
#' produce spurious results:
#' \itemize{
#'   \item The entries of \code{Theta} are \emph{not} in \eqn{[0, 2\pi)}.
#'   \item \code{Theta} is \emph{not} sorted increasingly when
#'   \code{data_sorted = TRUE}.
#' }
#' @examples
#' Theta <- cbind(c(pi, 0, 3 * pi / 2), c(0, 3 * pi / 2, pi), c(5, 3, 1))
#' cir_gaps(Theta)
#' @export
cir_gaps <- function(Theta, sorted = FALSE) {
    .Call('_sphunif_cir_gaps', PACKAGE = 'sphunif', Theta, sorted)
}

#' @title Efficient evaluation of the empirical cumulative distribution
#' function
#'
#' @description Evaluates the empirical cumulative distribution function
#' (ecdf) of a sample \code{data} at the evaluation points \code{sorted_x}.
#' This is done through binary search.
#'
#' @param data a vector or column matrix containing the sample.
#' @param sorted_x a vector or column matrix with the evaluation points
#' \bold{sorted increasingly}.
#' @param data_sorted is \code{data} is already sorted increasingly?
#' This avoids sorting the data internally.
#' @param efic use the more efficient version of the ecdf evaluation? Set to
#' \code{FALSE} only for debugging purposes.
#' @param divide_n if \code{FALSE}, returns the absolute frequencies instead
#' of the relative frequencies. Defaults to \code{TRUE}.
#' @return The ecdf evaluated at \code{sorted_x}.
#' @author Original code from Douglas Bates'
#' \url{https://github.com/dmbates/ecdfExample}. Minor adaptations by Eduardo
#' García-Portugués.
#' @section Warning:
#' Be careful on avoiding the next bad usages of the function, which will
#' produce spurious results:
#' \itemize{
#'   \item \code{sorted_x} is not sorted increasingly.
#'   \item \code{data} is not sorted increasingly when
#'   \code{data_sorted = TRUE}-
#' }
#' @keywords internal
ecdf_bin <- function(data, sorted_x, data_sorted = FALSE, efic = TRUE, divide_n = TRUE) {
    .Call('_sphunif_ecdf_bin', PACKAGE = 'sphunif', data, sorted_x, data_sorted, efic, divide_n)
}

#' @title The incomplete beta function and its inverse
#'
#' @description Computes the incomplete beta function
#' \deqn{I_x(a,b):=\int_0^x u^{a-1}(1-u)^{b-1}\,d\mathrm{u},\quad a,b>0}{
#' I_x(a,b):=\int_0^x u^{a-1}(1-u)^{b-1}du, a,b>0}
#' and its inverse function.
#'
#' @inheritParams cir_stat_distr
#' @param u a vector of probabilities of size \code{nu} or a matrix of size
#' \code{c(nu, 1)}.
#' @param a,b scalars giving the parameters of the beta function.
#' @param lower_tail accumulate the probability from the lower tail? If
#' \code{FALSE}, the probability is accumulated from the \emph{upper} tail.
#' Defaults to \code{FALSE}.
#' @param log use log-scale? If \code{TRUE}, returns the logarithm of the
#' incomplete beta function and uses log-scale for \code{u} in \code{beta_inc}.
#' Defaults to \code{FALSE}.
#' @return
#' \itemize{
#'   \item \code{beta_inc}: a matrix of size \code{c(nx, 1)} with the
#'   evaluation of the incomplete beta function at \code{x}.
#'   \item \code{beta_inc_inv}: a matrix of size \code{c(nu, 1)} with the
#'   evaluation of the inverse incomplete beta function at \code{u}.
#' }
#' @details
#' The functions are mere wrappers to R's internal \code{pbeta} and
#' \code{qbeta} functions.
#' @keywords internal
beta_inc <- function(x, a, b, lower_tail = TRUE, log = FALSE) {
    .Call('_sphunif_beta_inc', PACKAGE = 'sphunif', x, a, b, lower_tail, log)
}

#' @rdname beta_inc
beta_inc_inv <- function(u, a, b, lower_tail = TRUE, log = FALSE) {
    .Call('_sphunif_beta_inc_inv', PACKAGE = 'sphunif', u, a, b, lower_tail, log)
}

#' @rdname utils
#' @keywords internal
n_from_dist_vector <- function(n_dist) {
    .Call('_sphunif_n_from_dist_vector', PACKAGE = 'sphunif', n_dist)
}

#' @rdname utils
#' @keywords internal
t_inv_sqrt_one <- function(t) {
    .Call('_sphunif_t_inv_sqrt_one', PACKAGE = 'sphunif', t)
}

#' @rdname cir_stat
#' @export
cir_stat_Kuiper <- function(Theta, sorted = FALSE, KS = FALSE, Stephens = FALSE) {
    .Call('_sphunif_cir_stat_Kuiper', PACKAGE = 'sphunif', Theta, sorted, KS, Stephens)
}

#' @rdname cir_stat
#' @export
cir_stat_Watson <- function(Theta, sorted = FALSE, CvM = FALSE, Stephens = FALSE) {
    .Call('_sphunif_cir_stat_Watson', PACKAGE = 'sphunif', Theta, sorted, CvM, Stephens)
}

#' @rdname cir_stat
#' @export
cir_stat_Watson_1976 <- function(Theta, sorted = FALSE, minus = FALSE) {
    .Call('_sphunif_cir_stat_Watson_1976', PACKAGE = 'sphunif', Theta, sorted, minus)
}

#' @rdname cir_stat
#' @export
cir_stat_Range <- function(Theta, sorted = FALSE, gaps_in_Theta = FALSE, max_gap = TRUE) {
    .Call('_sphunif_cir_stat_Range', PACKAGE = 'sphunif', Theta, sorted, gaps_in_Theta, max_gap)
}

#' @rdname cir_stat
#' @export
cir_stat_Rao <- function(Theta, sorted = FALSE, gaps_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Rao', PACKAGE = 'sphunif', Theta, sorted, gaps_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Greenwood <- function(Theta, sorted = FALSE, gaps_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Greenwood', PACKAGE = 'sphunif', Theta, sorted, gaps_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Log_gaps <- function(Theta, sorted = FALSE, gaps_in_Theta = FALSE, abs_val = TRUE) {
    .Call('_sphunif_cir_stat_Log_gaps', PACKAGE = 'sphunif', Theta, sorted, gaps_in_Theta, abs_val)
}

#' @rdname cir_stat
#' @export
cir_stat_Vacancy <- function(Theta, a = 2 * pi, sorted = FALSE, gaps_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Vacancy', PACKAGE = 'sphunif', Theta, a, sorted, gaps_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Max_uncover <- function(Theta, a = 2 * pi, sorted = FALSE, gaps_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Max_uncover', PACKAGE = 'sphunif', Theta, a, sorted, gaps_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Num_uncover <- function(Theta, a = 2 * pi, sorted = FALSE, gaps_in_Theta = FALSE, minus_val = TRUE) {
    .Call('_sphunif_cir_stat_Num_uncover', PACKAGE = 'sphunif', Theta, a, sorted, gaps_in_Theta, minus_val)
}

#' @rdname cir_stat
#' @export
cir_stat_Gini <- function(Theta, sorted = FALSE, gaps_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Gini', PACKAGE = 'sphunif', Theta, sorted, gaps_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Gini_squared <- function(Theta, sorted = FALSE, gaps_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Gini_squared', PACKAGE = 'sphunif', Theta, sorted, gaps_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Ajne <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Ajne', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @keywords internal
cir_stat_An_Psi <- function(Psi, n) {
    .Call('_sphunif_cir_stat_An_Psi', PACKAGE = 'sphunif', Psi, n)
}

#' @rdname cir_stat
#' @export
cir_stat_Rothman <- function(Theta, t = 1.0 / 3.0, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Rothman', PACKAGE = 'sphunif', Theta, t, Psi_in_Theta)
}

#' @keywords internal
cir_stat_Rothman_Psi <- function(Psi, t_m2, t_min2, n) {
    .Call('_sphunif_cir_stat_Rothman_Psi', PACKAGE = 'sphunif', Psi, t_m2, t_min2, n)
}

#' @rdname cir_stat
#' @export
cir_stat_Hodges_Ajne <- function(Theta, asymp_std = FALSE, sorted = FALSE, use_Cressie = TRUE) {
    .Call('_sphunif_cir_stat_Hodges_Ajne', PACKAGE = 'sphunif', Theta, asymp_std, sorted, use_Cressie)
}

#' @rdname cir_stat
#' @export
cir_stat_Cressie <- function(Theta, t = 1.0 / 3.0, sorted = FALSE) {
    .Call('_sphunif_cir_stat_Cressie', PACKAGE = 'sphunif', Theta, t, sorted)
}

#' @rdname cir_stat
#' @export
cir_stat_FG01 <- function(Theta, sorted = FALSE) {
    .Call('_sphunif_cir_stat_FG01', PACKAGE = 'sphunif', Theta, sorted)
}

#' @rdname cir_stat
#' @export
cir_stat_Rayleigh <- function(Theta, m = 1L) {
    .Call('_sphunif_cir_stat_Rayleigh', PACKAGE = 'sphunif', Theta, m)
}

#' @rdname cir_stat
#' @export
cir_stat_Bingham <- function(Theta) {
    .Call('_sphunif_cir_stat_Bingham', PACKAGE = 'sphunif', Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Hermans_Rasson <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Hermans_Rasson', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @keywords internal
cir_stat_Hermans_Rasson_Psi <- function(Psi, n) {
    .Call('_sphunif_cir_stat_Hermans_Rasson_Psi', PACKAGE = 'sphunif', Psi, n)
}

#' @rdname cir_stat
#' @export
cir_stat_Gine_Gn <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Gine_Gn', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Gine_Fn <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Gine_Fn', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Pycke <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Pycke', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @keywords internal
cir_stat_Pycke_Psi <- function(Psi, n) {
    .Call('_sphunif_cir_stat_Pycke_Psi', PACKAGE = 'sphunif', Psi, n)
}

#' @rdname cir_stat
#' @export
cir_stat_Pycke_q <- function(Theta, Psi_in_Theta = FALSE, q = 0.5) {
    .Call('_sphunif_cir_stat_Pycke_q', PACKAGE = 'sphunif', Theta, Psi_in_Theta, q)
}

#' @keywords internal
cir_stat_Pycke_q_Psi <- function(Psi, n, q = 0.5) {
    .Call('_sphunif_cir_stat_Pycke_q_Psi', PACKAGE = 'sphunif', Psi, n, q)
}

#' @rdname cir_stat
#' @export
cir_stat_Bakshaev <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_Bakshaev', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_Riesz <- function(Theta, Psi_in_Theta = FALSE, s = 1.0) {
    .Call('_sphunif_cir_stat_Riesz', PACKAGE = 'sphunif', Theta, Psi_in_Theta, s)
}

#' @rdname cir_stat
#' @export
cir_stat_PCvM <- function(Theta, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_PCvM', PACKAGE = 'sphunif', Theta, Psi_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_PRt <- function(Theta, t = 1.0 / 3.0, Psi_in_Theta = FALSE) {
    .Call('_sphunif_cir_stat_PRt', PACKAGE = 'sphunif', Theta, t, Psi_in_Theta)
}

#' @rdname cir_stat
#' @export
cir_stat_PAD <- function(Theta, Psi_in_Theta = FALSE, AD = FALSE, sorted = FALSE) {
    .Call('_sphunif_cir_stat_PAD', PACKAGE = 'sphunif', Theta, Psi_in_Theta, AD, sorted)
}

#' @rdname cir_stat
#' @export
cir_stat_CCF09 <- function(Theta, dirs, K_CCF09 = 25L, original = FALSE) {
    .Call('_sphunif_cir_stat_CCF09', PACKAGE = 'sphunif', Theta, dirs, K_CCF09, original)
}

#' @rdname cir_stat_distr
#' @export
p_Kolmogorov <- function(x, K_Kolmogorov = 25L, alternating = TRUE) {
    .Call('_sphunif_p_Kolmogorov', PACKAGE = 'sphunif', x, K_Kolmogorov, alternating)
}

#' @rdname cir_stat_distr
#' @export
d_Kolmogorov <- function(x, K_Kolmogorov = 25L, alternating = TRUE) {
    .Call('_sphunif_d_Kolmogorov', PACKAGE = 'sphunif', x, K_Kolmogorov, alternating)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Ajne <- function(x, K_Ajne = 15L) {
    .Call('_sphunif_p_cir_stat_Ajne', PACKAGE = 'sphunif', x, K_Ajne)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Ajne <- function(x, K_Ajne = 15L) {
    .Call('_sphunif_d_cir_stat_Ajne', PACKAGE = 'sphunif', x, K_Ajne)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Bingham <- function(x) {
    .Call('_sphunif_p_cir_stat_Bingham', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Bingham <- function(x) {
    .Call('_sphunif_d_cir_stat_Bingham', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Greenwood <- function(x) {
    .Call('_sphunif_p_cir_stat_Greenwood', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Greenwood <- function(x) {
    .Call('_sphunif_d_cir_stat_Greenwood', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Gini <- function(x) {
    .Call('_sphunif_p_cir_stat_Gini', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Gini <- function(x) {
    .Call('_sphunif_d_cir_stat_Gini', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Gini_squared <- function(x) {
    .Call('_sphunif_p_cir_stat_Gini_squared', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Gini_squared <- function(x) {
    .Call('_sphunif_d_cir_stat_Gini_squared', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Hodges_Ajne2 <- function(x, n, asymp_std = FALSE) {
    .Call('_sphunif_p_cir_stat_Hodges_Ajne2', PACKAGE = 'sphunif', x, n, asymp_std)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Hodges_Ajne <- function(x, n, exact = TRUE, asymp_std = FALSE) {
    .Call('_sphunif_p_cir_stat_Hodges_Ajne', PACKAGE = 'sphunif', x, n, exact, asymp_std)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Hodges_Ajne <- function(x, n, exact = TRUE, asymp_std = FALSE) {
    .Call('_sphunif_d_cir_stat_Hodges_Ajne', PACKAGE = 'sphunif', x, n, exact, asymp_std)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Kuiper <- function(x, n, K_Kuiper = 12L, second_term = TRUE, Stephens = FALSE) {
    .Call('_sphunif_p_cir_stat_Kuiper', PACKAGE = 'sphunif', x, n, K_Kuiper, second_term, Stephens)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Kuiper <- function(x, n, K_Kuiper = 12L, second_term = TRUE, Stephens = FALSE) {
    .Call('_sphunif_d_cir_stat_Kuiper', PACKAGE = 'sphunif', x, n, K_Kuiper, second_term, Stephens)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Log_gaps <- function(x, abs_val = TRUE) {
    .Call('_sphunif_p_cir_stat_Log_gaps', PACKAGE = 'sphunif', x, abs_val)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Log_gaps <- function(x, abs_val = TRUE) {
    .Call('_sphunif_d_cir_stat_Log_gaps', PACKAGE = 'sphunif', x, abs_val)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Max_uncover <- function(x) {
    .Call('_sphunif_p_cir_stat_Max_uncover', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Max_uncover <- function(x) {
    .Call('_sphunif_d_cir_stat_Max_uncover', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Num_uncover <- function(x) {
    .Call('_sphunif_p_cir_stat_Num_uncover', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Num_uncover <- function(x) {
    .Call('_sphunif_d_cir_stat_Num_uncover', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Pycke <- function(x) {
    .Call('_sphunif_p_cir_stat_Pycke', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Pycke <- function(x) {
    .Call('_sphunif_d_cir_stat_Pycke', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Vacancy <- function(x) {
    .Call('_sphunif_p_cir_stat_Vacancy', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Vacancy <- function(x) {
    .Call('_sphunif_d_cir_stat_Vacancy', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Watson <- function(x, n = 0L, K_Watson = 25L, Stephens = FALSE) {
    .Call('_sphunif_p_cir_stat_Watson', PACKAGE = 'sphunif', x, n, K_Watson, Stephens)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Watson <- function(x, n = 0L, K_Watson = 25L, Stephens = FALSE) {
    .Call('_sphunif_d_cir_stat_Watson', PACKAGE = 'sphunif', x, n, K_Watson, Stephens)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Watson_1976 <- function(x, K_Watson_1976 = 8L, N = 40L) {
    .Call('_sphunif_p_cir_stat_Watson_1976', PACKAGE = 'sphunif', x, K_Watson_1976, N)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Watson_1976 <- function(x, K_Watson_1976 = 8L) {
    .Call('_sphunif_d_cir_stat_Watson_1976', PACKAGE = 'sphunif', x, K_Watson_1976)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Range <- function(x, n, max_gap = TRUE) {
    .Call('_sphunif_p_cir_stat_Range', PACKAGE = 'sphunif', x, n, max_gap)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Range <- function(x, n, max_gap = TRUE) {
    .Call('_sphunif_d_cir_stat_Range', PACKAGE = 'sphunif', x, n, max_gap)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Rao <- function(x) {
    .Call('_sphunif_p_cir_stat_Rao', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Rao <- function(x) {
    .Call('_sphunif_d_cir_stat_Rao', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
p_cir_stat_Rayleigh <- function(x) {
    .Call('_sphunif_p_cir_stat_Rayleigh', PACKAGE = 'sphunif', x)
}

#' @rdname cir_stat_distr
#' @export
d_cir_stat_Rayleigh <- function(x) {
    .Call('_sphunif_d_cir_stat_Rayleigh', PACKAGE = 'sphunif', x)
}

#' @title Projection of the spherical uniform distribution
#'
#' @description Density, distribution, and quantile functions of the
#' projection of the spherical uniform random variable on an arbitrary
#' direction, that is, the random variable
#' \eqn{\boldsymbol{\gamma}'{\bf X}}{\gamma'X}, where \eqn{{\bf X}}{X}
#' is uniformly distributed on the (hyper)sphere
#' \eqn{S^{p-1}:=\{{\bf x}\in R^p:||{\bf x}||=1\}}{S^{p-1}:=
#' \{x\in R^p:||x||=1\}}, \eqn{p\ge 2}, and
#' \eqn{\boldsymbol{\gamma}\in S^{p-1}}{\gamma\in S^{p-1}} is an
#' \emph{arbitrary} projection direction. Note that the distribution is
#' invariant to the choice of \eqn{\boldsymbol{\gamma}}{\gamma}. Also,
#' efficient simulation of \eqn{\boldsymbol{\gamma}'{\bf X}}{\gamma'X}.
#'
#' @inheritParams cir_stat_distr
#' @inheritParams r_unif
#' @param u vector of probabilities.
#' @param log compute the logarithm of the density or distribution?
#' @return A matrix of size \code{c(nx, 1)} with the evaluation of the
#' density, distribution, or quantile function at \code{x} or \code{u}.
#' For \code{r_proj_unif}, a random vector of size \code{n}.
#' @author Eduardo García-Portugués and Paula Navarro-Esteban.
#' @examples
#' # Density function
#' curve(d_proj_unif(x, p = 2), from = -2, to = 2, n = 2e2, ylim = c(0, 2))
#' curve(d_proj_unif(x, p = 3), n = 2e2, col = 2, add = TRUE)
#' curve(d_proj_unif(x, p = 4), n = 2e2, col = 3, add = TRUE)
#' curve(d_proj_unif(x, p = 5), n = 2e2, col = 4, add = TRUE)
#' curve(d_proj_unif(x, p = 6), n = 2e2, col = 5, add = TRUE)
#'
#' # Distribution function
#' curve(p_proj_unif(x, p = 2), from = -2, to = 2, n = 2e2, ylim = c(0, 1))
#' curve(p_proj_unif(x, p = 3), n = 2e2, col = 2, add = TRUE)
#' curve(p_proj_unif(x, p = 4), n = 2e2, col = 3, add = TRUE)
#' curve(p_proj_unif(x, p = 5), n = 2e2, col = 4, add = TRUE)
#' curve(p_proj_unif(x, p = 6), n = 2e2, col = 5, add = TRUE)
#'
#' # Quantile function
#' curve(q_proj_unif(u = x, p = 2), from = 0, to = 1, n = 2e2, ylim = c(-1, 1))
#' curve(q_proj_unif(u = x, p = 3), n = 2e2, col = 2, add = TRUE)
#' curve(q_proj_unif(u = x, p = 4), n = 2e2, col = 3, add = TRUE)
#' curve(q_proj_unif(u = x, p = 5), n = 2e2, col = 4, add = TRUE)
#' curve(q_proj_unif(u = x, p = 6), n = 2e2, col = 5, add = TRUE)
#'
#' # Sampling
#' hist(r_proj_unif(n = 1e4, p = 4), freq = FALSE, breaks = 50)
#' curve(d_proj_unif(x, p = 4), n = 2e2, col = 3, add = TRUE)
#' @name proj_unif
NULL

#' @title Sample uniformly distributed circular and spherical data
#'
#' @description Simulation of the uniform distribution on \eqn{[0, 2\pi)} and
#' \eqn{S^{p-1}:=\{{\bf x}\in R^p:||{\bf x}||=1\}}{
#' S^{p-1}:=\{x\in R^p:||x||=1\}}, \eqn{p\ge 2}.
#'
#' @param n sample size.
#' @param M number of samples of size \code{n}. Defaults to \code{1}.
#' @param p integer giving the dimension of the ambient space \eqn{R^p} that
#' contains \eqn{S^{p-1}}.
#' @param sorted return each circular sample sorted? Defaults to \code{FALSE}.
#' @return
#' \itemize{
#'   \item \code{r_unif_cir}: a \bold{matrix} of size \code{c(n, M)} with
#'   \code{M} random samples of size \code{n} of uniformly-generated circular
#'   data on \eqn{[0, 2\pi)}.
#'   \item \code{r_unif_sph}: an \bold{array} of size \code{c(n, p, M)} with
#'   \code{M} random samples of size \code{n} of uniformly-generated
#'   directions on \eqn{S^{p-1}}.
#' }
#' @examples
#' # A sample on [0, 2*pi)
#' n <- 5
#' r_unif_cir(n = n)
#'
#' # A sample on S^1
#' p <- 2
#' samp <- r_unif_sph(n = n, p = p)
#' samp
#' rowSums(samp^2)
#'
#' # A sample on S^2
#' p <- 3
#' samp <- r_unif_sph(n = n, p = p)
#' samp
#' rowSums(samp^2)
#' @name r_unif
NULL

#' @title Utilities for weighted sums of non-central chi squared random
#' variables
#'
#' @description Simulation from a weighted sum of non-central chi squared
#' random variables and Monte Carlo approximation of its distribution function.
#'
#' @inheritParams r_unif
#' @inheritParams wschisq
#' @param ncps non-negative non-centrality parameters. A vector with the same
#' length as \code{weights}.
#' @param M number of Monte Carlo samples for approximating the distribution.
#' Defaults to \code{1e4}.
#' @param sample if \code{use_sample = TRUE}, the Monte Carlo sample to
#' approximate the distribution. If not, it is computed internally. Defaults
#' to \code{1e4}.
#' @param use_sample use the already computed \code{sample}? If \code{FALSE}
#' (default), \code{sample} is computed internally.
#' @return
#' \itemize{
#'   \item \code{r_wschisq_Cpp}: a matrix of size \code{c(n, 1)}
#'   containing a random sample.
#'   \item \code{p_wschisq_MC}: a matrix of size \code{c(nx, 1)}
#'   with the evaluation of the distribution function at \code{x}.
#' }
#' @name wschisq_utils
NULL

#' @title Density and distribution of a chi squared
#'
#' @description Computation of the density and distribution functions of a chi
#' squared.
#'
#' @inheritParams cir_stat_distr
#' @param df degrees of freedom.
#' @param ncp non-centrality parameter.
#' @return A matrix of size \code{c(nx, 1)} with the evaluation of the
#' density or distribution function at \code{x}.
#' @name chisq
NULL

#' @rdname proj_unif
#' @export
d_proj_unif <- function(x, p, log = FALSE) {
    .Call('_sphunif_d_proj_unif', PACKAGE = 'sphunif', x, p, log)
}

#' @rdname proj_unif
#' @export
p_proj_unif <- function(x, p, log = FALSE) {
    .Call('_sphunif_p_proj_unif', PACKAGE = 'sphunif', x, p, log)
}

#' @rdname proj_unif
#' @export
q_proj_unif <- function(u, p) {
    .Call('_sphunif_q_proj_unif', PACKAGE = 'sphunif', u, p)
}

#' @rdname r_unif
#' @export
r_unif_cir <- function(n, M = 1L, sorted = FALSE) {
    .Call('_sphunif_r_unif_cir', PACKAGE = 'sphunif', n, M, sorted)
}

#' @rdname r_unif
#' @export
r_unif_sph <- function(n, p, M = 1L) {
    .Call('_sphunif_r_unif_sph', PACKAGE = 'sphunif', n, p, M)
}

#' @rdname wschisq_utils
#' @keywords internal
r_wschisq_Cpp <- function(n, weights, dfs, ncps) {
    .Call('_sphunif_r_wschisq_Cpp', PACKAGE = 'sphunif', n, weights, dfs, ncps)
}

#' @rdname wschisq_utils
#' @keywords internal
p_wschisq_MC <- function(x, weights = 0L, dfs = 0L, ncps = 0L, M = 1e4L, sample = 0L, use_sample = FALSE, x_sorted = FALSE) {
    .Call('_sphunif_p_wschisq_MC', PACKAGE = 'sphunif', x, weights, dfs, ncps, M, sample, use_sample, x_sorted)
}

#' @rdname chisq
#' @keywords internal
d_chisq <- function(x, df, ncp = 0L) {
    .Call('_sphunif_d_chisq', PACKAGE = 'sphunif', x, df, ncp)
}

#' @rdname chisq
#' @keywords internal
p_chisq <- function(x, df, ncp = 0L) {
    .Call('_sphunif_p_chisq', PACKAGE = 'sphunif', x, df, ncp)
}

#' @title Gauss--Legendre quadrature
#'
#' @description Convenience for computing the nodes \eqn{x_k} and weights
#' \eqn{w_k} of the \emph{Gauss--Legendre} quadrature formula
#' in \eqn{(a, b)}:
#' \deqn{\int_a^b f(x) w(x)\,\mathrm{d}x\approx\sum_{k=1}^N w_k f(x_k).}{
#' \int_a^b f(x) dx\approx\sum_{k=1}^N w_k f(x_k)}.
#'
#' @param a,b scalars giving the interval \eqn{(a, b)}. Defaults to
#' \eqn{(-1, 1)}.
#' @param N number of points used in the Gauss--Legendre quadrature. The
#' following choices are supported: \code{5}, \code{10}, \code{20},
#' \code{40}, \code{80}, \code{160}, \code{320}, \code{640}, \code{1280},
#' \code{2560}, and \code{5120}. Defaults to \code{40}.
#' @return A matrix of size \code{c(N, 1)} with the nodes \eqn{x_k}
#' (\code{Gauss_Legen_nodes}) or the corresponding weights \eqn{w_k}
#' (\code{Gauss_Legen_weights}).
#' @details For \eqn{C^\infty} functions, Gauss--Legendre quadrature
#' can be very efficient. It is exact for polynomials up to degree
#' \eqn{2N - 1}.
#'
#' The nodes and weights up to \eqn{N = 80} were retrieved from
#' \href{https://dlmf.nist.gov/3.5#v}{NIST} and have \eqn{10^{-21}} precision.
#' For \eqn{N = 160} onwards, the nodes and weights were computed with the
#' \code{gauss.quad} function from the \href{https://CRAN.R-project.org/package=statmod}{
#' \code{statmod}} package (Smyth, 1998), and have \eqn{10^{-15}} precision.
#' @references
#' \emph{NIST Digital Library of Mathematical Functions}. Release
#' 1.0.20 of 2018-09-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier,
#' B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller,
#' and B. V. Saunders, eds. \url{https://dlmf.nist.gov/}
#'
#' Smyth, G. K. (1998). Numerical integration. In: \emph{Encyclopedia of
#' Biostatistics}, P. Armitage and T. Colton (eds.), Wiley, London,
#' pp. 3088-3095.
#' @examples
#' ## Integration of a smooth function in (1, 10)
#'
#' # Weights and nodes for integrating
#' x_k <- Gauss_Legen_nodes(a = 1, b = 10, N = 40)
#' w_k <- Gauss_Legen_weights(a = 1, b = 10, N = 40)
#'
#' # Check quadrature
#' f <- function(x) sin(x) * x^2 - log(x + 1)
#' integrate(f, lower = 1, upper = 10, rel.tol = 1e-12)
#' sum(w_k * f(x_k))
#'
#' # Exact for polynomials up to degree 2 * N - 1
#' f <- function(x) (((x + 0.5) / 1e3)^5 - ((x - 0.5)/ 5)^4 +
#'   ((x - 0.25) / 10)^2 + 1)^20
#' sum(w_k * f(x_k))
#' integrate(f, lower = -1, upper = 1, rel.tol = 1e-12)
#'
#' ## Integration on (0, pi)
#'
#' # Weights and nodes for integrating
#' th_k <- Gauss_Legen_nodes(a = 0, b = pi, N = 40)
#' w_k <- Gauss_Legen_weights(a = 0, b = pi, N = 40)
#'
#' # Check quadrature
#' p <- 4
#' psi <- function(th) -sin(th / 2)
#' w <- function(th) sin(th)^(p - 2)
#' integrate(function(th) psi(th) * w(th), lower = 0, upper = pi,
#'           rel.tol = 1e-12)
#' sum(w_k * psi(th_k) * w(th_k))
#'
#' # Integral with Gegenbauer polynomial
#' k <- 3
#' C_k <- function(th) drop(Gegen_polyn(theta = th, k = k, p = p))
#' integrate(function(th) psi(th) * C_k(th) * w(th), lower = 0, upper = pi,
#'           rel.tol = 1e-12)
#' th_k <- drop(Gauss_Legen_nodes(a = 0, b = pi, N = 80))
#' w_k <- drop(Gauss_Legen_weights(a = 0, b = pi, N = 80))
#' sum(w_k * psi(th_k) * C_k(th_k) * w(th_k))
#' @name Gauss_Legen
NULL

#' @rdname Gauss_Legen
#' @export
Gauss_Legen_nodes <- function(a = -1, b = 1, N = 40L) {
    .Call('_sphunif_Gauss_Legen_nodes', PACKAGE = 'sphunif', a, b, N)
}

#' @rdname Gauss_Legen
#' @export
Gauss_Legen_weights <- function(a = -1, b = 1, N = 40L) {
    .Call('_sphunif_Gauss_Legen_weights', PACKAGE = 'sphunif', a, b, N)
}

#' @rdname sph_stat_distr
#' @export
p_sph_stat_Bingham <- function(x, p) {
    .Call('_sphunif_p_sph_stat_Bingham', PACKAGE = 'sphunif', x, p)
}

#' @rdname sph_stat_distr
#' @export
d_sph_stat_Bingham <- function(x, p) {
    .Call('_sphunif_d_sph_stat_Bingham', PACKAGE = 'sphunif', x, p)
}

#' @rdname sph_stat_distr
#' @export
p_sph_stat_CJ12 <- function(x, regime = 1L, beta = 0) {
    .Call('_sphunif_p_sph_stat_CJ12', PACKAGE = 'sphunif', x, regime, beta)
}

#' @rdname sph_stat_distr
#' @export
d_sph_stat_CJ12 <- function(x, regime = 3L, beta = 0) {
    .Call('_sphunif_d_sph_stat_CJ12', PACKAGE = 'sphunif', x, regime, beta)
}

#' @rdname sph_stat_distr
#' @export
p_sph_stat_Rayleigh <- function(x, p) {
    .Call('_sphunif_p_sph_stat_Rayleigh', PACKAGE = 'sphunif', x, p)
}

#' @rdname sph_stat_distr
#' @export
d_sph_stat_Rayleigh <- function(x, p) {
    .Call('_sphunif_d_sph_stat_Rayleigh', PACKAGE = 'sphunif', x, p)
}

#' @rdname sph_stat_distr
#' @export
p_sph_stat_Rayleigh_HD <- function(x, p) {
    .Call('_sphunif_p_sph_stat_Rayleigh_HD', PACKAGE = 'sphunif', x, p)
}

#' @rdname sph_stat_distr
#' @export
d_sph_stat_Rayleigh_HD <- function(x, p) {
    .Call('_sphunif_d_sph_stat_Rayleigh_HD', PACKAGE = 'sphunif', x, p)
}

#' @rdname sph_stat
#' @export
sph_stat_Rayleigh <- function(X) {
    .Call('_sphunif_sph_stat_Rayleigh', PACKAGE = 'sphunif', X)
}

#' @rdname sph_stat
#' @export
sph_stat_Bingham <- function(X) {
    .Call('_sphunif_sph_stat_Bingham', PACKAGE = 'sphunif', X)
}

#' @rdname sph_stat
#' @export
sph_stat_Ajne <- function(X, Psi_in_X = FALSE) {
    .Call('_sphunif_sph_stat_Ajne', PACKAGE = 'sphunif', X, Psi_in_X)
}

#' @rdname sph_stat
#' @export
sph_stat_Gine_Gn <- function(X, Psi_in_X = FALSE, p = 0L) {
    .Call('_sphunif_sph_stat_Gine_Gn', PACKAGE = 'sphunif', X, Psi_in_X, p)
}

#' @keywords internal
sph_stat_Gine_Gn_Psi <- function(Psi, n, p) {
    .Call('_sphunif_sph_stat_Gine_Gn_Psi', PACKAGE = 'sphunif', Psi, n, p)
}

#' @rdname sph_stat
#' @export
sph_stat_Gine_Fn <- function(X, Psi_in_X = FALSE, p = 0L) {
    .Call('_sphunif_sph_stat_Gine_Fn', PACKAGE = 'sphunif', X, Psi_in_X, p)
}

#' @keywords internal
sph_stat_Gine_Fn_Psi <- function(Psi, n, p) {
    .Call('_sphunif_sph_stat_Gine_Fn_Psi', PACKAGE = 'sphunif', Psi, n, p)
}

#' @rdname sph_stat
#' @export
sph_stat_Pycke <- function(X, Psi_in_X = FALSE, p = 0L) {
    .Call('_sphunif_sph_stat_Pycke', PACKAGE = 'sphunif', X, Psi_in_X, p)
}

#' @keywords internal
sph_stat_Pycke_Psi <- function(Psi, n, p) {
    .Call('_sphunif_sph_stat_Pycke_Psi', PACKAGE = 'sphunif', Psi, n, p)
}

#' @rdname sph_stat
#' @export
sph_stat_Bakshaev <- function(X, Psi_in_X = FALSE, p = 0L) {
    .Call('_sphunif_sph_stat_Bakshaev', PACKAGE = 'sphunif', X, Psi_in_X, p)
}

#' @rdname sph_stat
#' @export
sph_stat_Riesz <- function(X, Psi_in_X = FALSE, p = 0L, s = 1.0) {
    .Call('_sphunif_sph_stat_Riesz', PACKAGE = 'sphunif', X, Psi_in_X, p, s)
}

#' @keywords internal
sph_stat_Riesz_Psi <- function(Psi, n, s) {
    .Call('_sphunif_sph_stat_Riesz_Psi', PACKAGE = 'sphunif', Psi, n, s)
}

#' @rdname sph_stat
#' @export
sph_stat_PCvM <- function(X, Psi_in_X = FALSE, p = 0L, N = 160L, L = 1e3L) {
    .Call('_sphunif_sph_stat_PCvM', PACKAGE = 'sphunif', X, Psi_in_X, p, N, L)
}

#' @keywords internal
sph_stat_PCvM_Psi <- function(Psi, n, p, th_grid, int_grid) {
    .Call('_sphunif_sph_stat_PCvM_Psi', PACKAGE = 'sphunif', Psi, n, p, th_grid, int_grid)
}

#' @rdname sph_stat
#' @export
sph_stat_PRt <- function(X, t = 1.0 / 3.0, Psi_in_X = FALSE, p = 0L, N = 160L, L = 1e3L) {
    .Call('_sphunif_sph_stat_PRt', PACKAGE = 'sphunif', X, t, Psi_in_X, p, N, L)
}

#' @keywords internal
sph_stat_PRt_Psi <- function(Psi, t_m, theta_t_m, n, p, th_grid, int_grid) {
    .Call('_sphunif_sph_stat_PRt_Psi', PACKAGE = 'sphunif', Psi, t_m, theta_t_m, n, p, th_grid, int_grid)
}

#' @rdname sph_stat
#' @export
sph_stat_PAD <- function(X, Psi_in_X = FALSE, p = 0L, N = 160L, L = 1e3L) {
    .Call('_sphunif_sph_stat_PAD', PACKAGE = 'sphunif', X, Psi_in_X, p, N, L)
}

#' @keywords internal
sph_stat_PAD_Psi <- function(Psi, n, p, th_grid, int_grid) {
    .Call('_sphunif_sph_stat_PAD_Psi', PACKAGE = 'sphunif', Psi, n, p, th_grid, int_grid)
}

#' @rdname sph_stat
#' @export
sph_stat_CCF09 <- function(X, dirs, K_CCF09 = 25L, original = FALSE) {
    .Call('_sphunif_sph_stat_CCF09', PACKAGE = 'sphunif', X, dirs, K_CCF09, original)
}

#' @rdname sph_stat
#' @export
sph_stat_Rayleigh_HD <- function(X) {
    .Call('_sphunif_sph_stat_Rayleigh_HD', PACKAGE = 'sphunif', X)
}

#' @rdname sph_stat
#' @export
sph_stat_CJ12 <- function(X, regime = 3L, Psi_in_X = FALSE, p = 0L) {
    .Call('_sphunif_sph_stat_CJ12', PACKAGE = 'sphunif', X, regime, Psi_in_X, p)
}

#' @keywords internal
sph_stat_CJ12_Psi <- function(Psi, n, p) {
    .Call('_sphunif_sph_stat_CJ12_Psi', PACKAGE = 'sphunif', Psi, n, p)
}