### ----------------------------------------------------
### -- Math::Libgsl::Polynomial
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw-polynomial.rakutest
### ----------------------------------------------------
#!/usr/bin/env raku

use Test;
use lib 'lib';
use Math::Libgsl::Raw::Polynomial :ALL;
use Math::Libgsl::Raw::Complex :ALL;
use NativeCall;

subtest {
  my $c = CArray[num64].new: (1, 2, 3)».Num;
  ok gsl_poly_eval($c, $c.elems, 10e0) == 321, 'evaluate real polynomial for real x';
  my $res = CArray[num64].allocate: 3;
  gsl_poly_eval_derivs($c, $c.elems, 10e0, $res, $res.elems);
  is-deeply $res.list, (321e0, 62e0, 6e0), 'evaluate real polynomial and its derivatives for real x';
}, 'raw - evaluate real polynomial';

subtest {
  my gsl_complex $z = alloc_gsl_complex;
  mgsl_complex_rect(1e0, 2e0, $z);
  my $c = CArray[num64].new: (1, 2, 3)».Num;
  my gsl_complex $res = alloc_gsl_complex;
  lives-ok { mgsl_poly_complex_eval($c, $c.elems, $z, $res) }, 'eval real polynomial for complex variable';
  ok $res.dat[0] == -6e0 && $res.dat[1] == 16e0, 'complex result';
  my $z1 = alloc_gsl_complex;
  mgsl_complex_rect(.674e0, -1.423e0, $z1);
  my CArray[gsl_complex] $zc .= new: $z1,;
  mgsl_complex_rect(-1.44e0, 9.55e0, $z);
  lives-ok { mgsl_complex_poly_complex_eval($zc, $zc.elems, $z, $res) }, 'eval complex polynomial for complex variable';
  ok $res.dat[0] == 0.674 && $res.dat[1] == -1.423, 'complex result';
}, 'raw - evaluate complex polynomial';

subtest {
  my $xa = CArray[num64].new: (1e0, 2e0, 3e0, 4e0, 5e0, 6e0);
  my $ya = CArray[num64].new: (1e0, 2e0, 3e0, 2e0, 0e0, 2e0);
  my $dd = CArray[num64].allocate: 6;
  ok gsl_poly_dd_init($dd, $xa, $ya, $xa.elems) == 0, 'gsl_poly_dd_init: no error reported';
  is-deeply $dd.list».round(10⁻¹⁴),
    (1e0, 1e0, 0e0, -0.3333333333333333e0, 0.125e0, 0.008333333333333337e0)».round(10⁻¹⁴),
    'stored divided differences';
  ok gsl_poly_dd_eval($dd, $xa, $xa.elems, 7e0) == 18e0, 'evaluate polynomial in stored divided diferences';
  my $c = CArray[num64].allocate: 6;
  my $w = CArray[num64].allocate: 6;
  ok gsl_poly_dd_taylor($c, 2.5e0, $dd, $xa, $xa.elems, $w) == 0, 'gsl_poly_dd_taylor: no error reported';
  is-deeply $c.list».round(10⁻¹⁴),
    (2.68359375e0, 1.0880208333333332e0, -0.7604166666666666e0, -0.35416666666666663e0, 0.10416666666666666e0,
     0.008333333333333337e0)».round(10⁻¹⁴),
    'Taylor expansion';
  my $xxa  = CArray[num64].new: (1.3e0, 1.6e0, 1.9e0);
  my $yya  = CArray[num64].new: (0.6200860e0, 0.4554022e0, 0.2818186e0);
  my $dya  = CArray[num64].new: (-0.5220232e0, -0.5698959e0, -0.5811571e0);
  my $zza  = CArray[num64].allocate: 6;
  my $ddd  = CArray[num64].allocate: 6;
  ok gsl_poly_dd_hermite_init($ddd, $zza, $xxa, $yya, $dya, $xxa.elems) == 0,
    'gsl_poly_dd_hermite_init: no error reported';
  is-approx gsl_poly_dd_eval($ddd, $zza, $zza.elems, 7e0), -4.5079, 10⁻⁴,
    'evaluate the Hermite polynomial stored as divided diferences';
}, 'raw - divided difference';

subtest {
  my num64 $a = 1e0;
  my num64 $b = -2e0;
  my num64 $c = 1e0;
  my num64 $x0;
  my num64 $x1;
  my $ret = gsl_poly_solve_quadratic($a, $b, $c, $x0, $x1);
  ok $ret == 2, 'two roots found';
  ok $x0  == 1, 'first root';
  ok $x1  == 1, 'second root';
  $b = 1e0;
  my gsl_complex $z0 = alloc_gsl_complex;
  my gsl_complex $z1 = alloc_gsl_complex;
  $ret = gsl_poly_complex_solve_quadratic($a, $b, $c, $z0, $z1);
  ok $ret == 2, 'two roots found';
  ok $z0.dat[0] == -0.5, 'first root, re';
  is-approx $z0.dat[1], -0.8660254037844386, 10⁻¹⁴, 'first root, im';
  ok $z1.dat[0] == -0.5, 'second root, re';
  is-approx $z1.dat[1], 0.8660254037844386, 10⁻¹⁴, 'second root, im';
}, 'raw - quadratic equation';

subtest {
  my num64 $a = 1e0;
  my num64 $b = -2e0;
  my num64 $c = 1e0;
  my num64 $x0;
  my num64 $x1;
  my num64 $x2;
  my $ret = gsl_poly_solve_cubic($a, $b, $c, $x0, $x1, $x2);
  ok $ret == 1, 'one root found';
  is-approx $x0, -2.1478990357047874, 10⁻¹⁴, 'real root';
  my gsl_complex $z0 = alloc_gsl_complex;
  my gsl_complex $z1 = alloc_gsl_complex;
  my gsl_complex $z2 = alloc_gsl_complex;
  $ret = gsl_poly_complex_solve_cubic($a, $b, $c, $z0, $z1, $z2);
  ok $ret == 3, 'three roots found';
  is-approx $z0.dat[0], -2.1478990357047874, 10⁻¹⁴, 'first root, re';
  ok $z0.dat[1].round(10⁻¹⁴) == 0, 'first root, im';
  is-approx $z1.dat[0],  0.5739495178523937, 10⁻¹⁴, 'second root, re';
  is-approx $z1.dat[1], -0.3689894074818041, 10⁻¹⁴, 'second root, im';
  is-approx $z2.dat[0],  0.5739495178523937, 10⁻¹⁴, 'third root, re';
  is-approx $z2.dat[1],  0.3689894074818041, 10⁻¹⁴, 'third root, im';
}, 'raw - cubic equation';

subtest {
  my $a = CArray[num64].new: (-1, 0, 0, 0, 0, 1)».Num;
  my $z = CArray[num64].allocate(($a.elems - 1) * 2);
  my gsl_poly_complex_workspace $w = gsl_poly_complex_workspace_alloc($a.elems);
  isa-ok $w, gsl_poly_complex_workspace, 'allocate workspace';
  my $ret = gsl_poly_complex_solve($a, $a.elems, $w, $z);
  lives-ok { gsl_poly_complex_workspace_free($w) }, 'free workspace';
  ok $ret == 0, 'no error reported';
  is-deeply $z.list».round(10⁻¹⁴),
    (-0.8090169943749477,  0.5877852522924734,
     -0.8090169943749477, -0.5877852522924734,
      0.3090169943749475,  0.951056516295153,
      0.3090169943749475, -0.951056516295153,
      0.9999999999999999,  0)».round(10⁻¹⁴),
    'complex polynomial';
}, 'raw - solve complex polynomial';

done-testing;


