tidy_triangular <- function(.n = 50, .min = 0, .max = 1,
                             .mode = 1/2, .num_sims = 1, .return_tibble = TRUE){
  
  # Arguments
  n <- as.integer(.n)
  num_sims <- as.integer(.num_sims)
  mn <- as.numeric(.min)
  mx <- as.numeric(.max)
  md <- as.numeric(.mode)
  ret_tbl <- as.logical(.return_tibble)
  
  # Checks ----
  if (!is.integer(n) | n < 0) {
    rlang::abort(
      message = "The parameters '.n' must be of class integer. Please pass a whole
            number like 50 or 100. It must be greater than 0.",
      use_cli_format = TRUE
    )
  }
  
  if (!is.integer(num_sims) | num_sims < 0) {
    rlang::abort(
      message = "The parameter `.num_sims' must be of class integer. Please pass a
            whole number like 50 or 100. It must be greater than 0.",
      use_cli_format = TRUE
    )
  }
  
  if (mn > mx){
    rlang::abort(
      message = "The parameters .min and .max must satisfy .min < .max",
      use_cli_format = TRUE
    )
  }
  
  if (md < mn || md > mx){
    rlang::abort(
      message = "The parameters must follow .min <= .mode <= .max",
      use_cli_format = TRUE
    )
  }
  
  # Create a data.table with one row per simulation
  df <- data.table::CJ(sim_number = factor(1:num_sims), x = 1:n)
  
  # Group the data by sim_number and add columns for x and y
  df[, y := EnvStats::rtri(n = .N, min = mn, max = mx, mode = md)]
  
  # Compute the density of the y values and add columns for dx and dy
  df[, c("dx", "dy") := density(y, n = n)[c("x", "y")], by = sim_number]
  
  # Compute the p-values for the y values and add a column for p
  df[, p := EnvStats::ptri(y, min = mn, max = mx, mode = md)]
  
  # Compute the q-values for the p-values and add a column for q
  df[, q := EnvStats::qtri(p, min = mn, max = mx, mode = md)]
  
  if(.return_tibble){
    df <- dplyr::as_tibble(df)
  } else {
    data.table::setkey(df, NULL)
  }
  
  param_grid <- dplyr::tibble(mn, mx, md)
  
  # Attach descriptive attributes to tibble
  attr(df, "distribution_family_type") <- "continuous"
  attr(df, ".min") <- .min
  attr(df, ".max") <- .max
  attr(df, ".mode") <- .mode
  attr(df, ".n") <- .n
  attr(df, ".num_sims") <- .num_sims
  attr(df, ".ret_tbl") <- .return_tibble
  attr(df, "tibble_type") <- "tidy_triangular"
  attr(df, "param_grid") <- param_grid
  attr(df, "param_grid_txt") <- paste0(
    "c(",
    paste(param_grid[, names(param_grid)], collapse = ", "),
    ")"
  )
  attr(df, "dist_with_params") <- paste0(
    "Triangular",
    " ",
    paste0(
      "c(",
      paste(param_grid[, names(param_grid)], collapse = ", "),
      ")"
    )
  )
  
  return(df)
}

library(data.table)
library(EnvStats)

> set.seed(12)
> tidy_triangular()$y
 [1] 0.18622690 0.69815169 0.83062133 0.36700264 0.29098808 0.13018376 0.29898579 0.57671839
 [9] 0.10695266 0.06451677 0.44311240 0.69494309 0.43373290 0.43635547 0.36394943 0.46868663
[17] 0.47833417 0.52078582 0.59114784 0.23738043 0.33042939 0.67429795 0.22119340 0.61910006
[25] 0.33001746 0.36602158 0.50238969 0.30707241 0.46873731 0.59368688 0.34704699 0.76898582
[33] 0.75788059 0.69509285 0.57178546 0.82837183 0.60879125 0.72044879 0.43854558 0.44222729
[41] 0.54639325 0.51547833 0.89425477 0.31436412 0.72091756 0.21797697 0.43654954 0.15146637
[49] 0.27624612 0.66437982
> 
> set.seed(12)
> rtri(50)
 [1] 0.18622690 0.69815169 0.83062133 0.36700264 0.29098808 0.13018376 0.29898579 0.57671839
 [9] 0.10695266 0.06451677 0.44311240 0.69494309 0.43373290 0.43635547 0.36394943 0.46868663
[17] 0.47833417 0.52078582 0.59114784 0.23738043 0.33042939 0.67429795 0.22119340 0.61910006
[25] 0.33001746 0.36602158 0.50238969 0.30707241 0.46873731 0.59368688 0.34704699 0.76898582
[33] 0.75788059 0.69509285 0.57178546 0.82837183 0.60879125 0.72044879 0.43854558 0.44222729
[41] 0.54639325 0.51547833 0.89425477 0.31436412 0.72091756 0.21797697 0.43654954 0.15146637
[49] 0.27624612 0.66437982