%perlcode %{ @EXPORT_OK = qw/
               gsl_root_fsolver_alloc
               gsl_root_fsolver_free
               gsl_root_fsolver_set
               gsl_root_fsolver_iterate
               gsl_root_fsolver_name
               gsl_root_fsolver_root
               gsl_root_fsolver_x_lower
               gsl_root_fsolver_x_upper
               gsl_root_fdfsolver_alloc
               gsl_root_fdfsolver_set
               gsl_root_fdfsolver_iterate
               gsl_root_fdfsolver_free
               gsl_root_fdfsolver_name
               gsl_root_fdfsolver_root
               gsl_root_test_interval
               gsl_root_test_residual
               gsl_root_test_delta
               $gsl_root_fsolver_bisection
               $gsl_root_fsolver_brent
               $gsl_root_fsolver_falsepos
               $gsl_root_fdfsolver_newton
               $gsl_root_fdfsolver_secant
               $gsl_root_fdfsolver_steffenson
             /;
%EXPORT_TAGS = ( all => [ @EXPORT_OK ] );

__END__

=encoding utf8

=head1 NAME

Math::GSL::Roots - Find roots of arbitrary 1-D functions

=head1 SYNOPSIS

    use Math::GSL::Roots qw/:all/;

=head1 DESCRIPTION

=over

=item * C<gsl_root_fsolver_alloc($T)> -

This function returns a pointer to a newly allocated instance of a solver of
type C<$T>. C<$T> must be one of the constant included with this module. If there is
insufficient memory to create the solver then the function returns a null
pointer and the error handler is invoked with an error code of $GSL_ENOMEM.

=item * C<gsl_root_fsolver_free($s)> -

Don't call this function explicitly. It will be called automatically in DESTROY for fsolver.

=item * C<gsl_root_fsolver_set($s, $fspec, $x_lower, $x_upper)> -

This function initializes, or reinitializes, an existing solver C<$s>
to use the function described by C<$fspec> and the initial search interval
[C<$x_lower>, C<$x_upper>]. C<$fspec> may either be

=over

=item *

a I<coderef>, e.g.

  $fspec = sub { ... };

or

  sub f { ... };
  $fspec = \&f;

=item *

an arrayref with elements [ I<$coderef>, I<$params> ]

=back

The coderef is called as

      &$coderef( $x, $params );

and should return the function evaluated at C<$x, $params>. For example, to
find the root of a quadratic with run-time specified coefficients C<3, 2, 22>,

  $f = sub {
             my ( $x, $params ) = @_;
             return $params->[0] + $x * $params->[1] + $x**2 * $params->[2];
           };

  $fspec = [ $f, [ 3, 2, 22 ];

  gsl_root_fsolver_set( $s, $fspec, $x_lower, $x_upper );

If there are no extra parameters, set C<$fspec> to the function to be evaluated:

  $fspec = sub {
             my ( $x ) = shift;
             return  $x + $x**2;
           };

  gsl_root_fsolver_set( $s, $fspec, $x_lower, $x_upper );

Don't apply C<gsl_root_fsolver_set> twice to the same fsolver.
It will cause a memory leak. Instead of this you should create new fsolver.

=item * C<gsl_root_fsolver_iterate($s)> -

This function performs a single iteration of the solver $s. If the iteration
encounters an unexpected problem then an error code will be returned (the
Math::GSL::Errno has to be included),

$GSL_EBADFUNC - The iteration encountered a singular point where the function
or its derivative evaluated to Inf or NaN.

$GSL_EZERODIV - The derivative of the function vanished at the iteration point,
preventing the algorithm from continuing without a division by zero.

=item * C<gsl_root_fsolver_name($s)> -

This function returns the name of the solver use within the $s solver.

=item * C<gsl_root_fsolver_root($s)> -

This function returns the current estimate of the root for the solver $s.

=item * C<gsl_root_fsolver_x_lower($s)> -

This function returns the current lower value of the bracketing interval for the solver $s.

=item * C<gsl_root_fsolver_x_upper($s)> -

This function returns the current lower value of the bracketing interval for the solver $s.

=item * C<gsl_root_fdfsolver_alloc($T)> -

This function returns a pointer to a newly allocated instance of a
derivative-based solver of type $T. If there is insufficient memory to create
the solver then the function returns a null pointer and the error handler is
invoked with an error code of $GSL_ENOMEM.

=item * C<gsl_root_fdfsolver_set($s, $fspec, $root)> -

This function initializes, or reinitializes, an existing fdfsolver C<$s> to use the
function and its derivatives specified by C<$fspec> and the initial guess C<$root.> 

C<$fspec> may either be:


=over

=item *

A hashref with elements C<f>, C<df>, C<fdf>.


=item *

An arrayref with elements C<[ $hashref, $params ]>

where C<$hashref> has elements  C<f>, C<df>, C<fdf>;

=back

The hashref elements are

=over

=item *

C<f>

A coderef returning the value of the function at a given C<x>. It is called as C<&$f($x, $params)>.

=item *

C<df>

A coderef returning the value of the derivative of the function with
respect to C<x>. It is called as C<&$df($x, $params)>.

=item *

C<fdf>

A coderef returning the value of the function and its derivative with
respect to C<x>. It is called as C<&$fdf($x, $params)>.

=back


For example, to find the root of a quadratic with run-time specified coefficients C<3, 2, 22>,

  $fdf = {

      f => sub {
          my ( $x, $params ) = @_;
          return $params->[0] + $x * $params->[1] + $x**2 * $params->[2];
      },

      df => sub {
          my ( $x, $params ) = @_;

          $params->[1] + 2 * $x * $params->[2];
      },

      fdf => sub {
          my ( $x, $params ) = @_;

          return
            $params->[0] + $x * $params->[1] + $x**2 * $params->[2],
            $params->[1] + 2 * $x * $params->[2];
      },
  };

  $fspec = [ $fdf, [ 3, 2, 22 ];

  gsl_root_fdsolver_set( $s, $fspec );

If there are no extra parameters, set C<$fspec> to C<$fdf>:

  $fdf = {

      f => sub {
          my $x = shift;
          return $x + $x**2;
      },

      df => sub {
          my $x = shift;

          1 + 2 * $x;
      },

      fdf => sub {
          my $x = shift;

          return
            $x + $x**2,
            1 + 2 * $x;
      },
  };

  gsl_root_fdfsolver_set( $s, $fdf );

Don't apply C<gsl_root_fdffsolver_set> twice to the same fdfsolver.  It
will cause a memory leak. Instead of this you should create new
fdfsolver.


=item * C<gsl_root_fdfsolver_iterate($s)> -

This function performs a single iteration of the solver $s. If the iteration
encounters an unexpected problem then an error code will be returned (the
Math::GSL::Errno has to be included),

$GSL_EBADFUNC - The iteration encountered a singular point where the function or its derivative evaluated to Inf or NaN.
$GSL_EZERODIV - The derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.

=item * C<gsl_root_fdfsolver_free($s)> -

Don't call this function explicitly. It will be called automatically in DESTROY for fdfsolver.

=item * C<gsl_root_fdfsolver_name($s)> -

This function returns the name of the solver use within the $s solver.

=item * C<gsl_root_fdfsolver_root($s)> -

This function returns the current estimate of the root for the solver $s.

=item * C<gsl_root_test_interval($x_lower, $x_upper, $epsabs, $epsrel)> -

This function tests for the convergence of the interval [$x_lower, $x_upper]
with absolute error epsabs and relative error $epsrel. The test returns
$GSL_SUCCESS if the following condition is achieved,

    |a - b| < epsabs + epsrel min(|a|,|b|)

 when the interval x = [a,b] does not include the origin. If the interval
 includes the origin then \min(|a|,|b|) is replaced by zero (which is the
 minimum value of |x| over the interval). This ensures that the relative error
 is accurately estimated for roots close to the origin.  This condition on the
 interval also implies that any estimate of the root r in the interval
 satisfies the same condition with respect to the true root r^*,

    |r - r^*| < epsabs + epsrel r^*

  assuming that the true root r^* is contained within the interval.

=item * C<gsl_root_test_residual($f, $epsabs)> -

This function tests the residual value $f against the absolute error bound
$epsabs. The test returns $GSL_SUCCESS if the following condition is achieved,

    |$f| < $epsabs

and returns $GSL_CONTINUE otherwise. This criterion is suitable for situations
where the precise location of the root, x, is unimportant provided a value can
be found where the residual, |f(x)|, is small enough.

=item * C<gsl_root_test_delta($x1, $x0, $epsabs, $epsrel)> -

This function tests for the convergence of the sequence ..., $x0, $x1 with
absolute error $epsabs and relative error $epsrel. The test returns
$GSL_SUCCESS if the following condition is achieved,

    |x_1 - x_0| < epsabs + epsrel |x_1|

and returns $GSL_CONTINUE otherwise.

=back

This module also includes the following constants :

=over

=item * C<$gsl_root_fsolver_bisection> -

The bisection algorithm is the simplest method of bracketing the roots of a
function. It is the slowest algorithm provided by the library, with linear
convergence. On each iteration, the interval is bisected and the value of the
function at the midpoint is calculated. The sign of this value is used to
determine which half of the interval does not contain a root. That half is
discarded to give a new, smaller interval containing the root. This procedure
can be continued indefinitely until the interval is sufficiently small. At any
time the current estimate of the root is taken as the midpoint of the interval.

=item * C<$gsl_root_fsolver_brent> -

The Brent-Dekker method (referred to here as Brent's method) combines an
interpolation strategy with the bisection algorithm. This produces a fast
algorithm which is still robust. On each iteration Brent's method approximates
the function using an interpolating curve. On the first iteration this is a
linear interpolation of the two endpoints. For subsequent iterations the
algorithm uses an inverse quadratic fit to the last three points, for higher
accuracy. The intercept of the interpolating curve with the x-axis is taken as
a guess for the root. If it lies within the bounds of the current interval then
the interpolating point is accepted, and used to generate a smaller interval.
If the interpolating point is not accepted then the algorithm falls back to an
ordinary bisection step. The best estimate of the root is taken from the most
recent interpolation or bisection.

=item * C<$gsl_root_fsolver_falsepos> -

The false position algorithm is a method of finding roots based on linear
interpolation. Its convergence is linear, but it is usually faster than
bisection. On each iteration a line is drawn between the endpoints (a,f(a)) and
(b,f(b)) and the point where this line crosses the x-axis taken as a
"midpoint". The value of the function at this point is calculated and its sign
is used to determine which side of the interval does not contain a root. That
side is discarded to give a new, smaller interval containing the root. This
procedure can be continued indefinitely until the interval is sufficiently
small. The best estimate of the root is taken from the linear interpolation of
the interval on the current iteration.

=item * C<$gsl_root_fdfsolver_newton> -

Newton's Method is the standard root-polishing algorithm. The algorithm begins
with an initial guess for the location of the root. On each iteration, a line
tangent to the function f is drawn at that position. The point where this line
crosses the x-axis becomes the new guess. The iteration is defined by the
following sequence, x_{i+1} = x_i - f(x_i)/f'(x_i) Newton's method converges
quadratically for single roots, and linearly for multiple roots.

=item * C<$gsl_root_fdfsolver_secant> -

The secant method is a simplified version of Newton's method which does not
require the computation of the derivative on every step.  On its first
iteration the algorithm begins with Newton's method, using the derivative to
compute a first step,

    x_1 = x_0 - f(x_0)/f'(x_0)

Subsequent iterations avoid the evaluation of the derivative by replacing it
with a numerical estimate, the slope of the line through the previous two
points,

    x_{i+1} = x_i f(x_i) / f'_{est}

where

    f'_{est} = (f(x_i) - f(x_{i-1})/(x_i - x_{i-1})

When the derivative does not change significantly in the vicinity of the root
the secant method gives a useful saving. Asymptotically the secant method is
faster than Newton's method whenever the cost of evaluating the derivative is
more than 0.44 times the cost of evaluating the function itself. As with all
methods of computing a numerical derivative the estimate can suffer from
cancellation errors if the separation of the points becomes too small.

On single roots, the method has a convergence of order (1 + \sqrt 5)/2
(approximately 1.62). It converges linearly for multiple roots.

=item * C<$gsl_root_fdfsolver_steffenson> -

The Steffenson Method provides the fastest convergence of all the routines. It
combines the basic Newton algorithm with an Aitken “delta-squared”
acceleration. If the Newton iterates are x_i then the acceleration procedure
generates a new sequence R_i:

    R_i = x_i - (x_{i+1} - x_i)^2 / (x_{i+2} - 2 x_{i+1} + x_{i})

which converges faster than the original sequence under reasonable conditions.
The new sequence requires three terms before it can produce its first value so
the method returns accelerated values on the second and subsequent iterations.
On the first iteration it returns the ordinary Newton estimate. The Newton
iterate is also returned if the denominator of the acceleration term ever
becomes zero.

As with all acceleration procedures this method can become unstable if the
function is not well-behaved.

=back

For more information about these functions, we refer you to the official GSL
documentation: L<http://www.gnu.org/software/gsl/manual/html_node/>

=head1 AUTHORS

Jonathan "Duke" Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2008-2024 Jonathan "Duke" Leto and Thierry Moisan

This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=cut

%}