### ----------------------------------------------------
### -- Math::Libgsl::Function
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-function.rakutest
### ----------------------------------------------------
### -- Chunk 4 of 5
F₁(aᵣ+iaᵢ,aR-iaI,c,x) with error check';
  is-approx $val, 3.352857095662929028, :rel-tol(10⁻¹⁵), '₂F₁(aᵣ+iaᵢ,aR-iaI,c,x) with error check: value';
  is-approx F21-renorm(1, 1, 1, .5), 2, :rel-tol(10⁻¹⁵), '₂F₁(a,b,c,x)/Γ(c)';
  lives-ok { ($val, $err) = F21-renorm-e(1, 1, 1, .5) }, '₂F₁(a,b,c,x)/Γ(c) with error check';
  is-approx $val, 2, :rel-tol(10⁻¹⁵), '₂F₁(a,b,c,x)/Γ(c) with error check: value';
  is-approx F21-conj-renorm(9, 9, -1.5, -.99), 0.10834020229476124874, :rel-tol(10⁻¹⁵), '₂F₁(aᵣ+iaᵢ,aR-iaI,c,x)/Γ(c)';
  lives-ok { ($val, $err) = F21-conj-renorm-e(9, 9, -1.5, -.99) }, '₂F₁(aᵣ+iaᵢ,aR-iaI,c,x)/Γ(c) with error check';
  is-approx $val, 0.10834020229476124874, :rel-tol(10⁻¹⁵), '₂F₁(aᵣ+iaᵢ,aR-iaI,c,x)/Γ(c) with error check: value';
  is-approx F20(.1, .5, -.02), .99901595171179281891589794, :rel-tol(10⁻¹⁵), '₂F₀(a,b,x)';
  lives-ok { ($val, $err) = F20-e(.1, .5, -.02) }, '₂F₀(a,b,x) with error check';
  is-approx $val, .99901595171179281891589794, :rel-tol(10⁻¹⁵), '₂F₀(a,b,x) with error check: value';
}, 'Hypergeometric functions';

subtest {
  my ($val, $err);
  is-approx laguerre1(1, 1), 1, :rel-tol(10⁻¹⁵), 'Lʲ₁(x)';
  lives-ok { ($val, $err) = laguerre1-e(1, 1) }, 'Lʲ₁(x) with error check';
  is-approx $val, 1, :rel-tol(10⁻¹⁵), 'Lʲ₁(x) with error check: value';
  is-approx laguerre2(1, 1), .5, :rel-tol(10⁻¹⁵), 'Lʲ₂(x)';
  lives-ok { ($val, $err) = laguerre2-e(1, 1) }, 'Lʲ₂(x) with error check';
  is-approx $val, .5, :rel-tol(10⁻¹⁵), 'Lʲ₂(x) with error check: value';
  is-approx laguerre3(2, 1), 2.3333333333333333333, :rel-tol(10⁻¹⁵), 'Lʲ₃(x)';
  lives-ok { ($val, $err) = laguerre3-e(2, 1) }, 'Lʲ₃(x) with error check';
  is-approx $val, 2.3333333333333333333, :rel-tol(10⁻¹⁵), 'Lʲ₃(x) with error check: value';
  is-approx laguerre-n(1, .5, 1), .5, :rel-tol(10⁻¹⁵), 'Lʲₙ(x)';
  lives-ok { ($val, $err) = laguerre-n-e(1, .5, 1) }, 'Lʲₙ(x) with error check';
  is-approx $val, .5, :rel-tol(10⁻¹⁵), 'Lʲₙ(x) with error check: value';
}, 'Laguerre functions';

subtest {
  my ($val, $err);
  is-approx lambert-W0(0), 0, :rel-tol(10⁻¹⁵), 'W₀(x)';
  lives-ok { ($val, $err) = lambert-W0-e(0) }, 'W₀(x) with error check';
  is-approx $val, 0, :rel-tol(10⁻¹⁵), 'W₀(x) with error check: value';
  is-approx lambert-Wm1(0), 0, :rel-tol(10⁻¹⁵), 'W₋₁(x)';
  lives-ok { ($val, $err) = lambert-Wm1-e(0) }, 'W₋₁(x) with error check';
  is-approx $val, 0, :rel-tol(10⁻¹⁵), 'W₋₁(x) with error check: value';
}, 'Lambert functions';

subtest {
  my ($val, $err);
  is-approx legendre-P1(0.5), 0.5, :rel-tol(10⁻¹⁵), 'P₁(x)';
  lives-ok { ($val, $err) = legendre-P1-e(.5) }, 'P₁(x) with error check';
  is-approx $val, 0.5, :rel-tol(10⁻¹⁵), 'P₁(x) with error check: value';
  is-approx legendre-P2(.5), -0.125, :rel-tol(10⁻¹⁵), 'P₂(x)';
  lives-ok { ($val, $err) = legendre-P2-e(.5) }, 'P₂(x) with error check';
  is-approx $val, -0.125, :rel-tol(10⁻¹⁵), 'P₂(x) with error check: value';
  is-approx legendre-P3(.5), -0.4375, :rel-tol(10⁻¹⁵), 'P₃(x)';
  lives-ok { ($val, $err) = legendre-P3-e(.5) }, 'P₃(x) with error check';
  is-approx $val, -0.4375, :rel-tol(10⁻¹⁵), 'P₃(x) with error check: value';
  is-approx legendre-Pl(1, .5), .5, :rel-tol(10⁻¹⁵), 'Pₗ(x)';
  lives-ok { ($val, $err) = legendre-Pl-e(1, .5) }, 'Pₗ(x) with error check';
  is-approx $val, .5, :rel-tol(10⁻¹⁵), 'Pₗ(x) with error check: value';
  my @out;
  lives-ok { @out = legendre-Pl-array(10, .5) }, 'Pₗ(x) array';
  is-approx @out[0], 1, :rel-tol(10⁻¹⁵), 'Pₗ(x) array [0]: value';
  my %res;
  lives-ok { %res = legendre-Pl-deriv-array(10, .5) }, 'Pₗ(x) and dPₗ(x)/dx array';
  is-approx %res<out>[0], 1, :rel-tol(10⁻¹⁵), 'Pₗ(x) and dPₗ(x)/dx array Pₗ(x)[0]: value';
  is-approx %res<derout>[0], 0, :rel-tol(10⁻¹⁵), 'Pₗ(x) and dPₗ(x)/dx array dPₗ(x)/dx[0]: value';
  is-approx legendre-Q0(0), 0, :rel-tol(10⁻¹⁵), 'Q₀(x)';
  lives-ok { ($val, $err) = legendre-Q0-e(0) }, 'Q₀(x) with error check';
  is-approx $val, 0, :rel-tol(10⁻¹⁵), 'Q₀(x) with error check: value';
  is-approx legendre-Q1(0), -1, :rel-tol(10⁻¹⁵), 'Q₁(x)';
  lives-ok { ($val, $err) = legendre-Q1-e(0) }, 'Q₁(x) with error check';
  is-approx $val, -1, :rel-tol(10⁻¹⁵), 'Q₁(x) with error check: value';
  is-approx legendre-Ql(10, .5), 0.29165813966586752393, :rel-tol(10⁻¹⁵), 'Qₗ(x)';
  lives-ok { ($val, $err) = legendre-Ql-e(10, .5) }, 'Qₗ(x) with error check';
  is-approx $val, 0.29165813966586752393, :rel-tol(10⁻¹⁵), 'Qₗ(x) with error check: value';
  ok legendre-array-index(3, 3) == 9, 'position in a Legendre array';
  lives-ok { @out = legendre-array(GSL_SF_LEGENDRE_SCHMIDT, 10, .5) }, 'Pₗⁿ(x) array';
  is-approx @out[0], 1, :rel-tol(10⁻¹⁵), 'Pₗⁿ(x) array: value';
  lives-ok { @out = legendre-array-e(GSL_SF_LEGENDRE_SCHMIDT, 10, .5, 1) }, 'Pₗⁿ(x) array with error check';
  is-approx @out[legendre-array-index(3, 3)], 0.513489897661093, :rel-tol(10⁻¹⁵), 'Pₗⁿ(x) array with error check: value';
  dies-ok { @out = legendre-array-e(GSL_SF_LEGENDRE_SCHMIDT, 10, .5, 3) }, 'Pₗⁿ(x) array with error check dies with csphase ≠ 1|-1';
  lives-ok { %res = legendre-deriv-array(GSL_SF_LEGENDRE_SCHMIDT, 10, .15) }, 'dPₗⁿ(x)/dx array';
  is-approx %res<derout>[legendre-array-index(1, 0)], 1, :rel-tol(10⁻¹⁵), 'dPₗⁿ(x)/dx array: value';
  lives-ok { %res = legendre-deriv-array-e(GSL_SF_LEGENDRE_SCHMIDT, 10, .15, 1) }, 'dPₗⁿ(x)/dx array with error check';
  is-approx %res<derout>[legendre-array-index(1, 0)], 1, :rel-tol(10⁻¹⁵), 'dPₗⁿ(x)/dx array with error check: value';
  lives-ok { %res = legendre-deriv-alt-array(GSL_SF_LEGENDRE_SCHMIDT, 10, .15) }, 'alt dPₗⁿ(x)/dx array';
  is-approx %res<derout>[legendre-array-index(1, 0)], -0.9886859966642594, :rel-tol(10⁻¹⁵), 'alt dPₗⁿ(x)/dx array: value';
  lives-ok { %res = legendre-deriv-alt-array-e(GSL_SF_LEGENDRE_SCHMIDT, 10, .15, 1) }, 'alt dPₗⁿ(x)/dx array with error check';
  is-approx %res<derout>[legendre-array-index(1, 0)], -0.9886859966642594, :rel-tol(10⁻¹⁵), 'alt dPₗⁿ(x)/dx array with error check: value';
  lives-ok { %res = legendre-deriv2-array(GSL_SF_LEGENDRE_SCHMIDT, 10, .15) }, 'dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array';
  is-approx %res<derout2>[legendre-array-index(1, 0)], 0, 10⁻¹⁵, 'dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array: value';
  lives-ok { %res = legendre-deriv2-array-e(GSL_SF_LEGENDRE_SCHMIDT, 10, .15, 1) }, 'dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array with error check';
  is-approx %res<derout2>[legendre-array-index(1, 0)], 0, 10⁻¹⁵, 'dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array with error check: value';
  lives-ok { %res = legendre-deriv2-alt-array(GSL_SF_LEGENDRE_SCHMIDT, 10, .15) }, 'alt dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array';
  is-approx %res<derout2>[legendre-array-index(1, 0)], -.15, 10⁻¹⁵, 'alt dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array: value';
  lives-ok { %res = legendre-deriv2-alt-array-e(GSL_SF_LEGENDRE_SCHMIDT, 10, .15, 1) }, 'alt dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array with error check';
  is-approx %res<derout2>[legendre-array-index(1, 0)], -.15, 10⁻¹⁵, 'alt dPₗⁿ(x)/dx and d²Pₗⁿ(x)/dx² array with error check: value';
  is-approx legendre-Plm(10, 0, -.5), -0.18822860717773437500, :rel-tol(10⁻¹⁵), 'Pₗⁿ(x)';
  dies-ok { legendre-Plm(10, 20, -.5) }, 'Pₗⁿ(x) dies when l < m';
  lives-ok { ($val, $err) =legendre-Plm-e(10, 0, -.5) }, 'Pₗⁿ(x) with error check';
  is-approx $val, -0.18822860717773437500, :rel-tol(10⁻¹⁵), 'Pₗⁿ(x) with error check: value';
  is-approx legendre-sphPlm(10, 0, -.5e0), -0.24332702369300135, :rel-tol(10⁻¹⁵), 'spherical harmonics';
  lives-ok { ($val, $err) = legendre-sphPlm-e(10, 0, -.5) }, 'spherical harmonics with error check';
  is-approx $val, -0.24332702369300135, :rel-tol(10⁻¹⁵), 'spherical harmonics with error check: value';
  is-approx conicalP-half(0, .5), 0.8573827581049917129, :rel-tol(10⁻¹⁵), 'irregular Spherical Conical Function';
  lives-ok { ($val, $err) = conicalP-half-e(0, .5) }, 'irregular Spherical Conical Function with error check';
  is-approx $val, 0.8573827581049917129, :rel-tol(10⁻¹⁵), 'irregular Spherical Conical Function with error check: value';
  is-approx conicalP-mhalf(0, .5), 0.8978491247257322404, :rel-tol(10⁻¹⁵), 'regular Spherical Conical Function';
  lives-ok { ($val, $err) = conicalP-mhalf-e(0, .5) }, 'regular Spherical Conical Function with error check';
  is-approx $val, 0.8978491247257322404, :rel-tol(10⁻¹⁵), 'regular Spherical Conical Function with error check: value';
  is-approx conicalP0(0, .5), 1.0731820071493643751, :rel-tol(10⁻¹⁵), 'P⁰₋₁ ₂₊ᵢₗ(x) Function';
  lives-ok { ($val, $err) = conicalP0-e(0, .5) }, 'P⁰₋₁ ₂₊ᵢₗ(x) Function with error check';
  is-approx $val, 1.0731820071493643751, :rel-tol(10⁻¹⁵), 'P⁰₋₁ ₂₊ᵢₗ(x) Function with error check: value';
  is-approx conicalP1(0, .5), 0.14933621085538265636, :rel-tol(10⁻¹⁵), 'P¹₋₁ ₂₊ᵢₗ(x) Function';
  lives-ok { ($val, $err) = conicalP1-e(0, .5) }, 'P¹₋₁ ₂₊ᵢₗ(x) Function with error check';
  is-approx $val, 0.14933621085538265636, :rel-tol(10⁻¹⁵), 'P¹₋₁ ₂₊ᵢₗ(x) Function with error check: value';
  is-approx conicalP-sph-reg(2, 1, -.5), 1.6406279287008789526, :rel-tol(10⁻¹⁵), 'Regular Spherical Conical Function';
  lives-ok { ($val, $err) = conicalP-sph-reg-e(2, 1, -.5) }, 'Regular Spherical Conical Function with error check';
  is-approx $val, 1.6406279287008789526, :rel-tol(10⁻¹⁵), 'Regular Spherical Conical Function with error check: value';
  is-approx conicalP-cyl-reg(2, 1, -.5), 2.2048510472375258708, :rel-tol(10⁻¹⁵), 'Regular Cylindrical Conical Function';
  lives-ok { ($val, $err) = conicalP-cyl-reg-e(2, 1, -.5) }, 'Regular Cylindrical Conical Function with error check';
  is-approx $val, 2.2048510472375258708, :rel-tol(10⁻¹⁵), 'Regular Cylindrical Conical Function with error check: value';
  is-approx legendre-H3d0(1, 1), 0.7160229153604338713, :rel-tol(10⁻¹⁵), '0 radial eigenfunction of the Laplacian on 3D hyperbolic space';
  lives-ok { ($val, $err) = legendre-H3d0-e(1, 1) }, '0 radial eigenfunction of the Laplacian on 3D hyperbolic space with error check';
  is-approx $val, 0.7160229153604338713, :rel-tol(10⁻¹⁵), '0 radial eigenfunction of the Laplacian on 3D hyperbolic space with error check: value';
  is-approx legendre-H3d1(1, 1), 0.3397013994799344639, :rel-tol(10⁻¹⁵), '1 radial eigenfunction of the Laplacian on 3D hyperbolic space';
  lives-ok { ($val, $err) = legendre-H3d1-e(1, 1) }, '1 radial eigenfunction of the Laplacian on 3D hyperbolic space with error check';
  is-approx $val, 0.3397013994799344639, :rel-tol(10⁻¹⁵), '1 radial eigenfunction of the Laplacian on 3D hyperbolic space with error check: value';
  is-approx legendre-H3d(5, 1, 1), 0.011498635037491577728, 10⁻¹⁵, 'l radial eigenfunction of the Laplacian on 3D hyperbolic space';
  lives-ok { ($val, $err) = legendre-H3d-e(5, 1, 1) }, 'l radial eigenfunction of the Laplacian on 3D hyperbolic space with error check';
  is-approx $val, 0.011498635037491577728, 10⁻¹⁵, 'l radial eigenfunction of the Laplacian on 3D hyperbolic space with error check: value';
  lives-ok { @out = legendre-H3d-array(10, 1, 3) }, 'l radial eigenfunction of the Laplacian on 3D hyperbolic space array';
  is-approx @out[0], 0.014086820716211607, 10⁻¹⁵, 'l radial eigenfunction of the Laplacian on 3D hyperbolic space array: value';
}, 'Legendre functions';

subtest {
  my ($val, $err);
  is-approx gsl-log(1.1), 0.09531017980432486004, :rel-tol(10⁻¹⁵), 'log(x)';
  lives-ok { ($val, $err) = gsl-log-e(1.1) }, 'log(x) with error check';
  is-approx $val, 0.09531017980432486004, :rel-tol(10⁻¹⁵), 'log(x) with error check: value';
  is-approx gsl-log-abs(-1.1), 0.09531017980432486004, :rel-tol(10⁻¹⁵), 'log(|x|)';
  lives-ok { ($val, $err) = gsl-log-abs-e(-1.1) }, 'log(|x|) with error check';
  is-approx $val, 0.09531017980432486004, :rel-tol(10⁻¹⁵), 'log(|x|) with error check: value';
  my ($lnrval, $lnrerr, $thetaval, $thetaerr);
  lives-ok { ($lnrval, $lnrerr, $thetaval, $thetaerr) = gsl-log-complex-e(1, 1) }, 'log(z) with error check';
  is-approx $lnrval, 0.3465735902799726547, :rel-tol(10⁻¹⁵), 'log(z) with error check: lnr';
  is-approx $thetaval, 0.7853981633974483096, :rel-tol(10⁻¹⁵), 'log(z) with error check: theta';
  is-approx gsl-log1plusx(1e-4), 0.00009999500033330833533, :rel-tol(10⁻¹⁵), 'log(1+x)';
  lives-ok { ($val, $err) = gsl-log1plusx-e(1e-4) }, 'log(1+x) with error check';
  is-approx $val, 0.00009999500033330833533, :rel-tol(10⁻¹⁵), 'log(1+x) with error check: value';
  is-approx gsl-log1plusx-mx(1e-4), -4.999666691664666833e-09, :rel-tol(10⁻¹⁵), 'log(1+x)-x';
  lives-ok { ($val, $err) = gsl-log1plusx-mx-e(1e-4) }, 'log(1+x)-x with error check';
  is-approx $val, -4.999666691664666833e-09, :rel-tol(10⁻¹⁵), 'log(1+x)-x with error check: value';
}, 'logarithms';

subtest {
  my ($val, $err, @out);
  is-approx mathieu-a(2, 1.1), 4.440022293558291, :rel-tol(10⁻¹⁵), 'aₙ(q)';
  lives-ok { ($val, $err) = mathieu-a-e(2, 1.1) }, 'aₙ(q) with error check';
  is-approx $val, 4.440022293558291, :rel-tol(10⁻¹⁵), 'aₙ(q) with error check: value';
  is-approx mathieu-b(2, 1.1e0), 3.8996898851752686, :rel-tol(10⁻¹⁵), 'bₙ(q)';
  lives-ok { ($val, $err) = mathieu-b-e(2, 1.1) }, 'bₙ(q) with error check';
  is-approx $val, 3.8996898851752686, :rel-tol(10⁻¹⁵), 'bₙ(q) with error check: value';
  lives-ok { @out = mathieu-a-array(0, 5, 0) }, 'aₙ(q) array';
  is-deeply @out, [0e0, 1e0, 4e0, 9e0, 16e0, 25e0], 'aₙ(q) array: values';
  lives-ok { @out = mathieu-b-array(0, 5, 0) }, 'bₙ(q) array';
  is-deeply @out, [0e0, 1e0, 4e0, 9e0, 16e0, 25e0], 'bₙ(q) array: values';
  is-approx mathieu-ce(0, 0, 0), 0.7071067811865475, :rel-tol(10⁻¹⁵), 'ceₙ(q,x)';
  lives-ok { ($val, $err) = mathieu-ce-e(0, 0, 0) }, 'ceₙ(q,x) with error check';
  is-approx $val, 0.7071067811865475, :rel-tol(10⁻¹⁵), 'ceₙ(q,x) with error check: value';
  is-approx mathieu-se(15, 0, π/2), -1, :rel-tol(10⁻¹⁵), 'seₙ(q,x)';
  lives-ok { ($val, $err) = mathieu-se-e(15, 0, π/2) }, 'seₙ(q,x) with error check';
  is-approx $val, -1, :rel-tol(10⁻¹⁵), 'seₙ(q,x) with error check: value';
  lives-ok { @out = mathieu-ce-array(0, 5, 0, π/2) }, 'ceₙ(q,x) array';
  is-deeply @out, [0.7071067811865475e0, 6.123233995736766e-17, -1e0, -1.8369701987210297e-16, 1e0, 3.061616997868383e-16], 'ceₙ(q,x) array: values';
  lives-ok { @out = mathieu-se-array(0, 5, 0, π/2) }, 'seₙ(q,x) array';
  is-deeply @out, [0e0, 1e0, 1.2246467991473532e-16, -1e0, -2.4492935982947064e-16, 1e0], 'seₙ(q,x) array: values';
  is-approx mathieu-Mc(1, 0, 1, 0), 0.8447696705150902, :rel-tol(10⁻¹⁵), 'Mcʲₙ(q,x)';
  lives-ok { ($val, $err) = mathieu-Mc-e(1, 0, 1, 0) }, 'Mcʲₙ(q,x) with error check';
  is-approx $val, 0.8447696705150902, :rel-tol(10⁻¹⁵), 'Mcʲₙ(q,x) with error check: value';
  is-approx mathieu-Ms(1, 0, 1, 0), 0, :rel-tol(10⁻¹⁵), 'Msʲₙ(q,x)';
  lives-ok { ($val, $err) = mathieu-Ms-e(1, 0, 1, 0) }, 'Msʲₙ(q,x) with error check';
  is-approx $val, 0, :rel-tol(10⁻¹⁵), 'Msʲₙ(q,x) with error check: value';
  lives-ok { @out = mathieu-Mc-array(1, 0, 5, 1, π/2) }, 'Mcʲₙ(q,x) array';
  is-deeply @out, [-0.3118830825472924e0, -0.4013382117693621e0, 0.12436859192138153e0, 0.4930738160216241e0, 0.47545235504042044e0, 0.2964329279875498e0], 'Mcʲₙ(q,x) array: values';
  lives-ok { @out = mathieu-Ms-array(1, 0, 5, 1, π/2) }, 'Msʲₙ(q,x) array';
  is-deeply @out, [0e0, -0.3397911837218603e0, 0.1452203874466516e0, 0.4932413901187376e0, 0.47541486595587257e0, 0.2964318355087123e0], 'Msʲₙ(q,x) array: values';
}, 'Mathieu functions';

subtest {
  my ($val, $err);
  is-approx pow-int(2, 3), 8, :rel-tol(10⁻¹⁵), 'xⁿ';
  lives-ok { ($val, $err) = pow-int-e(2, 3) }, 'xⁿ with error check';
  is-approx $val, 8, :rel-tol(10⁻¹⁵), 'xⁿ with error check: value';
}, 'Power functions';

subtest {
  my ($val, $err);
  is-approx psi-int(1), -0.57721566490153286060, :rel-tol(10⁻¹⁵), 'ψ(n)';
  lives-ok { ($val, $err) = psi-int-e(1) }, 'ψ(n) with error check';
  is-approx $val, -0.57721566490153286060, :rel-tol(10⁻¹⁵), 'ψ(n) with error check: value';
  is-approx psi(1), -0.57721566490153286060, :rel-tol(10⁻¹⁵), 'ψ(x)';
  lives-ok { ($val, $err) = psi-e(1) }, 'ψ(x) with error check';
  is-approx $val, -0.57721566490153286