# include <math.h>
# include <stdlib.h>
# include <stdio.h>

# include "asa299.h"

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

void simplex_lattice_point_next ( int n, int t, int *more, int x[] )

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

    simplex_lattice_point_next() generates lattice points in a simplex.

  Discussion:

    The simplex is defined by N-dimensional points X such that:

        0 <= X(1:N)

    and

      sum ( X(1:N) ) <= T

    where T is an integer.

    Lattice points are points X which satisfy the simplex conditions and
    for which all the components are integers.

    This routine generates all the lattice points in a given simplex, one at 
    a time, in a reverse lexicographic order.

    To use the routine, initialize by setting N and T to appropriate values, 
    and MORE to FALSE.  The initial value of X is not important.

    Call the routine. On return, X will contain the first lattice point in 
    the simplex.  If MORE is TRUE, then the routine may be called again to 
    get the next point.  In fact, as long as the output value of MORE is 
    TRUE, there is at least one more lattice point that can be found by 
    making another call.  When MORE is returned as FALSE, then there are no 
    more lattice points; the value of X returned at that time is the 
    "last" such point.

    During the computation of a sequence of lattice points, the user should 
    not change the values of N, T, MORE or X.  

    The output for N = 3, T = 4 would be:

       1    4  0  0
       2    3  1  0
       3    3  0  1
       4    2  2  0
       5    2  1  1
       6    2  0  2
       7    1  3  0
       8    1  2  1
       9    1  1  2
      10    1  0  3
      11    0  4  0
      12    0  3  1
      13    0  2  2
      14    0  1  3
      15    0  0  4

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    10 January 2007

  Author:

    John Burkardt

  Reference:

    Scott Chasalow, Richard Brand,
    Algorithm AS 299:
    Generation of Simplex Lattice Points,
    Applied Statistics,
    Volume 44, Number 4, 1995, pages 534-545.

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms for Computers and Calculators,
    Second Edition,
    Academic Press, 1978,
    ISBN: 0-12-519260-6,
    LC: QA164.N54.

  Parameters:

    Input, int N, the spatial dimension.
    N must be positive.

    Input, int T, the characteristic of the simplex.
    T must be nonnegative.

    Input/output, int *MORE, initialized to FALSE by the user to
    begin a sequence of calculations, returned by the routine as TRUE,
    if there are more values of X that can be calculated, or FALSE
    if the accompanying value of X is the last one for this sequence.

    Input/output, int X[N], not initialized by the user, but not
    changed by the user on subsequent calls.  The routine returns
    a new point on each call, and on subsequent calls uses the input
    value (old point) to compute the output value (next point).
*/
{

  int i;
  int j;

  if ( !(*more) )
  {
    if ( n < 1 )
    {
      printf ( "\n" );
      printf ( "SIMPLEX_LATTICE_POINT_NEXT - Fatal error!\n" );
      printf ( "  N < 1.\n" );
      exit ( 1 );
    }

    if ( t < 0 )
    {
      printf ( "\n" );
      printf ( "SIMPLEX_LATTICE_POINT_NEXT - Fatal error!\n" );
      printf ( "  T < 0.\n" );
      exit ( 1 );
    }

    *more = 1;
    j = 1;

    x[0] = t;
    for ( i = 1; i < n; i++ )
    {
      x[i] = 0;
    }
/*
  The first point can actually also be the last!
*/
    if ( n == 1 )
    {
      *more = 0;
    }
  }
  else
  {
/*
  Search X(N-1 down to 1) for the first nonzero element.
  If none, then terminate.  (This should not happen!)
  Otherwise, set J to this index.
  Decrement X(J) by 1.
  Set X(J+1:N) to (T-X(1:J),0,0,...0).
*/
    j = n - 1;

    for ( i = n - 2; 0 <= i; i-- )
    {
      if ( 0 < x[i] )
      {
        j = i;
        break;
      }
    }

    if ( j == n - 1 )
    {
      printf ( "\n" );
      printf ( "SIMPLEX_LATTICE_POINT_NEXT - Fatal error!\n" );
      printf ( "  The input X vector is nonpositive in all entries\n" );
      printf ( "  except possibly the last one.\n" );
      printf ( "\n" );
      printf ( "  Perhaps the user has miscalled the routine\n" );
      printf ( "  or altered data between calls.\n" );
      printf ( "\n" );
      printf ( "ABNORMAL TERMINATION.\n" );
      exit ( 1 );
    }

    x[j] = x[j] - 1;
    x[j+1] = t;
    for ( i = 0; i <= j; i++ )
    {
      x[j+1] = x[j+1] - x[i];
    }
    for ( i = j+2; i < n; i++ )
    {
      x[i] = 0;
    } 
/*
  Is this the last point?
*/
    if ( x[n-1] == t )
    {
      *more = 0;
    }
  }
  return;
}