/* glpmpl03.c */

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

*  This code is part of GLPK (GNU Linear Programming Kit).

*

*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

*  E-mail: <mao@gnu.org>.

*

*  GLPK is free software: you can redistribute it and/or modify it

*  under the terms of the GNU General Public License as published by

*  the Free Software Foundation, either version 3 of the License, or

*  (at your option) any later version.

*

*  GLPK is distributed in the hope that it will be useful, but WITHOUT

*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

*  License for more details.

*

*  You should have received a copy of the GNU General Public License

*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

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

#define _GLPSTD_ERRNO

#define _GLPSTD_STDIO

#include "glpenv.h"

#include "glpmpl.h"

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

/* * *                   FLOATING-POINT NUMBERS                   * * */

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

/*----------------------------------------------------------------------

-- fp_add - floating-point addition.

--

-- This routine computes the sum x + y. */

double fp_add(MPL *mpl, double x, double y)

{     if (x > 0.0 && y > 0.0 && x > + 0.999 * DBL_MAX - y ||

          x < 0.0 && y < 0.0 && x < - 0.999 * DBL_MAX - y)

         error(mpl, "%.*g + %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      return x + y;

}

/*----------------------------------------------------------------------

-- fp_sub - floating-point subtraction.

--

-- This routine computes the difference x - y. */

double fp_sub(MPL *mpl, double x, double y)

{     if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y ||

          x < 0.0 && y > 0.0 && x < - 0.999 * DBL_MAX + y)

         error(mpl, "%.*g - %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      return x - y;

}

/*----------------------------------------------------------------------

-- fp_less - floating-point non-negative subtraction.

--

-- This routine computes the non-negative difference max(0, x - y). */

double fp_less(MPL *mpl, double x, double y)

{     if (x < y) return 0.0;

      if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y)

         error(mpl, "%.*g less %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      return x - y;

}

/*----------------------------------------------------------------------

-- fp_mul - floating-point multiplication.

--

-- This routine computes the product x * y. */

double fp_mul(MPL *mpl, double x, double y)

{     if (fabs(y) > 1.0 && fabs(x) > (0.999 * DBL_MAX) / fabs(y))

         error(mpl, "%.*g * %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      return x * y;

}

/*----------------------------------------------------------------------

-- fp_div - floating-point division.

--

-- This routine computes the quotient x / y. */

double fp_div(MPL *mpl, double x, double y)

{     if (fabs(y) < DBL_MIN)

         error(mpl, "%.*g / %.*g; floating-point zero divide",

            DBL_DIG, x, DBL_DIG, y);

      if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))

         error(mpl, "%.*g / %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      return x / y;

}

/*----------------------------------------------------------------------

-- fp_idiv - floating-point quotient of exact division.

--

-- This routine computes the quotient of exact division x div y. */

double fp_idiv(MPL *mpl, double x, double y)

{     if (fabs(y) < DBL_MIN)

         error(mpl, "%.*g div %.*g; floating-point zero divide",

            DBL_DIG, x, DBL_DIG, y);

      if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))

         error(mpl, "%.*g div %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      x /= y;

      return x > 0.0 ? floor(x) : x < 0.0 ? ceil(x) : 0.0;

}

/*----------------------------------------------------------------------

-- fp_mod - floating-point remainder of exact division.

--

-- This routine computes the remainder of exact division x mod y.

--

-- NOTE: By definition x mod y = x - y * floor(x / y). */

double fp_mod(MPL *mpl, double x, double y)

{     double r;

      xassert(mpl == mpl);

      if (x == 0.0)

         r = 0.0;

      else if (y == 0.0)

         r = x;

      else

      {  r = fmod(fabs(x), fabs(y));

         if (r != 0.0)

         {  if (x < 0.0) r = - r;

            if (x > 0.0 && y < 0.0 || x < 0.0 && y > 0.0) r += y;

         }

      }

      return r;

}

/*----------------------------------------------------------------------

-- fp_power - floating-point exponentiation (raise to power).

--

-- This routine computes the exponentiation x ** y. */

double fp_power(MPL *mpl, double x, double y)

{     double r;

      if (x == 0.0 && y <= 0.0 || x < 0.0 && y != floor(y))

         error(mpl, "%.*g ** %.*g; result undefined",

            DBL_DIG, x, DBL_DIG, y);

      if (x == 0.0) goto eval;

      if (fabs(x) > 1.0 && y > +1.0 &&

            +log(fabs(x)) > (0.999 * log(DBL_MAX)) / y ||

          fabs(x) < 1.0 && y < -1.0 &&

            +log(fabs(x)) < (0.999 * log(DBL_MAX)) / y)

         error(mpl, "%.*g ** %.*g; floating-point overflow",

            DBL_DIG, x, DBL_DIG, y);

      if (fabs(x) > 1.0 && y < -1.0 &&

            -log(fabs(x)) < (0.999 * log(DBL_MAX)) / y ||

          fabs(x) < 1.0 && y > +1.0 &&

            -log(fabs(x)) > (0.999 * log(DBL_MAX)) / y)

         r = 0.0;

      else

eval:    r = pow(x, y);

      return r;

}

/*----------------------------------------------------------------------

-- fp_exp - floating-point base-e exponential.

--

-- This routine computes the base-e exponential e ** x. */

double fp_exp(MPL *mpl, double x)

{     if (x > 0.999 * log(DBL_MAX))

         error(mpl, "exp(%.*g); floating-point overflow", DBL_DIG, x);

      return exp(x);

}

/*----------------------------------------------------------------------

-- fp_log - floating-point natural logarithm.

--

-- This routine computes the natural logarithm log x. */

double fp_log(MPL *mpl, double x)

{     if (x <= 0.0)

         error(mpl, "log(%.*g); non-positive argument", DBL_DIG, x);

      return log(x);

}

/*----------------------------------------------------------------------

-- fp_log10 - floating-point common (decimal) logarithm.

--

-- This routine computes the common (decimal) logarithm lg x. */

double fp_log10(MPL *mpl, double x)

{     if (x <= 0.0)

         error(mpl, "log10(%.*g); non-positive argument", DBL_DIG, x);

      return log10(x);

}

/*----------------------------------------------------------------------

-- fp_sqrt - floating-point square root.

--

-- This routine computes the square root x ** 0.5. */

double fp_sqrt(MPL *mpl, double x)

{     if (x < 0.0)

         error(mpl, "sqrt(%.*g); negative argument", DBL_DIG, x);

      return sqrt(x);

}

/*----------------------------------------------------------------------

-- fp_sin - floating-point trigonometric sine.

--

-- This routine computes the trigonometric sine sin(x). */

double fp_sin(MPL *mpl, double x)

{     if (!(-1e6 <= x && x <= +1e6))

         error(mpl, "sin(%.*g); argument too large", DBL_DIG, x);

      return sin(x);

}

/*----------------------------------------------------------------------

-- fp_cos - floating-point trigonometric cosine.

--

-- This routine computes the trigonometric cosine cos(x). */

double fp_cos(MPL *mpl, double x)

{     if (!(-1e6 <= x && x <= +1e6))

         error(mpl, "cos(%.*g); argument too large", DBL_DIG, x);

      return cos(x);

}

/*----------------------------------------------------------------------

-- fp_atan - floating-point trigonometric arctangent.

--

-- This routine computes the trigonometric arctangent atan(x). */

double fp_atan(MPL *mpl, double x)

{     xassert(mpl == mpl);

      return atan(x);

}

/*----------------------------------------------------------------------

-- fp_atan2 - floating-point trigonometric arctangent.

--

-- This routine computes the trigonometric arctangent atan(y / x). */

double fp_atan2(MPL *mpl, double y, double x)

{     xassert(mpl == mpl);

      return atan2(y, x);

}

/*----------------------------------------------------------------------

-- fp_round - round floating-point value to n fractional digits.

--

-- This routine rounds given floating-point value x to n fractional

-- digits with the formula:

--

--    round(x, n) = floor(x * 10^n + 0.5) / 10^n.

--

-- The parameter n is assumed to be integer. */

double fp_round(MPL *mpl, double x, double n)

{     double ten_to_n;

      if (n != floor(n))

         error(mpl, "round(%.*g, %.*g); non-integer second argument",

            DBL_DIG, x, DBL_DIG, n);

      if (n <= DBL_DIG + 2)

      {  ten_to_n = pow(10.0, n);

         if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)

         {  x = floor(x * ten_to_n + 0.5);

            if (x != 0.0) x /= ten_to_n;

         }

      }

      return x;

}

/*----------------------------------------------------------------------

-- fp_trunc - truncate floating-point value to n fractional digits.

--

-- This routine truncates given floating-point value x to n fractional

-- digits with the formula:

--

--                  ( floor(x * 10^n) / 10^n,  if x >= 0

--    trunc(x, n) = <

--                  ( ceil(x * 10^n) / 10^n,   if x < 0

--

-- The parameter n is assumed to be integer. */

double fp_trunc(MPL *mpl, double x, double n)

{     double ten_to_n;

      if (n != floor(n))

         error(mpl, "trunc(%.*g, %.*g); non-integer second argument",

            DBL_DIG, x, DBL_DIG, n);

      if (n <= DBL_DIG + 2)

      {  ten_to_n = pow(10.0, n);

         if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)

         {  x = (x >= 0.0 ? floor(x * ten_to_n) : ceil(x * ten_to_n));

            if (x != 0.0) x /= ten_to_n;

         }

      }

      return x;

}

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

/* * *              PSEUDO-RANDOM NUMBER GENERATORS               * * */

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

/*----------------------------------------------------------------------

-- fp_irand224 - pseudo-random integer in the range [0, 2^24).

--

-- This routine returns a next pseudo-random integer (converted to

-- floating-point) which is uniformly distributed between 0 and 2^24-1,

-- inclusive. */

#define two_to_the_24 0x1000000

double fp_irand224(MPL *mpl)

{     return

         (double)rng_unif_rand(mpl->rand, two_to_the_24);

}

/*----------------------------------------------------------------------

-- fp_uniform01 - pseudo-random number in the range [0, 1).

--

-- This routine returns a next pseudo-random number which is uniformly

-- distributed in the range [0, 1). */

#define two_to_the_31 ((unsigned int)0x80000000)

double fp_uniform01(MPL *mpl)

{     return

         (double)rng_next_rand(mpl->rand) / (double)two_to_the_31;

}

/*----------------------------------------------------------------------

-- fp_uniform - pseudo-random number in the range [a, b).

--

-- This routine returns a next pseudo-random number which is uniformly

-- distributed in the range [a, b). */

double fp_uniform(MPL *mpl, double a, double b)

{     double x;

      if (a >= b)

         error(mpl, "Uniform(%.*g, %.*g); invalid range",

            DBL_DIG, a, DBL_DIG, b);

      x = fp_uniform01(mpl);

#if 0

      x = a * (1.0 - x) + b * x;

#else

      x = fp_add(mpl, a * (1.0 - x), b * x);

#endif

      return x;

}

/*----------------------------------------------------------------------

-- fp_normal01 - Gaussian random variate with mu = 0 and sigma = 1.

--

-- This routine returns a Gaussian random variate with zero mean and

-- unit standard deviation. The polar (Box-Mueller) method is used.

--

-- This code is a modified version of the routine gsl_ran_gaussian from

-- the GNU Scientific Library Version 1.0. */

double fp_normal01(MPL *mpl)

{     double x, y, r2;

      do

      {  /* choose x, y in uniform square (-1,-1) to (+1,+1) */

         x = -1.0 + 2.0 * fp_uniform01(mpl);

         y = -1.0 + 2.0 * fp_uniform01(mpl);

         /* see if it is in the unit circle */

         r2 = x * x + y * y;

      } while (r2 > 1.0 || r2 == 0.0);

      /* Box-Muller transform */

      return y * sqrt(-2.0 * log (r2) / r2);

}

/*----------------------------------------------------------------------

-- fp_normal - Gaussian random variate with specified mu and sigma.

--

-- This routine returns a Gaussian random variate with mean mu and

-- standard deviation sigma. */

double fp_normal(MPL *mpl, double mu, double sigma)

{     double x;

#if 0

      x = mu + sigma * fp_normal01(mpl);

#else

      x = fp_add(mpl, mu, fp_mul(mpl, sigma, fp_normal01(mpl)));

#endif

      return x;

}

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

/* * *                SEGMENTED CHARACTER STRINGS                 * * */

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

/*----------------------------------------------------------------------

-- create_string - create character string.

--

-- This routine creates a segmented character string, which is exactly

-- equivalent to specified character string. */

STRING *create_string

(     MPL *mpl,

      char buf[MAX_LENGTH+1]  /* not changed */

)

#if 0

{     STRING *head, *tail;

      int i, j;

      xassert(buf != NULL);

      xassert(strlen(buf) <= MAX_LENGTH);

      head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));

      for (i = j = 0; ; i++)

      {  if ((tail->seg[j++] = buf[i]) == '\0') break;

         if (j == STRSEG_SIZE)

tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING))), j = 0;

      }

      tail->next = NULL;

      return head;

}

#else

{     STRING *str;

      xassert(strlen(buf) <= MAX_LENGTH);

      str = dmp_get_atom(mpl->strings, strlen(buf)+1);

      strcpy(str, buf);

      return str;

}

#endif

/*----------------------------------------------------------------------

-- copy_string - make copy of character string.

--

-- This routine returns an exact copy of segmented character string. */

STRING *copy_string

(     MPL *mpl,

      STRING *str             /* not changed */

)

#if 0

{     STRING *head, *tail;

      xassert(str != NULL);

      head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));

      for (; str != NULL; str = str->next)

      {  memcpy(tail->seg, str->seg, STRSEG_SIZE);

         if (str->next != NULL)

tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING)));

      }

      tail->next = NULL;

      return head;

}

#else

{     xassert(mpl == mpl);

      return create_string(mpl, str);

}

#endif

/*----------------------------------------------------------------------

-- compare_strings - compare one character string with another.

--

-- This routine compares one segmented character strings with another

-- and returns the result of comparison as follows:

--

-- = 0 - both strings are identical;

-- < 0 - the first string precedes the second one;

-- > 0 - the first string follows the second one. */

int compare_strings

(     MPL *mpl,

      STRING *str1,           /* not changed */

      STRING *str2            /* not changed */

)

#if 0

{     int j, c1, c2;

      xassert(mpl == mpl);

      for (;; str1 = str1->next, str2 = str2->next)

      {  xassert(str1 != NULL);

         xassert(str2 != NULL);

         for (j = 0; j < STRSEG_SIZE; j++)

         {  c1 = (unsigned char)str1->seg[j];

            c2 = (unsigned char)str2->seg[j];

            if (c1 < c2) return -1;

            if (c1 > c2) return +1;

            if (c1 == '\0') goto done;

         }

      }

done: return 0;

}

#else

{     xassert(mpl == mpl);

      return strcmp(str1, str2);

}

#endif

/*----------------------------------------------------------------------

-- fetch_string - extract content of character string.

--

-- This routine returns a character string, which is exactly equivalent

-- to specified segmented character string. */

char *fetch_string

(     MPL *mpl,

      STRING *str,            /* not changed */

      char buf[MAX_LENGTH+1]  /* modified */

)

#if 0

{     int i, j;

      xassert(mpl == mpl);

      xassert(buf != NULL);

      for (i = 0; ; str = str->next)

      {  xassert(str != NULL);

         for (j = 0; j < STRSEG_SIZE; j++)

            if ((buf[i++] = str->seg[j]) == '\0') goto done;

      }

done: xassert(strlen(buf) <= MAX_LENGTH);

      return buf;

}

#else

{     xassert(mpl == mpl);

      return strcpy(buf, str);

}

#endif

/*----------------------------------------------------------------------

-- delete_string - delete character string.

--

-- This routine deletes specified segmented character string. */

void delete_string

(     MPL *mpl,

      STRING *str             /* destroyed */

)

#if 0

{     STRING *temp;

      xassert(str != NULL);

      while (str != NULL)

      {  temp = str;

         str = str->next;

         dmp_free_atom(mpl->strings, temp, sizeof(STRING));

      }

      return;

}

#else

{     dmp_free_atom(mpl->strings, str, strlen(str)+1);

      return;

}

#endif

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

/* * *                          SYMBOLS                           * * */

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

/*----------------------------------------------------------------------

-- create_symbol_num - create symbol of numeric type.

--

-- This routine creates a symbol, which has a numeric value specified

-- as floating-point number. */

SYMBOL *create_symbol_num(MPL *mpl, double num)

{     SYMBOL *sym;

      sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));

      sym->num = num;

      sym->str = NULL;

      return sym;

}

/*----------------------------------------------------------------------

-- create_symbol_str - create symbol of abstract type.

--

-- This routine creates a symbol, which has an abstract value specified

-- as segmented character string. */

SYMBOL *create_symbol_str

(     MPL *mpl,

      STRING *str             /* destroyed */

)

{     SYMBOL *sym;

      xassert(str != NULL);

      sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));

      sym->num = 0.0;

      sym->str = str;

      return sym;

}

/*----------------------------------------------------------------------

-- copy_symbol - make copy of symbol.

--

-- This routine returns an exact copy of symbol. */

SYMBOL *copy_symbol

(     MPL *mpl,

      SYMBOL *sym             /* not changed */

)

{     SYMBOL *copy;

      xassert(sym != NULL);

      copy = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));

      if (sym->str == NULL)

      {  copy->num = sym->num;

         copy->str = NULL;

      }

      else

      {  copy->num = 0.0;

         copy->str = copy_string(mpl, sym->str);

      }

      return copy;

}

/*----------------------------------------------------------------------

-- compare_symbols - compare one symbol with another.

--

-- This routine compares one symbol with another and returns the result

-- of comparison as follows:

--

-- = 0 - both symbols are identical;

-- < 0 - the first symbol precedes the second one;

-- > 0 - the first symbol follows the second one.

--

-- Note that the linear order, in which symbols follow each other, is

-- implementation-dependent. It may be not an alphabetical order. */

int compare_symbols

(     MPL *mpl,

      SYMBOL *sym1,           /* not changed */

      SYMBOL *sym2            /* not changed */

)

{     xassert(sym1 != NULL);

      xassert(sym2 != NULL);

      /* let all numeric quantities precede all symbolic quantities */

      if (sym1->str == NULL && sym2->str == NULL)

      {  if (sym1->num < sym2->num) return -1;

         if (sym1->num > sym2->num) return +1;

         return 0;

      }

      if (sym1->str == NULL) return -1;

      if (sym2->str == NULL) return +1;

      return compare_strings(mpl, sym1->str, sym2->str);

}

/*----------------------------------------------------------------------

-- delete_symbol - delete symbol.

--

-- This routine deletes specified symbol. */

void delete_symbol

(     MPL *mpl,

      SYMBOL *sym             /* destroyed */

)

{     xassert(sym != NULL);

      if (sym->str != NULL) delete_string(mpl, sym->str);

      dmp_free_atom(mpl->symbols, sym, sizeof(SYMBOL));

      return;

}

/*----------------------------------------------------------------------

-- format_symbol - format symbol for displaying or printing.

--

-- This routine converts specified symbol to a charater string, which

-- is suitable for displaying or printing.

--

-- The resultant string is never longer than 255 characters. If it gets

-- longer, it is truncated from the right and appended by dots. */

char *format_symbol

(     MPL *mpl,

      SYMBOL *sym             /* not changed */

)

{     char *buf = mpl->sym_buf;

      xassert(sym != NULL);

      if (sym->str == NULL)

         sprintf(buf, "%.*g", DBL_DIG, sym->num);

      else

      {  char str[MAX_LENGTH+1];

         int quoted, j, len;

         fetch_string(mpl, sym->str, str);

         if (!(isalpha((unsigned char)str[0]) || str[0] == '_'))

            quoted = 1;

         else

         {  quoted = 0;

            for (j = 1; str[j] != '\0'; j++)

            {  if (!(isalnum((unsigned char)str[j]) ||

                     strchr("+-._", (unsigned char)str[j]) != NULL))

               {  quoted = 1;

                  break;

               }

            }

         }

#        define safe_append(c) \

            (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)

         buf[0] = '\0', len = 0;

         if (quoted) safe_append('\'');

         for (j = 0; str[j] != '\0'; j++)

         {  if (quoted && str[j] == '\'') safe_append('\'');

            safe_append(str[j]);

         }

         if (quoted) safe_append('\'');

#        undef safe_append

         buf[len] = '\0';

         if (len == 255) strcpy(buf+252, "...");

      }

      xassert(strlen(buf) <= 255);

      return buf;

}

/*----------------------------------------------------------------------

-- concat_symbols - concatenate one symbol with another.

--

-- This routine concatenates values of two given symbols and assigns

-- the resultant character string to a new symbol, which is returned on

-- exit. Both original symbols are destroyed. */

SYMBOL *concat_symbols

(     MPL *mpl,

      SYMBOL *sym1,           /* destroyed */

      SYMBOL *sym2            /* destroyed */

)

{     char str1[MAX_LENGTH+1], str2[MAX_LENGTH+1];

      xassert(MAX_LENGTH >= DBL_DIG + DBL_DIG);

      if (sym1->str == NULL)

         sprintf(str1, "%.*g", DBL_DIG, sym1->num);

      else

         fetch_string(mpl, sym1->str, str1);

      if (sym2->str == NULL)

         sprintf(str2, "%.*g", DBL_DIG, sym2->num);

      else

         fetch_string(mpl, sym2->str, str2);

      if (strlen(str1) + strlen(str2) > MAX_LENGTH)

      {  char buf[255+1];

         strcpy(buf, format_symbol(mpl, sym1));

         xassert(strlen(buf) < sizeof(buf));

         error(mpl, "%s & %s; resultant symbol exceeds %d characters",

            buf, format_symbol(mpl, sym2), MAX_LENGTH);

      }

      delete_symbol(mpl, sym1);

      delete_symbol(mpl, sym2);

      return create_symbol_str(mpl, create_string(mpl, strcat(str1,

         str2)));

}

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

/* * *                          N-TUPLES                          * * */

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

/*----------------------------------------------------------------------

-- create_tuple - create n-tuple.

--

-- This routine creates a n-tuple, which initially has no components,

-- i.e. which is 0-tuple. */

TUPLE *create_tuple(MPL *mpl)

{     TUPLE *tuple;

      xassert(mpl == mpl);

      tuple = NULL;

      return tuple;

}

/*----------------------------------------------------------------------

-- expand_tuple - append symbol to n-tuple.

--

-- This routine expands n-tuple appending to it a given symbol, which

-- becomes its new last component. */

TUPLE *expand_tuple

(     MPL *mpl,

      TUPLE *tuple,           /* destroyed */

      SYMBOL *sym             /* destroyed */

)

{     TUPLE *tail, *temp;

      xassert(sym != NULL);

      /* create a new component */

      tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));

      tail->sym = sym;

      tail->next = NULL;

      /* and append it to the component list */

      if (tuple == NULL)

         tuple = tail;

      else

      {  for (temp = tuple; temp->next != NULL; temp = temp->next);

         temp->next = tail;

      }

      return tuple;

}

/*----------------------------------------------------------------------

-- tuple_dimen - determine dimension of n-tuple.

--

-- This routine returns dimension of n-tuple, i.e. number of components

-- in the n-tuple. */

int tuple_dimen

(     MPL *mpl,

      TUPLE *tuple            /* not changed */

)

{     TUPLE *temp;

      int dim = 0;

      xassert(mpl == mpl);

      for (temp = tuple; temp != NULL; temp = temp->next) dim++;

      return dim;

}

/*----------------------------------------------------------------------

-- copy_tuple - make copy of n-tuple.

--

-- This routine returns an exact copy of n-tuple. */

TUPLE *copy_tuple

(     MPL *mpl,

      TUPLE *tuple            /* not changed */

)

{     TUPLE *head, *tail;

      if (tuple == NULL)

         head = NULL;

      else

      {  head = tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));

         for (; tuple != NULL; tuple = tuple->next)

         {  xassert(tuple->sym != NULL);

            tail->sym = copy_symbol(mpl, tuple->sym);

            if (tuple->next != NULL)

tail = (tail->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));

         }

         tail->next = NULL;

      }

      return head;

}

/*----------------------------------------------------------------------

-- compare_tuples - compare one n-tuple with another.

--

-- This routine compares two given n-tuples, which must have the same

-- dimension (not checked for the sake of efficiency), and returns one

-- of the following codes:

--

-- = 0 - both n-tuples are identical;

-- < 0 - the first n-tuple precedes the second one;

-- > 0 - the first n-tuple follows the second one.

--

-- Note that the linear order, in which n-tuples follow each other, is

-- implementation-dependent. It may be not an alphabetical order. */

int compare_tuples

(     MPL *mpl,

      TUPLE *tuple1,          /* not changed */

      TUPLE *tuple2           /* not changed */

)

{     TUPLE *item1, *item2;

      int ret;

      xassert(mpl == mpl);

      for (item1 = tuple1, item2 = tuple2; item1 != NULL;

           item1 = item1->next, item2 = item2->next)

      {  xassert(item2 != NULL);

         xassert(item1->sym != NULL);

         xassert(item2->sym != NULL);

         ret = compare_symbols(mpl, item1->sym, item2->sym);

         if (ret != 0) return ret;

      }

      xassert(item2 == NULL);

      return 0;

}

/*----------------------------------------------------------------------

-- build_subtuple - build subtuple of given n-tuple.

--

-- This routine builds subtuple, which consists of first dim components

-- of given n-tuple. */

TUPLE *build_subtuple

(     MPL *mpl,

      TUPLE *tuple,           /* not changed */

      int dim

)

{     TUPLE *head, *temp;

      int j;

      head = create_tuple(mpl);

      for (j = 1, temp = tuple; j <= dim; j++, temp = temp->next)

      {  xassert(temp != NULL);

         head = expand_tuple(mpl, head, copy_symbol(mpl, temp->sym));

      }

      return head;

}

/*----------------------------------------------------------------------

-- delete_tuple - delete n-tuple.

--

-- This routine deletes specified n-tuple. */

void delete_tuple

(     MPL *mpl,

      TUPLE *tuple            /* destroyed */

)

{     TUPLE *temp;

      while (tuple != NULL)

      {  temp = tuple;

         tuple = temp->next;

         xassert(temp->sym != NULL);

         delete_symbol(mpl, temp->sym);

         dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));

      }

      return;

}

/*----------------------------------------------------------------------

-- format_tuple - format n-tuple for displaying or printing.

--

-- This routine converts specified n-tuple to a character string, which

-- is suitable for displaying or printing.

--

-- The resultant string is never longer than 255 characters. If it gets

-- longer, it is truncated from the right and appended by dots. */

char *format_tuple

(     MPL *mpl,

      int c,

      TUPLE *tuple            /* not changed */

)

{     TUPLE *temp;

      int dim, j, len;

      char *buf = mpl->tup_buf, str[255+1], *save;

#     define safe_append(c) \

         (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)

      buf[0] = '\0', len = 0;

      dim = tuple_dimen(mpl, tuple);

      if (c == '[' && dim > 0) safe_append('[');

      if (c == '(' && dim > 1) safe_append('(');

      for (temp = tuple; temp != NULL; temp = temp->next)

      {  if (temp != tuple) safe_append(',');

         xassert(temp->sym != NULL);

         save = mpl->sym_buf;

         mpl->sym_buf = str;

         format_symbol(mpl, temp->sym);

         mpl->sym_buf = save;

         xassert(strlen(str) < sizeof(str));

         for (j = 0; str[j] != '\0'; j++) safe_append(str[j]);

      }

      if (c == '[' && dim > 0) safe_append(']');

      if (c == '(' && dim > 1) safe_append(')');

#     undef safe_append

      buf[len] = '\0';

      if (len == 255) strcpy(buf+252, "...");

      xassert(strlen(buf) <= 255);

      return buf;

}

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

/* * *                       ELEMENTAL SETS                       * * */

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

/*----------------------------------------------------------------------

-- create_elemset - create elemental set.

--

-- This routine creates an elemental set, whose members are n-tuples of

-- specified dimension. Being created the set is initially empty. */

ELEMSET *create_elemset(MPL *mpl, int dim)

{     ELEMSET *set;

      xassert(dim > 0);

      set = create_array(mpl, A_NONE, dim);

      return set;

}

/*----------------------------------------------------------------------

-- find_tuple - check if elemental set contains given n-tuple.

--

-- This routine finds given n-tuple in specified elemental set in order

-- to check if the set contains that n-tuple. If the n-tuple is found,

-- the routine returns pointer to corresponding array member. Otherwise

-- null pointer is returned. */

MEMBER *find_tuple

(     MPL *mpl,

      ELEMSET *set,           /* not changed */

      TUPLE *tuple            /* not changed */

)

{     xassert(set != NULL);

      xassert(set->type == A_NONE);

      xassert(set->dim == tuple_dimen(mpl, tuple));

      return find_member(mpl, set, tuple);

}

/*----------------------------------------------------------------------

-- add_tuple - add new n-tuple to elemental set.

--

-- This routine adds given n-tuple to specified elemental set.

--

-- For the sake of efficiency this routine doesn't check whether the

-- set already contains the same n-tuple or not. Therefore the calling

-- program should use the routine find_tuple (if necessary) in order to

-- make sure that the given n-tuple is not contained in the set, since

-- duplicate n-tuples within the same set are not allowed. */

MEMBER *add_tuple

(     MPL *mpl,

      ELEMSET *set,           /* modified */

      TUPLE *tuple            /* destroyed */

)

{     MEMBER *memb;