# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "r8utp.hpp"

//****************************************************************************80

int i4_log_10 ( int i )

//****************************************************************************80
//
//  Purpose:
//
//    i4_log_10() returns the integer part of the logarithm base 10 of an I4.
//
//  Example:
//
//        I  I4_LOG_10
//    -----  --------
//        0    0
//        1    0
//        2    0
//        9    0
//       10    1
//       11    1
//       99    1
//      100    2
//      101    2
//      999    2
//     1000    3
//     1001    3
//     9999    3
//    10000    4
//
//  Discussion:
//
//    I4_LOG_10 ( I ) + 1 is the number of decimal digits in I.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 January 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I, the number whose logarithm base 10 is desired.
//
//    Output, int I4_LOG_10, the integer part of the logarithm base 10 of
//    the absolute value of X.
//
{
  int i_abs;
  int ten_pow;
  int value;

  if ( i == 0 )
  {
    value = 0;
  }
  else
  {
    value = 0;
    ten_pow = 10;

    i_abs = abs ( i );

    while ( ten_pow <= i_abs )
    {
      value = value + 1;
      ten_pow = ten_pow * 10;
    }

  }

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

int i4_max ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    i4_max() returns the maximum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 October 1998
//
//  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;
}
//****************************************************************************80

int i4_min ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    i4_min() returns the minimum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 October 1998
//
//  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;
}
//****************************************************************************80

void r8ge_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8ge_print() prints an R8GE matrix.
//
//  Discussion:
//
//    The R8GE storage format is used for a "general" M by N matrix.  
//    A physical storage space is made for each logical entry.  The two 
//    dimensional logical array is mapped to a vector, in which storage is 
//    by columns.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    06 April 2006
//
//  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 R8GE matrix.
//
//    Input, string TITLE, a title.
//
{
  r8ge_print_some ( m, n, a, 1, 1, m, n, title );

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

void r8ge_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, 
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8ge_print_some() prints some of an R8GE matrix.
//
//  Discussion:
//
//    The R8GE storage format is used for a "general" M by N matrix.  
//    A physical storage space is made for each logical entry.  The two 
//    dimensional logical array is mapped to a vector, in which storage is 
//    by columns.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    06 April 2006
//
//  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 R8GE matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, string TITLE, a title.
//
{
# define INCX 5

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

  cout << "\n";
  cout << title << "\n";
//
//  Print the columns of the matrix, in strips of 5.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    j2hi = i4_min ( j2hi, n );
    j2hi = i4_min ( j2hi, jhi );

    cout << "\n";
//
//  For each column J in the current range...
//
//  Write the header.
//
    cout << "  Col:    ";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(7) << j << "       ";
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "  ---\n";
//
//  Determine the range of the rows in this strip.
//
    i2lo = i4_max ( ilo, 1 );
    i2hi = i4_min ( ihi, m );

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

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

double *r8ge_to_r8utp ( int m, int n, double a_ge[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8ge_to_r8utp() copies an R8GE matrix to an R8UTP matrix.
//
//  Discussion:
//
//    The R8GE storage format is used for a general M by N matrix.  A storage 
//    space is made for each entry.  The two dimensional logical
//    array can be thought of as a vector of M*N entries, starting with
//    the M entries in the column 1, then the M entries in column 2
//    and so on.  Considered as a vector, the entry A(I,J) is then stored
//    in vector location I+(J-1)*M.
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the number of rows and columns of 
//    the matrix.
//
//    double A_GE(N,N), the R8GE matrix.
//
//  Output:
//
//    double r8ge_to_r8utp[*], the R8UTP matrix.
//
{
  double *a_utp;
  int i;
  int j;
  int k;

  a_utp = r8utp_zeros ( m, n );

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < i4_min ( j + 1, m ); i++ )
    {
      a_utp[k] = a_ge[i+j*m];
      k = k + 1;
    }
  }

  return a_utp;
}
//****************************************************************************80

double r8utp_det ( int m, int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_det() computes the determinant of an R8UTP matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the order of the matrix.
//
//    double A[*], the R8UTP matrix.
//
//  Output:
//
//    double r8utp_det: the determinant of the matrix.
//
{
  double det;
  int j;
  int k;

  if ( m != n )
  {
    cout << "\n";
    cout << "r8utp_det(): Fatal error!\n";
    cout << "  m and n are not equal.\n";
    exit ( 1 );
  }

  det = 1.0;
  k = 0;
  for ( j = 0; j < n; j++ )
  {
    det = det * a[k];
    k = k + j + 2;
  }

  return det;
}
//****************************************************************************80

double *r8utp_indicator ( int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_indicator() sets up a R8UTP indicator matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the number of rows and columns of the matrix.
//
//  Output:
//
//    double r8utp_indicator[]: the R8UTP matrix.
//
{
  double *a;
  int fac;
  int i;
  int j;
  int k;
  int mn;

  mn = r8utp_size ( m, n );
  a = new double[mn];
  fac = pow ( 10, i4_log_10 ( n ) + 1 );

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < i4_min ( j + 1, m ); i++ )
    {
      a[k] = fac * ( i + 1 ) + ( j + 1 );
      k = k + 1;
    }
  }

  return a;
}
//****************************************************************************80

void r8utp_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_print() prints a R8UTP matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN). 
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the order of the matrix.
//
//    Input, double A[*], the matrix.
//
//    Input, string TITLE, a title.
//
{
  r8utp_print_some ( m, n, a, 1, 1, m, n, title );

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

void r8utp_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, 
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_print_some() prints some of an R8UTP matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN). 
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the order of the matrix.
//
//    Input, double A[*], the matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, string TITLE, a title.
//
{
# define INCX 5

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

  cout << "\n";
  cout << title << "\n";
//
//  Print the columns of the matrix, in strips of 5.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    j2hi = i4_min ( j2hi, n );
    j2hi = i4_min ( j2hi, jhi );

    cout << "\n";
    cout << "  Col: ";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(7) << j << "       ";
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "  ---\n";
//
//  Determine the range of the rows in this strip.
//
    i2lo = i4_max ( ilo, 1 );
    i2hi = i4_min ( ihi, m );

    for ( i = i2lo; i <= i2hi; i++ )
    {
      cout << setw(6) << i << "  ";
//
//  Print out (up to) 5 entries in row I, that lie in the current strip.
//
      for ( j = j2lo; j <= j2hi; j++ )
      {
        if ( i <= j )
        {
          aij = a[i-1+(j*(j-1))/2];
        }
        else
        {
          aij = 0.0;
        }

        cout << setw(12) << aij << "  ";
      }
      cout << "\n";
    }
  }

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

double *r8utp_random ( int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_random() randomizes an R8UTP matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the number of rows and columns of the matrix.
//    M and N must be positive.
//
//  Output:
//
//    double r8utp_random[*], the R8UTP matrix.
//
{
  double *a;
  int i;
  int mn;

  mn = r8utp_size ( m, n );
  a = new double[mn];

  for ( i = 0; i < mn; i++ )
  {
    a[i] = drand48 ( );
  }

  return a;
}
//****************************************************************************80

int r8utp_size ( int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_size() returns the size of an M x N R8UTP matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer M, N: the number of rows and columns.
//
//  Output:
//
//    integer r8utp_size: the length of the array needed to store the matrix.
//
{
  int value;

  if ( n < m )
  {
    value = ( n * ( n + 1 ) ) / 2;
  }
  else if ( m == n )
  {
    value = ( m * ( m + 1 ) ) / 2;
  }
  else
  {
    value = ( m * ( m + 1 ) / 2 ) + ( n - m ) * m;
  }

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

double *r8utp_to_r8ge ( int m, int n, double a_utp[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_to_r8ge() copies an R8UTP matrix to an R8GE matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//    The R8GE storage format is used for a general M by N matrix.  A storage 
//    space is made for each entry.  The two dimensional logical
//    array can be thought of as a vector of M*N entries, starting with
//    the M entries in the column 1, then the M entries in column 2
//    and so on.  Considered as a vector, the entry A(I,J) is then stored
//    in vector location I+(J-1)*M.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the order of the matrix.
//
//    double A_UTP(*), the R8UTP matrix.
//
//  Output:
//
//    double r8utp_to_r8ge[m,n]: the R8GE matrix.
//
{
  double *a_ge;
  int i;
  int j;
  int k;

  a_ge = new double[ m * n ];

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      if ( i < i4_min ( j + 1, m ) )
      {
        a_ge[i+j*m] = a_utp[k];
        k = k + 1;
      }
      else
      {
        a_ge[i+j*m] = 0.0;
      }
    }
  }

  return a_ge;
}
//****************************************************************************80

double *r8utp_zeros ( int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    r8utp_zeros() zeros an R8UTP matrix.
//
//  Discussion:
//
//    The R8UTP storage format is appropriate for an upper triangular
//    matrix.  Only the upper triangle of the matrix is stored,
//    by successive partial columns, in an array which contains
//    (A11,A12,A22,A13,A23,A33,A14,...,AMN).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 August 2022
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the number of rows and columns of the matrix.
//    M and N must be positive.
//
//  Output:
//
//    double r8utp_random[*], the R8UTP matrix.
//
{
  double *a;
  int i;
  int mn;

  mn = r8utp_size ( m, n );
  a = new double[mn];

  for ( i = 0; i < mn; i++ )
  {
    a[i] = 0.0;
  }

  return a;
}