lemon-project-template-glpk: deps/glpk/src/glpmpl03.c comparison

lemon-project-template-glpk

comparison deps/glpk/src/glpmpl03.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:77726e969f34
+/* 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, 2011 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,
+create_tuple(mpl),
+create_symbol_num
+(  mpl,
+arelset_member(mpl, t0, tf, dt, j)
+)
+)
+);
+}
+return set;
+}
+/*----------------------------------------------------------------------
+-- set_union - union of two elemental sets.
+--
+-- This routine computes the union:
+--
+--    X U Y = { j | (j in X) or (j in Y) },
+--
+-- where X and Y are given elemental sets (destroyed on exit). */
+ELEMSET *set_union
+(     MPL *mpl,
+ELEMSET *X,             /* destroyed */
+ELEMSET *Y              /* destroyed */
+)
+{     MEMBER *memb;
+xassert(X != NULL);
+xassert(X->type == A_NONE);
+xassert(X->dim > 0);
+xassert(Y != NULL);
+xassert(Y->type == A_NONE);
+xassert(Y->dim > 0);
+xassert(X->dim == Y->dim);
+for (memb = Y->head; memb != NULL; memb = memb->next)
+{  if (find_tuple(mpl, X, memb->tuple) == NULL)
+add_tuple(mpl, X, copy_tuple(mpl, memb->tuple));
+}
+delete_elemset(mpl, Y);
+return X;
+}
+/*----------------------------------------------------------------------
+-- set_diff - difference between two elemental sets.
+--
+-- This routine computes the difference:
+--
+--    X \ Y = { j | (j in X) and (j not in Y) },
+--
+-- where X and Y are given elemental sets (destroyed on exit). */
+ELEMSET *set_diff
+(     MPL *mpl,
+ELEMSET *X,             /* destroyed */
+ELEMSET *Y              /* destroyed */
+)
+{     ELEMSET *Z;
+MEMBER *memb;
+xassert(X != NULL);
+xassert(X->type == A_NONE);
+xassert(X->dim > 0);
+xassert(Y != NULL);
+xassert(Y->type == A_NONE);
+xassert(Y->dim > 0);
+xassert(X->dim == Y->dim);
+Z = create_elemset(mpl, X->dim);
+for (memb = X->head; memb != NULL; memb = memb->next)
+{  if (find_tuple(mpl, Y, memb->tuple) == NULL)
+add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
+}
+delete_elemset(mpl, X);
+delete_elemset(mpl, Y);
+return Z;
+}
+/*----------------------------------------------------------------------
+-- set_symdiff - symmetric difference between two elemental sets.
+--
+-- This routine computes the symmetric difference:
+--
+--    X (+) Y = (X \ Y) U (Y \ X),
+--
+-- where X and Y are given elemental sets (destroyed on exit). */
+ELEMSET *set_symdiff
+(     MPL *mpl,
+ELEMSET *X,             /* destroyed */
+ELEMSET *Y              /* destroyed */
+)
+{     ELEMSET *Z;
+MEMBER *memb;
+xassert(X != NULL);
+xassert(X->type == A_NONE);
+xassert(X->dim > 0);
+xassert(Y != NULL);
+xassert(Y->type == A_NONE);
+xassert(Y->dim > 0);
+xassert(X->dim == Y->dim);
+/* Z := X \ Y */
+Z = create_elemset(mpl, X->dim);
+for (memb = X->head; memb != NULL; memb = memb->next)
+{  if (find_tuple(mpl, Y, memb->tuple) == NULL)
+add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
+}
+/* Z := Z U (Y \ X) */
+for (memb = Y->head; memb != NULL; memb = memb->next)
+{  if (find_tuple(mpl, X, memb->tuple) == NULL)
+add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
+}
+delete_elemset(mpl, X);
+delete_elemset(mpl, Y);
+return Z;
+}
+/*----------------------------------------------------------------------
+-- set_inter - intersection of two elemental sets.
+--
+-- This routine computes the intersection:
+--
+--    X ^ Y = { j | (j in X) and (j in Y) },
+--
+-- where X and Y are given elemental sets (destroyed on exit). */
+ELEMSET *set_inter
+(     MPL *mpl,
+ELEMSET *X,             /* destroyed */
+ELEMSET *Y              /* destroyed */
+)
+{     ELEMSET *Z;
+MEMBER *memb;
+xassert(X != NULL);
+xassert(X->type == A_NONE);
+xassert(X->dim > 0);
+xassert(Y != NULL);
+xassert(Y->type == A_NONE);
+xassert(Y->dim > 0);
+xassert(X->dim == Y->dim);
+Z = create_elemset(mpl, X->dim);
+for (memb = X->head; memb != NULL; memb = memb->next)
+{  if (find_tuple(mpl, Y, memb->tuple) != NULL)
+add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
+}
+delete_elemset(mpl, X);
+delete_elemset(mpl, Y);
+return Z;
+}
+/*----------------------------------------------------------------------
+-- set_cross - cross (Cartesian) product of two elemental sets.
+--
+-- This routine computes the cross (Cartesian) product:
+--
+--    X x Y = { (i,j) | (i in X) and (j in Y) },
+--
+-- where X and Y are given elemental sets (destroyed on exit). */
+ELEMSET *set_cross
+(     MPL *mpl,
+ELEMSET *X,             /* destroyed */
+ELEMSET *Y              /* destroyed */
+)
+{     ELEMSET *Z;
+MEMBER *memx, *memy;
+TUPLE *tuple, *temp;
+xassert(X != NULL);
+xassert(X->type == A_NONE);
+xassert(X->dim > 0);
+xassert(Y != NULL);
+xassert(Y->type == A_NONE);
+xassert(Y->dim > 0);
+Z = create_elemset(mpl, X->dim + Y->dim);
+for (memx = X->head; memx != NULL; memx = memx->next)
+{  for (memy = Y->head; memy != NULL; memy = memy->next)
+{  tuple = copy_tuple(mpl, memx->tuple);
+for (temp = memy->tuple; temp != NULL; temp = temp->next)
+tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,
+temp->sym));
+add_tuple(mpl, Z, tuple);
+}
+}
+delete_elemset(mpl, X);
+delete_elemset(mpl, Y);
+return Z;
+}
+/**********************************************************************/
+/* * *                    ELEMENTAL VARIABLES                     * * */
+/**********************************************************************/
+/* (there are no specific routines for elemental variables) */
+/**********************************************************************/
+/* * *                        LINEAR FORMS                        * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- constant_term - create constant term.
+--
+-- This routine creates the linear form, which is a constant term. */
+FORMULA *constant_term(MPL *mpl, double coef)
+{     FORMULA *form;
+if (coef == 0.0)
+form = NULL;
+else
+{  form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
+form->coef = coef;
+form->var = NULL;
+form->next = NULL;
+}
+return form;
+}
+/*----------------------------------------------------------------------
+-- single_variable - create single variable.
+--
+-- This routine creates the linear form, which is a single elemental
+-- variable. */
+FORMULA *single_variable
+(     MPL *mpl,
+ELEMVAR *var            /* referenced */
+)
+{     FORMULA *form;
+xassert(var != NULL);
+form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
+form->coef = 1.0;
+form->var = var;
+form->next = NULL;
+return form;
+}
+/*----------------------------------------------------------------------
+-- copy_formula - make copy of linear form.
+--
+-- This routine returns an exact copy of linear form. */
+FORMULA *copy_formula
+(     MPL *mpl,
+FORMULA *form           /* not changed */
+)
+{     FORMULA *head, *tail;
+if (form == NULL)
+head = NULL;
+else
+{  head = tail = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
+for (; form != NULL; form = form->next)
+{  tail->coef = form->coef;
+tail->var = form->var;
+if (form->next != NULL)
+tail = (tail->next = dmp_get_atom(mpl->formulae, sizeof(FORMULA)));
+}
+tail->next = NULL;
+}
+return head;
+}
+/*----------------------------------------------------------------------
+-- delete_formula - delete linear form.
+--
+-- This routine deletes specified linear form. */
+void delete_formula
+(     MPL *mpl,
+FORMULA *form           /* destroyed */
+)
+{     FORMULA *temp;
+while (form != NULL)
+{  temp = form;
+form = form->next;
+dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- linear_comb - linear combination of two linear forms.
+--
+-- This routine computes the linear combination:
+--
+--    a * fx + b * fy,
+--
+-- where a and b are numeric coefficients, fx and fy are linear forms
+-- (destroyed on exit). */
+FORMULA *linear_comb
+(     MPL *mpl,
+double a, FORMULA *fx,  /* destroyed */
+double b, FORMULA *fy   /* destroyed */
+)
+{     FORMULA *form = NULL, *term, *temp;
+double c0 = 0.0;
+for (term = fx; term != NULL; term = term->next)
+{  if (term->var == NULL)
+c0 = fp_add(mpl, c0, fp_mul(mpl, a, term->coef));
+else
+term->var->temp =
+fp_add(mpl, term->var->temp, fp_mul(mpl, a, term->coef));
+}
+for (term = fy; term != NULL; term = term->next)
+{  if (term->var == NULL)
+c0 = fp_add(mpl, c0, fp_mul(mpl, b, term->coef));
+else
+term->var->temp =
+fp_add(mpl, term->var->temp, fp_mul(mpl, b, term->coef));
+}
+for (term = fx; term != NULL; term = term->next)
+{  if (term->var != NULL && term->var->temp != 0.0)
+{  temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
+temp->coef = term->var->temp, temp->var = term->var;
+temp->next = form, form = temp;
+term->var->temp = 0.0;
+}
+}
+for (term = fy; term != NULL; term = term->next)
+{  if (term->var != NULL && term->var->temp != 0.0)
+{  temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
+temp->coef = term->var->temp, temp->var = term->var;
+temp->next = form, form = temp;
+term->var->temp = 0.0;
+}
+}
+if (c0 != 0.0)
+{  temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
+temp->coef = c0, temp->var = NULL;
+temp->next = form, form = temp;
+}
+delete_formula(mpl, fx);
+delete_formula(mpl, fy);
+return form;
+}
+/*----------------------------------------------------------------------
+-- remove_constant - remove constant term from linear form.
+--
+-- This routine removes constant term from linear form and stores its
+-- value to given location. */
+FORMULA *remove_constant
+(     MPL *mpl,
+FORMULA *form,          /* destroyed */
+double *coef            /* modified */
+)
+{     FORMULA *head = NULL, *temp;
+*coef = 0.0;
+while (form != NULL)
+{  temp = form;
+form = form->next;
+if (temp->var == NULL)
+{  /* constant term */
+*coef = fp_add(mpl, *coef, temp->coef);
+dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));
+}
+else
+{  /* linear term */
+temp->next = head;
+head = temp;
+}
+}
+return head;
+}
+/*----------------------------------------------------------------------
+-- reduce_terms - reduce identical terms in linear form.
+--
+-- This routine reduces identical terms in specified linear form. */
+FORMULA *reduce_terms
+(     MPL *mpl,
+FORMULA *form           /* destroyed */
+)
+{     FORMULA *term, *next_term;
+double c0 = 0.0;
+for (term = form; term != NULL; term = term->next)
+{  if (term->var == NULL)
+c0 = fp_add(mpl, c0, term->coef);
+else
+term->var->temp = fp_add(mpl, term->var->temp, term->coef);
+}
+next_term = form, form = NULL;
+for (term = next_term; term != NULL; term = next_term)
+{  next_term = term->next;
+if (term->var == NULL && c0 != 0.0)
+{  term->coef = c0, c0 = 0.0;
+term->next = form, form = term;
+}
+else if (term->var != NULL && term->var->temp != 0.0)
+{  term->coef = term->var->temp, term->var->temp = 0.0;
+term->next = form, form = term;
+}
+else
+dmp_free_atom(mpl->formulae, term, sizeof(FORMULA));
+}
+return form;
+}
+/**********************************************************************/
+/* * *                   ELEMENTAL CONSTRAINTS                    * * */
+/**********************************************************************/
+/* (there are no specific routines for elemental constraints) */
+/**********************************************************************/
+/* * *                       GENERIC VALUES                       * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- delete_value - delete generic value.
+--
+-- This routine deletes specified generic value.
+--
+-- NOTE: The generic value to be deleted must be valid. */
+void delete_value
+(     MPL *mpl,
+int type,
+VALUE *value            /* content destroyed */
+)
+{     xassert(value != NULL);
+switch (type)
+{  case A_NONE:
+value->none = NULL;
+break;
+case A_NUMERIC:
+value->num = 0.0;
+break;
+case A_SYMBOLIC:
+delete_symbol(mpl, value->sym), value->sym = NULL;
+break;
+case A_LOGICAL:
+value->bit = 0;
+break;
+case A_TUPLE:
+delete_tuple(mpl, value->tuple), value->tuple = NULL;
+break;
+case A_ELEMSET:
+delete_elemset(mpl, value->set), value->set = NULL;
+break;
+case A_ELEMVAR:
+value->var = NULL;
+break;
+case A_FORMULA:
+delete_formula(mpl, value->form), value->form = NULL;
+break;
+case A_ELEMCON:
+value->con = NULL;
+break;
+default:
+xassert(type != type);
+}
+return;
+}
+/**********************************************************************/
+/* * *                SYMBOLICALLY INDEXED ARRAYS                 * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- create_array - create array.
+--
+-- This routine creates an array of specified type and dimension. Being
+-- created the array is initially empty.
+--
+-- The type indicator determines generic values, which can be assigned
+-- to the array members:
+--
+-- A_NONE     - none (members have no assigned values)
+-- A_NUMERIC  - floating-point numbers
+-- A_SYMBOLIC - symbols
+-- A_ELEMSET  - elemental sets
+-- A_ELEMVAR  - elemental variables
+-- A_ELEMCON  - elemental constraints
+--
+-- The dimension may be 0, in which case the array consists of the only
+-- member (such arrays represent 0-dimensional objects). */
+ARRAY *create_array(MPL *mpl, int type, int dim)
+{     ARRAY *array;
+xassert(type == A_NONE || type == A_NUMERIC ||
+type == A_SYMBOLIC || type == A_ELEMSET ||
+type == A_ELEMVAR || type == A_ELEMCON);
+xassert(dim >= 0);
+array = dmp_get_atom(mpl->arrays, sizeof(ARRAY));
+array->type = type;
+array->dim = dim;
+array->size = 0;
+array->head = NULL;
+array->tail = NULL;
+array->tree = NULL;
+array->prev = NULL;
+array->next = mpl->a_list;
+/* include the array in the global array list */
+if (array->next != NULL) array->next->prev = array;
+mpl->a_list = array;
+return array;
+}
+/*----------------------------------------------------------------------
+-- find_member - find array member with given n-tuple.
+--
+-- This routine finds an array member, which has given n-tuple. If the
+-- array is short, the linear search is used. Otherwise the routine
+-- autimatically creates the search tree (i.e. the array index) to find
+-- members for logarithmic time. */
+static int compare_member_tuples(void *info, const void *key1,
+const void *key2)
+{     /* this is an auxiliary routine used to compare keys, which are
+n-tuples assigned to array members */
+return compare_tuples((MPL *)info, (TUPLE *)key1, (TUPLE *)key2);
+}
+MEMBER *find_member
+(     MPL *mpl,
+ARRAY *array,           /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     MEMBER *memb;
+xassert(array != NULL);
+/* the n-tuple must have the same dimension as the array */
+xassert(tuple_dimen(mpl, tuple) == array->dim);
+/* if the array is large enough, create the search tree and index
+all existing members of the array */
+if (array->size > 30 && array->tree == NULL)
+{  array->tree = avl_create_tree(compare_member_tuples, mpl);
+for (memb = array->head; memb != NULL; memb = memb->next)
+avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
+(void *)memb);
+}
+/* find a member, which has the given tuple */
+if (array->tree == NULL)
+{  /* the search tree doesn't exist; use the linear search */
+for (memb = array->head; memb != NULL; memb = memb->next)
+if (compare_tuples(mpl, memb->tuple, tuple) == 0) break;
+}
+else
+{  /* the search tree exists; use the binary search */
+AVLNODE *node;
+node = avl_find_node(array->tree, tuple);
+memb = (MEMBER *)(node == NULL ? NULL : avl_get_node_link(node));
+}
+return memb;
+}
+/*----------------------------------------------------------------------
+-- add_member - add new member to array.
+--
+-- This routine creates a new member with given n-tuple and adds it to
+-- specified array.
+--
+-- For the sake of efficiency this routine doesn't check whether the
+-- array already contains a member with the given n-tuple or not. Thus,
+-- if necessary, the calling program should use the routine find_member
+-- in order to be sure that the array contains no member with the same
+-- n-tuple, because members with duplicate n-tuples are not allowed.
+--
+-- This routine assigns no generic value to the new member, because the
+-- calling program must do that. */
+MEMBER *add_member
+(     MPL *mpl,
+ARRAY *array,           /* modified */
+TUPLE *tuple            /* destroyed */
+)
+{     MEMBER *memb;
+xassert(array != NULL);
+/* the n-tuple must have the same dimension as the array */
+xassert(tuple_dimen(mpl, tuple) == array->dim);
+/* create new member */
+memb = dmp_get_atom(mpl->members, sizeof(MEMBER));
+memb->tuple = tuple;
+memb->next = NULL;
+memset(&memb->value, '?', sizeof(VALUE));
+/* and append it to the member list */
+array->size++;
+if (array->head == NULL)
+array->head = memb;
+else
+array->tail->next = memb;
+array->tail = memb;
+/* if the search tree exists, index the new member */
+if (array->tree != NULL)
+avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
+(void *)memb);
+return memb;
+}
+/*----------------------------------------------------------------------
+-- delete_array - delete array.
+--
+-- This routine deletes specified array.
+--
+-- Generic values assigned to the array members are not deleted by this
+-- routine. The calling program itself must delete all assigned generic
+-- values before deleting the array. */
+void delete_array
+(     MPL *mpl,
+ARRAY *array            /* destroyed */
+)
+{     MEMBER *memb;
+xassert(array != NULL);
+/* delete all existing array members */
+while (array->head != NULL)
+{  memb = array->head;
+array->head = memb->next;
+delete_tuple(mpl, memb->tuple);
+dmp_free_atom(mpl->members, memb, sizeof(MEMBER));
+}
+/* if the search tree exists, also delete it */
+if (array->tree != NULL) avl_delete_tree(array->tree);
+/* remove the array from the global array list */
+if (array->prev == NULL)
+mpl->a_list = array->next;
+else
+array->prev->next = array->next;
+if (array->next == NULL)
+;
+else
+array->next->prev = array->prev;
+/* delete the array descriptor */
+dmp_free_atom(mpl->arrays, array, sizeof(ARRAY));
+return;
+}
+/**********************************************************************/
+/* * *                 DOMAINS AND DUMMY INDICES                  * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- assign_dummy_index - assign new value to dummy index.
+--
+-- This routine assigns new value to specified dummy index and, that is
+-- important, invalidates all temporary resultant values, which depends
+-- on that dummy index. */
+void assign_dummy_index
+(     MPL *mpl,
+DOMAIN_SLOT *slot,      /* modified */
+SYMBOL *value           /* not changed */
+)
+{     CODE *leaf, *code;
+xassert(slot != NULL);
+xassert(value != NULL);
+/* delete the current value assigned to the dummy index */
+if (slot->value != NULL)
+{  /* if the current value and the new one are identical, actual
+assignment is not needed */
+if (compare_symbols(mpl, slot->value, value) == 0) goto done;
+/* delete a symbol, which is the current value */
+delete_symbol(mpl, slot->value), slot->value = NULL;
+}
+/* now walk through all the pseudo-codes with op = O_INDEX, which
+refer to the dummy index to be changed (these pseudo-codes are
+leaves in the forest of *all* expressions in the database) */
+for (leaf = slot->list; leaf != NULL; leaf = leaf->arg.index.
+next)
+{  xassert(leaf->op == O_INDEX);
+/* invalidate all resultant values, which depend on the dummy
+index, walking from the current leaf toward the root of the
+corresponding expression tree */
+for (code = leaf; code != NULL; code = code->up)
+{  if (code->valid)
+{  /* invalidate and delete resultant value */
+code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+}
+}
+/* assign new value to the dummy index */
+slot->value = copy_symbol(mpl, value);
+done: return;
+}
+/*----------------------------------------------------------------------
+-- update_dummy_indices - update current values of dummy indices.
+--
+-- This routine assigns components of "backup" n-tuple to dummy indices
+-- of specified domain block. If no "backup" n-tuple is defined for the
+-- domain block, values of the dummy indices remain untouched. */
+void update_dummy_indices
+(     MPL *mpl,
+DOMAIN_BLOCK *block     /* not changed */
+)
+{     DOMAIN_SLOT *slot;
+TUPLE *temp;
+if (block->backup != NULL)
+{  for (slot = block->list, temp = block->backup; slot != NULL;
+slot = slot->next, temp = temp->next)
+{  xassert(temp != NULL);
+xassert(temp->sym != NULL);
+assign_dummy_index(mpl, slot, temp->sym);
+}
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- enter_domain_block - enter domain block.
+--
+-- Let specified domain block have the form:
+--
+--    { ..., (j1, j2, ..., jn) in J, ... }
+--
+-- where j1, j2, ..., jn are dummy indices, J is a basic set.
+--
+-- This routine does the following:
+--
+-- 1. Checks if the given n-tuple is a member of the basic set J. Note
+--    that J being *out of the scope* of the domain block cannot depend
+--    on the dummy indices in the same and inner domain blocks, so it
+--    can be computed before the dummy indices are assigned new values.
+--    If this check fails, the routine returns with non-zero code.
+--
+-- 2. Saves current values of the dummy indices j1, j2, ..., jn.
+--
+-- 3. Assigns new values, which are components of the given n-tuple, to
+--    the dummy indices j1, j2, ..., jn. If dimension of the n-tuple is
+--    larger than n, its extra components n+1, n+2, ... are not used.
+--
+-- 4. Calls the formal routine func which either enters the next domain
+--    block or evaluates some code within the domain scope.
+--
+-- 5. Restores former values of the dummy indices j1, j2, ..., jn.
+--
+-- Since current values assigned to the dummy indices on entry to this
+-- routine are restored on exit, the formal routine func is allowed to
+-- call this routine recursively. */
+int enter_domain_block
+(     MPL *mpl,
+DOMAIN_BLOCK *block,    /* not changed */
+TUPLE *tuple,           /* not changed */
+void *info, void (*func)(MPL *mpl, void *info)
+)
+{     TUPLE *backup;
+int ret = 0;
+/* check if the given n-tuple is a member of the basic set */
+xassert(block->code != NULL);
+if (!is_member(mpl, block->code, tuple))
+{  ret = 1;
+goto done;
+}
+/* save reference to "backup" n-tuple, which was used to assign
+current values of the dummy indices (it is sufficient to save
+reference, not value, because that n-tuple is defined in some
+outer level of recursion and therefore cannot be changed on
+this and deeper recursive calls) */
+backup = block->backup;
+/* set up new "backup" n-tuple, which defines new values of the
+dummy indices */
+block->backup = tuple;
+/* assign new values to the dummy indices */
+update_dummy_indices(mpl, block);
+/* call the formal routine that does the rest part of the job */
+func(mpl, info);
+/* restore reference to the former "backup" n-tuple */
+block->backup = backup;
+/* restore former values of the dummy indices; note that if the
+domain block just escaped has no other active instances which
+may exist due to recursion (it is indicated by a null pointer
+to the former n-tuple), former values of the dummy indices are
+undefined; therefore in this case the routine keeps currently
+assigned values of the dummy indices that involves keeping all
+dependent temporary results and thereby, if this domain block
+is not used recursively, allows improving efficiency */
+update_dummy_indices(mpl, block);
+done: return ret;
+}
+/*----------------------------------------------------------------------
+-- eval_within_domain - perform evaluation within domain scope.
+--
+-- This routine assigns new values (symbols) to all dummy indices of
+-- specified domain and calls the formal routine func, which is used to
+-- evaluate some code in the domain scope. Each free dummy index in the
+-- domain is assigned a value specified in the corresponding component
+-- of given n-tuple. Non-free dummy indices are assigned values, which
+-- are computed by this routine.
+--
+-- Number of components in the given n-tuple must be the same as number
+-- of free indices in the domain.
+--
+-- If the given n-tuple is not a member of the domain set, the routine
+-- func is not called, and non-zero code is returned.
+--
+-- For the sake of convenience it is allowed to specify domain as NULL
+-- (then n-tuple also must be 0-tuple, i.e. empty), in which case this
+-- routine just calls the routine func and returns zero.
+--
+-- This routine allows recursive calls from the routine func providing
+-- correct values of dummy indices for each instance.
+--
+-- NOTE: The n-tuple passed to this routine must not be changed by any
+--       other routines called from the formal routine func until this
+--       routine has returned. */
+struct eval_domain_info
+{     /* working info used by the routine eval_within_domain */
+DOMAIN *domain;
+/* domain, which has to be entered */
+DOMAIN_BLOCK *block;
+/* domain block, which is currently processed */
+TUPLE *tuple;
+/* tail of original n-tuple, whose components have to be assigned
+to free dummy indices in the current domain block */
+void *info;
+/* transit pointer passed to the formal routine func */
+void (*func)(MPL *mpl, void *info);
+/* routine, which has to be executed in the domain scope */
+int failure;
+/* this flag indicates that given n-tuple is not a member of the
+domain set */
+};
+static void eval_domain_func(MPL *mpl, void *_my_info)
+{     /* this routine recursively enters into the domain scope and then
+calls the routine func */
+struct eval_domain_info *my_info = _my_info;
+if (my_info->block != NULL)
+{  /* the current domain block to be entered exists */
+DOMAIN_BLOCK *block;
+DOMAIN_SLOT *slot;
+TUPLE *tuple = NULL, *temp = NULL;
+/* save pointer to the current domain block */
+block = my_info->block;
+/* and get ready to enter the next block (if it exists) */
+my_info->block = block->next;
+/* construct temporary n-tuple, whose components correspond to
+dummy indices (slots) of the current domain; components of
+the temporary n-tuple that correspond to free dummy indices
+are assigned references (not values!) to symbols specified
+in the corresponding components of the given n-tuple, while
+other components that correspond to non-free dummy indices
+are assigned symbolic values computed here */
+for (slot = block->list; slot != NULL; slot = slot->next)
+{  /* create component that corresponds to the current slot */
+if (tuple == NULL)
+tuple = temp = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
+else
+temp = (temp->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));
+if (slot->code == NULL)
+{  /* dummy index is free; take reference to symbol, which
+is specified in the corresponding component of given
+n-tuple */
+xassert(my_info->tuple != NULL);
+temp->sym = my_info->tuple->sym;
+xassert(temp->sym != NULL);
+my_info->tuple = my_info->tuple->next;
+}
+else
+{  /* dummy index is non-free; compute symbolic value to be
+temporarily assigned to the dummy index */
+temp->sym = eval_symbolic(mpl, slot->code);
+}
+}
+temp->next = NULL;
+/* enter the current domain block */
+if (enter_domain_block(mpl, block, tuple, my_info,
+eval_domain_func)) my_info->failure = 1;
+/* delete temporary n-tuple as well as symbols that correspond
+to non-free dummy indices (they were computed here) */
+for (slot = block->list; slot != NULL; slot = slot->next)
+{  xassert(tuple != NULL);
+temp = tuple;
+tuple = tuple->next;
+if (slot->code != NULL)
+{  /* dummy index is non-free; delete symbolic value */
+delete_symbol(mpl, temp->sym);
+}
+/* delete component that corresponds to the current slot */
+dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));
+}
+}
+else
+{  /* there are no more domain blocks, i.e. we have reached the
+domain scope */
+xassert(my_info->tuple == NULL);
+/* check optional predicate specified for the domain */
+if (my_info->domain->code != NULL && !eval_logical(mpl,
+my_info->domain->code))
+{  /* the predicate is false */
+my_info->failure = 2;
+}
+else
+{  /* the predicate is true; do the job */
+my_info->func(mpl, my_info->info);
+}
+}
+return;
+}
+int eval_within_domain
+(     MPL *mpl,
+DOMAIN *domain,         /* not changed */
+TUPLE *tuple,           /* not changed */
+void *info, void (*func)(MPL *mpl, void *info)
+)
+{     /* this routine performs evaluation within domain scope */
+struct eval_domain_info _my_info, *my_info = &_my_info;
+if (domain == NULL)
+{  xassert(tuple == NULL);
+func(mpl, info);
+my_info->failure = 0;
+}
+else
+{  xassert(tuple != NULL);
+my_info->domain = domain;
+my_info->block = domain->list;
+my_info->tuple = tuple;
+my_info->info = info;
+my_info->func = func;
+my_info->failure = 0;
+/* enter the very first domain block */
+eval_domain_func(mpl, my_info);
+}
+return my_info->failure;
+}
+/*----------------------------------------------------------------------
+-- loop_within_domain - perform iterations within domain scope.
+--
+-- This routine iteratively assigns new values (symbols) to the dummy
+-- indices of specified domain by enumerating all n-tuples, which are
+-- members of the domain set, and for every n-tuple it calls the formal
+-- routine func to evaluate some code within the domain scope.
+--
+-- If the routine func returns non-zero, enumeration within the domain
+-- is prematurely terminated.
+--
+-- For the sake of convenience it is allowed to specify domain as NULL,
+-- in which case this routine just calls the routine func only once and
+-- returns zero.
+--
+-- This routine allows recursive calls from the routine func providing
+-- correct values of dummy indices for each instance. */
+struct loop_domain_info
+{     /* working info used by the routine loop_within_domain */
+DOMAIN *domain;
+/* domain, which has to be entered */
+DOMAIN_BLOCK *block;
+/* domain block, which is currently processed */
+int looping;
+/* clearing this flag leads to terminating enumeration */
+void *info;
+/* transit pointer passed to the formal routine func */
+int (*func)(MPL *mpl, void *info);
+/* routine, which needs to be executed in the domain scope */
+};
+static void loop_domain_func(MPL *mpl, void *_my_info)
+{     /* this routine enumerates all n-tuples in the basic set of the
+current domain block, enters recursively into the domain scope
+for every n-tuple, and then calls the routine func */
+struct loop_domain_info *my_info = _my_info;
+if (my_info->block != NULL)
+{  /* the current domain block to be entered exists */
+DOMAIN_BLOCK *block;
+DOMAIN_SLOT *slot;
+TUPLE *bound;
+/* save pointer to the current domain block */
+block = my_info->block;
+/* and get ready to enter the next block (if it exists) */
+my_info->block = block->next;
+/* compute symbolic values, at which non-free dummy indices of
+the current domain block are bound; since that values don't
+depend on free dummy indices of the current block, they can
+be computed once out of the enumeration loop */
+bound = create_tuple(mpl);
+for (slot = block->list; slot != NULL; slot = slot->next)
+{  if (slot->code != NULL)
+bound = expand_tuple(mpl, bound, eval_symbolic(mpl,
+slot->code));
+}
+/* start enumeration */
+xassert(block->code != NULL);
+if (block->code->op == O_DOTS)
+{  /* the basic set is "arithmetic", in which case it doesn't
+need to be computed explicitly */
+TUPLE *tuple;
+int n, j;
+double t0, tf, dt;
+/* compute "parameters" of the basic set */
+t0 = eval_numeric(mpl, block->code->arg.arg.x);
+tf = eval_numeric(mpl, block->code->arg.arg.y);
+if (block->code->arg.arg.z == NULL)
+dt = 1.0;
+else
+dt = eval_numeric(mpl, block->code->arg.arg.z);
+/* determine cardinality of the basic set */
+n = arelset_size(mpl, t0, tf, dt);
+/* create dummy 1-tuple for members of the basic set */
+tuple = expand_tuple(mpl, create_tuple(mpl),
+create_symbol_num(mpl, 0.0));
+/* in case of "arithmetic" set there is exactly one dummy
+index, which cannot be non-free */
+xassert(bound == NULL);
+/* walk through 1-tuples of the basic set */
+for (j = 1; j <= n && my_info->looping; j++)
+{  /* construct dummy 1-tuple for the current member */
+tuple->sym->num = arelset_member(mpl, t0, tf, dt, j);
+/* enter the current domain block */
+enter_domain_block(mpl, block, tuple, my_info,
+loop_domain_func);
+}
+/* delete dummy 1-tuple */
+delete_tuple(mpl, tuple);
+}
+else
+{  /* the basic set is of general kind, in which case it needs
+to be explicitly computed */
+ELEMSET *set;
+MEMBER *memb;
+TUPLE *temp1, *temp2;
+/* compute the basic set */
+set = eval_elemset(mpl, block->code);
+/* walk through all n-tuples of the basic set */
+for (memb = set->head; memb != NULL && my_info->looping;
+memb = memb->next)
+{  /* all components of the current n-tuple that correspond
+to non-free dummy indices must be feasible; otherwise
+the n-tuple is not in the basic set */
+temp1 = memb->tuple;
+temp2 = bound;
+for (slot = block->list; slot != NULL; slot = slot->next)
+{  xassert(temp1 != NULL);
+if (slot->code != NULL)
+{  /* non-free dummy index */
+xassert(temp2 != NULL);
+if (compare_symbols(mpl, temp1->sym, temp2->sym)
+!= 0)
+{  /* the n-tuple is not in the basic set */
+goto skip;
+}
+temp2 = temp2->next;
+}
+temp1 = temp1->next;
+}
+xassert(temp1 == NULL);
+xassert(temp2 == NULL);
+/* enter the current domain block */
+enter_domain_block(mpl, block, memb->tuple, my_info,
+loop_domain_func);
+skip:          ;
+}
+/* delete the basic set */
+delete_elemset(mpl, set);
+}
+/* delete symbolic values binding non-free dummy indices */
+delete_tuple(mpl, bound);
+/* restore pointer to the current domain block */
+my_info->block = block;
+}
+else
+{  /* there are no more domain blocks, i.e. we have reached the
+domain scope */
+/* check optional predicate specified for the domain */
+if (my_info->domain->code != NULL && !eval_logical(mpl,
+my_info->domain->code))
+{  /* the predicate is false */
+/* nop */;
+}
+else
+{  /* the predicate is true; do the job */
+my_info->looping = !my_info->func(mpl, my_info->info);
+}
+}
+return;
+}
+void loop_within_domain
+(     MPL *mpl,
+DOMAIN *domain,         /* not changed */
+void *info, int (*func)(MPL *mpl, void *info)
+)
+{     /* this routine performs iterations within domain scope */
+struct loop_domain_info _my_info, *my_info = &_my_info;
+if (domain == NULL)
+func(mpl, info);
+else
+{  my_info->domain = domain;
+my_info->block = domain->list;
+my_info->looping = 1;
+my_info->info = info;
+my_info->func = func;
+/* enter the very first domain block */
+loop_domain_func(mpl, my_info);
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- out_of_domain - raise domain exception.
+--
+-- This routine is called when a reference is made to a member of some
+-- model object, but its n-tuple is out of the object domain. */
+void out_of_domain
+(     MPL *mpl,
+char *name,             /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     xassert(name != NULL);
+xassert(tuple != NULL);
+error(mpl, "%s%s out of domain", name, format_tuple(mpl, '[',
+tuple));
+/* no return */
+}
+/*----------------------------------------------------------------------
+-- get_domain_tuple - obtain current n-tuple from domain.
+--
+-- This routine constructs n-tuple, whose components are current values
+-- assigned to *free* dummy indices of specified domain.
+--
+-- For the sake of convenience it is allowed to specify domain as NULL,
+-- in which case this routine returns 0-tuple.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+TUPLE *get_domain_tuple
+(     MPL *mpl,
+DOMAIN *domain          /* not changed */
+)
+{     DOMAIN_BLOCK *block;
+DOMAIN_SLOT *slot;
+TUPLE *tuple;
+tuple = create_tuple(mpl);
+if (domain != NULL)
+{  for (block = domain->list; block != NULL; block = block->next)
+{  for (slot = block->list; slot != NULL; slot = slot->next)
+{  if (slot->code == NULL)
+{  xassert(slot->value != NULL);
+tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,
+slot->value));
+}
+}
+}
+}
+return tuple;
+}
+/*----------------------------------------------------------------------
+-- clean_domain - clean domain.
+--
+-- This routine cleans specified domain that assumes deleting all stuff
+-- dynamically allocated during the generation phase. */
+void clean_domain(MPL *mpl, DOMAIN *domain)
+{     DOMAIN_BLOCK *block;
+DOMAIN_SLOT *slot;
+/* if no domain is specified, do nothing */
+if (domain == NULL) goto done;
+/* clean all domain blocks */
+for (block = domain->list; block != NULL; block = block->next)
+{  /* clean all domain slots */
+for (slot = block->list; slot != NULL; slot = slot->next)
+{  /* clean pseudo-code for computing bound value */
+clean_code(mpl, slot->code);
+/* delete symbolic value assigned to dummy index */
+if (slot->value != NULL)
+delete_symbol(mpl, slot->value), slot->value = NULL;
+}
+/* clean pseudo-code for computing basic set */
+clean_code(mpl, block->code);
+}
+/* clean pseudo-code for computing domain predicate */
+clean_code(mpl, domain->code);
+done: return;
+}
+/**********************************************************************/
+/* * *                         MODEL SETS                         * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- check_elem_set - check elemental set assigned to set member.
+--
+-- This routine checks if given elemental set being assigned to member
+-- of specified model set satisfies to all restrictions.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+void check_elem_set
+(     MPL *mpl,
+SET *set,               /* not changed */
+TUPLE *tuple,           /* not changed */
+ELEMSET *refer          /* not changed */
+)
+{     WITHIN *within;
+MEMBER *memb;
+int eqno;
+/* elemental set must be within all specified supersets */
+for (within = set->within, eqno = 1; within != NULL; within =
+within->next, eqno++)
+{  xassert(within->code != NULL);
+for (memb = refer->head; memb != NULL; memb = memb->next)
+{  if (!is_member(mpl, within->code, memb->tuple))
+{  char buf[255+1];
+strcpy(buf, format_tuple(mpl, '(', memb->tuple));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s contains %s which not within specified "
+"set; see (%d)", set->name, format_tuple(mpl, '[',
+tuple), buf, eqno);
+}
+}
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- take_member_set - obtain elemental set assigned to set member.
+--
+-- This routine obtains a reference to elemental set assigned to given
+-- member of specified model set and returns it on exit.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+ELEMSET *take_member_set      /* returns reference, not value */
+(     MPL *mpl,
+SET *set,               /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     MEMBER *memb;
+ELEMSET *refer;
+/* find member in the set array */
+memb = find_member(mpl, set->array, tuple);
+if (memb != NULL)
+{  /* member exists, so just take the reference */
+refer = memb->value.set;
+}
+else if (set->assign != NULL)
+{  /* compute value using assignment expression */
+refer = eval_elemset(mpl, set->assign);
+add:     /* check that the elemental set satisfies to all restrictions,
+assign it to new member, and add the member to the array */
+check_elem_set(mpl, set, tuple, refer);
+memb = add_member(mpl, set->array, copy_tuple(mpl, tuple));
+memb->value.set = refer;
+}
+else if (set->option != NULL)
+{  /* compute default elemental set */
+refer = eval_elemset(mpl, set->option);
+goto add;
+}
+else
+{  /* no value (elemental set) is provided */
+error(mpl, "no value for %s%s", set->name, format_tuple(mpl,
+'[', tuple));
+}
+return refer;
+}
+/*----------------------------------------------------------------------
+-- eval_member_set - evaluate elemental set assigned to set member.
+--
+-- This routine evaluates a reference to elemental set assigned to given
+-- member of specified model set and returns it on exit. */
+struct eval_set_info
+{     /* working info used by the routine eval_member_set */
+SET *set;
+/* model set */
+TUPLE *tuple;
+/* n-tuple, which defines set member */
+MEMBER *memb;
+/* normally this pointer is NULL; the routine uses this pointer
+to check data provided in the data section, in which case it
+points to a member currently checked; this check is performed
+automatically only once when a reference to any member occurs
+for the first time */
+ELEMSET *refer;
+/* evaluated reference to elemental set */
+};
+static void eval_set_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine to work within domain scope */
+struct eval_set_info *info = _info;
+if (info->memb != NULL)
+{  /* checking call; check elemental set being assigned */
+check_elem_set(mpl, info->set, info->memb->tuple,
+info->memb->value.set);
+}
+else
+{  /* normal call; evaluate member, which has given n-tuple */
+info->refer = take_member_set(mpl, info->set, info->tuple);
+}
+return;
+}
+#if 1 /* 12/XII-2008 */
+static void saturate_set(MPL *mpl, SET *set)
+{     GADGET *gadget = set->gadget;
+ELEMSET *data;
+MEMBER *elem, *memb;
+TUPLE *tuple, *work[20];
+int i;
+xprintf("Generating %s...\n", set->name);
+eval_whole_set(mpl, gadget->set);
+/* gadget set must have exactly one member */
+xassert(gadget->set->array != NULL);
+xassert(gadget->set->array->head != NULL);
+xassert(gadget->set->array->head == gadget->set->array->tail);
+data = gadget->set->array->head->value.set;
+xassert(data->type == A_NONE);
+xassert(data->dim == gadget->set->dimen);
+/* walk thru all elements of the plain set */
+for (elem = data->head; elem != NULL; elem = elem->next)
+{  /* create a copy of n-tuple */
+tuple = copy_tuple(mpl, elem->tuple);
+/* rearrange component of the n-tuple */
+for (i = 0; i < gadget->set->dimen; i++)
+work[i] = NULL;
+for (i = 0; tuple != NULL; tuple = tuple->next)
+work[gadget->ind[i++]-1] = tuple;
+xassert(i == gadget->set->dimen);
+for (i = 0; i < gadget->set->dimen; i++)
+{  xassert(work[i] != NULL);
+work[i]->next = work[i+1];
+}
+/* construct subscript list from first set->dim components */
+if (set->dim == 0)
+tuple = NULL;
+else
+tuple = work[0], work[set->dim-1]->next = NULL;
+/* find corresponding member of the set to be initialized */
+memb = find_member(mpl, set->array, tuple);
+if (memb == NULL)
+{  /* not found; add new member to the set and assign it empty
+elemental set */
+memb = add_member(mpl, set->array, tuple);
+memb->value.set = create_elemset(mpl, set->dimen);
+}
+else
+{  /* found; free subscript list */
+delete_tuple(mpl, tuple);
+}
+/* construct new n-tuple from rest set->dimen components */
+tuple = work[set->dim];
+xassert(set->dim + set->dimen == gadget->set->dimen);
+work[gadget->set->dimen-1]->next = NULL;
+/* and add it to the elemental set assigned to the member
+(no check for duplicates is needed) */
+add_tuple(mpl, memb->value.set, tuple);
+}
+/* the set has been saturated with data */
+set->data = 1;
+return;
+}
+#endif
+ELEMSET *eval_member_set      /* returns reference, not value */
+(     MPL *mpl,
+SET *set,               /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     /* this routine evaluates set member */
+struct eval_set_info _info, *info = &_info;
+xassert(set->dim == tuple_dimen(mpl, tuple));
+info->set = set;
+info->tuple = tuple;
+#if 1 /* 12/XII-2008 */
+if (set->gadget != NULL && set->data == 0)
+{  /* initialize the set with data from a plain set */
+saturate_set(mpl, set);
+}
+#endif
+if (set->data == 1)
+{  /* check data, which are provided in the data section, but not
+checked yet */
+/* save pointer to the last array member; note that during the
+check new members may be added beyond the last member due to
+references to the same parameter from default expression as
+well as from expressions that define restricting supersets;
+however, values assigned to the new members will be checked
+by other routine, so we don't need to check them here */
+MEMBER *tail = set->array->tail;
+/* change the data status to prevent infinite recursive loop
+due to references to the same set during the check */
+set->data = 2;
+/* check elemental sets assigned to array members in the data
+section until the marked member has been reached */
+for (info->memb = set->array->head; info->memb != NULL;
+info->memb = info->memb->next)
+{  if (eval_within_domain(mpl, set->domain, info->memb->tuple,
+info, eval_set_func))
+out_of_domain(mpl, set->name, info->memb->tuple);
+if (info->memb == tail) break;
+}
+/* the check has been finished */
+}
+/* evaluate member, which has given n-tuple */
+info->memb = NULL;
+if (eval_within_domain(mpl, info->set->domain, info->tuple, info,
+eval_set_func))
+out_of_domain(mpl, set->name, info->tuple);
+/* bring evaluated reference to the calling program */
+return info->refer;
+}
+/*----------------------------------------------------------------------
+-- eval_whole_set - evaluate model set over entire domain.
+--
+-- This routine evaluates all members of specified model set over entire
+-- domain. */
+static int whole_set_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+SET *set = (SET *)info;
+TUPLE *tuple = get_domain_tuple(mpl, set->domain);
+eval_member_set(mpl, set, tuple);
+delete_tuple(mpl, tuple);
+return 0;
+}
+void eval_whole_set(MPL *mpl, SET *set)
+{     loop_within_domain(mpl, set->domain, set, whole_set_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean set - clean model set.
+--
+-- This routine cleans specified model set that assumes deleting all
+-- stuff dynamically allocated during the generation phase. */
+void clean_set(MPL *mpl, SET *set)
+{     WITHIN *within;
+MEMBER *memb;
+/* clean subscript domain */
+clean_domain(mpl, set->domain);
+/* clean pseudo-code for computing supersets */
+for (within = set->within; within != NULL; within = within->next)
+clean_code(mpl, within->code);
+/* clean pseudo-code for computing assigned value */
+clean_code(mpl, set->assign);
+/* clean pseudo-code for computing default value */
+clean_code(mpl, set->option);
+/* reset data status flag */
+set->data = 0;
+/* delete content array */
+for (memb = set->array->head; memb != NULL; memb = memb->next)
+delete_value(mpl, set->array->type, &memb->value);
+delete_array(mpl, set->array), set->array = NULL;
+return;
+}
+/**********************************************************************/
+/* * *                      MODEL PARAMETERS                      * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- check_value_num - check numeric value assigned to parameter member.
+--
+-- This routine checks if numeric value being assigned to some member
+-- of specified numeric model parameter satisfies to all restrictions.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+void check_value_num
+(     MPL *mpl,
+PARAMETER *par,         /* not changed */
+TUPLE *tuple,           /* not changed */
+double value
+)
+{     CONDITION *cond;
+WITHIN *in;
+int eqno;
+/* the value must satisfy to the parameter type */
+switch (par->type)
+{  case A_NUMERIC:
+break;
+case A_INTEGER:
+if (value != floor(value))
+error(mpl, "%s%s = %.*g not integer", par->name,
+format_tuple(mpl, '[', tuple), DBL_DIG, value);
+break;
+case A_BINARY:
+if (!(value == 0.0 || value == 1.0))
+error(mpl, "%s%s = %.*g not binary", par->name,
+format_tuple(mpl, '[', tuple), DBL_DIG, value);
+break;
+default:
+xassert(par != par);
+}
+/* the value must satisfy to all specified conditions */
+for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,
+eqno++)
+{  double bound;
+char *rho;
+xassert(cond->code != NULL);
+bound = eval_numeric(mpl, cond->code);
+switch (cond->rho)
+{  case O_LT:
+if (!(value < bound))
+{  rho = "<";
+err:              error(mpl, "%s%s = %.*g not %s %.*g; see (%d)",
+par->name, format_tuple(mpl, '[', tuple), DBL_DIG,
+value, rho, DBL_DIG, bound, eqno);
+}
+break;
+case O_LE:
+if (!(value <= bound)) { rho = "<="; goto err; }
+break;
+case O_EQ:
+if (!(value == bound)) { rho = "="; goto err; }
+break;
+case O_GE:
+if (!(value >= bound)) { rho = ">="; goto err; }
+break;
+case O_GT:
+if (!(value > bound)) { rho = ">"; goto err; }
+break;
+case O_NE:
+if (!(value != bound)) { rho = "<>"; goto err; }
+break;
+default:
+xassert(cond != cond);
+}
+}
+/* the value must be in all specified supersets */
+for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)
+{  TUPLE *dummy;
+xassert(in->code != NULL);
+xassert(in->code->dim == 1);
+dummy = expand_tuple(mpl, create_tuple(mpl),
+create_symbol_num(mpl, value));
+if (!is_member(mpl, in->code, dummy))
+error(mpl, "%s%s = %.*g not in specified set; see (%d)",
+par->name, format_tuple(mpl, '[', tuple), DBL_DIG,
+value, eqno);
+delete_tuple(mpl, dummy);
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- take_member_num - obtain num. value assigned to parameter member.
+--
+-- This routine obtains a numeric value assigned to member of specified
+-- numeric model parameter and returns it on exit.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+double take_member_num
+(     MPL *mpl,
+PARAMETER *par,         /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     MEMBER *memb;
+double value;
+/* find member in the parameter array */
+memb = find_member(mpl, par->array, tuple);
+if (memb != NULL)
+{  /* member exists, so just take its value */
+value = memb->value.num;
+}
+else if (par->assign != NULL)
+{  /* compute value using assignment expression */
+value = eval_numeric(mpl, par->assign);
+add:     /* check that the value satisfies to all restrictions, assign
+it to new member, and add the member to the array */
+check_value_num(mpl, par, tuple, value);
+memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));
+memb->value.num = value;
+}
+else if (par->option != NULL)
+{  /* compute default value */
+value = eval_numeric(mpl, par->option);
+goto add;
+}
+else if (par->defval != NULL)
+{  /* take default value provided in the data section */
+if (par->defval->str != NULL)
+error(mpl, "cannot convert %s to floating-point number",
+format_symbol(mpl, par->defval));
+value = par->defval->num;
+goto add;
+}
+else
+{  /* no value is provided */
+error(mpl, "no value for %s%s", par->name, format_tuple(mpl,
+'[', tuple));
+}
+return value;
+}
+/*----------------------------------------------------------------------
+-- eval_member_num - evaluate num. value assigned to parameter member.
+--
+-- This routine evaluates a numeric value assigned to given member of
+-- specified numeric model parameter and returns it on exit. */
+struct eval_num_info
+{     /* working info used by the routine eval_member_num */
+PARAMETER *par;
+/* model parameter  */
+TUPLE *tuple;
+/* n-tuple, which defines parameter member */
+MEMBER *memb;
+/* normally this pointer is NULL; the routine uses this pointer
+to check data provided in the data section, in which case it
+points to a member currently checked; this check is performed
+automatically only once when a reference to any member occurs
+for the first time */
+double value;
+/* evaluated numeric value */
+};
+static void eval_num_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine to work within domain scope */
+struct eval_num_info *info = _info;
+if (info->memb != NULL)
+{  /* checking call; check numeric value being assigned */
+check_value_num(mpl, info->par, info->memb->tuple,
+info->memb->value.num);
+}
+else
+{  /* normal call; evaluate member, which has given n-tuple */
+info->value = take_member_num(mpl, info->par, info->tuple);
+}
+return;
+}
+double eval_member_num
+(     MPL *mpl,
+PARAMETER *par,         /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     /* this routine evaluates numeric parameter member */
+struct eval_num_info _info, *info = &_info;
+xassert(par->type == A_NUMERIC || par->type == A_INTEGER ||
+par->type == A_BINARY);
+xassert(par->dim == tuple_dimen(mpl, tuple));
+info->par = par;
+info->tuple = tuple;
+if (par->data == 1)
+{  /* check data, which are provided in the data section, but not
+checked yet */
+/* save pointer to the last array member; note that during the
+check new members may be added beyond the last member due to
+references to the same parameter from default expression as
+well as from expressions that define restricting conditions;
+however, values assigned to the new members will be checked
+by other routine, so we don't need to check them here */
+MEMBER *tail = par->array->tail;
+/* change the data status to prevent infinite recursive loop
+due to references to the same parameter during the check */
+par->data = 2;
+/* check values assigned to array members in the data section
+until the marked member has been reached */
+for (info->memb = par->array->head; info->memb != NULL;
+info->memb = info->memb->next)
+{  if (eval_within_domain(mpl, par->domain, info->memb->tuple,
+info, eval_num_func))
+out_of_domain(mpl, par->name, info->memb->tuple);
+if (info->memb == tail) break;
+}
+/* the check has been finished */
+}
+/* evaluate member, which has given n-tuple */
+info->memb = NULL;
+if (eval_within_domain(mpl, info->par->domain, info->tuple, info,
+eval_num_func))
+out_of_domain(mpl, par->name, info->tuple);
+/* bring evaluated value to the calling program */
+return info->value;
+}
+/*----------------------------------------------------------------------
+-- check_value_sym - check symbolic value assigned to parameter member.
+--
+-- This routine checks if symbolic value being assigned to some member
+-- of specified symbolic model parameter satisfies to all restrictions.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+void check_value_sym
+(     MPL *mpl,
+PARAMETER *par,         /* not changed */
+TUPLE *tuple,           /* not changed */
+SYMBOL *value           /* not changed */
+)
+{     CONDITION *cond;
+WITHIN *in;
+int eqno;
+/* the value must satisfy to all specified conditions */
+for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,
+eqno++)
+{  SYMBOL *bound;
+char buf[255+1];
+xassert(cond->code != NULL);
+bound = eval_symbolic(mpl, cond->code);
+switch (cond->rho)
+{
+#if 1 /* 13/VIII-2008 */
+case O_LT:
+if (!(compare_symbols(mpl, value, bound) < 0))
+{  strcpy(buf, format_symbol(mpl, bound));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s = %s not < %s",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), buf, eqno);
+}
+break;
+case O_LE:
+if (!(compare_symbols(mpl, value, bound) <= 0))
+{  strcpy(buf, format_symbol(mpl, bound));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s = %s not <= %s",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), buf, eqno);
+}
+break;
+#endif
+case O_EQ:
+if (!(compare_symbols(mpl, value, bound) == 0))
+{  strcpy(buf, format_symbol(mpl, bound));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s = %s not = %s",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), buf, eqno);
+}
+break;
+#if 1 /* 13/VIII-2008 */
+case O_GE:
+if (!(compare_symbols(mpl, value, bound) >= 0))
+{  strcpy(buf, format_symbol(mpl, bound));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s = %s not >= %s",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), buf, eqno);
+}
+break;
+case O_GT:
+if (!(compare_symbols(mpl, value, bound) > 0))
+{  strcpy(buf, format_symbol(mpl, bound));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s = %s not > %s",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), buf, eqno);
+}
+break;
+#endif
+case O_NE:
+if (!(compare_symbols(mpl, value, bound) != 0))
+{  strcpy(buf, format_symbol(mpl, bound));
+xassert(strlen(buf) < sizeof(buf));
+error(mpl, "%s%s = %s not <> %s",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), buf, eqno);
+}
+break;
+default:
+xassert(cond != cond);
+}
+delete_symbol(mpl, bound);
+}
+/* the value must be in all specified supersets */
+for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)
+{  TUPLE *dummy;
+xassert(in->code != NULL);
+xassert(in->code->dim == 1);
+dummy = expand_tuple(mpl, create_tuple(mpl), copy_symbol(mpl,
+value));
+if (!is_member(mpl, in->code, dummy))
+error(mpl, "%s%s = %s not in specified set; see (%d)",
+par->name, format_tuple(mpl, '[', tuple),
+format_symbol(mpl, value), eqno);
+delete_tuple(mpl, dummy);
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- take_member_sym - obtain symb. value assigned to parameter member.
+--
+-- This routine obtains a symbolic value assigned to member of specified
+-- symbolic model parameter and returns it on exit.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+SYMBOL *take_member_sym       /* returns value, not reference */
+(     MPL *mpl,
+PARAMETER *par,         /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     MEMBER *memb;
+SYMBOL *value;
+/* find member in the parameter array */
+memb = find_member(mpl, par->array, tuple);
+if (memb != NULL)
+{  /* member exists, so just take its value */
+value = copy_symbol(mpl, memb->value.sym);
+}
+else if (par->assign != NULL)
+{  /* compute value using assignment expression */
+value = eval_symbolic(mpl, par->assign);
+add:     /* check that the value satisfies to all restrictions, assign
+it to new member, and add the member to the array */
+check_value_sym(mpl, par, tuple, value);
+memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));
+memb->value.sym = copy_symbol(mpl, value);
+}
+else if (par->option != NULL)
+{  /* compute default value */
+value = eval_symbolic(mpl, par->option);
+goto add;
+}
+else if (par->defval != NULL)
+{  /* take default value provided in the data section */
+value = copy_symbol(mpl, par->defval);
+goto add;
+}
+else
+{  /* no value is provided */
+error(mpl, "no value for %s%s", par->name, format_tuple(mpl,
+'[', tuple));
+}
+return value;
+}
+/*----------------------------------------------------------------------
+-- eval_member_sym - evaluate symb. value assigned to parameter member.
+--
+-- This routine evaluates a symbolic value assigned to given member of
+-- specified symbolic model parameter and returns it on exit. */
+struct eval_sym_info
+{     /* working info used by the routine eval_member_sym */
+PARAMETER *par;
+/* model parameter */
+TUPLE *tuple;
+/* n-tuple, which defines parameter member */
+MEMBER *memb;
+/* normally this pointer is NULL; the routine uses this pointer
+to check data provided in the data section, in which case it
+points to a member currently checked; this check is performed
+automatically only once when a reference to any member occurs
+for the first time */
+SYMBOL *value;
+/* evaluated symbolic value */
+};
+static void eval_sym_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine to work within domain scope */
+struct eval_sym_info *info = _info;
+if (info->memb != NULL)
+{  /* checking call; check symbolic value being assigned */
+check_value_sym(mpl, info->par, info->memb->tuple,
+info->memb->value.sym);
+}
+else
+{  /* normal call; evaluate member, which has given n-tuple */
+info->value = take_member_sym(mpl, info->par, info->tuple);
+}
+return;
+}
+SYMBOL *eval_member_sym       /* returns value, not reference */
+(     MPL *mpl,
+PARAMETER *par,         /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     /* this routine evaluates symbolic parameter member */
+struct eval_sym_info _info, *info = &_info;
+xassert(par->type == A_SYMBOLIC);
+xassert(par->dim == tuple_dimen(mpl, tuple));
+info->par = par;
+info->tuple = tuple;
+if (par->data == 1)
+{  /* check data, which are provided in the data section, but not
+checked yet */
+/* save pointer to the last array member; note that during the
+check new members may be added beyond the last member due to
+references to the same parameter from default expression as
+well as from expressions that define restricting conditions;
+however, values assigned to the new members will be checked
+by other routine, so we don't need to check them here */
+MEMBER *tail = par->array->tail;
+/* change the data status to prevent infinite recursive loop
+due to references to the same parameter during the check */
+par->data = 2;
+/* check values assigned to array members in the data section
+until the marked member has been reached */
+for (info->memb = par->array->head; info->memb != NULL;
+info->memb = info->memb->next)
+{  if (eval_within_domain(mpl, par->domain, info->memb->tuple,
+info, eval_sym_func))
+out_of_domain(mpl, par->name, info->memb->tuple);
+if (info->memb == tail) break;
+}
+/* the check has been finished */
+}
+/* evaluate member, which has given n-tuple */
+info->memb = NULL;
+if (eval_within_domain(mpl, info->par->domain, info->tuple, info,
+eval_sym_func))
+out_of_domain(mpl, par->name, info->tuple);
+/* bring evaluated value to the calling program */
+return info->value;
+}
+/*----------------------------------------------------------------------
+-- eval_whole_par - evaluate model parameter over entire domain.
+--
+-- This routine evaluates all members of specified model parameter over
+-- entire domain. */
+static int whole_par_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+PARAMETER *par = (PARAMETER *)info;
+TUPLE *tuple = get_domain_tuple(mpl, par->domain);
+switch (par->type)
+{  case A_NUMERIC:
+case A_INTEGER:
+case A_BINARY:
+eval_member_num(mpl, par, tuple);
+break;
+case A_SYMBOLIC:
+delete_symbol(mpl, eval_member_sym(mpl, par, tuple));
+break;
+default:
+xassert(par != par);
+}
+delete_tuple(mpl, tuple);
+return 0;
+}
+void eval_whole_par(MPL *mpl, PARAMETER *par)
+{     loop_within_domain(mpl, par->domain, par, whole_par_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_parameter - clean model parameter.
+--
+-- This routine cleans specified model parameter that assumes deleting
+-- all stuff dynamically allocated during the generation phase. */
+void clean_parameter(MPL *mpl, PARAMETER *par)
+{     CONDITION *cond;
+WITHIN *in;
+MEMBER *memb;
+/* clean subscript domain */
+clean_domain(mpl, par->domain);
+/* clean pseudo-code for computing restricting conditions */
+for (cond = par->cond; cond != NULL; cond = cond->next)
+clean_code(mpl, cond->code);
+/* clean pseudo-code for computing restricting supersets */
+for (in = par->in; in != NULL; in = in->next)
+clean_code(mpl, in->code);
+/* clean pseudo-code for computing assigned value */
+clean_code(mpl, par->assign);
+/* clean pseudo-code for computing default value */
+clean_code(mpl, par->option);
+/* reset data status flag */
+par->data = 0;
+/* delete default symbolic value */
+if (par->defval != NULL)
+delete_symbol(mpl, par->defval), par->defval = NULL;
+/* delete content array */
+for (memb = par->array->head; memb != NULL; memb = memb->next)
+delete_value(mpl, par->array->type, &memb->value);
+delete_array(mpl, par->array), par->array = NULL;
+return;
+}
+/**********************************************************************/
+/* * *                      MODEL VARIABLES                       * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- take_member_var - obtain reference to elemental variable.
+--
+-- This routine obtains a reference to elemental variable assigned to
+-- given member of specified model variable and returns it on exit. If
+-- necessary, new elemental variable is created.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+ELEMVAR *take_member_var      /* returns reference */
+(     MPL *mpl,
+VARIABLE *var,          /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     MEMBER *memb;
+ELEMVAR *refer;
+/* find member in the variable array */
+memb = find_member(mpl, var->array, tuple);
+if (memb != NULL)
+{  /* member exists, so just take the reference */
+refer = memb->value.var;
+}
+else
+{  /* member is referenced for the first time and therefore does
+not exist; create new elemental variable, assign it to new
+member, and add the member to the variable array */
+memb = add_member(mpl, var->array, copy_tuple(mpl, tuple));
+refer = (memb->value.var =
+dmp_get_atom(mpl->elemvars, sizeof(ELEMVAR)));
+refer->j = 0;
+refer->var = var;
+refer->memb = memb;
+/* compute lower bound */
+if (var->lbnd == NULL)
+refer->lbnd = 0.0;
+else
+refer->lbnd = eval_numeric(mpl, var->lbnd);
+/* compute upper bound */
+if (var->ubnd == NULL)
+refer->ubnd = 0.0;
+else if (var->ubnd == var->lbnd)
+refer->ubnd = refer->lbnd;
+else
+refer->ubnd = eval_numeric(mpl, var->ubnd);
+/* nullify working quantity */
+refer->temp = 0.0;
+#if 1 /* 15/V-2010 */
+/* solution has not been obtained by the solver yet */
+refer->stat = 0;
+refer->prim = refer->dual = 0.0;
+#endif
+}
+return refer;
+}
+/*----------------------------------------------------------------------
+-- eval_member_var - evaluate reference to elemental variable.
+--
+-- This routine evaluates a reference to elemental variable assigned to
+-- member of specified model variable and returns it on exit. */
+struct eval_var_info
+{     /* working info used by the routine eval_member_var */
+VARIABLE *var;
+/* model variable */
+TUPLE *tuple;
+/* n-tuple, which defines variable member */
+ELEMVAR *refer;
+/* evaluated reference to elemental variable */
+};
+static void eval_var_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine to work within domain scope */
+struct eval_var_info *info = _info;
+info->refer = take_member_var(mpl, info->var, info->tuple);
+return;
+}
+ELEMVAR *eval_member_var      /* returns reference */
+(     MPL *mpl,
+VARIABLE *var,          /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     /* this routine evaluates variable member */
+struct eval_var_info _info, *info = &_info;
+xassert(var->dim == tuple_dimen(mpl, tuple));
+info->var = var;
+info->tuple = tuple;
+/* evaluate member, which has given n-tuple */
+if (eval_within_domain(mpl, info->var->domain, info->tuple, info,
+eval_var_func))
+out_of_domain(mpl, var->name, info->tuple);
+/* bring evaluated reference to the calling program */
+return info->refer;
+}
+/*----------------------------------------------------------------------
+-- eval_whole_var - evaluate model variable over entire domain.
+--
+-- This routine evaluates all members of specified model variable over
+-- entire domain. */
+static int whole_var_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+VARIABLE *var = (VARIABLE *)info;
+TUPLE *tuple = get_domain_tuple(mpl, var->domain);
+eval_member_var(mpl, var, tuple);
+delete_tuple(mpl, tuple);
+return 0;
+}
+void eval_whole_var(MPL *mpl, VARIABLE *var)
+{     loop_within_domain(mpl, var->domain, var, whole_var_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_variable - clean model variable.
+--
+-- This routine cleans specified model variable that assumes deleting
+-- all stuff dynamically allocated during the generation phase. */
+void clean_variable(MPL *mpl, VARIABLE *var)
+{     MEMBER *memb;
+/* clean subscript domain */
+clean_domain(mpl, var->domain);
+/* clean code for computing lower bound */
+clean_code(mpl, var->lbnd);
+/* clean code for computing upper bound */
+if (var->ubnd != var->lbnd) clean_code(mpl, var->ubnd);
+/* delete content array */
+for (memb = var->array->head; memb != NULL; memb = memb->next)
+dmp_free_atom(mpl->elemvars, memb->value.var, sizeof(ELEMVAR));
+delete_array(mpl, var->array), var->array = NULL;
+return;
+}
+/**********************************************************************/
+/* * *              MODEL CONSTRAINTS AND OBJECTIVES              * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- take_member_con - obtain reference to elemental constraint.
+--
+-- This routine obtains a reference to elemental constraint assigned
+-- to given member of specified model constraint and returns it on exit.
+-- If necessary, new elemental constraint is created.
+--
+-- NOTE: This routine must not be called out of domain scope. */
+ELEMCON *take_member_con      /* returns reference */
+(     MPL *mpl,
+CONSTRAINT *con,        /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     MEMBER *memb;
+ELEMCON *refer;
+/* find member in the constraint array */
+memb = find_member(mpl, con->array, tuple);
+if (memb != NULL)
+{  /* member exists, so just take the reference */
+refer = memb->value.con;
+}
+else
+{  /* member is referenced for the first time and therefore does
+not exist; create new elemental constraint, assign it to new
+member, and add the member to the constraint array */
+memb = add_member(mpl, con->array, copy_tuple(mpl, tuple));
+refer = (memb->value.con =
+dmp_get_atom(mpl->elemcons, sizeof(ELEMCON)));
+refer->i = 0;
+refer->con = con;
+refer->memb = memb;
+/* compute linear form */
+xassert(con->code != NULL);
+refer->form = eval_formula(mpl, con->code);
+/* compute lower and upper bounds */
+if (con->lbnd == NULL && con->ubnd == NULL)
+{  /* objective has no bounds */
+double temp;
+xassert(con->type == A_MINIMIZE || con->type == A_MAXIMIZE);
+/* carry the constant term to the right-hand side */
+refer->form = remove_constant(mpl, refer->form, &temp);
+refer->lbnd = refer->ubnd = - temp;
+}
+else if (con->lbnd != NULL && con->ubnd == NULL)
+{  /* constraint a * x + b >= c * y + d is transformed to the
+standard form a * x - c * y >= d - b */
+double temp;
+xassert(con->type == A_CONSTRAINT);
+refer->form = linear_comb(mpl,
++1.0, refer->form,
+-1.0, eval_formula(mpl, con->lbnd));
+refer->form = remove_constant(mpl, refer->form, &temp);
+refer->lbnd = - temp;
+refer->ubnd = 0.0;
+}
+else if (con->lbnd == NULL && con->ubnd != NULL)
+{  /* constraint a * x + b <= c * y + d is transformed to the
+standard form a * x - c * y <= d - b */
+double temp;
+xassert(con->type == A_CONSTRAINT);
+refer->form = linear_comb(mpl,
++1.0, refer->form,
+-1.0, eval_formula(mpl, con->ubnd));
+refer->form = remove_constant(mpl, refer->form, &temp);
+refer->lbnd = 0.0;
+refer->ubnd = - temp;
+}
+else if (con->lbnd == con->ubnd)
+{  /* constraint a * x + b = c * y + d is transformed to the
+standard form a * x - c * y = d - b */
+double temp;
+xassert(con->type == A_CONSTRAINT);
+refer->form = linear_comb(mpl,
++1.0, refer->form,
+-1.0, eval_formula(mpl, con->lbnd));
+refer->form = remove_constant(mpl, refer->form, &temp);
+refer->lbnd = refer->ubnd = - temp;
+}
+else
+{  /* ranged constraint c <= a * x + b <= d is transformed to
+the standard form c - b <= a * x <= d - b */
+double temp, temp1, temp2;
+xassert(con->type == A_CONSTRAINT);
+refer->form = remove_constant(mpl, refer->form, &temp);
+xassert(remove_constant(mpl, eval_formula(mpl, con->lbnd),
+&temp1) == NULL);
+xassert(remove_constant(mpl, eval_formula(mpl, con->ubnd),
+&temp2) == NULL);
+refer->lbnd = fp_sub(mpl, temp1, temp);
+refer->ubnd = fp_sub(mpl, temp2, temp);
+}
+#if 1 /* 15/V-2010 */
+/* solution has not been obtained by the solver yet */
+refer->stat = 0;
+refer->prim = refer->dual = 0.0;
+#endif
+}
+return refer;
+}
+/*----------------------------------------------------------------------
+-- eval_member_con - evaluate reference to elemental constraint.
+--
+-- This routine evaluates a reference to elemental constraint assigned
+-- to member of specified model constraint and returns it on exit. */
+struct eval_con_info
+{     /* working info used by the routine eval_member_con */
+CONSTRAINT *con;
+/* model constraint */
+TUPLE *tuple;
+/* n-tuple, which defines constraint member */
+ELEMCON *refer;
+/* evaluated reference to elemental constraint */
+};
+static void eval_con_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine to work within domain scope */
+struct eval_con_info *info = _info;
+info->refer = take_member_con(mpl, info->con, info->tuple);
+return;
+}
+ELEMCON *eval_member_con      /* returns reference */
+(     MPL *mpl,
+CONSTRAINT *con,        /* not changed */
+TUPLE *tuple            /* not changed */
+)
+{     /* this routine evaluates constraint member */
+struct eval_con_info _info, *info = &_info;
+xassert(con->dim == tuple_dimen(mpl, tuple));
+info->con = con;
+info->tuple = tuple;
+/* evaluate member, which has given n-tuple */
+if (eval_within_domain(mpl, info->con->domain, info->tuple, info,
+eval_con_func))
+out_of_domain(mpl, con->name, info->tuple);
+/* bring evaluated reference to the calling program */
+return info->refer;
+}
+/*----------------------------------------------------------------------
+-- eval_whole_con - evaluate model constraint over entire domain.
+--
+-- This routine evaluates all members of specified model constraint over
+-- entire domain. */
+static int whole_con_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+CONSTRAINT *con = (CONSTRAINT *)info;
+TUPLE *tuple = get_domain_tuple(mpl, con->domain);
+eval_member_con(mpl, con, tuple);
+delete_tuple(mpl, tuple);
+return 0;
+}
+void eval_whole_con(MPL *mpl, CONSTRAINT *con)
+{     loop_within_domain(mpl, con->domain, con, whole_con_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_constraint - clean model constraint.
+--
+-- This routine cleans specified model constraint that assumes deleting
+-- all stuff dynamically allocated during the generation phase. */
+void clean_constraint(MPL *mpl, CONSTRAINT *con)
+{     MEMBER *memb;
+/* clean subscript domain */
+clean_domain(mpl, con->domain);
+/* clean code for computing main linear form */
+clean_code(mpl, con->code);
+/* clean code for computing lower bound */
+clean_code(mpl, con->lbnd);
+/* clean code for computing upper bound */
+if (con->ubnd != con->lbnd) clean_code(mpl, con->ubnd);
+/* delete content array */
+for (memb = con->array->head; memb != NULL; memb = memb->next)
+{  delete_formula(mpl, memb->value.con->form);
+dmp_free_atom(mpl->elemcons, memb->value.con, sizeof(ELEMCON));
+}
+delete_array(mpl, con->array), con->array = NULL;
+return;
+}
+/**********************************************************************/
+/* * *                        PSEUDO-CODE                         * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- eval_numeric - evaluate pseudo-code to determine numeric value.
+--
+-- This routine evaluates specified pseudo-code to determine resultant
+-- numeric value, which is returned on exit. */
+struct iter_num_info
+{     /* working info used by the routine iter_num_func */
+CODE *code;
+/* pseudo-code for iterated operation to be performed */
+double value;
+/* resultant value */
+};
+static int iter_num_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine used to perform iterated operation
+on numeric "integrand" within domain scope */
+struct iter_num_info *info = _info;
+double temp;
+temp = eval_numeric(mpl, info->code->arg.loop.x);
+switch (info->code->op)
+{  case O_SUM:
+/* summation over domain */
+info->value = fp_add(mpl, info->value, temp);
+break;
+case O_PROD:
+/* multiplication over domain */
+info->value = fp_mul(mpl, info->value, temp);
+break;
+case O_MINIMUM:
+/* minimum over domain */
+if (info->value > temp) info->value = temp;
+break;
+case O_MAXIMUM:
+/* maximum over domain */
+if (info->value < temp) info->value = temp;
+break;
+default:
+xassert(info != info);
+}
+return 0;
+}
+double eval_numeric(MPL *mpl, CODE *code)
+{     double value;
+xassert(code != NULL);
+xassert(code->type == A_NUMERIC);
+xassert(code->dim == 0);
+/* if the operation has a side effect, invalidate and delete the
+resultant value */
+if (code->vflag && code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* if resultant value is valid, no evaluation is needed */
+if (code->valid)
+{  value = code->value.num;
+goto done;
+}
+/* evaluate pseudo-code recursively */
+switch (code->op)
+{  case O_NUMBER:
+/* take floating-point number */
+value = code->arg.num;
+break;
+case O_MEMNUM:
+/* take member of numeric parameter */
+{  TUPLE *tuple;
+ARG_LIST *e;
+tuple = create_tuple(mpl);
+for (e = code->arg.par.list; e != NULL; e = e->next)
+tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
+e->x));
+value = eval_member_num(mpl, code->arg.par.par, tuple);
+delete_tuple(mpl, tuple);
+}
+break;
+case O_MEMVAR:
+/* take computed value of elemental variable */
+{  TUPLE *tuple;
+ARG_LIST *e;
+#if 1 /* 15/V-2010 */
+ELEMVAR *var;
+#endif
+tuple = create_tuple(mpl);
+for (e = code->arg.var.list; e != NULL; e = e->next)
+tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
+e->x));
+#if 0 /* 15/V-2010 */
+value = eval_member_var(mpl, code->arg.var.var, tuple)
+->value;
+#else
+var = eval_member_var(mpl, code->arg.var.var, tuple);
+switch (code->arg.var.suff)
+{  case DOT_LB:
+if (var->var->lbnd == NULL)
+value = -DBL_MAX;
+else
+value = var->lbnd;
+break;
+case DOT_UB:
+if (var->var->ubnd == NULL)
+value = +DBL_MAX;
+else
+value = var->ubnd;
+break;
+case DOT_STATUS:
+value = var->stat;
+break;
+case DOT_VAL:
+value = var->prim;
+break;
+case DOT_DUAL:
+value = var->dual;
+break;
+default:
+xassert(code != code);
+}
+#endif
+delete_tuple(mpl, tuple);
+}
+break;
+#if 1 /* 15/V-2010 */
+case O_MEMCON:
+/* take computed value of elemental constraint */
+{  TUPLE *tuple;
+ARG_LIST *e;
+ELEMCON *con;
+tuple = create_tuple(mpl);
+for (e = code->arg.con.list; e != NULL; e = e->next)
+tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
+e->x));
+con = eval_member_con(mpl, code->arg.con.con, tuple);
+switch (code->arg.con.suff)
+{  case DOT_LB:
+if (con->con->lbnd == NULL)
+value = -DBL_MAX;
+else
+value = con->lbnd;
+break;
+case DOT_UB:
+if (con->con->ubnd == NULL)
+value = +DBL_MAX;
+else
+value = con->ubnd;
+break;
+case DOT_STATUS:
+value = con->stat;
+break;
+case DOT_VAL:
+value = con->prim;
+break;
+case DOT_DUAL:
+value = con->dual;
+break;
+default:
+xassert(code != code);
+}
+delete_tuple(mpl, tuple);
+}
+break;
+#endif
+case O_IRAND224:
+/* pseudo-random in [0, 2^24-1] */
+value = fp_irand224(mpl);
+break;
+case O_UNIFORM01:
+/* pseudo-random in [0, 1) */
+value = fp_uniform01(mpl);
+break;
+case O_NORMAL01:
+/* gaussian random, mu = 0, sigma = 1 */
+value = fp_normal01(mpl);
+break;
+case O_GMTIME:
+/* current calendar time */
+value = fn_gmtime(mpl);
+break;
+case O_CVTNUM:
+/* conversion to numeric */
+{  SYMBOL *sym;
+sym = eval_symbolic(mpl, code->arg.arg.x);
+#if 0 /* 23/XI-2008 */
+if (sym->str != NULL)
+error(mpl, "cannot convert %s to floating-point numbe"
+"r", format_symbol(mpl, sym));
+value = sym->num;
+#else
+if (sym->str == NULL)
+value = sym->num;
+else
+{  if (str2num(sym->str, &value))
+error(mpl, "cannot convert %s to floating-point nu"
+"mber", format_symbol(mpl, sym));
+}
+#endif
+delete_symbol(mpl, sym);
+}
+break;
+case O_PLUS:
+/* unary plus */
+value = + eval_numeric(mpl, code->arg.arg.x);
+break;
+case O_MINUS:
+/* unary minus */
+value = - eval_numeric(mpl, code->arg.arg.x);
+break;
+case O_ABS:
+/* absolute value */
+value = fabs(eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_CEIL:
+/* round upward ("ceiling of x") */
+value = ceil(eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_FLOOR:
+/* round downward ("floor of x") */
+value = floor(eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_EXP:
+/* base-e exponential */
+value = fp_exp(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_LOG:
+/* natural logarithm */
+value = fp_log(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_LOG10:
+/* common (decimal) logarithm */
+value = fp_log10(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_SQRT:
+/* square root */
+value = fp_sqrt(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_SIN:
+/* trigonometric sine */
+value = fp_sin(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_COS:
+/* trigonometric cosine */
+value = fp_cos(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_ATAN:
+/* trigonometric arctangent (one argument) */
+value = fp_atan(mpl, eval_numeric(mpl, code->arg.arg.x));
+break;
+case O_ATAN2:
+/* trigonometric arctangent (two arguments) */
+value = fp_atan2(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_ROUND:
+/* round to nearest integer */
+value = fp_round(mpl,
+eval_numeric(mpl, code->arg.arg.x), 0.0);
+break;
+case O_ROUND2:
+/* round to n fractional digits */
+value = fp_round(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_TRUNC:
+/* truncate to nearest integer */
+value = fp_trunc(mpl,
+eval_numeric(mpl, code->arg.arg.x), 0.0);
+break;
+case O_TRUNC2:
+/* truncate to n fractional digits */
+value = fp_trunc(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_ADD:
+/* addition */
+value = fp_add(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_SUB:
+/* subtraction */
+value = fp_sub(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_LESS:
+/* non-negative subtraction */
+value = fp_less(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_MUL:
+/* multiplication */
+value = fp_mul(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_DIV:
+/* division */
+value = fp_div(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_IDIV:
+/* quotient of exact division */
+value = fp_idiv(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_MOD:
+/* remainder of exact division */
+value = fp_mod(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_POWER:
+/* exponentiation (raise to power) */
+value = fp_power(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_UNIFORM:
+/* pseudo-random in [a, b) */
+value = fp_uniform(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_NORMAL:
+/* gaussian random, given mu and sigma */
+value = fp_normal(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y));
+break;
+case O_CARD:
+{  ELEMSET *set;
+set = eval_elemset(mpl, code->arg.arg.x);
+value = set->size;
+delete_array(mpl, set);
+}
+break;
+case O_LENGTH:
+{  SYMBOL *sym;
+char str[MAX_LENGTH+1];
+sym = eval_symbolic(mpl, code->arg.arg.x);
+if (sym->str == NULL)
+sprintf(str, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, str);
+delete_symbol(mpl, sym);
+value = strlen(str);
+}
+break;
+case O_STR2TIME:
+{  SYMBOL *sym;
+char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];
+sym = eval_symbolic(mpl, code->arg.arg.x);
+if (sym->str == NULL)
+sprintf(str, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, str);
+delete_symbol(mpl, sym);
+sym = eval_symbolic(mpl, code->arg.arg.y);
+if (sym->str == NULL)
+sprintf(fmt, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, fmt);
+delete_symbol(mpl, sym);
+value = fn_str2time(mpl, str, fmt);
+}
+break;
+case O_FORK:
+/* if-then-else */
+if (eval_logical(mpl, code->arg.arg.x))
+value = eval_numeric(mpl, code->arg.arg.y);
+else if (code->arg.arg.z == NULL)
+value = 0.0;
+else
+value = eval_numeric(mpl, code->arg.arg.z);
+break;
+case O_MIN:
+/* minimal value (n-ary) */
+{  ARG_LIST *e;
+double temp;
+value = +DBL_MAX;
+for (e = code->arg.list; e != NULL; e = e->next)
+{  temp = eval_numeric(mpl, e->x);
+if (value > temp) value = temp;
+}
+}
+break;
+case O_MAX:
+/* maximal value (n-ary) */
+{  ARG_LIST *e;
+double temp;
+value = -DBL_MAX;
+for (e = code->arg.list; e != NULL; e = e->next)
+{  temp = eval_numeric(mpl, e->x);
+if (value < temp) value = temp;
+}
+}
+break;
+case O_SUM:
+/* summation over domain */
+{  struct iter_num_info _info, *info = &_info;
+info->code = code;
+info->value = 0.0;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_num_func);
+value = info->value;
+}
+break;
+case O_PROD:
+/* multiplication over domain */
+{  struct iter_num_info _info, *info = &_info;
+info->code = code;
+info->value = 1.0;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_num_func);
+value = info->value;
+}
+break;
+case O_MINIMUM:
+/* minimum over domain */
+{  struct iter_num_info _info, *info = &_info;
+info->code = code;
+info->value = +DBL_MAX;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_num_func);
+if (info->value == +DBL_MAX)
+error(mpl, "min{} over empty set; result undefined");
+value = info->value;
+}
+break;
+case O_MAXIMUM:
+/* maximum over domain */
+{  struct iter_num_info _info, *info = &_info;
+info->code = code;
+info->value = -DBL_MAX;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_num_func);
+if (info->value == -DBL_MAX)
+error(mpl, "max{} over empty set; result undefined");
+value = info->value;
+}
+break;
+default:
+xassert(code != code);
+}
+/* save resultant value */
+xassert(!code->valid);
+code->valid = 1;
+code->value.num = value;
+done: return value;
+}
+/*----------------------------------------------------------------------
+-- eval_symbolic - evaluate pseudo-code to determine symbolic value.
+--
+-- This routine evaluates specified pseudo-code to determine resultant
+-- symbolic value, which is returned on exit. */
+SYMBOL *eval_symbolic(MPL *mpl, CODE *code)
+{     SYMBOL *value;
+xassert(code != NULL);
+xassert(code->type == A_SYMBOLIC);
+xassert(code->dim == 0);
+/* if the operation has a side effect, invalidate and delete the
+resultant value */
+if (code->vflag && code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* if resultant value is valid, no evaluation is needed */
+if (code->valid)
+{  value = copy_symbol(mpl, code->value.sym);
+goto done;
+}
+/* evaluate pseudo-code recursively */
+switch (code->op)
+{  case O_STRING:
+/* take character string */
+value = create_symbol_str(mpl, create_string(mpl,
+code->arg.str));
+break;
+case O_INDEX:
+/* take dummy index */
+xassert(code->arg.index.slot->value != NULL);
+value = copy_symbol(mpl, code->arg.index.slot->value);
+break;
+case O_MEMSYM:
+/* take member of symbolic parameter */
+{  TUPLE *tuple;
+ARG_LIST *e;
+tuple = create_tuple(mpl);
+for (e = code->arg.par.list; e != NULL; e = e->next)
+tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
+e->x));
+value = eval_member_sym(mpl, code->arg.par.par, tuple);
+delete_tuple(mpl, tuple);
+}
+break;
+case O_CVTSYM:
+/* conversion to symbolic */
+value = create_symbol_num(mpl, eval_numeric(mpl,
+code->arg.arg.x));
+break;
+case O_CONCAT:
+/* concatenation */
+value = concat_symbols(mpl,
+eval_symbolic(mpl, code->arg.arg.x),
+eval_symbolic(mpl, code->arg.arg.y));
+break;
+case O_FORK:
+/* if-then-else */
+if (eval_logical(mpl, code->arg.arg.x))
+value = eval_symbolic(mpl, code->arg.arg.y);
+else if (code->arg.arg.z == NULL)
+value = create_symbol_num(mpl, 0.0);
+else
+value = eval_symbolic(mpl, code->arg.arg.z);
+break;
+case O_SUBSTR:
+case O_SUBSTR3:
+{  double pos, len;
+char str[MAX_LENGTH+1];
+value = eval_symbolic(mpl, code->arg.arg.x);
+if (value->str == NULL)
+sprintf(str, "%.*g", DBL_DIG, value->num);
+else
+fetch_string(mpl, value->str, str);
+delete_symbol(mpl, value);
+if (code->op == O_SUBSTR)
+{  pos = eval_numeric(mpl, code->arg.arg.y);
+if (pos != floor(pos))
+error(mpl, "substr('...', %.*g); non-integer secon"
+"d argument", DBL_DIG, pos);
+if (pos < 1 || pos > strlen(str) + 1)
+error(mpl, "substr('...', %.*g); substring out of "
+"range", DBL_DIG, pos);
+}
+else
+{  pos = eval_numeric(mpl, code->arg.arg.y);
+len = eval_numeric(mpl, code->arg.arg.z);
+if (pos != floor(pos) || len != floor(len))
+error(mpl, "substr('...', %.*g, %.*g); non-integer"
+" second and/or third argument", DBL_DIG, pos,
+DBL_DIG, len);
+if (pos < 1 || len < 0 || pos + len > strlen(str) + 1)
+error(mpl, "substr('...', %.*g, %.*g); substring o"
+"ut of range", DBL_DIG, pos, DBL_DIG, len);
+str[(int)pos + (int)len - 1] = '\0';
+}
+value = create_symbol_str(mpl, create_string(mpl, str +
+(int)pos - 1));
+}
+break;
+case O_TIME2STR:
+{  double num;
+SYMBOL *sym;
+char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];
+num = eval_numeric(mpl, code->arg.arg.x);
+sym = eval_symbolic(mpl, code->arg.arg.y);
+if (sym->str == NULL)
+sprintf(fmt, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, fmt);
+delete_symbol(mpl, sym);
+fn_time2str(mpl, str, num, fmt);
+value = create_symbol_str(mpl, create_string(mpl, str));
+}
+break;
+default:
+xassert(code != code);
+}
+/* save resultant value */
+xassert(!code->valid);
+code->valid = 1;
+code->value.sym = copy_symbol(mpl, value);
+done: return value;
+}
+/*----------------------------------------------------------------------
+-- eval_logical - evaluate pseudo-code to determine logical value.
+--
+-- This routine evaluates specified pseudo-code to determine resultant
+-- logical value, which is returned on exit. */
+struct iter_log_info
+{     /* working info used by the routine iter_log_func */
+CODE *code;
+/* pseudo-code for iterated operation to be performed */
+int value;
+/* resultant value */
+};
+static int iter_log_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine used to perform iterated operation
+on logical "integrand" within domain scope */
+struct iter_log_info *info = _info;
+int ret = 0;
+switch (info->code->op)
+{  case O_FORALL:
+/* conjunction over domain */
+info->value &= eval_logical(mpl, info->code->arg.loop.x);
+if (!info->value) ret = 1;
+break;
+case O_EXISTS:
+/* disjunction over domain */
+info->value |= eval_logical(mpl, info->code->arg.loop.x);
+if (info->value) ret = 1;
+break;
+default:
+xassert(info != info);
+}
+return ret;
+}
+int eval_logical(MPL *mpl, CODE *code)
+{     int value;
+xassert(code->type == A_LOGICAL);
+xassert(code->dim == 0);
+/* if the operation has a side effect, invalidate and delete the
+resultant value */
+if (code->vflag && code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* if resultant value is valid, no evaluation is needed */
+if (code->valid)
+{  value = code->value.bit;
+goto done;
+}
+/* evaluate pseudo-code recursively */
+switch (code->op)
+{  case O_CVTLOG:
+/* conversion to logical */
+value = (eval_numeric(mpl, code->arg.arg.x) != 0.0);
+break;
+case O_NOT:
+/* negation (logical "not") */
+value = !eval_logical(mpl, code->arg.arg.x);
+break;
+case O_LT:
+/* comparison on 'less than' */
+#if 0 /* 02/VIII-2008 */
+value = (eval_numeric(mpl, code->arg.arg.x) <
+eval_numeric(mpl, code->arg.arg.y));
+#else
+xassert(code->arg.arg.x != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+value = (eval_numeric(mpl, code->arg.arg.x) <
+eval_numeric(mpl, code->arg.arg.y));
+else
+{  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
+SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
+value = (compare_symbols(mpl, sym1, sym2) < 0);
+delete_symbol(mpl, sym1);
+delete_symbol(mpl, sym2);
+}
+#endif
+break;
+case O_LE:
+/* comparison on 'not greater than' */
+#if 0 /* 02/VIII-2008 */
+value = (eval_numeric(mpl, code->arg.arg.x) <=
+eval_numeric(mpl, code->arg.arg.y));
+#else
+xassert(code->arg.arg.x != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+value = (eval_numeric(mpl, code->arg.arg.x) <=
+eval_numeric(mpl, code->arg.arg.y));
+else
+{  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
+SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
+value = (compare_symbols(mpl, sym1, sym2) <= 0);
+delete_symbol(mpl, sym1);
+delete_symbol(mpl, sym2);
+}
+#endif
+break;
+case O_EQ:
+/* comparison on 'equal to' */
+xassert(code->arg.arg.x != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+value = (eval_numeric(mpl, code->arg.arg.x) ==
+eval_numeric(mpl, code->arg.arg.y));
+else
+{  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
+SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
+value = (compare_symbols(mpl, sym1, sym2) == 0);
+delete_symbol(mpl, sym1);
+delete_symbol(mpl, sym2);
+}
+break;
+case O_GE:
+/* comparison on 'not less than' */
+#if 0 /* 02/VIII-2008 */
+value = (eval_numeric(mpl, code->arg.arg.x) >=
+eval_numeric(mpl, code->arg.arg.y));
+#else
+xassert(code->arg.arg.x != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+value = (eval_numeric(mpl, code->arg.arg.x) >=
+eval_numeric(mpl, code->arg.arg.y));
+else
+{  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
+SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
+value = (compare_symbols(mpl, sym1, sym2) >= 0);
+delete_symbol(mpl, sym1);
+delete_symbol(mpl, sym2);
+}
+#endif
+break;
+case O_GT:
+/* comparison on 'greater than' */
+#if 0 /* 02/VIII-2008 */
+value = (eval_numeric(mpl, code->arg.arg.x) >
+eval_numeric(mpl, code->arg.arg.y));
+#else
+xassert(code->arg.arg.x != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+value = (eval_numeric(mpl, code->arg.arg.x) >
+eval_numeric(mpl, code->arg.arg.y));
+else
+{  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
+SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
+value = (compare_symbols(mpl, sym1, sym2) > 0);
+delete_symbol(mpl, sym1);
+delete_symbol(mpl, sym2);
+}
+#endif
+break;
+case O_NE:
+/* comparison on 'not equal to' */
+xassert(code->arg.arg.x != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+value = (eval_numeric(mpl, code->arg.arg.x) !=
+eval_numeric(mpl, code->arg.arg.y));
+else
+{  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
+SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
+value = (compare_symbols(mpl, sym1, sym2) != 0);
+delete_symbol(mpl, sym1);
+delete_symbol(mpl, sym2);
+}
+break;
+case O_AND:
+/* conjunction (logical "and") */
+value = eval_logical(mpl, code->arg.arg.x) &&
+eval_logical(mpl, code->arg.arg.y);
+break;
+case O_OR:
+/* disjunction (logical "or") */
+value = eval_logical(mpl, code->arg.arg.x) ||
+eval_logical(mpl, code->arg.arg.y);
+break;
+case O_IN:
+/* test on 'x in Y' */
+{  TUPLE *tuple;
+tuple = eval_tuple(mpl, code->arg.arg.x);
+value = is_member(mpl, code->arg.arg.y, tuple);
+delete_tuple(mpl, tuple);
+}
+break;
+case O_NOTIN:
+/* test on 'x not in Y' */
+{  TUPLE *tuple;
+tuple = eval_tuple(mpl, code->arg.arg.x);
+value = !is_member(mpl, code->arg.arg.y, tuple);
+delete_tuple(mpl, tuple);
+}
+break;
+case O_WITHIN:
+/* test on 'X within Y' */
+{  ELEMSET *set;
+MEMBER *memb;
+set = eval_elemset(mpl, code->arg.arg.x);
+value = 1;
+for (memb = set->head; memb != NULL; memb = memb->next)
+{  if (!is_member(mpl, code->arg.arg.y, memb->tuple))
+{  value = 0;
+break;
+}
+}
+delete_elemset(mpl, set);
+}
+break;
+case O_NOTWITHIN:
+/* test on 'X not within Y' */
+{  ELEMSET *set;
+MEMBER *memb;
+set = eval_elemset(mpl, code->arg.arg.x);
+value = 1;
+for (memb = set->head; memb != NULL; memb = memb->next)
+{  if (is_member(mpl, code->arg.arg.y, memb->tuple))
+{  value = 0;
+break;
+}
+}
+delete_elemset(mpl, set);
+}
+break;
+case O_FORALL:
+/* conjunction (A-quantification) */
+{  struct iter_log_info _info, *info = &_info;
+info->code = code;
+info->value = 1;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_log_func);
+value = info->value;
+}
+break;
+case O_EXISTS:
+/* disjunction (E-quantification) */
+{  struct iter_log_info _info, *info = &_info;
+info->code = code;
+info->value = 0;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_log_func);
+value = info->value;
+}
+break;
+default:
+xassert(code != code);
+}
+/* save resultant value */
+xassert(!code->valid);
+code->valid = 1;
+code->value.bit = value;
+done: return value;
+}
+/*----------------------------------------------------------------------
+-- eval_tuple - evaluate pseudo-code to construct n-tuple.
+--
+-- This routine evaluates specified pseudo-code to construct resultant
+-- n-tuple, which is returned on exit. */
+TUPLE *eval_tuple(MPL *mpl, CODE *code)
+{     TUPLE *value;
+xassert(code != NULL);
+xassert(code->type == A_TUPLE);
+xassert(code->dim > 0);
+/* if the operation has a side effect, invalidate and delete the
+resultant value */
+if (code->vflag && code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* if resultant value is valid, no evaluation is needed */
+if (code->valid)
+{  value = copy_tuple(mpl, code->value.tuple);
+goto done;
+}
+/* evaluate pseudo-code recursively */
+switch (code->op)
+{  case O_TUPLE:
+/* make n-tuple */
+{  ARG_LIST *e;
+value = create_tuple(mpl);
+for (e = code->arg.list; e != NULL; e = e->next)
+value = expand_tuple(mpl, value, eval_symbolic(mpl,
+e->x));
+}
+break;
+case O_CVTTUP:
+/* convert to 1-tuple */
+value = expand_tuple(mpl, create_tuple(mpl),
+eval_symbolic(mpl, code->arg.arg.x));
+break;
+default:
+xassert(code != code);
+}
+/* save resultant value */
+xassert(!code->valid);
+code->valid = 1;
+code->value.tuple = copy_tuple(mpl, value);
+done: return value;
+}
+/*----------------------------------------------------------------------
+-- eval_elemset - evaluate pseudo-code to construct elemental set.
+--
+-- This routine evaluates specified pseudo-code to construct resultant
+-- elemental set, which is returned on exit. */
+struct iter_set_info
+{     /* working info used by the routine iter_set_func */
+CODE *code;
+/* pseudo-code for iterated operation to be performed */
+ELEMSET *value;
+/* resultant value */
+};
+static int iter_set_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine used to perform iterated operation
+on n-tuple "integrand" within domain scope */
+struct iter_set_info *info = _info;
+TUPLE *tuple;
+switch (info->code->op)
+{  case O_SETOF:
+/* compute next n-tuple and add it to the set; in this case
+duplicate n-tuples are silently ignored */
+tuple = eval_tuple(mpl, info->code->arg.loop.x);
+if (find_tuple(mpl, info->value, tuple) == NULL)
+add_tuple(mpl, info->value, tuple);
+else
+delete_tuple(mpl, tuple);
+break;
+case O_BUILD:
+/* construct next n-tuple using current values assigned to
+*free* dummy indices as its components and add it to the
+set; in this case duplicate n-tuples cannot appear */
+add_tuple(mpl, info->value, get_domain_tuple(mpl,
+info->code->arg.loop.domain));
+break;
+default:
+xassert(info != info);
+}
+return 0;
+}
+ELEMSET *eval_elemset(MPL *mpl, CODE *code)
+{     ELEMSET *value;
+xassert(code != NULL);
+xassert(code->type == A_ELEMSET);
+xassert(code->dim > 0);
+/* if the operation has a side effect, invalidate and delete the
+resultant value */
+if (code->vflag && code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* if resultant value is valid, no evaluation is needed */
+if (code->valid)
+{  value = copy_elemset(mpl, code->value.set);
+goto done;
+}
+/* evaluate pseudo-code recursively */
+switch (code->op)
+{  case O_MEMSET:
+/* take member of set */
+{  TUPLE *tuple;
+ARG_LIST *e;
+tuple = create_tuple(mpl);
+for (e = code->arg.set.list; e != NULL; e = e->next)
+tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
+e->x));
+value = copy_elemset(mpl,
+eval_member_set(mpl, code->arg.set.set, tuple));
+delete_tuple(mpl, tuple);
+}
+break;
+case O_MAKE:
+/* make elemental set of n-tuples */
+{  ARG_LIST *e;
+value = create_elemset(mpl, code->dim);
+for (e = code->arg.list; e != NULL; e = e->next)
+check_then_add(mpl, value, eval_tuple(mpl, e->x));
+}
+break;
+case O_UNION:
+/* union of two elemental sets */
+value = set_union(mpl,
+eval_elemset(mpl, code->arg.arg.x),
+eval_elemset(mpl, code->arg.arg.y));
+break;
+case O_DIFF:
+/* difference between two elemental sets */
+value = set_diff(mpl,
+eval_elemset(mpl, code->arg.arg.x),
+eval_elemset(mpl, code->arg.arg.y));
+break;
+case O_SYMDIFF:
+/* symmetric difference between two elemental sets */
+value = set_symdiff(mpl,
+eval_elemset(mpl, code->arg.arg.x),
+eval_elemset(mpl, code->arg.arg.y));
+break;
+case O_INTER:
+/* intersection of two elemental sets */
+value = set_inter(mpl,
+eval_elemset(mpl, code->arg.arg.x),
+eval_elemset(mpl, code->arg.arg.y));
+break;
+case O_CROSS:
+/* cross (Cartesian) product of two elemental sets */
+value = set_cross(mpl,
+eval_elemset(mpl, code->arg.arg.x),
+eval_elemset(mpl, code->arg.arg.y));
+break;
+case O_DOTS:
+/* build "arithmetic" elemental set */
+value = create_arelset(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_numeric(mpl, code->arg.arg.y),
+code->arg.arg.z == NULL ? 1.0 : eval_numeric(mpl,
+code->arg.arg.z));
+break;
+case O_FORK:
+/* if-then-else */
+if (eval_logical(mpl, code->arg.arg.x))
+value = eval_elemset(mpl, code->arg.arg.y);
+else
+value = eval_elemset(mpl, code->arg.arg.z);
+break;
+case O_SETOF:
+/* compute elemental set */
+{  struct iter_set_info _info, *info = &_info;
+info->code = code;
+info->value = create_elemset(mpl, code->dim);
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_set_func);
+value = info->value;
+}
+break;
+case O_BUILD:
+/* build elemental set identical to domain set */
+{  struct iter_set_info _info, *info = &_info;
+info->code = code;
+info->value = create_elemset(mpl, code->dim);
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_set_func);
+value = info->value;
+}
+break;
+default:
+xassert(code != code);
+}
+/* save resultant value */
+xassert(!code->valid);
+code->valid = 1;
+code->value.set = copy_elemset(mpl, value);
+done: return value;
+}
+/*----------------------------------------------------------------------
+-- is_member - check if n-tuple is in set specified by pseudo-code.
+--
+-- This routine checks if given n-tuple is a member of elemental set
+-- specified in the form of pseudo-code (i.e. by expression).
+--
+-- The n-tuple may have more components that dimension of the elemental
+-- set, in which case the extra components are ignored. */
+static void null_func(MPL *mpl, void *info)
+{     /* this is dummy routine used to enter the domain scope */
+xassert(mpl == mpl);
+xassert(info == NULL);
+return;
+}
+int is_member(MPL *mpl, CODE *code, TUPLE *tuple)
+{     int value;
+xassert(code != NULL);
+xassert(code->type == A_ELEMSET);
+xassert(code->dim > 0);
+xassert(tuple != NULL);
+switch (code->op)
+{  case O_MEMSET:
+/* check if given n-tuple is member of elemental set, which
+is assigned to member of model set */
+{  ARG_LIST *e;
+TUPLE *temp;
+ELEMSET *set;
+/* evaluate reference to elemental set */
+temp = create_tuple(mpl);
+for (e = code->arg.set.list; e != NULL; e = e->next)
+temp = expand_tuple(mpl, temp, eval_symbolic(mpl,
+e->x));
+set = eval_member_set(mpl, code->arg.set.set, temp);
+delete_tuple(mpl, temp);
+/* check if the n-tuple is contained in the set array */
+temp = build_subtuple(mpl, tuple, set->dim);
+value = (find_tuple(mpl, set, temp) != NULL);
+delete_tuple(mpl, temp);
+}
+break;
+case O_MAKE:
+/* check if given n-tuple is member of literal set */
+{  ARG_LIST *e;
+TUPLE *temp, *that;
+value = 0;
+temp = build_subtuple(mpl, tuple, code->dim);
+for (e = code->arg.list; e != NULL; e = e->next)
+{  that = eval_tuple(mpl, e->x);
+value = (compare_tuples(mpl, temp, that) == 0);
+delete_tuple(mpl, that);
+if (value) break;
+}
+delete_tuple(mpl, temp);
+}
+break;
+case O_UNION:
+value = is_member(mpl, code->arg.arg.x, tuple) ||
+is_member(mpl, code->arg.arg.y, tuple);
+break;
+case O_DIFF:
+value = is_member(mpl, code->arg.arg.x, tuple) &&
+!is_member(mpl, code->arg.arg.y, tuple);
+break;
+case O_SYMDIFF:
+{  int in1 = is_member(mpl, code->arg.arg.x, tuple);
+int in2 = is_member(mpl, code->arg.arg.y, tuple);
+value = (in1 && !in2) || (!in1 && in2);
+}
+break;
+case O_INTER:
+value = is_member(mpl, code->arg.arg.x, tuple) &&
+is_member(mpl, code->arg.arg.y, tuple);
+break;
+case O_CROSS:
+{  int j;
+value = is_member(mpl, code->arg.arg.x, tuple);
+if (value)
+{  for (j = 1; j <= code->arg.arg.x->dim; j++)
+{  xassert(tuple != NULL);
+tuple = tuple->next;
+}
+value = is_member(mpl, code->arg.arg.y, tuple);
+}
+}
+break;
+case O_DOTS:
+/* check if given 1-tuple is member of "arithmetic" set */
+{  int j;
+double x, t0, tf, dt;
+xassert(code->dim == 1);
+/* compute "parameters" of the "arithmetic" set */
+t0 = eval_numeric(mpl, code->arg.arg.x);
+tf = eval_numeric(mpl, code->arg.arg.y);
+if (code->arg.arg.z == NULL)
+dt = 1.0;
+else
+dt = eval_numeric(mpl, code->arg.arg.z);
+/* make sure the parameters are correct */
+arelset_size(mpl, t0, tf, dt);
+/* if component of 1-tuple is symbolic, not numeric, the
+1-tuple cannot be member of "arithmetic" set */
+xassert(tuple->sym != NULL);
+if (tuple->sym->str != NULL)
+{  value = 0;
+break;
+}
+/* determine numeric value of the component */
+x = tuple->sym->num;
+/* if the component value is out of the set range, the
+1-tuple is not in the set */
+if (dt > 0.0 && !(t0 <= x && x <= tf) ||
+dt < 0.0 && !(tf <= x && x <= t0))
+{  value = 0;
+break;
+}
+/* estimate ordinal number of the 1-tuple in the set */
+j = (int)(((x - t0) / dt) + 0.5) + 1;
+/* perform the main check */
+value = (arelset_member(mpl, t0, tf, dt, j) == x);
+}
+break;
+case O_FORK:
+/* check if given n-tuple is member of conditional set */
+if (eval_logical(mpl, code->arg.arg.x))
+value = is_member(mpl, code->arg.arg.y, tuple);
+else
+value = is_member(mpl, code->arg.arg.z, tuple);
+break;
+case O_SETOF:
+/* check if given n-tuple is member of computed set */
+/* it is not clear how to efficiently perform the check not
+computing the entire elemental set :+( */
+error(mpl, "implementation restriction; in/within setof{} n"
+"ot allowed");
+break;
+case O_BUILD:
+/* check if given n-tuple is member of domain set */
+{  TUPLE *temp;
+temp = build_subtuple(mpl, tuple, code->dim);
+/* try to enter the domain scope; if it is successful,
+the n-tuple is in the domain set */
+value = (eval_within_domain(mpl, code->arg.loop.domain,
+temp, NULL, null_func) == 0);
+delete_tuple(mpl, temp);
+}
+break;
+default:
+xassert(code != code);
+}
+return value;
+}
+/*----------------------------------------------------------------------
+-- eval_formula - evaluate pseudo-code to construct linear form.
+--
+-- This routine evaluates specified pseudo-code to construct resultant
+-- linear form, which is returned on exit. */
+struct iter_form_info
+{     /* working info used by the routine iter_form_func */
+CODE *code;
+/* pseudo-code for iterated operation to be performed */
+FORMULA *value;
+/* resultant value */
+FORMULA *tail;
+/* pointer to the last term */
+};
+static int iter_form_func(MPL *mpl, void *_info)
+{     /* this is auxiliary routine used to perform iterated operation
+on linear form "integrand" within domain scope */
+struct iter_form_info *info = _info;
+switch (info->code->op)
+{  case O_SUM:
+/* summation over domain */
+#if 0
+info->value =
+linear_comb(mpl,
++1.0, info->value,
++1.0, eval_formula(mpl, info->code->arg.loop.x));
+#else
+/* the routine linear_comb needs to look through all terms
+of both linear forms to reduce identical terms, so using
+it here is not a good idea (for example, evaluation of
+sum{i in 1..n} x[i] required quadratic time); the better
+idea is to gather all terms of the integrand in one list
+and reduce identical terms only once after all terms of
+the resultant linear form have been evaluated */
+{  FORMULA *form, *term;
+form = eval_formula(mpl, info->code->arg.loop.x);
+if (info->value == NULL)
+{  xassert(info->tail == NULL);
+info->value = form;
+}
+else
+{  xassert(info->tail != NULL);
+info->tail->next = form;
+}
+for (term = form; term != NULL; term = term->next)
+info->tail = term;
+}
+#endif
+break;
+default:
+xassert(info != info);
+}
+return 0;
+}
+FORMULA *eval_formula(MPL *mpl, CODE *code)
+{     FORMULA *value;
+xassert(code != NULL);
+xassert(code->type == A_FORMULA);
+xassert(code->dim == 0);
+/* if the operation has a side effect, invalidate and delete the
+resultant value */
+if (code->vflag && code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* if resultant value is valid, no evaluation is needed */
+if (code->valid)
+{  value = copy_formula(mpl, code->value.form);
+goto done;
+}
+/* evaluate pseudo-code recursively */
+switch (code->op)
+{  case O_MEMVAR:
+/* take member of variable */
+{  TUPLE *tuple;
+ARG_LIST *e;
+tuple = create_tuple(mpl);
+for (e = code->arg.var.list; e != NULL; e = e->next)
+tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
+e->x));
+#if 1 /* 15/V-2010 */
+xassert(code->arg.var.suff == DOT_NONE);
+#endif
+value = single_variable(mpl,
+eval_member_var(mpl, code->arg.var.var, tuple));
+delete_tuple(mpl, tuple);
+}
+break;
+case O_CVTLFM:
+/* convert to linear form */
+value = constant_term(mpl, eval_numeric(mpl,
+code->arg.arg.x));
+break;
+case O_PLUS:
+/* unary plus */
+value = linear_comb(mpl,
+0.0, constant_term(mpl, 0.0),
++1.0, eval_formula(mpl, code->arg.arg.x));
+break;
+case O_MINUS:
+/* unary minus */
+value = linear_comb(mpl,
+0.0, constant_term(mpl, 0.0),
+-1.0, eval_formula(mpl, code->arg.arg.x));
+break;
+case O_ADD:
+/* addition */
+value = linear_comb(mpl,
++1.0, eval_formula(mpl, code->arg.arg.x),
++1.0, eval_formula(mpl, code->arg.arg.y));
+break;
+case O_SUB:
+/* subtraction */
+value = linear_comb(mpl,
++1.0, eval_formula(mpl, code->arg.arg.x),
+-1.0, eval_formula(mpl, code->arg.arg.y));
+break;
+case O_MUL:
+/* multiplication */
+xassert(code->arg.arg.x != NULL);
+xassert(code->arg.arg.y != NULL);
+if (code->arg.arg.x->type == A_NUMERIC)
+{  xassert(code->arg.arg.y->type == A_FORMULA);
+value = linear_comb(mpl,
+eval_numeric(mpl, code->arg.arg.x),
+eval_formula(mpl, code->arg.arg.y),
+0.0, constant_term(mpl, 0.0));
+}
+else
+{  xassert(code->arg.arg.x->type == A_FORMULA);
+xassert(code->arg.arg.y->type == A_NUMERIC);
+value = linear_comb(mpl,
+eval_numeric(mpl, code->arg.arg.y),
+eval_formula(mpl, code->arg.arg.x),
+0.0, constant_term(mpl, 0.0));
+}
+break;
+case O_DIV:
+/* division */
+value = linear_comb(mpl,
+fp_div(mpl, 1.0, eval_numeric(mpl, code->arg.arg.y)),
+eval_formula(mpl, code->arg.arg.x),
+0.0, constant_term(mpl, 0.0));
+break;
+case O_FORK:
+/* if-then-else */
+if (eval_logical(mpl, code->arg.arg.x))
+value = eval_formula(mpl, code->arg.arg.y);
+else if (code->arg.arg.z == NULL)
+value = constant_term(mpl, 0.0);
+else
+value = eval_formula(mpl, code->arg.arg.z);
+break;
+case O_SUM:
+/* summation over domain */
+{  struct iter_form_info _info, *info = &_info;
+info->code = code;
+info->value = constant_term(mpl, 0.0);
+info->tail = NULL;
+loop_within_domain(mpl, code->arg.loop.domain, info,
+iter_form_func);
+value = reduce_terms(mpl, info->value);
+}
+break;
+default:
+xassert(code != code);
+}
+/* save resultant value */
+xassert(!code->valid);
+code->valid = 1;
+code->value.form = copy_formula(mpl, value);
+done: return value;
+}
+/*----------------------------------------------------------------------
+-- clean_code - clean pseudo-code.
+--
+-- This routine recursively cleans specified pseudo-code that assumes
+-- deleting all temporary resultant values. */
+void clean_code(MPL *mpl, CODE *code)
+{     ARG_LIST *e;
+/* if no pseudo-code is specified, do nothing */
+if (code == NULL) goto done;
+/* if resultant value is valid (exists), delete it */
+if (code->valid)
+{  code->valid = 0;
+delete_value(mpl, code->type, &code->value);
+}
+/* recursively clean pseudo-code for operands */
+switch (code->op)
+{  case O_NUMBER:
+case O_STRING:
+case O_INDEX:
+break;
+case O_MEMNUM:
+case O_MEMSYM:
+for (e = code->arg.par.list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+break;
+case O_MEMSET:
+for (e = code->arg.set.list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+break;
+case O_MEMVAR:
+for (e = code->arg.var.list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+break;
+#if 1 /* 15/V-2010 */
+case O_MEMCON:
+for (e = code->arg.con.list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+break;
+#endif
+case O_TUPLE:
+case O_MAKE:
+for (e = code->arg.list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+break;
+case O_SLICE:
+xassert(code != code);
+case O_IRAND224:
+case O_UNIFORM01:
+case O_NORMAL01:
+case O_GMTIME:
+break;
+case O_CVTNUM:
+case O_CVTSYM:
+case O_CVTLOG:
+case O_CVTTUP:
+case O_CVTLFM:
+case O_PLUS:
+case O_MINUS:
+case O_NOT:
+case O_ABS:
+case O_CEIL:
+case O_FLOOR:
+case O_EXP:
+case O_LOG:
+case O_LOG10:
+case O_SQRT:
+case O_SIN:
+case O_COS:
+case O_ATAN:
+case O_ROUND:
+case O_TRUNC:
+case O_CARD:
+case O_LENGTH:
+/* unary operation */
+clean_code(mpl, code->arg.arg.x);
+break;
+case O_ADD:
+case O_SUB:
+case O_LESS:
+case O_MUL:
+case O_DIV:
+case O_IDIV:
+case O_MOD:
+case O_POWER:
+case O_ATAN2:
+case O_ROUND2:
+case O_TRUNC2:
+case O_UNIFORM:
+case O_NORMAL:
+case O_CONCAT:
+case O_LT:
+case O_LE:
+case O_EQ:
+case O_GE:
+case O_GT:
+case O_NE:
+case O_AND:
+case O_OR:
+case O_UNION:
+case O_DIFF:
+case O_SYMDIFF:
+case O_INTER:
+case O_CROSS:
+case O_IN:
+case O_NOTIN:
+case O_WITHIN:
+case O_NOTWITHIN:
+case O_SUBSTR:
+case O_STR2TIME:
+case O_TIME2STR:
+/* binary operation */
+clean_code(mpl, code->arg.arg.x);
+clean_code(mpl, code->arg.arg.y);
+break;
+case O_DOTS:
+case O_FORK:
+case O_SUBSTR3:
+/* ternary operation */
+clean_code(mpl, code->arg.arg.x);
+clean_code(mpl, code->arg.arg.y);
+clean_code(mpl, code->arg.arg.z);
+break;
+case O_MIN:
+case O_MAX:
+/* n-ary operation */
+for (e = code->arg.list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+break;
+case O_SUM:
+case O_PROD:
+case O_MINIMUM:
+case O_MAXIMUM:
+case O_FORALL:
+case O_EXISTS:
+case O_SETOF:
+case O_BUILD:
+/* iterated operation */
+clean_domain(mpl, code->arg.loop.domain);
+clean_code(mpl, code->arg.loop.x);
+break;
+default:
+xassert(code->op != code->op);
+}
+done: return;
+}
+#if 1 /* 11/II-2008 */
+/**********************************************************************/
+/* * *                        DATA TABLES                         * * */
+/**********************************************************************/
+int mpl_tab_num_args(TABDCA *dca)
+{     /* returns the number of arguments */
+return dca->na;
+}
+const char *mpl_tab_get_arg(TABDCA *dca, int k)
+{     /* returns pointer to k-th argument */
+xassert(1 <= k && k <= dca->na);
+return dca->arg[k];
+}
+int mpl_tab_num_flds(TABDCA *dca)
+{     /* returns the number of fields */
+return dca->nf;
+}
+const char *mpl_tab_get_name(TABDCA *dca, int k)
+{     /* returns pointer to name of k-th field */
+xassert(1 <= k && k <= dca->nf);
+return dca->name[k];
+}
+int mpl_tab_get_type(TABDCA *dca, int k)
+{     /* returns type of k-th field */
+xassert(1 <= k && k <= dca->nf);
+return dca->type[k];
+}
+double mpl_tab_get_num(TABDCA *dca, int k)
+{     /* returns numeric value of k-th field */
+xassert(1 <= k && k <= dca->nf);
+xassert(dca->type[k] == 'N');
+return dca->num[k];
+}
+const char *mpl_tab_get_str(TABDCA *dca, int k)
+{     /* returns pointer to string value of k-th field */
+xassert(1 <= k && k <= dca->nf);
+xassert(dca->type[k] == 'S');
+xassert(dca->str[k] != NULL);
+return dca->str[k];
+}
+void mpl_tab_set_num(TABDCA *dca, int k, double num)
+{     /* assign numeric value to k-th field */
+xassert(1 <= k && k <= dca->nf);
+xassert(dca->type[k] == '?');
+dca->type[k] = 'N';
+dca->num[k] = num;
+return;
+}
+void mpl_tab_set_str(TABDCA *dca, int k, const char *str)
+{     /* assign string value to k-th field */
+xassert(1 <= k && k <= dca->nf);
+xassert(dca->type[k] == '?');
+xassert(strlen(str) <= MAX_LENGTH);
+xassert(dca->str[k] != NULL);
+dca->type[k] = 'S';
+strcpy(dca->str[k], str);
+return;
+}
+static int write_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+TABLE *tab = info;
+TABDCA *dca = mpl->dca;
+TABOUT *out;
+SYMBOL *sym;
+int k;
+char buf[MAX_LENGTH+1];
+/* evaluate field values */
+k = 0;
+for (out = tab->u.out.list; out != NULL; out = out->next)
+{  k++;
+switch (out->code->type)
+{  case A_NUMERIC:
+dca->type[k] = 'N';
+dca->num[k] = eval_numeric(mpl, out->code);
+dca->str[k][0] = '\0';
+break;
+case A_SYMBOLIC:
+sym = eval_symbolic(mpl, out->code);
+if (sym->str == NULL)
+{  dca->type[k] = 'N';
+dca->num[k] = sym->num;
+dca->str[k][0] = '\0';
+}
+else
+{  dca->type[k] = 'S';
+dca->num[k] = 0.0;
+fetch_string(mpl, sym->str, buf);
+strcpy(dca->str[k], buf);
+}
+delete_symbol(mpl, sym);
+break;
+default:
+xassert(out != out);
+}
+}
+/* write record to output table */
+mpl_tab_drv_write(mpl);
+return 0;
+}
+void execute_table(MPL *mpl, TABLE *tab)
+{     /* execute table statement */
+TABARG *arg;
+TABFLD *fld;
+TABIN *in;
+TABOUT *out;
+TABDCA *dca;
+SET *set;
+int k;
+char buf[MAX_LENGTH+1];
+/* allocate table driver communication area */
+xassert(mpl->dca == NULL);
+mpl->dca = dca = xmalloc(sizeof(TABDCA));
+dca->id = 0;
+dca->link = NULL;
+dca->na = 0;
+dca->arg = NULL;
+dca->nf = 0;
+dca->name = NULL;
+dca->type = NULL;
+dca->num = NULL;
+dca->str = NULL;
+/* allocate arguments */
+xassert(dca->na == 0);
+for (arg = tab->arg; arg != NULL; arg = arg->next)
+dca->na++;
+dca->arg = xcalloc(1+dca->na, sizeof(char *));
+#if 1 /* 28/IX-2008 */
+for (k = 1; k <= dca->na; k++) dca->arg[k] = NULL;
+#endif
+/* evaluate argument values */
+k = 0;
+for (arg = tab->arg; arg != NULL; arg = arg->next)
+{  SYMBOL *sym;
+k++;
+xassert(arg->code->type == A_SYMBOLIC);
+sym = eval_symbolic(mpl, arg->code);
+if (sym->str == NULL)
+sprintf(buf, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, buf);
+delete_symbol(mpl, sym);
+dca->arg[k] = xmalloc(strlen(buf)+1);
+strcpy(dca->arg[k], buf);
+}
+/* perform table input/output */
+switch (tab->type)
+{  case A_INPUT:  goto read_table;
+case A_OUTPUT: goto write_table;
+default:       xassert(tab != tab);
+}
+read_table:
+/* read data from input table */
+/* add the only member to the control set and assign it empty
+elemental set */
+set = tab->u.in.set;
+if (set != NULL)
+{  if (set->data)
+error(mpl, "%s already provided with data", set->name);
+xassert(set->array->head == NULL);
+add_member(mpl, set->array, NULL)->value.set =
+create_elemset(mpl, set->dimen);
+set->data = 1;
+}
+/* check parameters specified in the input list */
+for (in = tab->u.in.list; in != NULL; in = in->next)
+{  if (in->par->data)
+error(mpl, "%s already provided with data", in->par->name);
+in->par->data = 1;
+}
+/* allocate and initialize fields */
+xassert(dca->nf == 0);
+for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
+dca->nf++;
+for (in = tab->u.in.list; in != NULL; in = in->next)
+dca->nf++;
+dca->name = xcalloc(1+dca->nf, sizeof(char *));
+dca->type = xcalloc(1+dca->nf, sizeof(int));
+dca->num = xcalloc(1+dca->nf, sizeof(double));
+dca->str = xcalloc(1+dca->nf, sizeof(char *));
+k = 0;
+for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
+{  k++;
+dca->name[k] = fld->name;
+dca->type[k] = '?';
+dca->num[k] = 0.0;
+dca->str[k] = xmalloc(MAX_LENGTH+1);
+dca->str[k][0] = '\0';
+}
+for (in = tab->u.in.list; in != NULL; in = in->next)
+{  k++;
+dca->name[k] = in->name;
+dca->type[k] = '?';
+dca->num[k] = 0.0;
+dca->str[k] = xmalloc(MAX_LENGTH+1);
+dca->str[k][0] = '\0';
+}
+/* open input table */
+mpl_tab_drv_open(mpl, 'R');
+/* read and process records */
+for (;;)
+{  TUPLE *tup;
+/* reset field types */
+for (k = 1; k <= dca->nf; k++)
+dca->type[k] = '?';
+/* read next record */
+if (mpl_tab_drv_read(mpl)) break;
+/* all fields must be set by the driver */
+for (k = 1; k <= dca->nf; k++)
+{  if (dca->type[k] == '?')
+error(mpl, "field %s missing in input table",
+dca->name[k]);
+}
+/* construct n-tuple */
+tup = create_tuple(mpl);
+k = 0;
+for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
+{  k++;
+xassert(k <= dca->nf);
+switch (dca->type[k])
+{  case 'N':
+tup = expand_tuple(mpl, tup, create_symbol_num(mpl,
+dca->num[k]));
+break;
+case 'S':
+xassert(strlen(dca->str[k]) <= MAX_LENGTH);
+tup = expand_tuple(mpl, tup, create_symbol_str(mpl,
+create_string(mpl, dca->str[k])));
+break;
+default:
+xassert(dca != dca);
+}
+}
+/* add n-tuple just read to the control set */
+if (tab->u.in.set != NULL)
+check_then_add(mpl, tab->u.in.set->array->head->value.set,
+copy_tuple(mpl, tup));
+/* assign values to the parameters in the input list */
+for (in = tab->u.in.list; in != NULL; in = in->next)
+{  MEMBER *memb;
+k++;
+xassert(k <= dca->nf);
+/* there must be no member with the same n-tuple */
+if (find_member(mpl, in->par->array, tup) != NULL)
+error(mpl, "%s%s already defined", in->par->name,
+format_tuple(mpl, '[', tup));
+/* create new parameter member with given n-tuple */
+memb = add_member(mpl, in->par->array, copy_tuple(mpl, tup))
+;
+/* assign value to the parameter member */
+switch (in->par->type)
+{  case A_NUMERIC:
+case A_INTEGER:
+case A_BINARY:
+if (dca->type[k] != 'N')
+error(mpl, "%s requires numeric data",
+in->par->name);
+memb->value.num = dca->num[k];
+break;
+case A_SYMBOLIC:
+switch (dca->type[k])
+{  case 'N':
+memb->value.sym = create_symbol_num(mpl,
+dca->num[k]);
+break;
+case 'S':
+xassert(strlen(dca->str[k]) <= MAX_LENGTH);
+memb->value.sym = create_symbol_str(mpl,
+create_string(mpl,dca->str[k]));
+break;
+default:
+xassert(dca != dca);
+}
+break;
+default:
+xassert(in != in);
+}
+}
+/* n-tuple is no more needed */
+delete_tuple(mpl, tup);
+}
+/* close input table */
+mpl_tab_drv_close(mpl);
+goto done;
+write_table:
+/* write data to output table */
+/* allocate and initialize fields */
+xassert(dca->nf == 0);
+for (out = tab->u.out.list; out != NULL; out = out->next)
+dca->nf++;
+dca->name = xcalloc(1+dca->nf, sizeof(char *));
+dca->type = xcalloc(1+dca->nf, sizeof(int));
+dca->num = xcalloc(1+dca->nf, sizeof(double));
+dca->str = xcalloc(1+dca->nf, sizeof(char *));
+k = 0;
+for (out = tab->u.out.list; out != NULL; out = out->next)
+{  k++;
+dca->name[k] = out->name;
+dca->type[k] = '?';
+dca->num[k] = 0.0;
+dca->str[k] = xmalloc(MAX_LENGTH+1);
+dca->str[k][0] = '\0';
+}
+/* open output table */
+mpl_tab_drv_open(mpl, 'W');
+/* evaluate fields and write records */
+loop_within_domain(mpl, tab->u.out.domain, tab, write_func);
+/* close output table */
+mpl_tab_drv_close(mpl);
+done: /* free table driver communication area */
+free_dca(mpl);
+return;
+}
+void free_dca(MPL *mpl)
+{     /* free table driver communucation area */
+TABDCA *dca = mpl->dca;
+int k;
+if (dca != NULL)
+{  if (dca->link != NULL)
+mpl_tab_drv_close(mpl);
+if (dca->arg != NULL)
+{  for (k = 1; k <= dca->na; k++)
+#if 1 /* 28/IX-2008 */
+if (dca->arg[k] != NULL)
+#endif
+xfree(dca->arg[k]);
+xfree(dca->arg);
+}
+if (dca->name != NULL) xfree(dca->name);
+if (dca->type != NULL) xfree(dca->type);
+if (dca->num != NULL) xfree(dca->num);
+if (dca->str != NULL)
+{  for (k = 1; k <= dca->nf; k++)
+xfree(dca->str[k]);
+xfree(dca->str);
+}
+xfree(dca), mpl->dca = NULL;
+}
+return;
+}
+void clean_table(MPL *mpl, TABLE *tab)
+{     /* clean table statement */
+TABARG *arg;
+TABOUT *out;
+/* clean string list */
+for (arg = tab->arg; arg != NULL; arg = arg->next)
+clean_code(mpl, arg->code);
+switch (tab->type)
+{  case A_INPUT:
+break;
+case A_OUTPUT:
+/* clean subscript domain */
+clean_domain(mpl, tab->u.out.domain);
+/* clean output list */
+for (out = tab->u.out.list; out != NULL; out = out->next)
+clean_code(mpl, out->code);
+break;
+default:
+xassert(tab != tab);
+}
+return;
+}
+#endif
+/**********************************************************************/
+/* * *                      MODEL STATEMENTS                      * * */
+/**********************************************************************/
+/*----------------------------------------------------------------------
+-- execute_check - execute check statement.
+--
+-- This routine executes specified check statement. */
+static int check_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+CHECK *chk = (CHECK *)info;
+if (!eval_logical(mpl, chk->code))
+error(mpl, "check%s failed", format_tuple(mpl, '[',
+get_domain_tuple(mpl, chk->domain)));
+return 0;
+}
+void execute_check(MPL *mpl, CHECK *chk)
+{     loop_within_domain(mpl, chk->domain, chk, check_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_check - clean check statement.
+--
+-- This routine cleans specified check statement that assumes deleting
+-- all stuff dynamically allocated on generating/postsolving phase. */
+void clean_check(MPL *mpl, CHECK *chk)
+{     /* clean subscript domain */
+clean_domain(mpl, chk->domain);
+/* clean pseudo-code for computing predicate */
+clean_code(mpl, chk->code);
+return;
+}
+/*----------------------------------------------------------------------
+-- execute_display - execute display statement.
+--
+-- This routine executes specified display statement. */
+static void display_set(MPL *mpl, SET *set, MEMBER *memb)
+{     /* display member of model set */
+ELEMSET *s = memb->value.set;
+MEMBER *m;
+write_text(mpl, "%s%s%s\n", set->name,
+format_tuple(mpl, '[', memb->tuple),
+s->head == NULL ? " is empty" : ":");
+for (m = s->head; m != NULL; m = m->next)
+write_text(mpl, "   %s\n", format_tuple(mpl, '(', m->tuple));
+return;
+}
+static void display_par(MPL *mpl, PARAMETER *par, MEMBER *memb)
+{     /* display member of model parameter */
+switch (par->type)
+{  case A_NUMERIC:
+case A_INTEGER:
+case A_BINARY:
+write_text(mpl, "%s%s = %.*g\n", par->name,
+format_tuple(mpl, '[', memb->tuple),
+DBL_DIG, memb->value.num);
+break;
+case A_SYMBOLIC:
+write_text(mpl, "%s%s = %s\n", par->name,
+format_tuple(mpl, '[', memb->tuple),
+format_symbol(mpl, memb->value.sym));
+break;
+default:
+xassert(par != par);
+}
+return;
+}
+#if 1 /* 15/V-2010 */
+static void display_var(MPL *mpl, VARIABLE *var, MEMBER *memb,
+int suff)
+{     /* display member of model variable */
+if (suff == DOT_NONE || suff == DOT_VAL)
+write_text(mpl, "%s%s.val = %.*g\n", var->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.var->prim);
+else if (suff == DOT_LB)
+write_text(mpl, "%s%s.lb = %.*g\n", var->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.var->var->lbnd == NULL ? -DBL_MAX :
+memb->value.var->lbnd);
+else if (suff == DOT_UB)
+write_text(mpl, "%s%s.ub = %.*g\n", var->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.var->var->ubnd == NULL ? +DBL_MAX :
+memb->value.var->ubnd);
+else if (suff == DOT_STATUS)
+write_text(mpl, "%s%s.status = %d\n", var->name, format_tuple
+(mpl, '[', memb->tuple), memb->value.var->stat);
+else if (suff == DOT_DUAL)
+write_text(mpl, "%s%s.dual = %.*g\n", var->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.var->dual);
+else
+xassert(suff != suff);
+return;
+}
+#endif
+#if 1 /* 15/V-2010 */
+static void display_con(MPL *mpl, CONSTRAINT *con, MEMBER *memb,
+int suff)
+{     /* display member of model constraint */
+if (suff == DOT_NONE || suff == DOT_VAL)
+write_text(mpl, "%s%s.val = %.*g\n", con->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.con->prim);
+else if (suff == DOT_LB)
+write_text(mpl, "%s%s.lb = %.*g\n", con->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.con->con->lbnd == NULL ? -DBL_MAX :
+memb->value.con->lbnd);
+else if (suff == DOT_UB)
+write_text(mpl, "%s%s.ub = %.*g\n", con->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.con->con->ubnd == NULL ? +DBL_MAX :
+memb->value.con->ubnd);
+else if (suff == DOT_STATUS)
+write_text(mpl, "%s%s.status = %d\n", con->name, format_tuple
+(mpl, '[', memb->tuple), memb->value.con->stat);
+else if (suff == DOT_DUAL)
+write_text(mpl, "%s%s.dual = %.*g\n", con->name,
+format_tuple(mpl, '[', memb->tuple), DBL_DIG,
+memb->value.con->dual);
+else
+xassert(suff != suff);
+return;
+}
+#endif
+static void display_memb(MPL *mpl, CODE *code)
+{     /* display member specified by pseudo-code */
+MEMBER memb;
+ARG_LIST *e;
+xassert(code->op == O_MEMNUM || code->op == O_MEMSYM
+|| code->op == O_MEMSET || code->op == O_MEMVAR
+|| code->op == O_MEMCON);
+memb.tuple = create_tuple(mpl);
+for (e = code->arg.par.list; e != NULL; e = e->next)
+memb.tuple = expand_tuple(mpl, memb.tuple, eval_symbolic(mpl,
+e->x));
+switch (code->op)
+{  case O_MEMNUM:
+memb.value.num = eval_member_num(mpl, code->arg.par.par,
+memb.tuple);
+display_par(mpl, code->arg.par.par, &memb);
+break;
+case O_MEMSYM:
+memb.value.sym = eval_member_sym(mpl, code->arg.par.par,
+memb.tuple);
+display_par(mpl, code->arg.par.par, &memb);
+delete_symbol(mpl, memb.value.sym);
+break;
+case O_MEMSET:
+memb.value.set = eval_member_set(mpl, code->arg.set.set,
+memb.tuple);
+display_set(mpl, code->arg.set.set, &memb);
+break;
+case O_MEMVAR:
+memb.value.var = eval_member_var(mpl, code->arg.var.var,
+memb.tuple);
+display_var
+(mpl, code->arg.var.var, &memb, code->arg.var.suff);
+break;
+case O_MEMCON:
+memb.value.con = eval_member_con(mpl, code->arg.con.con,
+memb.tuple);
+display_con
+(mpl, code->arg.con.con, &memb, code->arg.con.suff);
+break;
+default:
+xassert(code != code);
+}
+delete_tuple(mpl, memb.tuple);
+return;
+}
+static void display_code(MPL *mpl, CODE *code)
+{     /* display value of expression */
+switch (code->type)
+{  case A_NUMERIC:
+/* numeric value */
+{  double num;
+num = eval_numeric(mpl, code);
+write_text(mpl, "%.*g\n", DBL_DIG, num);
+}
+break;
+case A_SYMBOLIC:
+/* symbolic value */
+{  SYMBOL *sym;
+sym = eval_symbolic(mpl, code);
+write_text(mpl, "%s\n", format_symbol(mpl, sym));
+delete_symbol(mpl, sym);
+}
+break;
+case A_LOGICAL:
+/* logical value */
+{  int bit;
+bit = eval_logical(mpl, code);
+write_text(mpl, "%s\n", bit ? "true" : "false");
+}
+break;
+case A_TUPLE:
+/* n-tuple */
+{  TUPLE *tuple;
+tuple = eval_tuple(mpl, code);
+write_text(mpl, "%s\n", format_tuple(mpl, '(', tuple));
+delete_tuple(mpl, tuple);
+}
+break;
+case A_ELEMSET:
+/* elemental set */
+{  ELEMSET *set;
+MEMBER *memb;
+set = eval_elemset(mpl, code);
+if (set->head == 0)
+write_text(mpl, "set is empty\n");
+for (memb = set->head; memb != NULL; memb = memb->next)
+write_text(mpl, "   %s\n", format_tuple(mpl, '(',
+memb->tuple));
+delete_elemset(mpl, set);
+}
+break;
+case A_FORMULA:
+/* linear form */
+{  FORMULA *form, *term;
+form = eval_formula(mpl, code);
+if (form == NULL)
+write_text(mpl, "linear form is empty\n");
+for (term = form; term != NULL; term = term->next)
+{  if (term->var == NULL)
+write_text(mpl, "   %.*g\n", term->coef);
+else
+write_text(mpl, "   %.*g %s%s\n", DBL_DIG,
+term->coef, term->var->var->name,
+format_tuple(mpl, '[', term->var->memb->tuple));
+}
+delete_formula(mpl, form);
+}
+break;
+default:
+xassert(code != code);
+}
+return;
+}
+static int display_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+DISPLAY *dpy = (DISPLAY *)info;
+DISPLAY1 *entry;
+for (entry = dpy->list; entry != NULL; entry = entry->next)
+{  if (entry->type == A_INDEX)
+{  /* dummy index */
+DOMAIN_SLOT *slot = entry->u.slot;
+write_text(mpl, "%s = %s\n", slot->name,
+format_symbol(mpl, slot->value));
+}
+else if (entry->type == A_SET)
+{  /* model set */
+SET *set = entry->u.set;
+MEMBER *memb;
+if (set->assign != NULL)
+{  /* the set has assignment expression; evaluate all its
+members over entire domain */
+eval_whole_set(mpl, set);
+}
+else
+{  /* the set has no assignment expression; refer to its
+any existing member ignoring resultant value to check
+the data provided the data section */
+#if 1 /* 12/XII-2008 */
+if (set->gadget != NULL && set->data == 0)
+{  /* initialize the set with data from a plain set */
+saturate_set(mpl, set);
+}
+#endif
+if (set->array->head != NULL)
+eval_member_set(mpl, set, set->array->head->tuple);
+}
+/* display all members of the set array */
+if (set->array->head == NULL)
+write_text(mpl, "%s has empty content\n", set->name);
+for (memb = set->array->head; memb != NULL; memb =
+memb->next) display_set(mpl, set, memb);
+}
+else if (entry->type == A_PARAMETER)
+{  /* model parameter */
+PARAMETER *par = entry->u.par;
+MEMBER *memb;
+if (par->assign != NULL)
+{  /* the parameter has an assignment expression; evaluate
+all its member over entire domain */
+eval_whole_par(mpl, par);
+}
+else
+{  /* the parameter has no assignment expression; refer to
+its any existing member ignoring resultant value to
+check the data provided in the data section */
+if (par->array->head != NULL)
+{  if (par->type != A_SYMBOLIC)
+eval_member_num(mpl, par, par->array->head->tuple);
+else
+delete_symbol(mpl, eval_member_sym(mpl, par,
+par->array->head->tuple));
+}
+}
+/* display all members of the parameter array */
+if (par->array->head == NULL)
+write_text(mpl, "%s has empty content\n", par->name);
+for (memb = par->array->head; memb != NULL; memb =
+memb->next) display_par(mpl, par, memb);
+}
+else if (entry->type == A_VARIABLE)
+{  /* model variable */
+VARIABLE *var = entry->u.var;
+MEMBER *memb;
+xassert(mpl->flag_p);
+/* display all members of the variable array */
+if (var->array->head == NULL)
+write_text(mpl, "%s has empty content\n", var->name);
+for (memb = var->array->head; memb != NULL; memb =
+memb->next) display_var(mpl, var, memb, DOT_NONE);
+}
+else if (entry->type == A_CONSTRAINT)
+{  /* model constraint */
+CONSTRAINT *con = entry->u.con;
+MEMBER *memb;
+xassert(mpl->flag_p);
+/* display all members of the constraint array */
+if (con->array->head == NULL)
+write_text(mpl, "%s has empty content\n", con->name);
+for (memb = con->array->head; memb != NULL; memb =
+memb->next) display_con(mpl, con, memb, DOT_NONE);
+}
+else if (entry->type == A_EXPRESSION)
+{  /* expression */
+CODE *code = entry->u.code;
+if (code->op == O_MEMNUM || code->op == O_MEMSYM ||
+code->op == O_MEMSET || code->op == O_MEMVAR ||
+code->op == O_MEMCON)
+display_memb(mpl, code);
+else
+display_code(mpl, code);
+}
+else
+xassert(entry != entry);
+}
+return 0;
+}
+void execute_display(MPL *mpl, DISPLAY *dpy)
+{     loop_within_domain(mpl, dpy->domain, dpy, display_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_display - clean display statement.
+--
+-- This routine cleans specified display statement that assumes deleting
+-- all stuff dynamically allocated on generating/postsolving phase. */
+void clean_display(MPL *mpl, DISPLAY *dpy)
+{     DISPLAY1 *d;
+#if 0 /* 15/V-2010 */
+ARG_LIST *e;
+#endif
+/* clean subscript domain */
+clean_domain(mpl, dpy->domain);
+/* clean display list */
+for (d = dpy->list; d != NULL; d = d->next)
+{  /* clean pseudo-code for computing expression */
+if (d->type == A_EXPRESSION)
+clean_code(mpl, d->u.code);
+#if 0 /* 15/V-2010 */
+/* clean pseudo-code for computing subscripts */
+for (e = d->list; e != NULL; e = e->next)
+clean_code(mpl, e->x);
+#endif
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- execute_printf - execute printf statement.
+--
+-- This routine executes specified printf statement. */
+#if 1 /* 14/VII-2006 */
+static void print_char(MPL *mpl, int c)
+{     if (mpl->prt_fp == NULL)
+write_char(mpl, c);
+else
+xfputc(c, mpl->prt_fp);
+return;
+}
+static void print_text(MPL *mpl, char *fmt, ...)
+{     va_list arg;
+char buf[OUTBUF_SIZE], *c;
+va_start(arg, fmt);
+vsprintf(buf, fmt, arg);
+xassert(strlen(buf) < sizeof(buf));
+va_end(arg);
+for (c = buf; *c != '\0'; c++) print_char(mpl, *c);
+return;
+}
+#endif
+static int printf_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+PRINTF *prt = (PRINTF *)info;
+PRINTF1 *entry;
+SYMBOL *sym;
+char fmt[MAX_LENGTH+1], *c, *from, save;
+/* evaluate format control string */
+sym = eval_symbolic(mpl, prt->fmt);
+if (sym->str == NULL)
+sprintf(fmt, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, fmt);
+delete_symbol(mpl, sym);
+/* scan format control string and perform formatting output */
+entry = prt->list;
+for (c = fmt; *c != '\0'; c++)
+{  if (*c == '%')
+{  /* scan format specifier */
+from = c++;
+if (*c == '%')
+{  print_char(mpl, '%');
+continue;
+}
+if (entry == NULL) break;
+/* scan optional flags */
+while (*c == '-' || *c == '+' || *c == ' ' || *c == '#' ||
+*c == '0') c++;
+/* scan optional minimum field width */
+while (isdigit((unsigned char)*c)) c++;
+/* scan optional precision */
+if (*c == '.')
+{  c++;
+while (isdigit((unsigned char)*c)) c++;
+}
+/* scan conversion specifier and perform formatting */
+save = *(c+1), *(c+1) = '\0';
+if (*c == 'd' || *c == 'i' || *c == 'e' || *c == 'E' ||
+*c == 'f' || *c == 'F' || *c == 'g' || *c == 'G')
+{  /* the specifier requires numeric value */
+double value;
+xassert(entry != NULL);
+switch (entry->code->type)
+{  case A_NUMERIC:
+value = eval_numeric(mpl, entry->code);
+break;
+case A_SYMBOLIC:
+sym = eval_symbolic(mpl, entry->code);
+if (sym->str != NULL)
+error(mpl, "cannot convert %s to floating-point"
+" number", format_symbol(mpl, sym));
+value = sym->num;
+delete_symbol(mpl, sym);
+break;
+case A_LOGICAL:
+if (eval_logical(mpl, entry->code))
+value = 1.0;
+else
+value = 0.0;
+break;
+default:
+xassert(entry != entry);
+}
+if (*c == 'd' || *c == 'i')
+{  double int_max = (double)INT_MAX;
+if (!(-int_max <= value && value <= +int_max))
+error(mpl, "cannot convert %.*g to integer",
+DBL_DIG, value);
+print_text(mpl, from, (int)floor(value + 0.5));
+}
+else
+print_text(mpl, from, value);
+}
+else if (*c == 's')
+{  /* the specifier requires symbolic value */
+char value[MAX_LENGTH+1];
+switch (entry->code->type)
+{  case A_NUMERIC:
+sprintf(value, "%.*g", DBL_DIG, eval_numeric(mpl,
+entry->code));
+break;
+case A_LOGICAL:
+if (eval_logical(mpl, entry->code))
+strcpy(value, "T");
+else
+strcpy(value, "F");
+break;
+case A_SYMBOLIC:
+sym = eval_symbolic(mpl, entry->code);
+if (sym->str == NULL)
+sprintf(value, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, value);
+delete_symbol(mpl, sym);
+break;
+default:
+xassert(entry != entry);
+}
+print_text(mpl, from, value);
+}
+else
+error(mpl, "format specifier missing or invalid");
+*(c+1) = save;
+entry = entry->next;
+}
+else if (*c == '\\')
+{  /* write some control character */
+c++;
+if (*c == 't')
+print_char(mpl, '\t');
+else if (*c == 'n')
+print_char(mpl, '\n');
+#if 1 /* 28/X-2010 */
+else if (*c == '\0')
+{  /* format string ends with backslash */
+error(mpl, "invalid use of escape character \\ in format"
+" control string");
+}
+#endif
+else
+print_char(mpl, *c);
+}
+else
+{  /* write character without formatting */
+print_char(mpl, *c);
+}
+}
+return 0;
+}
+#if 0 /* 14/VII-2006 */
+void execute_printf(MPL *mpl, PRINTF *prt)
+{     loop_within_domain(mpl, prt->domain, prt, printf_func);
+return;
+}
+#else
+void execute_printf(MPL *mpl, PRINTF *prt)
+{     if (prt->fname == NULL)
+{  /* switch to the standard output */
+if (mpl->prt_fp != NULL)
+{  xfclose(mpl->prt_fp), mpl->prt_fp = NULL;
+xfree(mpl->prt_file), mpl->prt_file = NULL;
+}
+}
+else
+{  /* evaluate file name string */
+SYMBOL *sym;
+char fname[MAX_LENGTH+1];
+sym = eval_symbolic(mpl, prt->fname);
+if (sym->str == NULL)
+sprintf(fname, "%.*g", DBL_DIG, sym->num);
+else
+fetch_string(mpl, sym->str, fname);
+delete_symbol(mpl, sym);
+/* close the current print file, if necessary */
+if (mpl->prt_fp != NULL &&
+(!prt->app || strcmp(mpl->prt_file, fname) != 0))
+{  xfclose(mpl->prt_fp), mpl->prt_fp = NULL;
+xfree(mpl->prt_file), mpl->prt_file = NULL;
+}
+/* open the specified print file, if necessary */
+if (mpl->prt_fp == NULL)
+{  mpl->prt_fp = xfopen(fname, prt->app ? "a" : "w");
+if (mpl->prt_fp == NULL)
+error(mpl, "unable to open `%s' for writing - %s",
+fname, xerrmsg());
+mpl->prt_file = xmalloc(strlen(fname)+1);
+strcpy(mpl->prt_file, fname);
+}
+}
+loop_within_domain(mpl, prt->domain, prt, printf_func);
+if (mpl->prt_fp != NULL)
+{  xfflush(mpl->prt_fp);
+if (xferror(mpl->prt_fp))
+error(mpl, "writing error to `%s' - %s", mpl->prt_file,
+xerrmsg());
+}
+return;
+}
+#endif
+/*----------------------------------------------------------------------
+-- clean_printf - clean printf statement.
+--
+-- This routine cleans specified printf statement that assumes deleting
+-- all stuff dynamically allocated on generating/postsolving phase. */
+void clean_printf(MPL *mpl, PRINTF *prt)
+{     PRINTF1 *p;
+/* clean subscript domain */
+clean_domain(mpl, prt->domain);
+/* clean pseudo-code for computing format string */
+clean_code(mpl, prt->fmt);
+/* clean printf list */
+for (p = prt->list; p != NULL; p = p->next)
+{  /* clean pseudo-code for computing value to be printed */
+clean_code(mpl, p->code);
+}
+#if 1 /* 14/VII-2006 */
+/* clean pseudo-code for computing file name string */
+clean_code(mpl, prt->fname);
+#endif
+return;
+}
+/*----------------------------------------------------------------------
+-- execute_for - execute for statement.
+--
+-- This routine executes specified for statement. */
+static int for_func(MPL *mpl, void *info)
+{     /* this is auxiliary routine to work within domain scope */
+FOR *fur = (FOR *)info;
+STATEMENT *stmt, *save;
+save = mpl->stmt;
+for (stmt = fur->list; stmt != NULL; stmt = stmt->next)
+execute_statement(mpl, stmt);
+mpl->stmt = save;
+return 0;
+}
+void execute_for(MPL *mpl, FOR *fur)
+{     loop_within_domain(mpl, fur->domain, fur, for_func);
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_for - clean for statement.
+--
+-- This routine cleans specified for statement that assumes deleting all
+-- stuff dynamically allocated on generating/postsolving phase. */
+void clean_for(MPL *mpl, FOR *fur)
+{     STATEMENT *stmt;
+/* clean subscript domain */
+clean_domain(mpl, fur->domain);
+/* clean all sub-statements */
+for (stmt = fur->list; stmt != NULL; stmt = stmt->next)
+clean_statement(mpl, stmt);
+return;
+}
+/*----------------------------------------------------------------------
+-- execute_statement - execute specified model statement.
+--
+-- This routine executes specified model statement. */
+void execute_statement(MPL *mpl, STATEMENT *stmt)
+{     mpl->stmt = stmt;
+switch (stmt->type)
+{  case A_SET:
+case A_PARAMETER:
+case A_VARIABLE:
+break;
+case A_CONSTRAINT:
+xprintf("Generating %s...\n", stmt->u.con->name);
+eval_whole_con(mpl, stmt->u.con);
+break;
+case A_TABLE:
+switch (stmt->u.tab->type)
+{  case A_INPUT:
+xprintf("Reading %s...\n", stmt->u.tab->name);
+break;
+case A_OUTPUT:
+xprintf("Writing %s...\n", stmt->u.tab->name);
+break;
+default:
+xassert(stmt != stmt);
+}
+execute_table(mpl, stmt->u.tab);
+break;
+case A_SOLVE:
+break;
+case A_CHECK:
+xprintf("Checking (line %d)...\n", stmt->line);
+execute_check(mpl, stmt->u.chk);
+break;
+case A_DISPLAY:
+write_text(mpl, "Display statement at line %d\n",
+stmt->line);
+execute_display(mpl, stmt->u.dpy);
+break;
+case A_PRINTF:
+execute_printf(mpl, stmt->u.prt);
+break;
+case A_FOR:
+execute_for(mpl, stmt->u.fur);
+break;
+default:
+xassert(stmt != stmt);
+}
+return;
+}
+/*----------------------------------------------------------------------
+-- clean_statement - clean specified model statement.
+--
+-- This routine cleans specified model statement that assumes deleting
+-- all stuff dynamically allocated on generating/postsolving phase. */
+void clean_statement(MPL *mpl, STATEMENT *stmt)
+{     switch(stmt->type)
+{  case A_SET:
+clean_set(mpl, stmt->u.set); break;
+case A_PARAMETER:
+clean_parameter(mpl, stmt->u.par); break;
+case A_VARIABLE:
+clean_variable(mpl, stmt->u.var); break;
+case A_CONSTRAINT:
+clean_constraint(mpl, stmt->u.con); break;
+#if 1 /* 11/II-2008 */
+case A_TABLE:
+clean_table(mpl, stmt->u.tab); break;
+#endif
+case A_SOLVE:
+break;
+case A_CHECK:
+clean_check(mpl, stmt->u.chk); break;
+case A_DISPLAY:
+clean_display(mpl, stmt->u.dpy); break;
+case A_PRINTF:
+clean_printf(mpl, stmt->u.prt); break;
+case A_FOR:
+clean_for(mpl, stmt->u.fur); break;
+default:
+xassert(stmt != stmt);
+}
+return;
+}
+/* eof */