# include <float.h>
# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "laguerre_polynomial.h"

/******************************************************************************/

int i4_max ( int i1, int i2 )

/******************************************************************************/
/*
  Purpose:

    I4_MAX returns the maximum of two I4's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Input, int I1, I2, are two integers to be compared.

    Output, int I4_MAX, the larger of I1 and I2.
*/
{
  int value;

  if ( i2 < i1 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
/******************************************************************************/

int i4_min ( int i1, int i2 )

/******************************************************************************/
/*
  Purpose:

    I4_MIN returns the smaller of two I4's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Input, int I1, I2, two integers to be compared.

    Output, int I4_MIN, the smaller of I1 and I2.
*/
{
  int value;

  if ( i1 < i2 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
/******************************************************************************/

void imtqlx ( int n, double d[], double e[], double z[] )

/******************************************************************************/
/*
  Purpose:

    IMTQLX diagonalizes a symmetric tridiagonal matrix.

  Discussion:

    This routine is a slightly modified version of the EISPACK routine to
    perform the implicit QL algorithm on a symmetric tridiagonal matrix.

    The authors thank the authors of EISPACK for permission to use this
    routine.

    It has been modified to produce the product Q' * Z, where Z is an input
    vector and Q is the orthogonal matrix diagonalizing the input matrix.
    The changes consist (essentialy) of applying the orthogonal transformations
    directly to Z as they are generated.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 January 2010

  Author:

    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
    C version by John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

    Roger Martin, James Wilkinson,
    The Implicit QL Algorithm,
    Numerische Mathematik,
    Volume 12, Number 5, December 1968, pages 377-383.

  Parameters:

    Input, int N, the order of the matrix.

    Input/output, double D(N), the diagonal entries of the matrix.
    On output, the information in D has been overwritten.

    Input/output, double E(N), the subdiagonal entries of the
    matrix, in entries E(1) through E(N-1).  On output, the information in
    E has been overwritten.

    Input/output, double Z(N).  On input, a vector.  On output,
    the value of Q' * Z, where Q is the matrix that diagonalizes the
    input symmetric tridiagonal matrix.
*/
{
  double b;
  double c;
  double f;
  double g;
  int i;
  int ii;
  int itn = 30;
  int j;
  int k;
  int l;
  int m;
  int mml;
  double p;
  double prec;
  double r;
  double s;

  prec = DBL_EPSILON;

  if ( n == 1 )
  {
    return;
  }

  e[n-1] = 0.0;

  for ( l = 1; l <= n; l++ )
  {
    j = 0;
    for ( ; ; )
    {
      for ( m = l; m <= n; m++ )
      {
        if ( m == n )
        {
          break;
        }

        if ( fabs ( e[m-1] ) <= prec * ( fabs ( d[m-1] ) + fabs ( d[m] ) ) )
        {
          break;
        }
      }
      p = d[l-1];
      if ( m == l )
      {
        break;
      }
      if ( itn <= j )
      {
        printf ( "\n" );
        printf ( "IMTQLX - Fatal error!\n" );
        printf ( "  Iteration limit exceeded\n" );
        exit ( 1 );
      }
      j = j + 1;
      g = ( d[l] - p ) / ( 2.0 * e[l-1] );
      r =  sqrt ( g * g + 1.0 );
      g = d[m-1] - p + e[l-1] / ( g + fabs ( r ) * r8_sign ( g ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ( ii = 1; ii <= mml; ii++ )
      {
        i = m - ii;
        f = s * e[i-1];
        b = c * e[i-1];

        if ( fabs ( g ) <= fabs ( f ) )
        {
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e[i] = f * r;
          s = 1.0 / r;
          c = c * s;
        }
        else
        {
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e[i] = g * r;
          c = 1.0 / r;
          s = s * c;
        }
        g = d[i] - p;
        r = ( d[i-1] - g ) * s + 2.0 * c * b;
        p = s * r;
        d[i] = g + p;
        g = c * r - b;
        f = z[i];
        z[i] = s * z[i-1] + c * f;
        z[i-1] = c * z[i-1] - s * f;
      }
      d[l-1] = d[l-1] - p;
      e[l-1] = g;
      e[m-1] = 0.0;
    }
  }
/*
  Sorting.
*/
  for ( ii = 2; ii <= m; ii++ )
  {
    i = ii - 1;
    k = i;
    p = d[i-1];

    for ( j = ii; j <= n; j++ )
    {
      if ( d[j-1] < p )
      {
         k = j;
         p = d[j-1];
      }
    }

    if ( k != i )
    {
      d[k-1] = d[i-1];
      d[i-1] = p;
      p = z[i-1];
      z[i-1] = z[k-1];
      z[k-1] = p;
    }
  }
  return;
}
/******************************************************************************/

double *l_exponential_product ( int p, double b )

/******************************************************************************/
/*
  Purpose:

    L_EXPONENTIAL_PRODUCT: exponential product table for L(n,x).

  Discussion:

    Let L(n,x) represent the Laguerre polynomial of degree n.

    For polynomial chaos applications, it is of interest to know the
    value of the integrals of products of exp(B*X) with every possible pair
    of basis functions.  That is, we'd like to form

      Tij = Integral ( 0 <= X < +oo ) exp(b*x) * L(i,x) * L(j,x) * exp (-x) dx

    Because of the exponential factor, the quadrature will not be exact.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Parameters:

    Input, int P, the maximum degree of the polyonomial
    factors.  0 <= P.

    Input, double B, the coefficient of X in the exponential factor.

    Output, double L_EXPONENTIAL_PRODUCT[(P+1)*(P+1)], the table of integrals.
    TABLE(I,J) represents the weighted integral of exp(B*X) * L(i,x) * L(j,x).
*/
{
  int i;
  int j;
  int k;
  double *l_table;
  int order;
  double *table;
  double *w_table;
  double x;
  double *x_table;

  table = ( double * ) malloc ( (p+1)*(p+1) * sizeof ( double ) );

  for ( j = 0; j <= p; j++ )
  {
    for ( i = 0; i <= p; i++ )
    {
      table[i+j*(p+1)] = 0.0;
    }
  }

  order = ( 3 * p + 4 ) / 2;

  x_table = ( double * ) malloc ( order * sizeof ( double ) );
  w_table = ( double * ) malloc ( order * sizeof ( double ) );

  l_quadrature_rule ( order, x_table, w_table );

  for ( k = 0; k < order; k++ )
  {
    x = x_table[k];
    l_table = l_polynomial ( 1, p, x_table+k );
    for ( j = 0; j <= p; j++ )
    {
      for ( i = 0; i <= p; i++ )
      {
        table[i+j*(p+1)] = table[i+j*(p+1)]
          + w_table[k] * exp ( b * x ) * l_table[i] * l_table[j];
      }
    }
    free ( l_table );
  }
  free ( w_table );
  free ( x_table );

  return table;
}
/******************************************************************************/

double l_integral ( int n )

/******************************************************************************/
/*
  Purpose:

    L_INTEGRAL evaluates a monomial integral associated with L(n,x).

  Discussion:

    The integral:

      integral ( 0 <= x < +oo ) x^n * exp ( -x ) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the exponent.
    0 <= N.

    Output, double L_INTEGRAL, the value of the integral.
*/
{
  double value;

  value = r8_factorial ( n );

  return value;
}
/******************************************************************************/

double *l_polynomial ( int m, int n, double x[] )

/******************************************************************************/
/*
  Purpose:

    L_POLYNOMIAL evaluates the Laguerre polynomials L(n,x).

  Differential equation:

    X * Y'' + (1-X) * Y' + N * Y = 0

  First terms:

      1
     -X    +  1
   (  X^2 -  4 X     +  2 ) / 2
   ( -X^3 +  9 X^2 -  18 X    +    6 ) / 6
   (  X^4 - 16 X^3 +  72 X^2 -   96 X +      24 ) / 24
   ( -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -  600 x    +  120 ) / 120
   (  X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720 ) / 720
   ( -X^7 + 49 X^6 - 882 X^5 + 7350 X^4 - 29400 X^3
      + 52920 X^2 - 35280 X + 5040 ) / 5040

  Recursion:

    L(0,X) = 1,
    L(1,X) = 1-X,
    N * L(N,X) = (2*N-1-X) * L(N-1,X) - (N-1) * L(N-2,X)

  Orthogonality:

    Integral ( 0 <= X < +oo ) exp ( - X ) * L(N,X) * L(M,X) dX
    = 0 if N /= M
    = 1 if N == M

  Special values:

    L(N,0) = 1.

  Relations:

    L(N,X) = (-1)^N / N! * exp ( x ) * (d/dx)^n ( exp ( - x ) * x^n )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

  Parameters:

    Input, int M, the number of evaluation points.

    Input, int N, the highest order polynomial to compute.

    Input, double X[M], the evaluation points.

    Output, double L_POLYNOMIAL[M*(N+1)], the function values.
*/
{
  int i;
  int j;
  double *v;

  if ( n < 0 )
  {
    return NULL;
  }

  v = ( double * ) malloc ( m*(n+1) * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    v[i+0*m] = 1.0;
  }

  if ( n == 0 )
  {
    return v;
  }

  for ( i = 0; i < m; i++ )
  {
    v[i+1*m] = 1.0 - x[i];
  }

  for ( j = 2; j <= n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      v[i+j*m] = ( ( ( double ) ( 2 * j - 1 ) - x[i] ) * v[i+(j-1)*m]
                   + ( double ) (   - j + 1 )          * v[i+(j-2)*m] )
                   / ( double ) (     j     );
    }
  }

  return v;
}
/******************************************************************************/

double *l_polynomial_coefficients ( int n )

/******************************************************************************/
/*
  Purpose:

    L_POLYNOMIAL_COEFFICIENTS: coeffs for Laguerre polynomial L(n,x).

  First terms:

    0: 1
    1: 1  -1
    2: 1  -2  1/2
    3: 1  -3  3/2  1/6
    4: 1  -4  4   -2/3  1/24
    5: 1  -5  5   -5/3  5/24  -1/120

  Recursion:

    L(0,X) = ( 1,  0, 0, ..., 0 )
    L(1,X) = ( 1, -1, 0, ..., 0 )
    L(N,X) = (2*N-1-X) * L(N-1,X) - (N-1) * L(N-2,X) / N

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

  Parameters:

    Input, int N, the highest order polynomial to compute.
    Note that polynomials 0 through N will be computed.

    Output, double L_POLYNOMIAL_COEFFICIENTS[(N+1)*(N+1)], the coefficients
    of the Laguerre polynomials of degree 0 through N.
*/
{
  double *c;
  int i;
  int j;

  if ( n < 0 )
  {
    return NULL;
  }

  c = ( double * ) malloc ( (n+1)*(n+1) * sizeof ( double ) );

  for ( i = 0; i <= n; i++ )
  {
    for ( j = 0; j <= n; j++ )
    {
      c[i+j*(n+1)] = 0.0;
    }
  }

  for ( i = 0; i <= n; i++ )
  {
    c[i+0*(n+1)] = 1.0;
  }

  if ( n == 0 )
  {
    return c;
  }

  c[1+1*(n+1)] = -1.0;

  for ( i = 2; i <= n; i++ )
  {
    for ( j = 1; j <= n; j++ )
    {
      c[i+j*(n+1)] = (
          ( double ) ( 2 * i - 1 ) * c[i-1+j*(n+1)]
        + ( double ) (   - i + 1 ) * c[i-2+j*(n+1)]
        -                            c[i-1+(j-1)*(n+1)] )
        / ( double )       i;
    }
  }
  return c;
}
/******************************************************************************/

void l_polynomial_values ( int *n_data, int *n, double *x, double *fx )

/******************************************************************************/
/*
  Purpose:

    L_POLYNOMIAL_VALUES returns some values of the Laguerre polynomial L(n,x).

  Discussion:

    In Mathematica, the function can be evaluated by:

      LaguerreL[n,x]

  Differential equation:

    X * Y'' + (1-X) * Y' + N * Y = 0;

  First terms:

      1
     -X    +  1
   (  X^2 -  4 X     +  2 ) / 2
   ( -X^3 +  9 X^2 -  18 X    +    6 ) / 6
   (  X^4 - 16 X^3 +  72 X^2 -   96 X +      24 ) / 24
   ( -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -  600 x    +  120 ) / 120
   (  X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720 ) / 720
   ( -X^7 + 49 X^6 - 882 X^5 + 7350 X^4 - 29400 X^3
      + 52920 X^2 - 35280 X + 5040 ) / 5040

  Recursion:

    L(0,X) = 1,
    L(1,X) = 1-X,
    N * L(N,X) = (2*N-1-X) * L(N-1,X) - (N-1) * L(N-2,X)

  Orthogonality:

    Integral ( 0 <= X < oo ) exp ( - X ) * L(N,X) * L(M,X) dX
    = 0 if N /= M
    = 1 if N == M

  Special values:

    L(N,0) = 1.

  Relations:

    L(N,X) = (-1)^N / N! * exp ( x ) * (d/dx)^n ( exp ( - x ) * x^n )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Parameters:

    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
    first call.  On each call, the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    Output, int *N, the order of the polynomial.

    Output, double *X, the point where the polynomial is evaluated.

    Output, double *FX, the value of the function.
*/
{
# define N_MAX 17

  static double fx_vec[N_MAX] = {
      0.1000000000000000E+01,
      0.0000000000000000E+00,
     -0.5000000000000000E+00,
     -0.6666666666666667E+00,
     -0.6250000000000000E+00,
     -0.4666666666666667E+00,
     -0.2569444444444444E+00,
     -0.4047619047619048E-01,
      0.1539930555555556E+00,
      0.3097442680776014E+00,
      0.4189459325396825E+00,
      0.4801341790925124E+00,
      0.4962122235082305E+00,
     -0.4455729166666667E+00,
      0.8500000000000000E+00,
     -0.3166666666666667E+01,
      0.3433333333333333E+02  };

  static int n_vec[N_MAX] = {
     0,  1,  2,
     3,  4,  5,
     6,  7,  8,
     9, 10, 11,
    12,  5,  5,
     5,  5 };

  static double x_vec[N_MAX] = {
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     0.5E+00,
     3.0E+00,
     5.0E+00,
     1.0E+01 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

double *l_polynomial_zeros ( int n )

/******************************************************************************/
/*
  Purpose:

    L_POLYNOMIAL_ZEROS: zeros of the Laguerre polynomial L(n,x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int N, the order of the polynomial.

    Output, double L_POLYNOMIAL_ZEROS[N], the zeros.
*/
{
  double *bj;
  int i;
  double *w;
  double *x;
  double zemu;
/*
  Define the zero-th moment.
*/
  zemu = 1.0;
/*
  Define the Jacobi matrix.
*/
  bj = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    bj[i] = ( double ) ( i + 1 );
  }

  x = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    x[i] = ( double ) ( 2 * i + 1 );
  }

  w = ( double * ) malloc ( n * sizeof ( double ) );
  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
/*
  Diagonalize the Jacobi matrix.
*/
  imtqlx ( n, x, bj, w );

  free ( bj );
  free ( w );

  return x;
}
/******************************************************************************/

double *l_power_product ( int p, int e )

/******************************************************************************/
/*
  Purpose:

    L_POWER_PRODUCT: power product table for L(n,x).

  Discussion:

    Let L(n,x) represent the Laguerre polynomial of degree n.

    For polynomial chaos applications, it is of interest to know the
    value of the integrals of products of X^E with every possible pair
    of basis functions.  That is, we'd like to form

      Tij = Integral ( 0 <= X < +oo ) x^e * L(i,x) * L(j,x) * exp (-x) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Parameters:

    Input, int P, the maximum degree of the polyonomial
    factors.  0 <= P.

    Input, int E, the exponent of X.
    0 <= E.

    Output, double L_POWER_PRODUCT[(P+1)*(P+1)], the table of integrals.
    TABLE(I,J) represents the weighted integral of X^E * L(i,x) * L(j,x).
*/
{
  int i;
  int j;
  int k;
  double *l_table;
  int order;
  double *table;
  double *w_table;
  double x;
  double *x_table;

  table = ( double * ) malloc ( (p+1)*(p+1) * sizeof ( double ) );

  for ( j = 0; j <= p; j++ )
  {
    for ( i = 0; i <= p; i++ )
    {
      table[i+j*(p+1)] = 0.0;
    }
  }

  order = p + 1 + ( e + 1 ) / 2;

  x_table = ( double * ) malloc ( order * sizeof ( double ) );
  w_table = ( double * ) malloc ( order * sizeof ( double ) );

  l_quadrature_rule ( order, x_table, w_table );

  for ( k = 0; k < order; k++ )
  {
    x = x_table[k];
    l_table = l_polynomial ( 1, p, x_table+k );

    if ( e == 0 )
    {
      for ( j = 0; j <= p; j++ )
      {
        for ( i = 0; i <= p; i++ )
        {
          table[i+j*(p+1)] = table[i+j*(p+1)] + w_table[k] * l_table[i] * l_table[j];
        }
      }
    }
    else
    {
      for ( j = 0; j <= p; j++ )
      {
        for ( i = 0; i <= p; i++ )
        {
          table[i+j*(p+1)] = table[i+j*(p+1)]
            + w_table[k] * pow ( x, e ) * l_table[i] * l_table[j];
        }
      }
    }
    free ( l_table );
  }
  free ( w_table );
  free ( x_table );

  return table;
}
/******************************************************************************/

void l_quadrature_rule ( int n, double x[], double w[] )

/******************************************************************************/
/*
  Purpose:

    L_QUADRATURE_RULE: Gauss-Laguerre quadrature based on L(n,x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int N, the order.

    Output, double X[N], the abscissas.

    Output, double W[N], the weights.
*/
{
  double *bj;
  int i;
  double zemu;
/*
  Define the zero-th moment.
*/
  zemu = 1.0;
/*
  Define the Jacobi matrix.
*/
  bj = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    bj[i] = ( double ) ( i + 1 );
  }

  for ( i = 0; i < n; i++ )
  {
    x[i] = ( double ) ( 2 * i + 1 );
  }

  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
/*
  Diagonalize the Jacobi matrix.
*/
  imtqlx ( n, x, bj, w );

  for ( i = 0; i < n; i++ )
  {
    w[i] = w[i] * w[i];
  }

  free ( bj );

  return;
}
/******************************************************************************/

double *lf_function ( int mm, int n, double alpha, double x[] )

/******************************************************************************/
/*
  Purpose:

    LF_FUNCTION evaluates the Laguerre function Lf(n,alpha,x).

  Recursion:

    Lf(0,ALPHA,X) = 1
    Lf(1,ALPHA,X) = 1+ALPHA-X

    Lf(N,ALPHA,X) = (2*N-1+ALPHA-X)/N * Lf(N-1,ALPHA,X)
                      - (N-1+ALPHA)/N * Lf(N-2,ALPHA,X)

  Restrictions:

    -1 < ALPHA

  Special values:

    Lf(N,0,X) = L(N,X).
    Lf(N,ALPHA,X) = LM(N,ALPHA,X) for ALPHA integral.

  Norm:

    Integral ( 0 <= X < +oo ) exp ( - X ) * Lf(N,ALPHA,X)^2 dX
    = Gamma ( N + ALPHA + 1 ) / N!

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

  Parameters:

    Input, int MM, the number of evaluation points.

    Input, int N, the highest order polynomial to compute.

    Input, double ALPHA, the parameter.  -1.0 < ALPHA.

    Input, double X[MM], the evaluation points.

    Output, double LM_POLYNOMIAL[MM*(N+1)], the function values.
*/
{
  int i;
  int j;
  double *v;

  if ( n < 0 )
  {
    return NULL;
  }

  v = ( double * ) malloc ( mm*(n+1) * sizeof ( double ) );

  for ( i = 0; i < mm; i++ )
  {
    v[i+0*mm] = 1.0;
  }

  if ( n == 0 )
  {
    return v;
  }

  for ( i = 0; i < mm; i++ )
  {
    v[i+1*mm] = alpha + 1.0 - x[i];
  }

  for ( j = 2; j <= n; j++ )
  {
    for ( i = 0; i < mm; i++ )
    {
      v[i+j*mm] = ( (   alpha + ( double ) ( 2 * i - 1 ) - x[i] ) * v[i+(j-1)*mm]
                  + ( - alpha + ( double ) (   - i + 1 )        ) * v[i+(j-2)*mm] )
                    /           ( double ) (     i     );
    }
  }

  return v;
}
/******************************************************************************/

void lf_function_values ( int *n_data, int *n, double *a, double *x,
  double *fx )

/******************************************************************************/
/*
  Purpose:

    LF_FUNCTION_VALUES: some values of the Laguerre function Lf(n,alpha,x).

  Discussion:

    In Mathematica, the function can be evaluated by:

      LaguerreL[n,a,x]

    The functions satisfy the following differential equation:

      X * Y'' + (ALPHA+1-X) * Y' + N * Y = 0;

    Function values can be generated by the recursion:

      Lf(0,ALPHA,X) = 1
      Lf(1,ALPHA,X) = 1+ALPHA-X

      Lf(N,ALPHA,X) = ( (2*N-1+ALPHA-X) * Lf(N-1,ALPHA,X)
                     - (N-1+ALPHA) * Lf(N-2,ALPHA,X) ) / N

    The parameter ALPHA is required to be greater than -1.

    For ALPHA = 0, the generalized Laguerre function Lf(N,ALPHA,X)
    is equal to the Laguerre polynomial L(N,X).

    For ALPHA integral, the generalized Laguerre function
    Lf(N,ALPHA,X) equals the associated Laguerre polynomial Lm(N,ALPHA,X).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Parameters:

    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
    first call.  On each call, the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    Output, int *N, the order of the function.

    Output, double *A, the parameter.

    Output, double *X, the point where the function is evaluated.

    Output, double *FX, the value of the function.
*/
{
# define N_MAX 20

  static double a_vec[N_MAX] = {
     0.00E+00,
     0.25E+00,
     0.50E+00,
     0.75E+00,
     1.50E+00,
     2.50E+00,
     5.00E+00,
     1.20E+00,
     1.20E+00,
     1.20E+00,
     1.20E+00,
     1.20E+00,
     1.20E+00,
     5.20E+00,
     5.20E+00,
     5.20E+00,
     5.20E+00,
     5.20E+00,
     5.20E+00,
     5.20E+00 };

  static double fx_vec[N_MAX] = {
      0.3726399739583333E-01,
      0.3494791666666667E+00,
      0.8710042317708333E+00,
      0.1672395833333333E+01,
      0.6657625325520833E+01,
      0.2395726725260417E+02,
      0.2031344319661458E+03,
      0.1284193996800000E+02,
      0.5359924801587302E+01,
      0.9204589064126984E+00,
     -0.1341585114857143E+01,
     -0.2119726307555556E+01,
     -0.1959193658349206E+01,
      0.1000000000000000E+01,
      0.5450000000000000E+01,
      0.1720125000000000E+02,
      0.4110393750000000E+02,
      0.8239745859375000E+02,
      0.1460179186171875E+03,
      0.2359204608298828E+03 };

  static int n_vec[N_MAX] = {
     5,
     5,
     5,
     5,
     5,
     5,
     5,
     8,
     8,
     8,
     8,
     8,
     8,
     0,
     1,
     2,
     3,
     4,
     5,
     6 };

  static double x_vec[N_MAX] = {
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.00E+00,
     0.20E+00,
     0.40E+00,
     0.60E+00,
     0.80E+00,
     1.00E+00,
     0.75E+00,
     0.75E+00,
     0.75E+00,
     0.75E+00,
     0.75E+00,
     0.75E+00,
     0.75E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *a = 0.0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *a = a_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

double *lf_function_zeros ( int n, double alpha )

/******************************************************************************/
/*
  Purpose:

    LF_FUNCTION_ZEROS returns the zeros of Lf(n,alpha,x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int N, the order.

    Input, double ALPHA, the exponent of the X factor.
    ALPHA must be nonnegative.

    Output, double LF_FUNCTION_ZEROS[N], the zeros.
*/
{
  double *bj;
  int i;
  double i_r8;
  double *w;
  double *x;
  double zemu;
/*
  Define the zero-th moment.
*/
  zemu = tgamma ( alpha + 1.0 );
/*
  Define the Jacobi matrix.
*/
  bj = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    i_r8 = ( double ) ( i );
    bj[i] = ( i_r8 + 1.0 ) * ( i_r8 + 1.0 + alpha );
  }

  x = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    i_r8 = ( double ) ( i );
    x[i] = ( double ) ( 2.0 * i_r8 + 1.0 + alpha );
  }

  w = ( double * ) malloc ( n * sizeof ( double ) );
  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
/*
  Diagonalize the Jacobi matrix.
*/
  imtqlx ( n, x, bj, w );

  for ( i = 0; i < n; i++ )
  {
    w[i] = w[i] * w[i];
  }

  free ( bj );
  free ( w );

  return x;
}
/******************************************************************************/

double lf_integral ( int n, double alpha )

/******************************************************************************/
/*
  Purpose:

    LF_INTEGRAL evaluates a monomial integral associated with Lf(n,alpha,x).

  Discussion:

    The integral:

      integral ( 0 <= x < +oo ) x^n * x^alpha * exp ( -x ) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the exponent.
    0 <= N.

    Input, double ALPHA, the exponent of X in the weight function.

    Output, double LF_INTEGRAL, the value of the integral.
*/
{
  double arg;
  double value;

  arg = alpha + ( double ) ( n + 1 );

  value = tgamma ( arg );

  return value;
}
/******************************************************************************/

void lf_quadrature_rule ( int n, double alpha, double x[], double w[] )

/******************************************************************************/
/*
  Purpose:

    LF_QUADRATURE_RULE: Gauss-Laguerre quadrature rule for Lf(n,alpha,x);

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int N, the order.

    Input, double ALPHA, the exponent of the X factor.
    ALPHA must be nonnegative.

    Output, double X[N], the abscissas.

    Output, double W[N], the weights.
*/
{
  double *bj;
  int i;
  double i_r8;
  double zemu;
/*
  Define the zero-th moment.
*/
  zemu = tgamma ( alpha + 1.0 );
/*
  Define the Jacobi matrix.
*/
  bj = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    i_r8 = ( double ) ( i );
    bj[i] = ( i_r8 + 1.0 ) * ( i_r8 + 1.0 + alpha );
  }

  for ( i = 0; i < n; i++ )
  {
    i_r8 = ( double ) ( i );
    x[i] = ( double ) ( 2.0 * i_r8 + 1.0 + alpha );
  }

  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
/*
  Diagonalize the Jacobi matrix.
*/
  imtqlx ( n, x, bj, w );

  for ( i = 0; i < n; i++ )
  {
    w[i] = w[i] * w[i];
  }

  free ( bj );

  return;
}
/******************************************************************************/

double lm_integral ( int n, int m )

/******************************************************************************/
/*
  Purpose:

    LM_INTEGRAL evaluates a monomial integral associated with Lm(n,m,x).

  Discussion:

    The integral:

      integral ( 0 <= x < +oo ) x^n * x^m * exp ( -x ) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the exponent.
    0 <= N.

    Input, int M, the parameter.
    0 <= M.

    Output, double LM_INTEGRAL, the value of the integral.
*/
{
  double value;

  value = r8_factorial ( n + m );

  return value;
}
/******************************************************************************/

double *lm_polynomial ( int mm, int n, int m, double x[] )

/******************************************************************************/
/*
  Purpose:

    LM_POLYNOMIAL evaluates Laguerre polynomials Lm(n,m,x).

  First terms:

    M = 0

    Lm(0,0,X) =   1
    Lm(1,0,X) =  -X   +  1
    Lm(2,0,X) =   X^2 -  4 X   +  2
    Lm(3,0,X) =  -X^3 +  9 X^2 -  18 X   +    6
    Lm(4,0,X) =   X^4 - 16 X^3 +  72 X^2 -   96 X +     24
    Lm(5,0,X) =  -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -  600 x   +  120
    Lm(6,0,X) =   X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720

    M = 1

    Lm(0,1,X) =    0
    Lm(1,1,X) =   -1,
    Lm(2,1,X) =    2 X - 4,
    Lm(3,1,X) =   -3 X^2 + 18 X - 18,
    Lm(4,1,X) =    4 X^3 - 48 X^2 + 144 X - 96

    M = 2

    Lm(0,2,X) =    0
    Lm(1,2,X) =    0,
    Lm(2,2,X) =    2,
    Lm(3,2,X) =   -6 X + 18,
    Lm(4,2,X) =   12 X^2 - 96 X + 144

    M = 3

    Lm(0,3,X) =    0
    Lm(1,3,X) =    0,
    Lm(2,3,X) =    0,
    Lm(3,3,X) =   -6,
    Lm(4,3,X) =   24 X - 96

    M = 4

    Lm(0,4,X) =    0
    Lm(1,4,X) =    0
    Lm(2,4,X) =    0
    Lm(3,4,X) =    0
    Lm(4,4,X) =   24

  Recursion:

    Lm(0,M,X)   = 1
    Lm(1,M,X)   = (M+1-X)

    if 2 <= N:

      Lm(N,M,X)   = ( (M+2*N-1-X) * Lm(N-1,M,X)
                   +   (1-M-N)    * Lm(N-2,M,X) ) / N

  Special values:

    For M = 0, the associated Laguerre polynomials Lm(N,M,X) are equal
    to the Laguerre polynomials L(N,X).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

  Parameters:

    Input, int MM, the number of evaluation points.

    Input, int N, the highest order polynomial to compute.

    Input, int M, the parameter.  M must be nonnegative.

    Input, double X[MM], the evaluation points.

    Output, double LM_POLYNOMIAL[MM*(N+1)], the function values.
*/
{
  int i;
  int j;
  double *v;

  if ( n < 0 )
  {
    return NULL;
  }

  v = ( double * ) malloc ( mm*(n+1) * sizeof ( double ) );

  for ( j = 0; j <= n; j++ )
  {
    for ( i = 0; i < mm; i++ )
    {
      v[i+j*mm] = 0.0;
    }
  }

  for ( i = 0; i < mm; i++ )
  {
    v[i+0*mm] = 1.0;
  }

  if ( n == 0 )
  {
    return v;
  }

  for ( i = 0; i < mm; i++ )
  {
    v[i+1*mm] = ( double ) ( m + 1 ) - x[i];
  }

  for ( j = 2; j <= n; j++ )
  {
    for ( i = 0; i < mm; i++ )
    {
      v[i+j*mm] = ( ( ( double ) (   m + 2 * j - 1 ) - x[i] ) * v[i+(j-1)*mm]
                    + ( double ) ( - m     - j + 1 )          * v[i+(j-2)*mm] )
                    / ( double ) (           j     );
    }
  }

  return v;
}
/******************************************************************************/

double *lm_polynomial_coefficients ( int n, int m )

/******************************************************************************/
/*
  Purpose:

    LM_POLYNOMIAL_COEFFICIENTS: coefficients of Laguerre polynomial Lm(n,m,x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

  Parameters:

    Input, int N, the highest order polynomial to compute.
    Note that polynomials 0 through N will be computed.

    Input, int M, the parameter.

    Output, double LM_POLYNOMIAL_COEFFICIENTS[(N+1)*(N+1)], the coefficients
    of the Laguerre polynomials of degree 0 through N.
*/
{
  double *c;
  int i;
  int j;

  if ( n < 0 )
  {
    return NULL;
  }

  c = ( double * ) malloc ( (n+1)*(n+1) * sizeof ( double ) );

  for ( i = 0; i <= n; i++ )
  {
    for ( j = 0; j <= n; j++ )
    {
      c[i+j*(n+1)] = 0.0;
    }
  }

  c[0+0*(n+1)] = 1.0;

  if ( n == 0 )
  {
    return c;
  }

  c[1+0*(n+1)] = ( double ) ( m + 1 );
  c[1+1*(n+1)] = -1.0;

  for ( i = 2; i <= n; i++ )
  {
    for ( j = 0; j <= i; j++ )
    {
      c[i+j*(n+1)] = (
          ( double ) (   m + 2 * i - 1 ) * c[i-1+j*(n+1)]
        + ( double ) ( - m     - i + 1 ) * c[i-2+j*(n+1)] )
        / ( double )             i;
    }
    for ( j = 1; j <= i; j++ )
    {
      c[i+j*(n+1)] = c[i+j*(n+1)] - c[i-1+(j-1)*(n+1)] / ( double ) i;
    }
  }
  return c;
}
/******************************************************************************/

void lm_polynomial_values ( int *n_data, int *n, int *m, double *x,
  double *fx )

/******************************************************************************/
/*
  Purpose:

    LM_POLYNOMIAL_VALUES: some values of the Laguerre polynomial Lm(n,m,x).

  Discussion:

    In Mathematica, the function can be evaluated by:

      LaguerreL[n,m,x]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt

  Reference:

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Parameters:

    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
    first call.  On each call, the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    Output, int *N, the order of the function.

    Output, int *M, the parameter.

    Output, double *X, the point where the function is evaluated.

    Output, double *FX, the value of the function.
*/
{
# define N_MAX 20

  static double fx_vec[N_MAX] = {
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1500000000000000E+01,
     0.1625000000000000E+01,
     0.1479166666666667E+01,
     0.1148437500000000E+01,
     0.4586666666666667E+00,
     0.2878666666666667E+01,
     0.8098666666666667E+01,
     0.1711866666666667E+02,
     0.1045328776041667E+02,
     0.1329019368489583E+02,
     0.5622453647189670E+02,
     0.7484729341779436E+02,
     0.3238912982762806E+03,
     0.4426100000097533E+03,
     0.1936876572288250E+04 };

  static int m_vec[N_MAX] = {
    0, 0, 0, 0,
    0, 1, 1, 1,
    1, 0, 1, 2,
    3, 2, 2, 3,
    3, 4, 4, 5 };

  static int n_vec[N_MAX] = {
    1,  2,  3,  4,
    5,  1,  2,  3,
    4,  3,  3,  3,
    3,  4,  5,  6,
    7,  8,  9, 10 };

  static double x_vec[N_MAX] = {
     0.00E+00,
     0.00E+00,
     0.00E+00,
     0.00E+00,
     0.00E+00,
     0.50E+00,
     0.50E+00,
     0.50E+00,
     0.50E+00,
     0.20E+00,
     0.20E+00,
     0.20E+00,
     0.20E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00,
     0.25E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *m = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *m = m_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

double *lm_polynomial_zeros ( int n, int m )

/******************************************************************************/
/*
  Purpose:

    LM_POLYNOMIAL_ZEROS returns the zeros for Lm(n,m,x).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int N, the order.

    Input, int M, the parameter.
    0 <= M.

    Output, double X[N], the zeros.
*/
{
  double *bj;
  int i;
  double *w;
  double *x;
  double zemu;
/*
  Define the zero-th moment.
*/
  zemu = r8_factorial ( m );
/*
  Define the Jacobi matrix.
*/
  bj = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    bj[i] = ( double ) ( i + 1 ) * ( i + 1 + m );
  }

  x = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    x[i] = ( double ) ( 2 * i + 1 + m );
  }

  w = ( double * ) malloc ( n * sizeof ( double ) );
  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
/*
  Diagonalize the Jacobi matrix.
*/
  imtqlx ( n, x, bj, w );

  free ( bj );

  return x;
}
/******************************************************************************/

void lm_quadrature_rule ( int n, int m, double x[], double w[] )

/******************************************************************************/
/*
  Purpose:

    LM_QUADRATURE_RULE: Gauss-Laguerre quadrature rule for Lm(n,m,x);

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 August 2013

  Author:

    John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int N, the order.

    Input, int M, the parameter.
    0 <= M.

    Output, double X[N], the abscissas.

    Output, double W[N], the weights.
*/
{
  double *bj;
  int i;
  double zemu;
/*
  Define the zero-th moment.
*/
  zemu = r8_factorial ( m );
/*
  Define the Jacobi matrix.
*/
  bj = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    bj[i] = ( double ) ( i + 1 ) * ( i + 1 + m );
  }

  for ( i = 0; i < n; i++ )
  {
    x[i] = ( double ) ( 2 * i + 1 + m );
  }

  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
/*
  Diagonalize the Jacobi matrix.
*/
  imtqlx ( n, x, bj, w );

  for ( i = 0; i < n; i++ )
  {
    w[i] = w[i] * w[i];
  }

  free ( bj );

  return;
}
/******************************************************************************/

double r8_factorial ( int n )

/******************************************************************************/
/*
  Purpose:

    R8_FACTORIAL computes the factorial of N.

  Discussion:

    factorial ( N ) = product ( 1 <= I <= N ) I

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the argument of the factorial function.
    If N is less than 1, the function value is returned as 1.

    Output, double R8_FACTORIAL, the factorial of N.
*/
{
  int i;
  double value;

  value = 1.0;

  for ( i = 1; i <= n; i++ )
  {
    value = value * ( double ) ( i );
  }

  return value;
}
/******************************************************************************/

double r8_sign ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_SIGN returns the sign of an R8.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, the number whose sign is desired.

    Output, double R8_SIGN, the sign of X.
*/
{
  double value;

  if ( x < 0.0 )
  {
    value = - 1.0;
  }
  else
  {
    value = + 1.0;
  }
  return value;
}
/******************************************************************************/

void r8mat_print ( int m, int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8MAT_PRINT prints an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

    Entry A(I,J) is stored as A[I+J*M]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 May 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows in A.

    Input, int N, the number of columns in A.

    Input, double A[M*N], the M by N matrix.

    Input, char *TITLE, a title.
*/
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    R8MAT_PRINT_SOME prints some of an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Input, double A[M*N], the matrix.

    Input, int ILO, JLO, IHI, JHI, designate the first row and
    column, and the last row and column to be printed.

    Input, char *TITLE, a title.
*/
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }
/*
  Print the columns of the matrix, in strips of 5.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:  ");
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %7d     ", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( m < ihi )
    {
      i2hi = m;
    }
    else
    {
      i2hi = ihi;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
/*
  Print out (up to) 5 entries in row I, that lie in the current strip.
*/
      fprintf ( stdout, "%5d:", i - 1 );
      for ( j = j2lo; j <= j2hi; j++ )
      {
        fprintf ( stdout, "  %14f", a[i-1+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}
/******************************************************************************/

double r8vec_dot_product ( int n, double a1[], double a2[] )

/******************************************************************************/
/*
  Purpose:

    R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 July 2007

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vectors.

    Input, double A1[N], A2[N], the two vectors to be considered.

    Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors.
*/
{
  int i;
  double value;

  value = 0.0;
  for ( i = 0; i < n; i++ )
  {
    value = value + a1[i] * a2[i];
  }
  return value;
}
/******************************************************************************/

void r8vec_print ( int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8VEC_PRINT prints an R8VEC.

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 April 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of components of the vector.

    Input, double A[N], the vector to be printed.

    Input, char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );
  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %8d: %14f\n", i, a[i] );
  }

  return;
}
/******************************************************************************/

void r8vec2_print ( int n, double a1[], double a2[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8VEC2_PRINT prints an R8VEC2.

  Discussion:

    An R8VEC2 is a dataset consisting of N pairs of real values, stored
    as two separate vectors A1 and A2.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 March 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of components of the vector.

    Input, double A1[N], double A2[N], the vectors to be printed.

    Input, char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );
  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %4d: %14f  %14f\n", i, a1[i], a2[i] );
  }

  return;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}