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

       xassert(set != NULL);

       xassert(set->type == A_NONE);

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

       memb = add_member(mpl, set, tuple);

       memb->value.none = NULL;

       return memb;

 }

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

 -- check_then_add - check and add new n-tuple to elemental set.

 --

 -- This routine is equivalent to the routine add_tuple except that it

 -- does check for duplicate n-tuples. */

 MEMBER *check_then_add

 (     MPL *mpl,

       ELEMSET *set,           /* modified */

       TUPLE *tuple            /* destroyed */

 )

 {     if (find_tuple(mpl, set, tuple) != NULL)

          error(mpl, "duplicate tuple %s detected", format_tuple(mpl,

             '(', tuple));

       return add_tuple(mpl, set, tuple);

 }

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

 -- copy_elemset - make copy of elemental set.

 --

 -- This routine makes an exact copy of elemental set. */

 ELEMSET *copy_elemset

 (     MPL *mpl,

       ELEMSET *set            /* not changed */

 )

 {     ELEMSET *copy;

       MEMBER *memb;

       xassert(set != NULL);

       xassert(set->type == A_NONE);

       xassert(set->dim > 0);

       copy = create_elemset(mpl, set->dim);

       for (memb = set->head; memb != NULL; memb = memb->next)

          add_tuple(mpl, copy, copy_tuple(mpl, memb->tuple));

       return copy;

 }

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

 -- delete_elemset - delete elemental set.

 --

 -- This routine deletes specified elemental set. */

 void delete_elemset

 (     MPL *mpl,

       ELEMSET *set            /* destroyed */

 )

 {     xassert(set != NULL);

       xassert(set->type == A_NONE);

       delete_array(mpl, set);

       return;

 }

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

 -- arelset_size - compute size of "arithmetic" elemental set.

 --

 -- This routine computes the size of "arithmetic" elemental set, which

 -- is specified in the form of arithmetic progression:

 --

 --    { t0 .. tf by dt }.

 --

 -- The size is computed using the formula:

 --

 --    n = max(0, floor((tf - t0) / dt) + 1). */

 int arelset_size(MPL *mpl, double t0, double tf, double dt)

 {     double temp;

       if (dt == 0.0)

          error(mpl, "%.*g .. %.*g by %.*g; zero stride not allowed",

             DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);

       if (tf > 0.0 && t0 < 0.0 && tf > + 0.999 * DBL_MAX + t0)

          temp = +DBL_MAX;

       else if (tf < 0.0 && t0 > 0.0 && tf < - 0.999 * DBL_MAX + t0)

          temp = -DBL_MAX;

       else

          temp = tf - t0;

       if (fabs(dt) < 1.0 && fabs(temp) > (0.999 * DBL_MAX) * fabs(dt))

       {  if (temp > 0.0 && dt > 0.0 || temp < 0.0 && dt < 0.0)

             temp = +DBL_MAX;

          else

             temp = 0.0;

       }

       else

       {  temp = floor(temp / dt) + 1.0;

          if (temp < 0.0) temp = 0.0;

       }

       xassert(temp >= 0.0);

       if (temp > (double)(INT_MAX - 1))

          error(mpl, "%.*g .. %.*g by %.*g; set too large",

             DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);

       return (int)(temp + 0.5);

 }

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

 -- arelset_member - compute member of "arithmetic" elemental set.

 --

 -- This routine returns a numeric value of symbol, which is equivalent

 -- to j-th member of given "arithmetic" elemental set specified in the

 -- form of arithmetic progression:

 --

 --    { t0 .. tf by dt }.

 --

 -- The symbol value is computed with the formula:

 --

 --    j-th member = t0 + (j - 1) * dt,

 --

 -- The number j must satisfy to the restriction 1 <= j <= n, where n is

 -- the set size computed by the routine arelset_size. */

 double arelset_member(MPL *mpl, double t0, double tf, double dt, int j)

 {     xassert(1 <= j && j <= arelset_size(mpl, t0, tf, dt));

       return t0 + (double)(j - 1) * dt;

 }

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

 -- create_arelset - create "arithmetic" elemental set.

 --

 -- This routine creates "arithmetic" elemental set, which is specified

 -- in the form of arithmetic progression:

 --

 --    { t0 .. tf by dt }.

 --

 -- Components of this set are 1-tuples. */

 ELEMSET *create_arelset(MPL *mpl, double t0, double tf, double dt)

 {     ELEMSET *set;

       int j, n;

       set = create_elemset(mpl, 1);

       n = arelset_size(mpl, t0, tf, dt);

       for (j = 1; j <= n; j++)

       {  add_tuple

          (  mpl,

             set,

             expand_tuple

             (  mpl,