### ----------------------------------------------------
### -- Math::Libgsl::RandomDistribution
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw-randist.rakutest
### ----------------------------------------------------
### -- Chunk 1 of 2
#!/usr/bin/env raku

use Test;
use NativeCall;
use NativeHelpers::Blob;
use lib 'lib';
use Math::Libgsl::Raw::RandomDistribution :ALL;
use Math::Libgsl::Raw::Random;
use Math::Libgsl::Raw::Matrix :ALL;
use Math::Libgsl::Raw::BLAS :ALL;
use Math::Libgsl::Raw::LinearAlgebra :ALL;
use Math::Libgsl::Constants;

subtest 'gaussian distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gaussian($r, 3e0), 0.4017558243560277, :rel-tol(10⁻¹²), 'gaussian random variate';
  is-approx gsl_ran_gaussian_pdf(1e0, 3e0), 0.12579440923099774, :rel-tol(10⁻¹²), 'gaussian probability density';
  is-approx gsl_ran_gaussian_ziggurat($r, 3e0), 1.3570084962645521, :rel-tol(10⁻¹²), 'gaussian ziggurat';
  is-approx gsl_ran_gaussian_ratio_method($r, 3e0), 4.571692762640136, :rel-tol(10⁻¹²), 'gaussian ratio method';
  is-approx gsl_ran_ugaussian($r), 0.2712258328755056, :rel-tol(10⁻¹²), 'unit gaussian';
  is-approx gsl_ran_ugaussian_pdf(3e0), 0.0044318484119380075, :rel-tol(10⁻¹²), 'unit gaussian probability density';
  is-approx gsl_ran_ugaussian_ratio_method($r), 1.7149190276802055, :rel-tol(10⁻¹²), 'unit gaussian ratio method';
  is-approx gsl_cdf_gaussian_P(1e0, 3e0), 0.6305586598182363, :rel-tol(10⁻¹²), 'cumulative distribution P(x)';
  is-approx gsl_cdf_gaussian_Q(1e0, 3e0), 0.36944134018176367, :rel-tol(10⁻¹²), 'cumulative distribution Q(x)';
  is-approx gsl_cdf_gaussian_Pinv(0.6305586598182363e0, 3e0), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution P(x)';
  is-approx gsl_cdf_gaussian_Qinv(0.36944134018176367e0, 3e0), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution Q(x)';
  is-approx gsl_cdf_ugaussian_P(1e0), 0.8413447460685429, :rel-tol(10⁻¹²), 'unit cumulative distribution P(x)';
  is-approx gsl_cdf_ugaussian_Q(1e0), 0.15865525393145705, :rel-tol(10⁻¹²), 'unit cumulative distribution Q(x)';
  is-approx gsl_cdf_ugaussian_Pinv(0.8413447460685429e0), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution P(x)';
  is-approx gsl_cdf_ugaussian_Qinv(0.15865525393145705e0), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution Q(x)';
  gsl_rng_free($r);
}

subtest 'gaussian tail distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gaussian_tail($r, 1e0, 3e0), 5.023225218761322, :rel-tol(10⁻¹²), 'gaussian tail';
  is-approx gsl_ran_gaussian_tail_pdf(1e0, 1e0, 3e0), 0.3404990063350988, :rel-tol(10⁻¹²), 'gaussian tail probability density';
  is-approx gsl_ran_ugaussian_tail($r, 1e0), 1.198044133837379, :rel-tol(10⁻¹²), 'unit gaussian tail';
  is-approx gsl_ran_ugaussian_tail_pdf(1e0, 1e0), 1.525135276160981, :rel-tol(10⁻¹²), 'unit gaussian tail probability density';
  gsl_rng_free($r);
}

subtest 'bivariate gaussian distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my num64 ($x, $y) = (0e0, 0e0);
  gsl_ran_bivariate_gaussian($r, 3e0, 2e0, .3e0, $x, $y);
  ok ([&&] (($x, $y) Z≅ (-0.1952914837346669e0, 0.21644212342026375e0))), 'bivariate gaussian';
  is-approx gsl_ran_bivariate_gaussian_pdf(0e0, 0e0, 3e0, 2e0, .3e0), 0.027806618922095617, :rel-tol(10⁻¹²), 'bivariate gaussian probability density';
  gsl_rng_free($r);
}

subtest 'multivariate gaussian distribution' => {
  my $N = 100_000;
  my $d = 2;
  my num64 $t2        = 0e0;
  my num64 $threshold = 0e0;
  my num64 $alpha     = 0.05e0;
  my num64 $pvalue    = 0e0;
  my gsl_vector $mu            = gsl_vector_calloc($d);
  my gsl_matrix $sigma         = gsl_matrix_calloc($d, $d);
  my gsl_matrix $L             = gsl_matrix_calloc($d, $d);
  my gsl_vector $sample        = gsl_vector_calloc($d);
  my gsl_matrix $samples       = gsl_matrix_calloc($N, $d);
  my gsl_vector $mu_hat        = gsl_vector_calloc($d);
  my gsl_matrix $sigma_hat     = gsl_matrix_calloc($d, $d);
  my gsl_vector $mu_hat_ctr    = gsl_vector_calloc($d);
  my gsl_matrix $sigma_hat_inv = gsl_matrix_calloc($d, $d);
  my gsl_vector $tmp           = gsl_vector_calloc($d);
  my gsl_rng    $r             = mgsl_rng_setup(DEFAULT);

  gsl_vector_set($mu, 0, 1e0);
  gsl_vector_set($mu, 1, 2e0);
  gsl_matrix_set($sigma, 0, 0, 4e0);
  gsl_matrix_set($sigma, 0, 1, 2e0);
  gsl_matrix_set($sigma, 1, 0, 2e0);
  gsl_matrix_set($sigma, 1, 1, 3e0);

  gsl_matrix_memcpy($L, $sigma);
  gsl_linalg_cholesky_decomp1($L);

  for ^$N -> $i {
    gsl_ran_multivariate_gaussian($r, $mu, $L, $sample);
    gsl_matrix_set_row($samples, $i, $sample);
  }

  gsl_ran_multivariate_gaussian_mean($samples, $mu_hat);
  gsl_ran_multivariate_gaussian_vcov($samples, $sigma_hat);

  gsl_vector_memcpy($mu_hat_ctr, $mu_hat);
  gsl_vector_sub($mu_hat_ctr, $mu);
  gsl_matrix_memcpy($sigma_hat_inv, $sigma_hat);
  gsl_linalg_cholesky_decomp1($sigma_hat_inv);
  gsl_linalg_cholesky_invert($sigma_hat_inv);
  gsl_blas_dgemv(CblasNoTrans, 1e0, $sigma_hat_inv, $mu_hat_ctr, 0e0, $tmp);
  gsl_blas_ddot($mu_hat_ctr, $tmp, $t2);
  $t2 *= $N;

  my num64 $diff = ($N - $d).Num;
  my num64 $dd   = $d.Num;
  $threshold = ($N - 1e0) * $d / $diff * gsl_cdf_fdist_Pinv(1e0 - $alpha, $dd, $diff);

  ok $t2 > $threshold && $pvalue < $alpha, 'multivariate gaussian, mean, vcov';

  my num64 $obs_res = 0e0;
  my gsl_vector $x  = gsl_vector_calloc($d);

  gsl_ran_multivariate_gaussian_pdf($x, $mu, $L, $obs_res, $tmp);

  is-approx $obs_res, 0.028294217120391e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian with pdf';

  gsl_ran_multivariate_gaussian_log_pdf($x, $mu, $L, $obs_res, $tmp);

  is-approx $obs_res, -3.565097837249263e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian log with pdf';

  gsl_vector_free($mu);
  gsl_matrix_free($sigma);
  gsl_matrix_free($L);
  gsl_vector_free($sample);
  gsl_matrix_free($samples);
  gsl_vector_free($mu_hat);
  gsl_matrix_free($sigma_hat);
  gsl_vector_free($mu_hat_ctr);
  gsl_matrix_free($sigma_hat_inv);
  gsl_vector_free($tmp);
  gsl_rng_free($r);
}

subtest 'exponential distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_exponential($r, 2e0), 16.523156432740787e0, :abs-tol(10⁻¹⁰), 'exponential';
  is-approx gsl_ran_exponential_pdf(2e0, 2e0), 0.18393972058572117e0, :abs-tol(10⁻¹⁰), 'exponential with pdf';
  is-approx gsl_cdf_exponential_P(2e0, 2e0), 0.6321205588285577e0, :abs-tol(10⁻¹⁰), 'P cumulative exponential';
  is-approx gsl_cdf_exponential_Q(2e0, 2e0), 0.36787944117144233e0, :abs-tol(10⁻¹⁰), 'Q cumulative exponential';
  is-approx gsl_cdf_exponential_Pinv(0.6321205588285577e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative exponential';
  is-approx gsl_cdf_exponential_Qinv(0.36787944117144233e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative exponential';
  gsl_rng_free($r);
}

subtest 'Laplace distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_laplace($r, 2.75e0), 0.001420747954609721e0, :abs-tol(10⁻¹⁰), 'laplace';
  is-approx gsl_ran_laplace_pdf(2e0, 2.75e0), 0.08785910567087736e0, :abs-tol(10⁻¹⁰), 'laplace with pdf';
  is-approx gsl_cdf_laplace_P(2e0, 2.75e0), 0.7583874594050872e0, :abs-tol(10⁻¹⁰), 'P cumulative laplace';
  is-approx gsl_cdf_laplace_Q(2e0, 2.75e0), 0.24161254059491272e0, :abs-tol(10⁻¹⁰), 'Q cumulative laplace';
  is-approx gsl_cdf_laplace_Pinv(0.7583874594050872e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative laplace';
  is-approx gsl_cdf_laplace_Qinv(0.24161254059491272e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative laplace';
  gsl_rng_free($r);
}

subtest 'exponential power distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_exppow($r, 3.7e0, 0.3e0), 96.18418180115754e0, :abs-tol(10⁻¹⁰), 'exppow';
  is-approx gsl_ran_exppow_pdf(2e0, 3.7e0, 0.3e0), 0.006353720526045815e0, :abs-tol(10⁻¹⁰), 'exppow with pdf';
  is-approx gsl_cdf_exppow_P(2e0, 3.7e0, 0.3e0), 0.5155820414218526e0, :abs-tol(10⁻¹⁰), 'P cumulative expow';
  is-approx gsl_cdf_exppow_Q(2e0, 3.7e0, 0.3e0), 0.4844179585781474e0, :abs-tol(10⁻¹⁰), 'Q cumulative expow';
  gsl_rng_free($r);
}

subtest 'Cauchy distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_cauchy($r, 2e0), -0.001622639831190868e0, :abs-tol(10⁻¹⁰), 'cauchy';
  is-approx gsl_ran_cauchy_pdf(2e0, 2e0), 0.07957747154594767e0, :abs-tol(10⁻¹⁰), 'cauchy with pdf';
  is-approx gsl_cdf_cauchy_P(2e0, 2e0), 0.75e0, :abs-tol(10⁻¹⁰), 'P cumulative cauchy';
  is-approx gsl_cdf_cauchy_Q(2e0, 2e0), 0.25e0, :abs-tol(10⁻¹⁰), 'Q cumulative cauchy';
  is-approx gsl_cdf_cauchy_Pinv(.75e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative cauchy';
  is-approx gsl_cdf_cauchy_Qinv(.25e0, 2e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative cauchy';
  gsl_rng_free($r);
}

subtest 'Rayleigh distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_rayleigh($r, 1.9e0), 0.04318348873449323e0, :abs-tol(10⁻¹⁰), 'rayleigh';
  is-approx gsl_ran_rayleigh_pdf(2e0, 1.9e0), 0.3183584873991547e0, :abs-tol(10⁻¹⁰), 'rayleigh with pdf';
  is-approx gsl_cdf_rayleigh_P(2e0, 1.9e0), 0.4253629302445258e0, :abs-tol(10⁻¹⁰), 'P cumulative rayleigh';
  is-approx gsl_cdf_rayleigh_Q(2e0, 1.9e0), 0.5746370697554741e0, :abs-tol(10⁻¹⁰), 'Q cumulative rayleigh';
  is-approx gsl_cdf_rayleigh_Pinv(0.4253629302445258e0, 1.9e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative rayleigh';
  is-approx gsl_cdf_rayleigh_Qinv(0.5746370697554741e0, 1.9e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative rayleigh';
  gsl_rng_free($r);
}

subtest 'Rayleigh tail distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_rayleigh_tail($r, 2.7e0, 1.9e0), 2.7003453137884574e0, :abs-tol(10⁻¹⁰), 'rayleigh tail';
  is-approx gsl_ran_rayleigh_tail_pdf(2e0, 2.7e0, 1.9e0), 0e0, :abs-tol(10⁻¹⁰), 'rayleigh tail with pdf';
  gsl_rng_free($r);
}

subtest 'Landau distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_landau($r), 3880.037426254597e0, :abs-tol(10⁻¹⁰), 'landau';
  is-approx gsl_ran_landau_pdf(2e0), 0.10491298949080555e0, :abs-tol(10⁻¹⁰), 'landau with pdf';
  gsl_rng_free($r);
}

subtest 'Levy alpha-Stable distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_levy($r, 5e0, 1e0), 6162.79707164862e0, :abs-tol(10⁻¹⁰), 'levy';
  gsl_rng_free($r);
}

subtest 'Levy skew alpha-Stable distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_levy_skew($r, 5e0, 1e0, 0e0), 6162.79707164862e0, :abs-tol(10⁻¹⁰), 'levy skew';
  gsl_rng_free($r);
}

subtest 'gamma distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gamma($r, 2.5e0, 2.17e0), 7.592395030767021e0, :abs-tol(10⁻¹⁰), 'gamma';
  is-approx gsl_ran_gamma_knuth($r, 2.5e0, 2.17e0), 7.554777177884447e0, :abs-tol(10⁻¹⁰), 'gamma knuth';
  is-approx gsl_ran_gamma_pdf(2e0, 2.5e0, 2.17e0), 0.12203602332361416e0, :abs-tol(10⁻¹⁰), 'gamma with pdf';
  is-approx gsl_cdf_gamma_P(2e0, 2.5e0, 2.17e0), 0.12962769871792748e0, :abs-tol(10⁻¹⁰), 'P cumulative gamma';
  is-approx gsl_cdf_gamma_Q(2e0, 2.5e0, 2.17e0), 0.8703723012820725e0, :abs-tol(10⁻¹⁰), 'Q cumulative gamma';
  is-approx gsl_cdf_gamma_Pinv(0.12962769871792748e0, 2.5e0, 2.17e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative gamma';
  is-approx gsl_cdf_gamma_Qinv(0.8703723012820725e0, 2.5e0, 2.17e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gamma';
  gsl_rng_free($r);
}

subtest 'flat distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_flat($r, 3e0, 4e0), 3.999741748906672e0, :abs-tol(10⁻¹⁰), 'flat';
  is-approx gsl_ran_flat_pdf(2e0, 3e0, 4e0), 0e0, :abs-tol(10⁻¹⁰), 'flat with pdf';
  is-approx gsl_cdf_flat_P(2e0, 3e0, 4e0), 0e0, :abs-tol(10⁻¹⁰), 'P cumulative flat';
  is-approx gsl_cdf_flat_Q(2e0, 3e0, 4e0), 1e0, :abs-tol(10⁻¹⁰), 'Q cumulative flat';
  is-approx gsl_cdf_flat_Pinv(0e0, 3e0, 4e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative flat';
  is-approx gsl_cdf_flat_Qinv(1e0, 3e0, 4e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative flat';
  gsl_rng_free($r);
}

subtest 'lognormal distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_lognormal($r, 2.7e0, 3e0), 12.239990687580857e0, :abs-tol(10⁻¹⁰), 'lognormal';
  is-approx gsl_ran_lognormal_pdf(2e0, 2.7e0, 3e0), 0.053160178764929525e0, :abs-tol(10⁻¹⁰), 'lognormal with pdf';
  is-approx gsl_cdf_lognormal_P(2e0, 2.7e0, 3e0), 0.25176338701440004e0, :abs-tol(10⁻¹⁰), 'P cumulative lognormal';
  is-approx gsl_cdf_lognormal_Q(2e0, 2.7e0, 3e0), 0.7482366129856e0, :abs-tol(10⁻¹⁰), 'Q cumulative lognormal';
  is-approx gsl_cdf_lognormal_Pinv(0.25176338701440004e0, 2.7e0, 3e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative lognormal';
  is-approx gsl_cdf_lognormal_Qinv(0.7482366129856e0, 2.7e0, 3e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative lognormal';
  gsl_rng_free($r);
}

subtest 'chi-squared distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_chisq($r, 13e0), 16.53548783112717e0, :abs-tol(10⁻¹⁰), 'chisq';
  is-approx gsl_ran_chisq_pdf(2e0, 13e0), 0.0006389340989638759e0, :abs-tol(10⁻¹⁰), 'chisq with pdf';
  is-approx gsl_cdf_chisq_P(2e0, 13e0), 0.00022625008465604695e0, :abs-tol(10⁻¹⁰), 'P cumulative chisq';
  is-approx gsl_cdf_chisq_Q(2e0, 13e0), 0.999773749915344e0, :abs-tol(10⁻¹⁰), 'Q cumulative chisq';
  is-approx gsl_cdf_chisq_Pinv(0.00022625008465604695e0, 13e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative chisq';
  is-approx gsl_cdf_chisq_Qinv(0.999773749915344e0, 13e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative chisq';
  gsl_rng_free($r);
}

subtest 'f-distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_fdist($r, 3e0, 4e0), 1.265383001140288e0, :abs-tol(10⁻¹⁰), 'fdist';
  is-approx gsl_ran_fdist_pdf(2e0, 3e0, 4e0), 0.13942740046346694e0, :abs-tol(10⁻¹⁰), 'fdist with pdf';
  is-approx gsl_cdf_fdist_P(2e0, 3e0, 4e0), 0.743612802471824e0, :abs-tol(10⁻¹⁰), 'P cumulative fdist';
  is-approx gsl_cdf_fdist_Q(2e0, 3e0, 4e0), 0.2563871975281759e0, :abs-tol(10⁻¹⁰), 'Q cumulative fdist';
  is-approx gsl_cdf_fdist_Pinv(0.743612802471824e0, 3e0, 4e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative fdist';
  is-approx gsl_cdf_fdist_Qinv(0.2563871975281759e0, 3e0, 4e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative fdist';
  gsl_rng_free($r);
}

subtest 't-distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_tdist($r, 1.75e0), 0.2825086397013947e0, :abs-tol(10⁻¹⁰), 'tdist';
  is-approx gsl_ran_tdist_pdf(2e0, 1.75e0), 0.06778157292044809e0, :abs-tol(10⁻¹⁰), 'tdist with pdf';
  is-approx gsl_cdf_tdist_P(2e0, 1.75e0), 0.8993101012459747e0, :abs-tol(10⁻¹⁰), 'P cumulative tdist';
  is-approx gsl_cdf_tdist_Q(2e0, 1.75e0), 0.1006898987540253e0, :abs-tol(10⁻¹⁰), 'Q cumulative tdist';
  is-approx gsl_cdf_tdist_Pinv(0.8993101012459747e0, 1.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative tdist';
  is-approx gsl_cdf_tdist_Qinv(0.1006898987540253e0, 1.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative tdist';
  gsl_rng_free($r);
}

subtest 'beta distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_beta($r, 2e0, 3e0), 0.45158505658193976e0, :abs-tol(10⁻¹⁰), 'beta';
  is-approx gsl_ran_beta_pdf(2e0, 2e0, 3e0), 0e0, :abs-tol(10⁻¹⁰), 'beta with pdf';
  is-approx gsl_cdf_beta_P(2e0, 2e0, 3e0), 1e0, :abs-tol(10⁻¹⁰), 'P cumulative beta';
  is-approx gsl_cdf_beta_Q(2e0, 2e0, 3e0), 0e0, :abs-tol(10⁻¹⁰), 'Q cumulative beta';
  is-approx gsl_cdf_beta_Pinv(1e0, 2e0, 3e0), 1e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative beta';
  is-approx gsl_cdf_beta_Qinv(0e0, 2e0, 3e0), 1e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative beta';
  gsl_rng_free($r);
}

subtest 'logistic distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_logistic($r, 3.1e0), 25.61009178896598e0, :abs-tol(10⁻¹⁰), 'logistic';
  is-approx gsl_ran_logistic_pdf(2e0, 3.1e0), 0.07280296434630028e0, :abs-tol(10⁻¹⁰), 'logistic with pdf';
  is-approx gsl_cdf_logistic_P(2e0, 3.1e0), 0.6559192436053651e0, :abs-tol(10⁻¹⁰), 'P cumulative logistic';
  is-approx gsl_cdf_logistic_Q(2e0, 3.1e0), 0.344080756394635e0, :abs-tol(10⁻¹⁰), 'Q cumulative logistic';
  is-approx gsl_cdf_logistic_Pinv(0.6559192436053651e0, 3.1e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative logistic';
  is-approx gsl_cdf_logistic_Qinv(0.344080756394635e0, 3.1e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative logistic';
  gsl_rng_free($r);
}

subtest 'Pareto distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_pareto($r, 1.9e0, 2.75e0), 2.7503738581610317e0, :abs-tol(10⁻¹⁰), 'pareto';
  is-approx gsl_ran_pareto_pdf(3e0, 1.9e0, 2.75e0), 0.5368266659600216e0, :abs