// Copyright (C) 2017 Michele Ginesi
//
// This file is part of Octave.
//
// Octave is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 3 of the License, or
// (at your option) any later version.
//
// Octave is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Octave; see the file COPYING.  If not, see
// <http://www.gnu.org/licenses/>.

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include "defun.h"
#include <complex>

DEFUN (__expint_lentz_cmplx__, args, , "Continued fraction for exponential integral function, complex case")
{
  int nargin = args.length ();

  if (nargin != 1)
    print_usage ();
  else
    {
      octave_value x_arg = args(0);
      if (x_arg.is_single_type ())
        {
          const std::complex<float> x = args(0).float_value ();
          static const std::complex<float> tiny = pow (2, -50);
          static const std::complex<float> eps = std::numeric_limits<float>::epsilon();
          std::complex<float> f = tiny;
          std::complex<float> C = f;
          std::complex<float> D = 0;
          std::complex<float> alpha_m = 1;
          std::complex<float> beta_m = x;
          std::complex<float> Delta = 0;
          int m = 1;
          int maxit = 100;
          while((std::abs(Delta - 1) > eps) & (m < maxit))
            {
               D = beta_m + alpha_m * D;
               if (D == 0)
                 D = tiny;
               C = beta_m + alpha_m / C;
               if (C == 0)
                 C = tiny;
               D = 1 / D;
               Delta = C * D;
               f *= Delta;
               alpha_m = floor ((m + 1) / 2);
               if (beta_m == x)
                 beta_m = 1;
               else
                 beta_m = x;
               m++;
             }
           if (! error_state)
            return octave_value (f);
        }
      else
        {
          const std::complex<double> x = args(0).std::double_value ();
          static const std::complex<double> tiny = pow (2, -100);
          static const std::complex<double> eps = std::numeric_limits<double>::epsilon();
          std::complex<double> f = tiny;
          std::complex<double> C = f;
          std::complex<double> D = 0;
          std::complex<double> alpha_m = 1;
          std::complex<double> beta_m = x;
          std::complex<double> Delta = 0;
          int m = 1;
          int maxit = 200;
          while((std::abs(Delta - 1) > eps) & (m < maxit))
            {
               D = beta_m + alpha_m * D;
               if (D == 0)
                 D = tiny;
               C = beta_m + alpha_m / C;
               if (C == 0)
                 C = tiny;
               D = 1 / D;
               Delta = C * D;
               f *= Delta;
               alpha_m = floor ((m + 1) / 2);
               if (beta_m == x)
                 beta_m = 1;
               else
                 beta_m = x;
               m++;
             }
           if (! error_state)
            return octave_value (f);
        }
      }
  return octave_value_list ();
}