Overhaul special functions modified by GSOC2018 project. * NEWS: Add note about new functions added. Add note explaining the changes done to existing functions. * __betainc__.cc: Renamed from __betainc_lentz__.cc. Use standard GPL v3 copyright block. Add missing #include "dNDArray.h". Add one-line texinfo documentation for internal function. Remove fourth input argument to function. Use is_single_type() to decide whether to operate with FloatNDArray or NDArray. Delete temporary variables x_arg_s, a_arg_s, b_arg_s used in input validation. Rename len_x, len_a, len_b to numel_[xab] for clarity. Remove input validation that numel match (internal function, no need). Rename function output to retval according to Octave coding conventions. Use more spacing and newlines for readability of code. * __expint__.cc: Renamed from __expint_lentz__.cc. Use standard GPL v3 copyright block. Cull list of #includes to just the necessary ones. Use Complex and FloatComplex typedefs defined by Octave. Add one-line texinfo documentation for internal function. Remove second input argument to function. Use is_single_type() to decide whether to operate with FloatComplexNDArray or ComplexNDArray. Delete temporary variables x_arg_s used in input validation. Rename len_x to numel_x for clarity. Use constructor with dim_vector and scalar value rather than fill() after creating array. Rename function output to retval according to Octave coding conventions. Use more spacing and newlines for readability of code. * __gammainc__.cc: Renamed from __gammainc_lentz__.cc. Use standard GPL v3 copyright block. Add one-line texinfo documentation for internal function. Remove third input argument to function. Remove input validation that numel match (internal function, no need). Use is_single_type() to decide whether to operate with FloatNDArray or NDArray. Delete temporary variables x_arg, a_arg used in input validation. Use constructor with dim_vector and scalar value rather than fill() after creating array. Rename function output to retval according to Octave coding conventions. Use more spacing and newlines for readability of code. * __betainc_lentz__.cc, __expint_lentz__.cc, __gammainc_lentz__.cc: Removed. * libinterp/corefcn/module.mk: Add renamed functions to build system. * betainc.m: Use Octave standard GPL block. Rewrite parts of docstring. Don't use array brackets around single output of function. Remove isscalar checks on inputs because common_size() function will already handle it. Use capital variable names in error messages to match documentation as displayed in terminal. Reshape all inputs in to column vectors quickly so that input validation tests that depend on all/any will pass with N-D arrays. Add comments to code. Check for specific error messages in input validation BIST tests. * betaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring. Don't use array brackets around single output of function. Remove isscalar checks on inputs because common_size() function will already handle it. Use capital variable names in error messages to match documentation as displayed in terminal. Reshape all inputs in to column vectors quickly so that input validation tests that depend on all/any will pass with N-D arrays. Put most common case of tail ("lower") first in if/elseif trees. Call functions directly with function handle rather than using unnecessary feval() call. Use numel in preference to length. Rename variable i_miss to todo for clarity. Add comments to code. Check for specific error messages in input validation BIST tests. * cosint.m: Use Octave standard GPL block. Rewrite parts of docstring. Don't use array brackets around single output of function. Remove isscalar checks on inputs because common_size() function will already handle it. Add input validation check for isnumeric value. Convert integer classes to double before proceeding. Rename i_miss to todo for clarity. Use isinf to detect both -Inf and +Inf rather than separate tests. Use ++it in while loop conditional to shorten loop blocks. Add Input validation BIST tests. * expint.m: Remove Sylvain who was not actually an author on this file. Rewrite parts of docstring. Rename variable sparse_x to orig_sparse. Eliminate temporary variables res_tmp, x_s_tmp, ssum_tmp. Rename i_miss to todo. Use Octave coding conventions throughout. Add comments to code. * gammainc.m: Use Octave standard GPL block. Rewrite parts of docstring. Remove isscalar checks on inputs because common_size() function will already handle it. Use capital variable names in error messages to match documentation as displayed in terminal. Reshape all inputs in to column vectors quickly so that input validation tests that depend on all/any will pass with N-D arrays. Put most common case of tail ("lower") first in if/elseif trees. Add input validation of tail. Rename variable ii to idx for clarity. Rename variable i_done to todo and switch polarity so that the '!' operator is not required every time the variable is updated. Use indexing and direct assignment to update todo rather than logical operator '&' which is slower. Use tolower on tail variable and then switch strcmpi calls to strcmp. Reformat %!test blocks in to %!assert blocks to be more compact. Check for specific error messages in input validation BIST tests. * gammaincinv.m: Use Octave standard GPL block. Rewrite parts of docstring. Don't use array brackets around single output of function. Remove isscalar checks on inputs because common_size() function will already handle it. Add input validation check for iscomplex value. Use capital variable names in error messages to match documentation as displayed in terminal. Reshape all inputs in to column vectors quickly so that input validation tests that depend on all/any will pass with N-D arrays. Rename i_miss to todo. Use numel in preference to length. Call functions directly with function handle rather than using unnecessary feval() call.Rename variable i_miss to todo for clarity. Use it++ in while loop conditional to shorten loop blocks. Add comments to code. Check for specific error messages in input validation BIST tests. Add input validation BIST tests for all error messages. * sinint.m: Use Octave standard GPL block. Rewrite parts of docstring. Add input validation for isnumeric. Convert integers to double for calculation. Reshape input to column vector. Rename variable sz to orig_sz for clarity. rename i_miss to todo. Reformat BIST tests to mak them more compact. Add input validation BIST tests.

/*

Copyright (C) 2017 Nir Krakauer
Copyright (C) 2018 Michele Ginesi
Copyright (C) 2018 Stefan Schlögl

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 "fNDArray.h"

DEFUN (__gammainc__, args, , 
       doc: /* -*- texinfo -*-
@deftypefn {} {@var{y} =} __gammainc__ (@var{x}, @var{a})
Continued fraction for incomplete gamma function.
@end deftypefn */)
{
  int nargin = args.length ();

  if (nargin != 2)
    print_usage ();

  bool is_single = args(0).is_single_type () || args(1).is_single_type ();

  // Total number of scenarios: get maximum of length of all vectors
  int numel_x = args(0).numel ();
  int numel_a = args(1).numel ();
  int len = std::max (numel_x, numel_a);

  octave_value_list retval;
  // Initialize output dimension vector
  dim_vector output_dv (len, 1);

  // Lentz's algorithm in two cases: single and double precision
  if (is_single)
    {
      // Initialize output and inputs
      FloatColumnVector output (output_dv);
      FloatNDArray x, a;

      if (numel_x == 1)
        x = FloatNDArray (output_dv, args(0).float_scalar_value ());
      else
        x = args(0).float_array_value ();

      if (numel_a == 1)
        a = FloatNDArray (output_dv, args(1).float_scalar_value ());
      else
        a = args(1).float_array_value ();

      // Initialize variables used in algorithm
      static const float tiny = pow (2, -50);
      static const float eps = std::numeric_limits<float>::epsilon();
      float y, Cj, Dj, bj, aj, Deltaj;
      int j, maxit;
      maxit = 200;

      // Loop over all elements
      for (octave_idx_type i = 0; i < len; ++i)
        {
          // Catch Ctrl+C
          OCTAVE_QUIT;

          // Variable initialization for the current element
          y = tiny;
          Cj = y;
          Dj = 0;
          bj = x(i) - a(i) + 1;
          aj = a(i);
          Deltaj = 0;
          j = 1;

          // Lentz's algorithm
          while ((std::abs ((Deltaj - 1) / y) > eps) && (j < maxit))
            {
              Cj = bj + aj/Cj;
              Dj = 1 / (bj + aj*Dj);
              Deltaj = Cj * Dj;
              y *= Deltaj;
              bj += 2;
              aj = j * (a(i) - j);
              j++;
            }

          output(i) = y;
        }

      retval(0) = output;
    }
  else
    {
      // Initialize output and inputs
      ColumnVector output (output_dv);
      NDArray x, a;

      if (numel_x == 1)
        x = NDArray (output_dv, args(0).scalar_value ());
      else
        x = args(0).array_value ();

      if (numel_a == 1)
        a = NDArray (output_dv, args(1).scalar_value ());
      else
        a = args(1).array_value ();

      // Initialize variables used in algorithm
      static const double tiny = pow (2, -100);
      static const double eps = std::numeric_limits<double>::epsilon();
      double y, Cj, Dj, bj, aj, Deltaj;
      int j, maxit;
      maxit = 200;

      // Loop over all elements
      for (octave_idx_type i = 0; i < len; ++i)
        {
          // Catch Ctrl+C
          OCTAVE_QUIT;

          // Variable initialization for the current element
          y = tiny;
          Cj = y;
          Dj = 0;
          bj = x(i) - a(i) + 1;
          aj = a(i);
          Deltaj = 0;
          j = 1;

          // Lentz's algorithm
          while ((std::abs ((Deltaj - 1) / y) > eps) && (j < maxit))
            {
              Cj = bj + aj/Cj;
              Dj = 1 / (bj + aj*Dj);
              Deltaj = Cj * Dj;
              y *= Deltaj;
              bj += 2;
              aj = j * (a(i) - j);
              j++;
            }

          output(i) = y;
        }

      retval(0) = output;
    }

  return retval;
}