# include <limits.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>

# include "matrix_chain_brute.h"

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

int catalan_number ( int n )

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

    catalan_number() computes the N-th Catalan number.

  First values:

     C(0)     1
     C(1)     1
     C(2)     2
     C(3)     5
     C(4)    14
     C(5)    42
     C(6)   132
     C(7)   429
     C(8)  1430
     C(9)  4862
    C(10) 16796

  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 )
  {
    return 0;
  }

  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;
}
/******************************************************************************/

int i4_factorial ( int n )

/******************************************************************************/
/*
  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;
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "i4_factorial(): Fatal error!\n" );
    fprintf ( stderr, "  i4_factorial(n) cannot be computed as an integer\n" );
    fprintf ( stderr, "  for 13 < N.\n" );
    fprintf ( stderr, "  Input value N = %d\n", n );
    exit ( 1 );
  }

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

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

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

/******************************************************************************/
/*
  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 = ( int * ) malloc ( n_mults * sizeof ( int ) );
  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];
      }
    }
  }

  free ( this_p );

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

bool pivot_sequence_check ( int n, int t[] )

/******************************************************************************/
/*
  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 )
      {
        printf ( "\n" );
        printf ( "pivot_sequence_check(): Fatal error!\n" );
        printf ( "  t[%d] = %d <= 0\n", i, t[i] );
      }
      check = false;
      return check;
    }
    else if ( n - i < t[i] )
    {
      if ( verbose )
      {
        printf ( "\n" );
        printf ( "pivot_seq_check - Fatal error!\n" );
        printf ( "  n-i = %d < t[%d] = %d\n", n - i, i, t[i] );
      }
      check = false;
      return false;
    }

  }

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

int pivot_sequence_enum ( int n )

/******************************************************************************/
/*
  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=linspace(n,1,n).

  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 )
  {
    printf ( "\n" );
    printf ( "pivot_sequence_successor(): Fatal error!\n" );
    printf ( "  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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 = ( int * ) malloc ( n_dims2 * sizeof ( int ) );
  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;
  }

  free ( dims2 );

  return;
}