# include <climits>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "matrix_chain_brute.hpp"

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

int catalan_number ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    catalan_number() computes the N-th Catalan number.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    23 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the index of the Catalan number.
//
//  Output:
//
//    int C, the value of the Catalan number.
//
{
  int c;
  int i;

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

  c = 1;
//
//  The extra parentheses ensure that the integer division is
//  done AFTER the integer multiplication.
//
  for ( i = 1; i <= n; i++ )
  {
    c = ( c * 2 * ( 2 * i - 1 ) ) / ( i + 1 );
  }

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

int i4_factorial ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    i4_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
//
//  Input:
//
//    int N, the argument of the factorial function.
//    If N is less than 1, the function value is returned as 1.
//    0 <= N <= 13 is required.
//
//  Output:
//
//    int I4_FACTORIAL, the factorial of N.
//
{
  int i;
  int value;

  value = 1;

  if ( 13 < n ) 
  {
    value = - 1;
    cerr << "i4_factorial(): Fatal error!\n";
    cerr << "  i4_factorial(n) cannot be computed as an integer\n";
    cerr << "  for 13 < N.\n";
    cerr << "  Input value N = " << n << "\n";
    exit ( 1 );
  }

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

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

void matrix_chain_brute ( int n_mats, int dims[], int &cost, int p[] )

//****************************************************************************80
//
//  Purpose:
//
//    matrix_chain_brute() finds the lowest cost to form a multiple matrix product.
//
//  Discussion:
//
//    This code represents a brute force approach.
//
//    An "efficient" brute force approach would only go through every
//    possible parenthesization of the multiplication, and a cost of
//    catalan(n_mats-1).  But it's not clear to me how to generate a sequence
//    of parentheses and properly interpret them as a multiplication ordering.
//
//    Instead, we rely on the fact that, to multiply N matrices, there are
//    N-1 multiplications to carry out, in any order.  This means that,
//    if we do some careful bookkeeping, we have N-1 choices for the first
//    multiplication, N-2 for the second, and a total of factorial(N-1)    
//    distinct cases to consider.  Some of these cases collapse to the
//    same paretheses representation, but we don't care.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    27 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n_mats: the number of matrices in the product.
//
//    integer dims(n_mats+1): matrix dimension information.  Matrix A(i)
//    has dimensions dims(i) x dims(i+1).  All entries must be positive.
//
//  Output:
//
//    integer cost: the minimal cost, in terms of scalar multiplications,
//    for the optimal ordering of the matrix multiplications.
//
//    integer p(n_mats-1): defines the order of the multiplications.
//
{
  int i;
  int j;
  int n_mults;
  int pivot_sequence_num;
  int this_cost;
  int *this_p;

  n_mults = n_mats - 1;
//
//  Deal with stupidity.
//
  if ( n_mats == 1 )
  {
    cost = 0;
    return;
  }
//
//  Initialize the output.
//
  cost = INT_MAX;
  for ( i = 0; i < n_mults; i++ )
  {
    p[i] = n_mults - i;
  }
//
//  Prepare to loop over all orderings.
//
  this_p = new int[n_mults];
  for ( i = 0; i < n_mults; i++ )
  {
    this_p[i] = p[i];
  }
  pivot_sequence_num = pivot_sequence_enum ( n_mults );

  for ( j = 0; j < pivot_sequence_num; j++ )
  {
    pivot_sequence_successor ( n_mults, this_p );
    pivot_sequence_to_matrix_chain_cost ( n_mats, this_p, dims, this_cost );

    if ( this_cost < cost )
    {
      cost = this_cost;
      for ( i = 0; i < n_mults; i++ )
      {
        p[i] = this_p[i];
      }
    }
  }

  delete [] this_p;

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

bool pivot_sequence_check ( int n, int t[] )

//****************************************************************************80
//
//  Purpose:
//
//    pivot_sequence_check() checks a pivot sequence.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of pivot steps.
//
//    integer t(n): a pivot sequence.
//
//  Output:
//
//    logical CHECK, error flag.
//    true, T is a pivot sequence.
//    false, T is not a legal pivot sequence.
//
{
  bool check;
  int i;
  bool verbose;

  verbose = false;
  check = true;

  for ( i = 0; i < n; i++ )
  {
    if ( t[i] <= 0 )
    {
      if ( verbose )
      {
        cout << "\n";
        cout << "pivot_sequence_check(): Fatal error!\n";
        cout << "  t[" << i << "] = " << t[i] << " <= 0\n";
      }
      check = false;
      return check;
    }
    else if ( n - i < t[i] )
    {
      if ( verbose )
      {
        cout << "\n";
        cout << "pivot_seq_check(): Fatal error!\n";
        cout << "  n-i = " << n - i << " < t[" << i << "] = " << t[i] << "\n";
      }
      check = false;
      return false;
    }

  }

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

int pivot_sequence_enum ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    pivot_sequence_enum() enumerates pivot sequences.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of pivot steps.
//
//  Output:
//
//    integer pivot_sequence_num: the number of pivot sequences of n steps.
//
{
  int pivot_sequence_num;

  pivot_sequence_num = i4_factorial ( n );

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

void pivot_sequence_successor ( int n, int t[] )

/******************************************************************************/
//
//  Purpose:
//
//    pivot_sequence_successor() computes the pivot sequence successor.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of pivot steps.
//
//    integer t(n): the previous pivot sequence.
//    To initiate the routine, call with t={n,n-1,...,1}.
//
//  Output:
//
//    integer t(n): the lexical successor of the input.
//
{
  bool check;
  int i;
  int last;
//
//  Check.
//
  check = pivot_sequence_check ( n, t );

  if ( ! check )
  {
    cout << "\n";
    cout << "pivot_sequence_successor(): Fatal error!\n";
    cout << "  The input array is illegal.\n";
  }
//
//  Find last entry that is less than its maximum value.
//
  last = -1;
  for ( i = n - 1; 0 <= i; i-- )
  {
    if ( t[i] < n - i )
    {
      last = i;
      break;
    }
  }

  if ( last == -1 )
  {
    for ( i = 0; i < n; i++ )
    {
      t[i] = 1;
    }
  }
  else
  {
    t[last] = t[last] + 1;
    for ( i = last + 1; i < n; i++ )
    {
      t[i] = 1;
    }
  }

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

void pivot_sequence_to_matrix_chain_cost ( int n_mats, int p[], int dims[],
  int &cost )

//****************************************************************************80
//
//  Purpose:
//
//    pivot_sequence_to_matrix_chain_cost() evaluates a matrix chain product.
//
//  Discussion:
//
//    There are n_mats matrices to multiply.
//    There are n_dims = n_mats+1 dimensions to consider.
//    There are n_mults = n_mats-1 matrix multiplications to carry out.
//
//    The cost, in terms of multiplications of pairs of real numbers,
//    of a multiplying a single pair of matrices A and B,
//    of dimensions LxM and MxN, is L*M*N.
//
//    The cost of multiplying a sequence of N matrices A*B*C*...*Z
//    will in general vary, depending on the order in which the N-1
//    multiplications are carried out.
//
//    This function assumes that a particular multiplication ordering 
//    has been specified in the pivot sequence p(), and returns the
//    corresponding cost.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    27 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n_mats: the number of matrices to multiply.
//
//    integer p(n_mats-1): the pivot sequence.
//    On the i-th step of the multiplication procedure, there are
//    n-i multiplications to choose from, and we choose the p(i) one.
//    It must be the case that, for each i, 1 <= p(i) <= n - i
//
//    integer dims(n_mats+1): the matrix dimensions.
//    For 1 <= i < n, the i-th matrix has dimension dims(i) x dims(i+1).
//
//  Output:
//
//    integer cost: the matrix multiplication cost.
//
{
  int d1;
  int d2;
  int d3;
  int *dims2;
  int i;
  int j;
  int k;
  int n_dims;
  int n_dims2;
  int n_mults;

  cost = 0;
  n_dims = n_mats + 1;
  n_mults = n_mats - 1;
  n_dims2 = n_dims;
  dims2 = new int[n_dims2];
  for ( i = 0; i < n_dims2; i++ )
  {
    dims2[i] = dims[i];
  }
//
//  Carry out n_mults multiplications.
//  On step I, we carry out multiplication J.
//  Multiplication J involves 
//    dims(j)dims(j+1) * dims(j+1)dims(j+2) => dims(j)dims(j+2)
//  at an additional cost of dims(j) * dims(j+1) * dims(j+2),
//  while removing dims(j+1) from the list of dimensions.
//
  for ( i = 0; i < n_mults; i++ )
  {
    j = p[i] - 1;
    d1 = dims2[j];
    d2 = dims2[j+1];
    d3 = dims2[j+2];
    cost = cost + d1 * d2 * d3;
    for ( k = j + 1; k < n_dims2 - 1; k++ )
    {
      dims2[k] = dims2[k+1];
    }
    n_dims2 = n_dims2 - 1;
  }

  delete [] dims2;

  return;
}