1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpmpl03.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,6078 @@
1.4 +/* glpmpl03.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#define _GLPSTD_ERRNO
1.29 +#define _GLPSTD_STDIO
1.30 +#include "glpenv.h"
1.31 +#include "glpmpl.h"
1.32 +
1.33 +/**********************************************************************/
1.34 +/* * * FLOATING-POINT NUMBERS * * */
1.35 +/**********************************************************************/
1.36 +
1.37 +/*----------------------------------------------------------------------
1.38 +-- fp_add - floating-point addition.
1.39 +--
1.40 +-- This routine computes the sum x + y. */
1.41 +
1.42 +double fp_add(MPL *mpl, double x, double y)
1.43 +{ if (x > 0.0 && y > 0.0 && x > + 0.999 * DBL_MAX - y ||
1.44 + x < 0.0 && y < 0.0 && x < - 0.999 * DBL_MAX - y)
1.45 + error(mpl, "%.*g + %.*g; floating-point overflow",
1.46 + DBL_DIG, x, DBL_DIG, y);
1.47 + return x + y;
1.48 +}
1.49 +
1.50 +/*----------------------------------------------------------------------
1.51 +-- fp_sub - floating-point subtraction.
1.52 +--
1.53 +-- This routine computes the difference x - y. */
1.54 +
1.55 +double fp_sub(MPL *mpl, double x, double y)
1.56 +{ if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y ||
1.57 + x < 0.0 && y > 0.0 && x < - 0.999 * DBL_MAX + y)
1.58 + error(mpl, "%.*g - %.*g; floating-point overflow",
1.59 + DBL_DIG, x, DBL_DIG, y);
1.60 + return x - y;
1.61 +}
1.62 +
1.63 +/*----------------------------------------------------------------------
1.64 +-- fp_less - floating-point non-negative subtraction.
1.65 +--
1.66 +-- This routine computes the non-negative difference max(0, x - y). */
1.67 +
1.68 +double fp_less(MPL *mpl, double x, double y)
1.69 +{ if (x < y) return 0.0;
1.70 + if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y)
1.71 + error(mpl, "%.*g less %.*g; floating-point overflow",
1.72 + DBL_DIG, x, DBL_DIG, y);
1.73 + return x - y;
1.74 +}
1.75 +
1.76 +/*----------------------------------------------------------------------
1.77 +-- fp_mul - floating-point multiplication.
1.78 +--
1.79 +-- This routine computes the product x * y. */
1.80 +
1.81 +double fp_mul(MPL *mpl, double x, double y)
1.82 +{ if (fabs(y) > 1.0 && fabs(x) > (0.999 * DBL_MAX) / fabs(y))
1.83 + error(mpl, "%.*g * %.*g; floating-point overflow",
1.84 + DBL_DIG, x, DBL_DIG, y);
1.85 + return x * y;
1.86 +}
1.87 +
1.88 +/*----------------------------------------------------------------------
1.89 +-- fp_div - floating-point division.
1.90 +--
1.91 +-- This routine computes the quotient x / y. */
1.92 +
1.93 +double fp_div(MPL *mpl, double x, double y)
1.94 +{ if (fabs(y) < DBL_MIN)
1.95 + error(mpl, "%.*g / %.*g; floating-point zero divide",
1.96 + DBL_DIG, x, DBL_DIG, y);
1.97 + if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))
1.98 + error(mpl, "%.*g / %.*g; floating-point overflow",
1.99 + DBL_DIG, x, DBL_DIG, y);
1.100 + return x / y;
1.101 +}
1.102 +
1.103 +/*----------------------------------------------------------------------
1.104 +-- fp_idiv - floating-point quotient of exact division.
1.105 +--
1.106 +-- This routine computes the quotient of exact division x div y. */
1.107 +
1.108 +double fp_idiv(MPL *mpl, double x, double y)
1.109 +{ if (fabs(y) < DBL_MIN)
1.110 + error(mpl, "%.*g div %.*g; floating-point zero divide",
1.111 + DBL_DIG, x, DBL_DIG, y);
1.112 + if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))
1.113 + error(mpl, "%.*g div %.*g; floating-point overflow",
1.114 + DBL_DIG, x, DBL_DIG, y);
1.115 + x /= y;
1.116 + return x > 0.0 ? floor(x) : x < 0.0 ? ceil(x) : 0.0;
1.117 +}
1.118 +
1.119 +/*----------------------------------------------------------------------
1.120 +-- fp_mod - floating-point remainder of exact division.
1.121 +--
1.122 +-- This routine computes the remainder of exact division x mod y.
1.123 +--
1.124 +-- NOTE: By definition x mod y = x - y * floor(x / y). */
1.125 +
1.126 +double fp_mod(MPL *mpl, double x, double y)
1.127 +{ double r;
1.128 + xassert(mpl == mpl);
1.129 + if (x == 0.0)
1.130 + r = 0.0;
1.131 + else if (y == 0.0)
1.132 + r = x;
1.133 + else
1.134 + { r = fmod(fabs(x), fabs(y));
1.135 + if (r != 0.0)
1.136 + { if (x < 0.0) r = - r;
1.137 + if (x > 0.0 && y < 0.0 || x < 0.0 && y > 0.0) r += y;
1.138 + }
1.139 + }
1.140 + return r;
1.141 +}
1.142 +
1.143 +/*----------------------------------------------------------------------
1.144 +-- fp_power - floating-point exponentiation (raise to power).
1.145 +--
1.146 +-- This routine computes the exponentiation x ** y. */
1.147 +
1.148 +double fp_power(MPL *mpl, double x, double y)
1.149 +{ double r;
1.150 + if (x == 0.0 && y <= 0.0 || x < 0.0 && y != floor(y))
1.151 + error(mpl, "%.*g ** %.*g; result undefined",
1.152 + DBL_DIG, x, DBL_DIG, y);
1.153 + if (x == 0.0) goto eval;
1.154 + if (fabs(x) > 1.0 && y > +1.0 &&
1.155 + +log(fabs(x)) > (0.999 * log(DBL_MAX)) / y ||
1.156 + fabs(x) < 1.0 && y < -1.0 &&
1.157 + +log(fabs(x)) < (0.999 * log(DBL_MAX)) / y)
1.158 + error(mpl, "%.*g ** %.*g; floating-point overflow",
1.159 + DBL_DIG, x, DBL_DIG, y);
1.160 + if (fabs(x) > 1.0 && y < -1.0 &&
1.161 + -log(fabs(x)) < (0.999 * log(DBL_MAX)) / y ||
1.162 + fabs(x) < 1.0 && y > +1.0 &&
1.163 + -log(fabs(x)) > (0.999 * log(DBL_MAX)) / y)
1.164 + r = 0.0;
1.165 + else
1.166 +eval: r = pow(x, y);
1.167 + return r;
1.168 +}
1.169 +
1.170 +/*----------------------------------------------------------------------
1.171 +-- fp_exp - floating-point base-e exponential.
1.172 +--
1.173 +-- This routine computes the base-e exponential e ** x. */
1.174 +
1.175 +double fp_exp(MPL *mpl, double x)
1.176 +{ if (x > 0.999 * log(DBL_MAX))
1.177 + error(mpl, "exp(%.*g); floating-point overflow", DBL_DIG, x);
1.178 + return exp(x);
1.179 +}
1.180 +
1.181 +/*----------------------------------------------------------------------
1.182 +-- fp_log - floating-point natural logarithm.
1.183 +--
1.184 +-- This routine computes the natural logarithm log x. */
1.185 +
1.186 +double fp_log(MPL *mpl, double x)
1.187 +{ if (x <= 0.0)
1.188 + error(mpl, "log(%.*g); non-positive argument", DBL_DIG, x);
1.189 + return log(x);
1.190 +}
1.191 +
1.192 +/*----------------------------------------------------------------------
1.193 +-- fp_log10 - floating-point common (decimal) logarithm.
1.194 +--
1.195 +-- This routine computes the common (decimal) logarithm lg x. */
1.196 +
1.197 +double fp_log10(MPL *mpl, double x)
1.198 +{ if (x <= 0.0)
1.199 + error(mpl, "log10(%.*g); non-positive argument", DBL_DIG, x);
1.200 + return log10(x);
1.201 +}
1.202 +
1.203 +/*----------------------------------------------------------------------
1.204 +-- fp_sqrt - floating-point square root.
1.205 +--
1.206 +-- This routine computes the square root x ** 0.5. */
1.207 +
1.208 +double fp_sqrt(MPL *mpl, double x)
1.209 +{ if (x < 0.0)
1.210 + error(mpl, "sqrt(%.*g); negative argument", DBL_DIG, x);
1.211 + return sqrt(x);
1.212 +}
1.213 +
1.214 +/*----------------------------------------------------------------------
1.215 +-- fp_sin - floating-point trigonometric sine.
1.216 +--
1.217 +-- This routine computes the trigonometric sine sin(x). */
1.218 +
1.219 +double fp_sin(MPL *mpl, double x)
1.220 +{ if (!(-1e6 <= x && x <= +1e6))
1.221 + error(mpl, "sin(%.*g); argument too large", DBL_DIG, x);
1.222 + return sin(x);
1.223 +}
1.224 +
1.225 +/*----------------------------------------------------------------------
1.226 +-- fp_cos - floating-point trigonometric cosine.
1.227 +--
1.228 +-- This routine computes the trigonometric cosine cos(x). */
1.229 +
1.230 +double fp_cos(MPL *mpl, double x)
1.231 +{ if (!(-1e6 <= x && x <= +1e6))
1.232 + error(mpl, "cos(%.*g); argument too large", DBL_DIG, x);
1.233 + return cos(x);
1.234 +}
1.235 +
1.236 +/*----------------------------------------------------------------------
1.237 +-- fp_atan - floating-point trigonometric arctangent.
1.238 +--
1.239 +-- This routine computes the trigonometric arctangent atan(x). */
1.240 +
1.241 +double fp_atan(MPL *mpl, double x)
1.242 +{ xassert(mpl == mpl);
1.243 + return atan(x);
1.244 +}
1.245 +
1.246 +/*----------------------------------------------------------------------
1.247 +-- fp_atan2 - floating-point trigonometric arctangent.
1.248 +--
1.249 +-- This routine computes the trigonometric arctangent atan(y / x). */
1.250 +
1.251 +double fp_atan2(MPL *mpl, double y, double x)
1.252 +{ xassert(mpl == mpl);
1.253 + return atan2(y, x);
1.254 +}
1.255 +
1.256 +/*----------------------------------------------------------------------
1.257 +-- fp_round - round floating-point value to n fractional digits.
1.258 +--
1.259 +-- This routine rounds given floating-point value x to n fractional
1.260 +-- digits with the formula:
1.261 +--
1.262 +-- round(x, n) = floor(x * 10^n + 0.5) / 10^n.
1.263 +--
1.264 +-- The parameter n is assumed to be integer. */
1.265 +
1.266 +double fp_round(MPL *mpl, double x, double n)
1.267 +{ double ten_to_n;
1.268 + if (n != floor(n))
1.269 + error(mpl, "round(%.*g, %.*g); non-integer second argument",
1.270 + DBL_DIG, x, DBL_DIG, n);
1.271 + if (n <= DBL_DIG + 2)
1.272 + { ten_to_n = pow(10.0, n);
1.273 + if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)
1.274 + { x = floor(x * ten_to_n + 0.5);
1.275 + if (x != 0.0) x /= ten_to_n;
1.276 + }
1.277 + }
1.278 + return x;
1.279 +}
1.280 +
1.281 +/*----------------------------------------------------------------------
1.282 +-- fp_trunc - truncate floating-point value to n fractional digits.
1.283 +--
1.284 +-- This routine truncates given floating-point value x to n fractional
1.285 +-- digits with the formula:
1.286 +--
1.287 +-- ( floor(x * 10^n) / 10^n, if x >= 0
1.288 +-- trunc(x, n) = <
1.289 +-- ( ceil(x * 10^n) / 10^n, if x < 0
1.290 +--
1.291 +-- The parameter n is assumed to be integer. */
1.292 +
1.293 +double fp_trunc(MPL *mpl, double x, double n)
1.294 +{ double ten_to_n;
1.295 + if (n != floor(n))
1.296 + error(mpl, "trunc(%.*g, %.*g); non-integer second argument",
1.297 + DBL_DIG, x, DBL_DIG, n);
1.298 + if (n <= DBL_DIG + 2)
1.299 + { ten_to_n = pow(10.0, n);
1.300 + if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)
1.301 + { x = (x >= 0.0 ? floor(x * ten_to_n) : ceil(x * ten_to_n));
1.302 + if (x != 0.0) x /= ten_to_n;
1.303 + }
1.304 + }
1.305 + return x;
1.306 +}
1.307 +
1.308 +/**********************************************************************/
1.309 +/* * * PSEUDO-RANDOM NUMBER GENERATORS * * */
1.310 +/**********************************************************************/
1.311 +
1.312 +/*----------------------------------------------------------------------
1.313 +-- fp_irand224 - pseudo-random integer in the range [0, 2^24).
1.314 +--
1.315 +-- This routine returns a next pseudo-random integer (converted to
1.316 +-- floating-point) which is uniformly distributed between 0 and 2^24-1,
1.317 +-- inclusive. */
1.318 +
1.319 +#define two_to_the_24 0x1000000
1.320 +
1.321 +double fp_irand224(MPL *mpl)
1.322 +{ return
1.323 + (double)rng_unif_rand(mpl->rand, two_to_the_24);
1.324 +}
1.325 +
1.326 +/*----------------------------------------------------------------------
1.327 +-- fp_uniform01 - pseudo-random number in the range [0, 1).
1.328 +--
1.329 +-- This routine returns a next pseudo-random number which is uniformly
1.330 +-- distributed in the range [0, 1). */
1.331 +
1.332 +#define two_to_the_31 ((unsigned int)0x80000000)
1.333 +
1.334 +double fp_uniform01(MPL *mpl)
1.335 +{ return
1.336 + (double)rng_next_rand(mpl->rand) / (double)two_to_the_31;
1.337 +}
1.338 +
1.339 +/*----------------------------------------------------------------------
1.340 +-- fp_uniform - pseudo-random number in the range [a, b).
1.341 +--
1.342 +-- This routine returns a next pseudo-random number which is uniformly
1.343 +-- distributed in the range [a, b). */
1.344 +
1.345 +double fp_uniform(MPL *mpl, double a, double b)
1.346 +{ double x;
1.347 + if (a >= b)
1.348 + error(mpl, "Uniform(%.*g, %.*g); invalid range",
1.349 + DBL_DIG, a, DBL_DIG, b);
1.350 + x = fp_uniform01(mpl);
1.351 +#if 0
1.352 + x = a * (1.0 - x) + b * x;
1.353 +#else
1.354 + x = fp_add(mpl, a * (1.0 - x), b * x);
1.355 +#endif
1.356 + return x;
1.357 +}
1.358 +
1.359 +/*----------------------------------------------------------------------
1.360 +-- fp_normal01 - Gaussian random variate with mu = 0 and sigma = 1.
1.361 +--
1.362 +-- This routine returns a Gaussian random variate with zero mean and
1.363 +-- unit standard deviation. The polar (Box-Mueller) method is used.
1.364 +--
1.365 +-- This code is a modified version of the routine gsl_ran_gaussian from
1.366 +-- the GNU Scientific Library Version 1.0. */
1.367 +
1.368 +double fp_normal01(MPL *mpl)
1.369 +{ double x, y, r2;
1.370 + do
1.371 + { /* choose x, y in uniform square (-1,-1) to (+1,+1) */
1.372 + x = -1.0 + 2.0 * fp_uniform01(mpl);
1.373 + y = -1.0 + 2.0 * fp_uniform01(mpl);
1.374 + /* see if it is in the unit circle */
1.375 + r2 = x * x + y * y;
1.376 + } while (r2 > 1.0 || r2 == 0.0);
1.377 + /* Box-Muller transform */
1.378 + return y * sqrt(-2.0 * log (r2) / r2);
1.379 +}
1.380 +
1.381 +/*----------------------------------------------------------------------
1.382 +-- fp_normal - Gaussian random variate with specified mu and sigma.
1.383 +--
1.384 +-- This routine returns a Gaussian random variate with mean mu and
1.385 +-- standard deviation sigma. */
1.386 +
1.387 +double fp_normal(MPL *mpl, double mu, double sigma)
1.388 +{ double x;
1.389 +#if 0
1.390 + x = mu + sigma * fp_normal01(mpl);
1.391 +#else
1.392 + x = fp_add(mpl, mu, fp_mul(mpl, sigma, fp_normal01(mpl)));
1.393 +#endif
1.394 + return x;
1.395 +}
1.396 +
1.397 +/**********************************************************************/
1.398 +/* * * SEGMENTED CHARACTER STRINGS * * */
1.399 +/**********************************************************************/
1.400 +
1.401 +/*----------------------------------------------------------------------
1.402 +-- create_string - create character string.
1.403 +--
1.404 +-- This routine creates a segmented character string, which is exactly
1.405 +-- equivalent to specified character string. */
1.406 +
1.407 +STRING *create_string
1.408 +( MPL *mpl,
1.409 + char buf[MAX_LENGTH+1] /* not changed */
1.410 +)
1.411 +#if 0
1.412 +{ STRING *head, *tail;
1.413 + int i, j;
1.414 + xassert(buf != NULL);
1.415 + xassert(strlen(buf) <= MAX_LENGTH);
1.416 + head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));
1.417 + for (i = j = 0; ; i++)
1.418 + { if ((tail->seg[j++] = buf[i]) == '\0') break;
1.419 + if (j == STRSEG_SIZE)
1.420 +tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING))), j = 0;
1.421 + }
1.422 + tail->next = NULL;
1.423 + return head;
1.424 +}
1.425 +#else
1.426 +{ STRING *str;
1.427 + xassert(strlen(buf) <= MAX_LENGTH);
1.428 + str = dmp_get_atom(mpl->strings, strlen(buf)+1);
1.429 + strcpy(str, buf);
1.430 + return str;
1.431 +}
1.432 +#endif
1.433 +
1.434 +/*----------------------------------------------------------------------
1.435 +-- copy_string - make copy of character string.
1.436 +--
1.437 +-- This routine returns an exact copy of segmented character string. */
1.438 +
1.439 +STRING *copy_string
1.440 +( MPL *mpl,
1.441 + STRING *str /* not changed */
1.442 +)
1.443 +#if 0
1.444 +{ STRING *head, *tail;
1.445 + xassert(str != NULL);
1.446 + head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));
1.447 + for (; str != NULL; str = str->next)
1.448 + { memcpy(tail->seg, str->seg, STRSEG_SIZE);
1.449 + if (str->next != NULL)
1.450 +tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING)));
1.451 + }
1.452 + tail->next = NULL;
1.453 + return head;
1.454 +}
1.455 +#else
1.456 +{ xassert(mpl == mpl);
1.457 + return create_string(mpl, str);
1.458 +}
1.459 +#endif
1.460 +
1.461 +/*----------------------------------------------------------------------
1.462 +-- compare_strings - compare one character string with another.
1.463 +--
1.464 +-- This routine compares one segmented character strings with another
1.465 +-- and returns the result of comparison as follows:
1.466 +--
1.467 +-- = 0 - both strings are identical;
1.468 +-- < 0 - the first string precedes the second one;
1.469 +-- > 0 - the first string follows the second one. */
1.470 +
1.471 +int compare_strings
1.472 +( MPL *mpl,
1.473 + STRING *str1, /* not changed */
1.474 + STRING *str2 /* not changed */
1.475 +)
1.476 +#if 0
1.477 +{ int j, c1, c2;
1.478 + xassert(mpl == mpl);
1.479 + for (;; str1 = str1->next, str2 = str2->next)
1.480 + { xassert(str1 != NULL);
1.481 + xassert(str2 != NULL);
1.482 + for (j = 0; j < STRSEG_SIZE; j++)
1.483 + { c1 = (unsigned char)str1->seg[j];
1.484 + c2 = (unsigned char)str2->seg[j];
1.485 + if (c1 < c2) return -1;
1.486 + if (c1 > c2) return +1;
1.487 + if (c1 == '\0') goto done;
1.488 + }
1.489 + }
1.490 +done: return 0;
1.491 +}
1.492 +#else
1.493 +{ xassert(mpl == mpl);
1.494 + return strcmp(str1, str2);
1.495 +}
1.496 +#endif
1.497 +
1.498 +/*----------------------------------------------------------------------
1.499 +-- fetch_string - extract content of character string.
1.500 +--
1.501 +-- This routine returns a character string, which is exactly equivalent
1.502 +-- to specified segmented character string. */
1.503 +
1.504 +char *fetch_string
1.505 +( MPL *mpl,
1.506 + STRING *str, /* not changed */
1.507 + char buf[MAX_LENGTH+1] /* modified */
1.508 +)
1.509 +#if 0
1.510 +{ int i, j;
1.511 + xassert(mpl == mpl);
1.512 + xassert(buf != NULL);
1.513 + for (i = 0; ; str = str->next)
1.514 + { xassert(str != NULL);
1.515 + for (j = 0; j < STRSEG_SIZE; j++)
1.516 + if ((buf[i++] = str->seg[j]) == '\0') goto done;
1.517 + }
1.518 +done: xassert(strlen(buf) <= MAX_LENGTH);
1.519 + return buf;
1.520 +}
1.521 +#else
1.522 +{ xassert(mpl == mpl);
1.523 + return strcpy(buf, str);
1.524 +}
1.525 +#endif
1.526 +
1.527 +/*----------------------------------------------------------------------
1.528 +-- delete_string - delete character string.
1.529 +--
1.530 +-- This routine deletes specified segmented character string. */
1.531 +
1.532 +void delete_string
1.533 +( MPL *mpl,
1.534 + STRING *str /* destroyed */
1.535 +)
1.536 +#if 0
1.537 +{ STRING *temp;
1.538 + xassert(str != NULL);
1.539 + while (str != NULL)
1.540 + { temp = str;
1.541 + str = str->next;
1.542 + dmp_free_atom(mpl->strings, temp, sizeof(STRING));
1.543 + }
1.544 + return;
1.545 +}
1.546 +#else
1.547 +{ dmp_free_atom(mpl->strings, str, strlen(str)+1);
1.548 + return;
1.549 +}
1.550 +#endif
1.551 +
1.552 +/**********************************************************************/
1.553 +/* * * SYMBOLS * * */
1.554 +/**********************************************************************/
1.555 +
1.556 +/*----------------------------------------------------------------------
1.557 +-- create_symbol_num - create symbol of numeric type.
1.558 +--
1.559 +-- This routine creates a symbol, which has a numeric value specified
1.560 +-- as floating-point number. */
1.561 +
1.562 +SYMBOL *create_symbol_num(MPL *mpl, double num)
1.563 +{ SYMBOL *sym;
1.564 + sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));
1.565 + sym->num = num;
1.566 + sym->str = NULL;
1.567 + return sym;
1.568 +}
1.569 +
1.570 +/*----------------------------------------------------------------------
1.571 +-- create_symbol_str - create symbol of abstract type.
1.572 +--
1.573 +-- This routine creates a symbol, which has an abstract value specified
1.574 +-- as segmented character string. */
1.575 +
1.576 +SYMBOL *create_symbol_str
1.577 +( MPL *mpl,
1.578 + STRING *str /* destroyed */
1.579 +)
1.580 +{ SYMBOL *sym;
1.581 + xassert(str != NULL);
1.582 + sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));
1.583 + sym->num = 0.0;
1.584 + sym->str = str;
1.585 + return sym;
1.586 +}
1.587 +
1.588 +/*----------------------------------------------------------------------
1.589 +-- copy_symbol - make copy of symbol.
1.590 +--
1.591 +-- This routine returns an exact copy of symbol. */
1.592 +
1.593 +SYMBOL *copy_symbol
1.594 +( MPL *mpl,
1.595 + SYMBOL *sym /* not changed */
1.596 +)
1.597 +{ SYMBOL *copy;
1.598 + xassert(sym != NULL);
1.599 + copy = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));
1.600 + if (sym->str == NULL)
1.601 + { copy->num = sym->num;
1.602 + copy->str = NULL;
1.603 + }
1.604 + else
1.605 + { copy->num = 0.0;
1.606 + copy->str = copy_string(mpl, sym->str);
1.607 + }
1.608 + return copy;
1.609 +}
1.610 +
1.611 +/*----------------------------------------------------------------------
1.612 +-- compare_symbols - compare one symbol with another.
1.613 +--
1.614 +-- This routine compares one symbol with another and returns the result
1.615 +-- of comparison as follows:
1.616 +--
1.617 +-- = 0 - both symbols are identical;
1.618 +-- < 0 - the first symbol precedes the second one;
1.619 +-- > 0 - the first symbol follows the second one.
1.620 +--
1.621 +-- Note that the linear order, in which symbols follow each other, is
1.622 +-- implementation-dependent. It may be not an alphabetical order. */
1.623 +
1.624 +int compare_symbols
1.625 +( MPL *mpl,
1.626 + SYMBOL *sym1, /* not changed */
1.627 + SYMBOL *sym2 /* not changed */
1.628 +)
1.629 +{ xassert(sym1 != NULL);
1.630 + xassert(sym2 != NULL);
1.631 + /* let all numeric quantities precede all symbolic quantities */
1.632 + if (sym1->str == NULL && sym2->str == NULL)
1.633 + { if (sym1->num < sym2->num) return -1;
1.634 + if (sym1->num > sym2->num) return +1;
1.635 + return 0;
1.636 + }
1.637 + if (sym1->str == NULL) return -1;
1.638 + if (sym2->str == NULL) return +1;
1.639 + return compare_strings(mpl, sym1->str, sym2->str);
1.640 +}
1.641 +
1.642 +/*----------------------------------------------------------------------
1.643 +-- delete_symbol - delete symbol.
1.644 +--
1.645 +-- This routine deletes specified symbol. */
1.646 +
1.647 +void delete_symbol
1.648 +( MPL *mpl,
1.649 + SYMBOL *sym /* destroyed */
1.650 +)
1.651 +{ xassert(sym != NULL);
1.652 + if (sym->str != NULL) delete_string(mpl, sym->str);
1.653 + dmp_free_atom(mpl->symbols, sym, sizeof(SYMBOL));
1.654 + return;
1.655 +}
1.656 +
1.657 +/*----------------------------------------------------------------------
1.658 +-- format_symbol - format symbol for displaying or printing.
1.659 +--
1.660 +-- This routine converts specified symbol to a charater string, which
1.661 +-- is suitable for displaying or printing.
1.662 +--
1.663 +-- The resultant string is never longer than 255 characters. If it gets
1.664 +-- longer, it is truncated from the right and appended by dots. */
1.665 +
1.666 +char *format_symbol
1.667 +( MPL *mpl,
1.668 + SYMBOL *sym /* not changed */
1.669 +)
1.670 +{ char *buf = mpl->sym_buf;
1.671 + xassert(sym != NULL);
1.672 + if (sym->str == NULL)
1.673 + sprintf(buf, "%.*g", DBL_DIG, sym->num);
1.674 + else
1.675 + { char str[MAX_LENGTH+1];
1.676 + int quoted, j, len;
1.677 + fetch_string(mpl, sym->str, str);
1.678 + if (!(isalpha((unsigned char)str[0]) || str[0] == '_'))
1.679 + quoted = 1;
1.680 + else
1.681 + { quoted = 0;
1.682 + for (j = 1; str[j] != '\0'; j++)
1.683 + { if (!(isalnum((unsigned char)str[j]) ||
1.684 + strchr("+-._", (unsigned char)str[j]) != NULL))
1.685 + { quoted = 1;
1.686 + break;
1.687 + }
1.688 + }
1.689 + }
1.690 +# define safe_append(c) \
1.691 + (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)
1.692 + buf[0] = '\0', len = 0;
1.693 + if (quoted) safe_append('\'');
1.694 + for (j = 0; str[j] != '\0'; j++)
1.695 + { if (quoted && str[j] == '\'') safe_append('\'');
1.696 + safe_append(str[j]);
1.697 + }
1.698 + if (quoted) safe_append('\'');
1.699 +# undef safe_append
1.700 + buf[len] = '\0';
1.701 + if (len == 255) strcpy(buf+252, "...");
1.702 + }
1.703 + xassert(strlen(buf) <= 255);
1.704 + return buf;
1.705 +}
1.706 +
1.707 +/*----------------------------------------------------------------------
1.708 +-- concat_symbols - concatenate one symbol with another.
1.709 +--
1.710 +-- This routine concatenates values of two given symbols and assigns
1.711 +-- the resultant character string to a new symbol, which is returned on
1.712 +-- exit. Both original symbols are destroyed. */
1.713 +
1.714 +SYMBOL *concat_symbols
1.715 +( MPL *mpl,
1.716 + SYMBOL *sym1, /* destroyed */
1.717 + SYMBOL *sym2 /* destroyed */
1.718 +)
1.719 +{ char str1[MAX_LENGTH+1], str2[MAX_LENGTH+1];
1.720 + xassert(MAX_LENGTH >= DBL_DIG + DBL_DIG);
1.721 + if (sym1->str == NULL)
1.722 + sprintf(str1, "%.*g", DBL_DIG, sym1->num);
1.723 + else
1.724 + fetch_string(mpl, sym1->str, str1);
1.725 + if (sym2->str == NULL)
1.726 + sprintf(str2, "%.*g", DBL_DIG, sym2->num);
1.727 + else
1.728 + fetch_string(mpl, sym2->str, str2);
1.729 + if (strlen(str1) + strlen(str2) > MAX_LENGTH)
1.730 + { char buf[255+1];
1.731 + strcpy(buf, format_symbol(mpl, sym1));
1.732 + xassert(strlen(buf) < sizeof(buf));
1.733 + error(mpl, "%s & %s; resultant symbol exceeds %d characters",
1.734 + buf, format_symbol(mpl, sym2), MAX_LENGTH);
1.735 + }
1.736 + delete_symbol(mpl, sym1);
1.737 + delete_symbol(mpl, sym2);
1.738 + return create_symbol_str(mpl, create_string(mpl, strcat(str1,
1.739 + str2)));
1.740 +}
1.741 +
1.742 +/**********************************************************************/
1.743 +/* * * N-TUPLES * * */
1.744 +/**********************************************************************/
1.745 +
1.746 +/*----------------------------------------------------------------------
1.747 +-- create_tuple - create n-tuple.
1.748 +--
1.749 +-- This routine creates a n-tuple, which initially has no components,
1.750 +-- i.e. which is 0-tuple. */
1.751 +
1.752 +TUPLE *create_tuple(MPL *mpl)
1.753 +{ TUPLE *tuple;
1.754 + xassert(mpl == mpl);
1.755 + tuple = NULL;
1.756 + return tuple;
1.757 +}
1.758 +
1.759 +/*----------------------------------------------------------------------
1.760 +-- expand_tuple - append symbol to n-tuple.
1.761 +--
1.762 +-- This routine expands n-tuple appending to it a given symbol, which
1.763 +-- becomes its new last component. */
1.764 +
1.765 +TUPLE *expand_tuple
1.766 +( MPL *mpl,
1.767 + TUPLE *tuple, /* destroyed */
1.768 + SYMBOL *sym /* destroyed */
1.769 +)
1.770 +{ TUPLE *tail, *temp;
1.771 + xassert(sym != NULL);
1.772 + /* create a new component */
1.773 + tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
1.774 + tail->sym = sym;
1.775 + tail->next = NULL;
1.776 + /* and append it to the component list */
1.777 + if (tuple == NULL)
1.778 + tuple = tail;
1.779 + else
1.780 + { for (temp = tuple; temp->next != NULL; temp = temp->next);
1.781 + temp->next = tail;
1.782 + }
1.783 + return tuple;
1.784 +}
1.785 +
1.786 +/*----------------------------------------------------------------------
1.787 +-- tuple_dimen - determine dimension of n-tuple.
1.788 +--
1.789 +-- This routine returns dimension of n-tuple, i.e. number of components
1.790 +-- in the n-tuple. */
1.791 +
1.792 +int tuple_dimen
1.793 +( MPL *mpl,
1.794 + TUPLE *tuple /* not changed */
1.795 +)
1.796 +{ TUPLE *temp;
1.797 + int dim = 0;
1.798 + xassert(mpl == mpl);
1.799 + for (temp = tuple; temp != NULL; temp = temp->next) dim++;
1.800 + return dim;
1.801 +}
1.802 +
1.803 +/*----------------------------------------------------------------------
1.804 +-- copy_tuple - make copy of n-tuple.
1.805 +--
1.806 +-- This routine returns an exact copy of n-tuple. */
1.807 +
1.808 +TUPLE *copy_tuple
1.809 +( MPL *mpl,
1.810 + TUPLE *tuple /* not changed */
1.811 +)
1.812 +{ TUPLE *head, *tail;
1.813 + if (tuple == NULL)
1.814 + head = NULL;
1.815 + else
1.816 + { head = tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
1.817 + for (; tuple != NULL; tuple = tuple->next)
1.818 + { xassert(tuple->sym != NULL);
1.819 + tail->sym = copy_symbol(mpl, tuple->sym);
1.820 + if (tuple->next != NULL)
1.821 +tail = (tail->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));
1.822 + }
1.823 + tail->next = NULL;
1.824 + }
1.825 + return head;
1.826 +}
1.827 +
1.828 +/*----------------------------------------------------------------------
1.829 +-- compare_tuples - compare one n-tuple with another.
1.830 +--
1.831 +-- This routine compares two given n-tuples, which must have the same
1.832 +-- dimension (not checked for the sake of efficiency), and returns one
1.833 +-- of the following codes:
1.834 +--
1.835 +-- = 0 - both n-tuples are identical;
1.836 +-- < 0 - the first n-tuple precedes the second one;
1.837 +-- > 0 - the first n-tuple follows the second one.
1.838 +--
1.839 +-- Note that the linear order, in which n-tuples follow each other, is
1.840 +-- implementation-dependent. It may be not an alphabetical order. */
1.841 +
1.842 +int compare_tuples
1.843 +( MPL *mpl,
1.844 + TUPLE *tuple1, /* not changed */
1.845 + TUPLE *tuple2 /* not changed */
1.846 +)
1.847 +{ TUPLE *item1, *item2;
1.848 + int ret;
1.849 + xassert(mpl == mpl);
1.850 + for (item1 = tuple1, item2 = tuple2; item1 != NULL;
1.851 + item1 = item1->next, item2 = item2->next)
1.852 + { xassert(item2 != NULL);
1.853 + xassert(item1->sym != NULL);
1.854 + xassert(item2->sym != NULL);
1.855 + ret = compare_symbols(mpl, item1->sym, item2->sym);
1.856 + if (ret != 0) return ret;
1.857 + }
1.858 + xassert(item2 == NULL);
1.859 + return 0;
1.860 +}
1.861 +
1.862 +/*----------------------------------------------------------------------
1.863 +-- build_subtuple - build subtuple of given n-tuple.
1.864 +--
1.865 +-- This routine builds subtuple, which consists of first dim components
1.866 +-- of given n-tuple. */
1.867 +
1.868 +TUPLE *build_subtuple
1.869 +( MPL *mpl,
1.870 + TUPLE *tuple, /* not changed */
1.871 + int dim
1.872 +)
1.873 +{ TUPLE *head, *temp;
1.874 + int j;
1.875 + head = create_tuple(mpl);
1.876 + for (j = 1, temp = tuple; j <= dim; j++, temp = temp->next)
1.877 + { xassert(temp != NULL);
1.878 + head = expand_tuple(mpl, head, copy_symbol(mpl, temp->sym));
1.879 + }
1.880 + return head;
1.881 +}
1.882 +
1.883 +/*----------------------------------------------------------------------
1.884 +-- delete_tuple - delete n-tuple.
1.885 +--
1.886 +-- This routine deletes specified n-tuple. */
1.887 +
1.888 +void delete_tuple
1.889 +( MPL *mpl,
1.890 + TUPLE *tuple /* destroyed */
1.891 +)
1.892 +{ TUPLE *temp;
1.893 + while (tuple != NULL)
1.894 + { temp = tuple;
1.895 + tuple = temp->next;
1.896 + xassert(temp->sym != NULL);
1.897 + delete_symbol(mpl, temp->sym);
1.898 + dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));
1.899 + }
1.900 + return;
1.901 +}
1.902 +
1.903 +/*----------------------------------------------------------------------
1.904 +-- format_tuple - format n-tuple for displaying or printing.
1.905 +--
1.906 +-- This routine converts specified n-tuple to a character string, which
1.907 +-- is suitable for displaying or printing.
1.908 +--
1.909 +-- The resultant string is never longer than 255 characters. If it gets
1.910 +-- longer, it is truncated from the right and appended by dots. */
1.911 +
1.912 +char *format_tuple
1.913 +( MPL *mpl,
1.914 + int c,
1.915 + TUPLE *tuple /* not changed */
1.916 +)
1.917 +{ TUPLE *temp;
1.918 + int dim, j, len;
1.919 + char *buf = mpl->tup_buf, str[255+1], *save;
1.920 +# define safe_append(c) \
1.921 + (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)
1.922 + buf[0] = '\0', len = 0;
1.923 + dim = tuple_dimen(mpl, tuple);
1.924 + if (c == '[' && dim > 0) safe_append('[');
1.925 + if (c == '(' && dim > 1) safe_append('(');
1.926 + for (temp = tuple; temp != NULL; temp = temp->next)
1.927 + { if (temp != tuple) safe_append(',');
1.928 + xassert(temp->sym != NULL);
1.929 + save = mpl->sym_buf;
1.930 + mpl->sym_buf = str;
1.931 + format_symbol(mpl, temp->sym);
1.932 + mpl->sym_buf = save;
1.933 + xassert(strlen(str) < sizeof(str));
1.934 + for (j = 0; str[j] != '\0'; j++) safe_append(str[j]);
1.935 + }
1.936 + if (c == '[' && dim > 0) safe_append(']');
1.937 + if (c == '(' && dim > 1) safe_append(')');
1.938 +# undef safe_append
1.939 + buf[len] = '\0';
1.940 + if (len == 255) strcpy(buf+252, "...");
1.941 + xassert(strlen(buf) <= 255);
1.942 + return buf;
1.943 +}
1.944 +
1.945 +/**********************************************************************/
1.946 +/* * * ELEMENTAL SETS * * */
1.947 +/**********************************************************************/
1.948 +
1.949 +/*----------------------------------------------------------------------
1.950 +-- create_elemset - create elemental set.
1.951 +--
1.952 +-- This routine creates an elemental set, whose members are n-tuples of
1.953 +-- specified dimension. Being created the set is initially empty. */
1.954 +
1.955 +ELEMSET *create_elemset(MPL *mpl, int dim)
1.956 +{ ELEMSET *set;
1.957 + xassert(dim > 0);
1.958 + set = create_array(mpl, A_NONE, dim);
1.959 + return set;
1.960 +}
1.961 +
1.962 +/*----------------------------------------------------------------------
1.963 +-- find_tuple - check if elemental set contains given n-tuple.
1.964 +--
1.965 +-- This routine finds given n-tuple in specified elemental set in order
1.966 +-- to check if the set contains that n-tuple. If the n-tuple is found,
1.967 +-- the routine returns pointer to corresponding array member. Otherwise
1.968 +-- null pointer is returned. */
1.969 +
1.970 +MEMBER *find_tuple
1.971 +( MPL *mpl,
1.972 + ELEMSET *set, /* not changed */
1.973 + TUPLE *tuple /* not changed */
1.974 +)
1.975 +{ xassert(set != NULL);
1.976 + xassert(set->type == A_NONE);
1.977 + xassert(set->dim == tuple_dimen(mpl, tuple));
1.978 + return find_member(mpl, set, tuple);
1.979 +}
1.980 +
1.981 +/*----------------------------------------------------------------------
1.982 +-- add_tuple - add new n-tuple to elemental set.
1.983 +--
1.984 +-- This routine adds given n-tuple to specified elemental set.
1.985 +--
1.986 +-- For the sake of efficiency this routine doesn't check whether the
1.987 +-- set already contains the same n-tuple or not. Therefore the calling
1.988 +-- program should use the routine find_tuple (if necessary) in order to
1.989 +-- make sure that the given n-tuple is not contained in the set, since
1.990 +-- duplicate n-tuples within the same set are not allowed. */
1.991 +
1.992 +MEMBER *add_tuple
1.993 +( MPL *mpl,
1.994 + ELEMSET *set, /* modified */
1.995 + TUPLE *tuple /* destroyed */
1.996 +)
1.997 +{ MEMBER *memb;
1.998 + xassert(set != NULL);
1.999 + xassert(set->type == A_NONE);
1.1000 + xassert(set->dim == tuple_dimen(mpl, tuple));
1.1001 + memb = add_member(mpl, set, tuple);
1.1002 + memb->value.none = NULL;
1.1003 + return memb;
1.1004 +}
1.1005 +
1.1006 +/*----------------------------------------------------------------------
1.1007 +-- check_then_add - check and add new n-tuple to elemental set.
1.1008 +--
1.1009 +-- This routine is equivalent to the routine add_tuple except that it
1.1010 +-- does check for duplicate n-tuples. */
1.1011 +
1.1012 +MEMBER *check_then_add
1.1013 +( MPL *mpl,
1.1014 + ELEMSET *set, /* modified */
1.1015 + TUPLE *tuple /* destroyed */
1.1016 +)
1.1017 +{ if (find_tuple(mpl, set, tuple) != NULL)
1.1018 + error(mpl, "duplicate tuple %s detected", format_tuple(mpl,
1.1019 + '(', tuple));
1.1020 + return add_tuple(mpl, set, tuple);
1.1021 +}
1.1022 +
1.1023 +/*----------------------------------------------------------------------
1.1024 +-- copy_elemset - make copy of elemental set.
1.1025 +--
1.1026 +-- This routine makes an exact copy of elemental set. */
1.1027 +
1.1028 +ELEMSET *copy_elemset
1.1029 +( MPL *mpl,
1.1030 + ELEMSET *set /* not changed */
1.1031 +)
1.1032 +{ ELEMSET *copy;
1.1033 + MEMBER *memb;
1.1034 + xassert(set != NULL);
1.1035 + xassert(set->type == A_NONE);
1.1036 + xassert(set->dim > 0);
1.1037 + copy = create_elemset(mpl, set->dim);
1.1038 + for (memb = set->head; memb != NULL; memb = memb->next)
1.1039 + add_tuple(mpl, copy, copy_tuple(mpl, memb->tuple));
1.1040 + return copy;
1.1041 +}
1.1042 +
1.1043 +/*----------------------------------------------------------------------
1.1044 +-- delete_elemset - delete elemental set.
1.1045 +--
1.1046 +-- This routine deletes specified elemental set. */
1.1047 +
1.1048 +void delete_elemset
1.1049 +( MPL *mpl,
1.1050 + ELEMSET *set /* destroyed */
1.1051 +)
1.1052 +{ xassert(set != NULL);
1.1053 + xassert(set->type == A_NONE);
1.1054 + delete_array(mpl, set);
1.1055 + return;
1.1056 +}
1.1057 +
1.1058 +/*----------------------------------------------------------------------
1.1059 +-- arelset_size - compute size of "arithmetic" elemental set.
1.1060 +--
1.1061 +-- This routine computes the size of "arithmetic" elemental set, which
1.1062 +-- is specified in the form of arithmetic progression:
1.1063 +--
1.1064 +-- { t0 .. tf by dt }.
1.1065 +--
1.1066 +-- The size is computed using the formula:
1.1067 +--
1.1068 +-- n = max(0, floor((tf - t0) / dt) + 1). */
1.1069 +
1.1070 +int arelset_size(MPL *mpl, double t0, double tf, double dt)
1.1071 +{ double temp;
1.1072 + if (dt == 0.0)
1.1073 + error(mpl, "%.*g .. %.*g by %.*g; zero stride not allowed",
1.1074 + DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);
1.1075 + if (tf > 0.0 && t0 < 0.0 && tf > + 0.999 * DBL_MAX + t0)
1.1076 + temp = +DBL_MAX;
1.1077 + else if (tf < 0.0 && t0 > 0.0 && tf < - 0.999 * DBL_MAX + t0)
1.1078 + temp = -DBL_MAX;
1.1079 + else
1.1080 + temp = tf - t0;
1.1081 + if (fabs(dt) < 1.0 && fabs(temp) > (0.999 * DBL_MAX) * fabs(dt))
1.1082 + { if (temp > 0.0 && dt > 0.0 || temp < 0.0 && dt < 0.0)
1.1083 + temp = +DBL_MAX;
1.1084 + else
1.1085 + temp = 0.0;
1.1086 + }
1.1087 + else
1.1088 + { temp = floor(temp / dt) + 1.0;
1.1089 + if (temp < 0.0) temp = 0.0;
1.1090 + }
1.1091 + xassert(temp >= 0.0);
1.1092 + if (temp > (double)(INT_MAX - 1))
1.1093 + error(mpl, "%.*g .. %.*g by %.*g; set too large",
1.1094 + DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);
1.1095 + return (int)(temp + 0.5);
1.1096 +}
1.1097 +
1.1098 +/*----------------------------------------------------------------------
1.1099 +-- arelset_member - compute member of "arithmetic" elemental set.
1.1100 +--
1.1101 +-- This routine returns a numeric value of symbol, which is equivalent
1.1102 +-- to j-th member of given "arithmetic" elemental set specified in the
1.1103 +-- form of arithmetic progression:
1.1104 +--
1.1105 +-- { t0 .. tf by dt }.
1.1106 +--
1.1107 +-- The symbol value is computed with the formula:
1.1108 +--
1.1109 +-- j-th member = t0 + (j - 1) * dt,
1.1110 +--
1.1111 +-- The number j must satisfy to the restriction 1 <= j <= n, where n is
1.1112 +-- the set size computed by the routine arelset_size. */
1.1113 +
1.1114 +double arelset_member(MPL *mpl, double t0, double tf, double dt, int j)
1.1115 +{ xassert(1 <= j && j <= arelset_size(mpl, t0, tf, dt));
1.1116 + return t0 + (double)(j - 1) * dt;
1.1117 +}
1.1118 +
1.1119 +/*----------------------------------------------------------------------
1.1120 +-- create_arelset - create "arithmetic" elemental set.
1.1121 +--
1.1122 +-- This routine creates "arithmetic" elemental set, which is specified
1.1123 +-- in the form of arithmetic progression:
1.1124 +--
1.1125 +-- { t0 .. tf by dt }.
1.1126 +--
1.1127 +-- Components of this set are 1-tuples. */
1.1128 +
1.1129 +ELEMSET *create_arelset(MPL *mpl, double t0, double tf, double dt)
1.1130 +{ ELEMSET *set;
1.1131 + int j, n;
1.1132 + set = create_elemset(mpl, 1);
1.1133 + n = arelset_size(mpl, t0, tf, dt);
1.1134 + for (j = 1; j <= n; j++)
1.1135 + { add_tuple
1.1136 + ( mpl,
1.1137 + set,
1.1138 + expand_tuple
1.1139 + ( mpl,
1.1140 + create_tuple(mpl),
1.1141 + create_symbol_num
1.1142 + ( mpl,
1.1143 + arelset_member(mpl, t0, tf, dt, j)
1.1144 + )
1.1145 + )
1.1146 + );
1.1147 + }
1.1148 + return set;
1.1149 +}
1.1150 +
1.1151 +/*----------------------------------------------------------------------
1.1152 +-- set_union - union of two elemental sets.
1.1153 +--
1.1154 +-- This routine computes the union:
1.1155 +--
1.1156 +-- X U Y = { j | (j in X) or (j in Y) },
1.1157 +--
1.1158 +-- where X and Y are given elemental sets (destroyed on exit). */
1.1159 +
1.1160 +ELEMSET *set_union
1.1161 +( MPL *mpl,
1.1162 + ELEMSET *X, /* destroyed */
1.1163 + ELEMSET *Y /* destroyed */
1.1164 +)
1.1165 +{ MEMBER *memb;
1.1166 + xassert(X != NULL);
1.1167 + xassert(X->type == A_NONE);
1.1168 + xassert(X->dim > 0);
1.1169 + xassert(Y != NULL);
1.1170 + xassert(Y->type == A_NONE);
1.1171 + xassert(Y->dim > 0);
1.1172 + xassert(X->dim == Y->dim);
1.1173 + for (memb = Y->head; memb != NULL; memb = memb->next)
1.1174 + { if (find_tuple(mpl, X, memb->tuple) == NULL)
1.1175 + add_tuple(mpl, X, copy_tuple(mpl, memb->tuple));
1.1176 + }
1.1177 + delete_elemset(mpl, Y);
1.1178 + return X;
1.1179 +}
1.1180 +
1.1181 +/*----------------------------------------------------------------------
1.1182 +-- set_diff - difference between two elemental sets.
1.1183 +--
1.1184 +-- This routine computes the difference:
1.1185 +--
1.1186 +-- X \ Y = { j | (j in X) and (j not in Y) },
1.1187 +--
1.1188 +-- where X and Y are given elemental sets (destroyed on exit). */
1.1189 +
1.1190 +ELEMSET *set_diff
1.1191 +( MPL *mpl,
1.1192 + ELEMSET *X, /* destroyed */
1.1193 + ELEMSET *Y /* destroyed */
1.1194 +)
1.1195 +{ ELEMSET *Z;
1.1196 + MEMBER *memb;
1.1197 + xassert(X != NULL);
1.1198 + xassert(X->type == A_NONE);
1.1199 + xassert(X->dim > 0);
1.1200 + xassert(Y != NULL);
1.1201 + xassert(Y->type == A_NONE);
1.1202 + xassert(Y->dim > 0);
1.1203 + xassert(X->dim == Y->dim);
1.1204 + Z = create_elemset(mpl, X->dim);
1.1205 + for (memb = X->head; memb != NULL; memb = memb->next)
1.1206 + { if (find_tuple(mpl, Y, memb->tuple) == NULL)
1.1207 + add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1.1208 + }
1.1209 + delete_elemset(mpl, X);
1.1210 + delete_elemset(mpl, Y);
1.1211 + return Z;
1.1212 +}
1.1213 +
1.1214 +/*----------------------------------------------------------------------
1.1215 +-- set_symdiff - symmetric difference between two elemental sets.
1.1216 +--
1.1217 +-- This routine computes the symmetric difference:
1.1218 +--
1.1219 +-- X (+) Y = (X \ Y) U (Y \ X),
1.1220 +--
1.1221 +-- where X and Y are given elemental sets (destroyed on exit). */
1.1222 +
1.1223 +ELEMSET *set_symdiff
1.1224 +( MPL *mpl,
1.1225 + ELEMSET *X, /* destroyed */
1.1226 + ELEMSET *Y /* destroyed */
1.1227 +)
1.1228 +{ ELEMSET *Z;
1.1229 + MEMBER *memb;
1.1230 + xassert(X != NULL);
1.1231 + xassert(X->type == A_NONE);
1.1232 + xassert(X->dim > 0);
1.1233 + xassert(Y != NULL);
1.1234 + xassert(Y->type == A_NONE);
1.1235 + xassert(Y->dim > 0);
1.1236 + xassert(X->dim == Y->dim);
1.1237 + /* Z := X \ Y */
1.1238 + Z = create_elemset(mpl, X->dim);
1.1239 + for (memb = X->head; memb != NULL; memb = memb->next)
1.1240 + { if (find_tuple(mpl, Y, memb->tuple) == NULL)
1.1241 + add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1.1242 + }
1.1243 + /* Z := Z U (Y \ X) */
1.1244 + for (memb = Y->head; memb != NULL; memb = memb->next)
1.1245 + { if (find_tuple(mpl, X, memb->tuple) == NULL)
1.1246 + add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1.1247 + }
1.1248 + delete_elemset(mpl, X);
1.1249 + delete_elemset(mpl, Y);
1.1250 + return Z;
1.1251 +}
1.1252 +
1.1253 +/*----------------------------------------------------------------------
1.1254 +-- set_inter - intersection of two elemental sets.
1.1255 +--
1.1256 +-- This routine computes the intersection:
1.1257 +--
1.1258 +-- X ^ Y = { j | (j in X) and (j in Y) },
1.1259 +--
1.1260 +-- where X and Y are given elemental sets (destroyed on exit). */
1.1261 +
1.1262 +ELEMSET *set_inter
1.1263 +( MPL *mpl,
1.1264 + ELEMSET *X, /* destroyed */
1.1265 + ELEMSET *Y /* destroyed */
1.1266 +)
1.1267 +{ ELEMSET *Z;
1.1268 + MEMBER *memb;
1.1269 + xassert(X != NULL);
1.1270 + xassert(X->type == A_NONE);
1.1271 + xassert(X->dim > 0);
1.1272 + xassert(Y != NULL);
1.1273 + xassert(Y->type == A_NONE);
1.1274 + xassert(Y->dim > 0);
1.1275 + xassert(X->dim == Y->dim);
1.1276 + Z = create_elemset(mpl, X->dim);
1.1277 + for (memb = X->head; memb != NULL; memb = memb->next)
1.1278 + { if (find_tuple(mpl, Y, memb->tuple) != NULL)
1.1279 + add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1.1280 + }
1.1281 + delete_elemset(mpl, X);
1.1282 + delete_elemset(mpl, Y);
1.1283 + return Z;
1.1284 +}
1.1285 +
1.1286 +/*----------------------------------------------------------------------
1.1287 +-- set_cross - cross (Cartesian) product of two elemental sets.
1.1288 +--
1.1289 +-- This routine computes the cross (Cartesian) product:
1.1290 +--
1.1291 +-- X x Y = { (i,j) | (i in X) and (j in Y) },
1.1292 +--
1.1293 +-- where X and Y are given elemental sets (destroyed on exit). */
1.1294 +
1.1295 +ELEMSET *set_cross
1.1296 +( MPL *mpl,
1.1297 + ELEMSET *X, /* destroyed */
1.1298 + ELEMSET *Y /* destroyed */
1.1299 +)
1.1300 +{ ELEMSET *Z;
1.1301 + MEMBER *memx, *memy;
1.1302 + TUPLE *tuple, *temp;
1.1303 + xassert(X != NULL);
1.1304 + xassert(X->type == A_NONE);
1.1305 + xassert(X->dim > 0);
1.1306 + xassert(Y != NULL);
1.1307 + xassert(Y->type == A_NONE);
1.1308 + xassert(Y->dim > 0);
1.1309 + Z = create_elemset(mpl, X->dim + Y->dim);
1.1310 + for (memx = X->head; memx != NULL; memx = memx->next)
1.1311 + { for (memy = Y->head; memy != NULL; memy = memy->next)
1.1312 + { tuple = copy_tuple(mpl, memx->tuple);
1.1313 + for (temp = memy->tuple; temp != NULL; temp = temp->next)
1.1314 + tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,
1.1315 + temp->sym));
1.1316 + add_tuple(mpl, Z, tuple);
1.1317 + }
1.1318 + }
1.1319 + delete_elemset(mpl, X);
1.1320 + delete_elemset(mpl, Y);
1.1321 + return Z;
1.1322 +}
1.1323 +
1.1324 +/**********************************************************************/
1.1325 +/* * * ELEMENTAL VARIABLES * * */
1.1326 +/**********************************************************************/
1.1327 +
1.1328 +/* (there are no specific routines for elemental variables) */
1.1329 +
1.1330 +/**********************************************************************/
1.1331 +/* * * LINEAR FORMS * * */
1.1332 +/**********************************************************************/
1.1333 +
1.1334 +/*----------------------------------------------------------------------
1.1335 +-- constant_term - create constant term.
1.1336 +--
1.1337 +-- This routine creates the linear form, which is a constant term. */
1.1338 +
1.1339 +FORMULA *constant_term(MPL *mpl, double coef)
1.1340 +{ FORMULA *form;
1.1341 + if (coef == 0.0)
1.1342 + form = NULL;
1.1343 + else
1.1344 + { form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1.1345 + form->coef = coef;
1.1346 + form->var = NULL;
1.1347 + form->next = NULL;
1.1348 + }
1.1349 + return form;
1.1350 +}
1.1351 +
1.1352 +/*----------------------------------------------------------------------
1.1353 +-- single_variable - create single variable.
1.1354 +--
1.1355 +-- This routine creates the linear form, which is a single elemental
1.1356 +-- variable. */
1.1357 +
1.1358 +FORMULA *single_variable
1.1359 +( MPL *mpl,
1.1360 + ELEMVAR *var /* referenced */
1.1361 +)
1.1362 +{ FORMULA *form;
1.1363 + xassert(var != NULL);
1.1364 + form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1.1365 + form->coef = 1.0;
1.1366 + form->var = var;
1.1367 + form->next = NULL;
1.1368 + return form;
1.1369 +}
1.1370 +
1.1371 +/*----------------------------------------------------------------------
1.1372 +-- copy_formula - make copy of linear form.
1.1373 +--
1.1374 +-- This routine returns an exact copy of linear form. */
1.1375 +
1.1376 +FORMULA *copy_formula
1.1377 +( MPL *mpl,
1.1378 + FORMULA *form /* not changed */
1.1379 +)
1.1380 +{ FORMULA *head, *tail;
1.1381 + if (form == NULL)
1.1382 + head = NULL;
1.1383 + else
1.1384 + { head = tail = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1.1385 + for (; form != NULL; form = form->next)
1.1386 + { tail->coef = form->coef;
1.1387 + tail->var = form->var;
1.1388 + if (form->next != NULL)
1.1389 +tail = (tail->next = dmp_get_atom(mpl->formulae, sizeof(FORMULA)));
1.1390 + }
1.1391 + tail->next = NULL;
1.1392 + }
1.1393 + return head;
1.1394 +}
1.1395 +
1.1396 +/*----------------------------------------------------------------------
1.1397 +-- delete_formula - delete linear form.
1.1398 +--
1.1399 +-- This routine deletes specified linear form. */
1.1400 +
1.1401 +void delete_formula
1.1402 +( MPL *mpl,
1.1403 + FORMULA *form /* destroyed */
1.1404 +)
1.1405 +{ FORMULA *temp;
1.1406 + while (form != NULL)
1.1407 + { temp = form;
1.1408 + form = form->next;
1.1409 + dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));
1.1410 + }
1.1411 + return;
1.1412 +}
1.1413 +
1.1414 +/*----------------------------------------------------------------------
1.1415 +-- linear_comb - linear combination of two linear forms.
1.1416 +--
1.1417 +-- This routine computes the linear combination:
1.1418 +--
1.1419 +-- a * fx + b * fy,
1.1420 +--
1.1421 +-- where a and b are numeric coefficients, fx and fy are linear forms
1.1422 +-- (destroyed on exit). */
1.1423 +
1.1424 +FORMULA *linear_comb
1.1425 +( MPL *mpl,
1.1426 + double a, FORMULA *fx, /* destroyed */
1.1427 + double b, FORMULA *fy /* destroyed */
1.1428 +)
1.1429 +{ FORMULA *form = NULL, *term, *temp;
1.1430 + double c0 = 0.0;
1.1431 + for (term = fx; term != NULL; term = term->next)
1.1432 + { if (term->var == NULL)
1.1433 + c0 = fp_add(mpl, c0, fp_mul(mpl, a, term->coef));
1.1434 + else
1.1435 + term->var->temp =
1.1436 + fp_add(mpl, term->var->temp, fp_mul(mpl, a, term->coef));
1.1437 + }
1.1438 + for (term = fy; term != NULL; term = term->next)
1.1439 + { if (term->var == NULL)
1.1440 + c0 = fp_add(mpl, c0, fp_mul(mpl, b, term->coef));
1.1441 + else
1.1442 + term->var->temp =
1.1443 + fp_add(mpl, term->var->temp, fp_mul(mpl, b, term->coef));
1.1444 + }
1.1445 + for (term = fx; term != NULL; term = term->next)
1.1446 + { if (term->var != NULL && term->var->temp != 0.0)
1.1447 + { temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1.1448 + temp->coef = term->var->temp, temp->var = term->var;
1.1449 + temp->next = form, form = temp;
1.1450 + term->var->temp = 0.0;
1.1451 + }
1.1452 + }
1.1453 + for (term = fy; term != NULL; term = term->next)
1.1454 + { if (term->var != NULL && term->var->temp != 0.0)
1.1455 + { temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1.1456 + temp->coef = term->var->temp, temp->var = term->var;
1.1457 + temp->next = form, form = temp;
1.1458 + term->var->temp = 0.0;
1.1459 + }
1.1460 + }
1.1461 + if (c0 != 0.0)
1.1462 + { temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1.1463 + temp->coef = c0, temp->var = NULL;
1.1464 + temp->next = form, form = temp;
1.1465 + }
1.1466 + delete_formula(mpl, fx);
1.1467 + delete_formula(mpl, fy);
1.1468 + return form;
1.1469 +}
1.1470 +
1.1471 +/*----------------------------------------------------------------------
1.1472 +-- remove_constant - remove constant term from linear form.
1.1473 +--
1.1474 +-- This routine removes constant term from linear form and stores its
1.1475 +-- value to given location. */
1.1476 +
1.1477 +FORMULA *remove_constant
1.1478 +( MPL *mpl,
1.1479 + FORMULA *form, /* destroyed */
1.1480 + double *coef /* modified */
1.1481 +)
1.1482 +{ FORMULA *head = NULL, *temp;
1.1483 + *coef = 0.0;
1.1484 + while (form != NULL)
1.1485 + { temp = form;
1.1486 + form = form->next;
1.1487 + if (temp->var == NULL)
1.1488 + { /* constant term */
1.1489 + *coef = fp_add(mpl, *coef, temp->coef);
1.1490 + dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));
1.1491 + }
1.1492 + else
1.1493 + { /* linear term */
1.1494 + temp->next = head;
1.1495 + head = temp;
1.1496 + }
1.1497 + }
1.1498 + return head;
1.1499 +}
1.1500 +
1.1501 +/*----------------------------------------------------------------------
1.1502 +-- reduce_terms - reduce identical terms in linear form.
1.1503 +--
1.1504 +-- This routine reduces identical terms in specified linear form. */
1.1505 +
1.1506 +FORMULA *reduce_terms
1.1507 +( MPL *mpl,
1.1508 + FORMULA *form /* destroyed */
1.1509 +)
1.1510 +{ FORMULA *term, *next_term;
1.1511 + double c0 = 0.0;
1.1512 + for (term = form; term != NULL; term = term->next)
1.1513 + { if (term->var == NULL)
1.1514 + c0 = fp_add(mpl, c0, term->coef);
1.1515 + else
1.1516 + term->var->temp = fp_add(mpl, term->var->temp, term->coef);
1.1517 + }
1.1518 + next_term = form, form = NULL;
1.1519 + for (term = next_term; term != NULL; term = next_term)
1.1520 + { next_term = term->next;
1.1521 + if (term->var == NULL && c0 != 0.0)
1.1522 + { term->coef = c0, c0 = 0.0;
1.1523 + term->next = form, form = term;
1.1524 + }
1.1525 + else if (term->var != NULL && term->var->temp != 0.0)
1.1526 + { term->coef = term->var->temp, term->var->temp = 0.0;
1.1527 + term->next = form, form = term;
1.1528 + }
1.1529 + else
1.1530 + dmp_free_atom(mpl->formulae, term, sizeof(FORMULA));
1.1531 + }
1.1532 + return form;
1.1533 +}
1.1534 +
1.1535 +/**********************************************************************/
1.1536 +/* * * ELEMENTAL CONSTRAINTS * * */
1.1537 +/**********************************************************************/
1.1538 +
1.1539 +/* (there are no specific routines for elemental constraints) */
1.1540 +
1.1541 +/**********************************************************************/
1.1542 +/* * * GENERIC VALUES * * */
1.1543 +/**********************************************************************/
1.1544 +
1.1545 +/*----------------------------------------------------------------------
1.1546 +-- delete_value - delete generic value.
1.1547 +--
1.1548 +-- This routine deletes specified generic value.
1.1549 +--
1.1550 +-- NOTE: The generic value to be deleted must be valid. */
1.1551 +
1.1552 +void delete_value
1.1553 +( MPL *mpl,
1.1554 + int type,
1.1555 + VALUE *value /* content destroyed */
1.1556 +)
1.1557 +{ xassert(value != NULL);
1.1558 + switch (type)
1.1559 + { case A_NONE:
1.1560 + value->none = NULL;
1.1561 + break;
1.1562 + case A_NUMERIC:
1.1563 + value->num = 0.0;
1.1564 + break;
1.1565 + case A_SYMBOLIC:
1.1566 + delete_symbol(mpl, value->sym), value->sym = NULL;
1.1567 + break;
1.1568 + case A_LOGICAL:
1.1569 + value->bit = 0;
1.1570 + break;
1.1571 + case A_TUPLE:
1.1572 + delete_tuple(mpl, value->tuple), value->tuple = NULL;
1.1573 + break;
1.1574 + case A_ELEMSET:
1.1575 + delete_elemset(mpl, value->set), value->set = NULL;
1.1576 + break;
1.1577 + case A_ELEMVAR:
1.1578 + value->var = NULL;
1.1579 + break;
1.1580 + case A_FORMULA:
1.1581 + delete_formula(mpl, value->form), value->form = NULL;
1.1582 + break;
1.1583 + case A_ELEMCON:
1.1584 + value->con = NULL;
1.1585 + break;
1.1586 + default:
1.1587 + xassert(type != type);
1.1588 + }
1.1589 + return;
1.1590 +}
1.1591 +
1.1592 +/**********************************************************************/
1.1593 +/* * * SYMBOLICALLY INDEXED ARRAYS * * */
1.1594 +/**********************************************************************/
1.1595 +
1.1596 +/*----------------------------------------------------------------------
1.1597 +-- create_array - create array.
1.1598 +--
1.1599 +-- This routine creates an array of specified type and dimension. Being
1.1600 +-- created the array is initially empty.
1.1601 +--
1.1602 +-- The type indicator determines generic values, which can be assigned
1.1603 +-- to the array members:
1.1604 +--
1.1605 +-- A_NONE - none (members have no assigned values)
1.1606 +-- A_NUMERIC - floating-point numbers
1.1607 +-- A_SYMBOLIC - symbols
1.1608 +-- A_ELEMSET - elemental sets
1.1609 +-- A_ELEMVAR - elemental variables
1.1610 +-- A_ELEMCON - elemental constraints
1.1611 +--
1.1612 +-- The dimension may be 0, in which case the array consists of the only
1.1613 +-- member (such arrays represent 0-dimensional objects). */
1.1614 +
1.1615 +ARRAY *create_array(MPL *mpl, int type, int dim)
1.1616 +{ ARRAY *array;
1.1617 + xassert(type == A_NONE || type == A_NUMERIC ||
1.1618 + type == A_SYMBOLIC || type == A_ELEMSET ||
1.1619 + type == A_ELEMVAR || type == A_ELEMCON);
1.1620 + xassert(dim >= 0);
1.1621 + array = dmp_get_atom(mpl->arrays, sizeof(ARRAY));
1.1622 + array->type = type;
1.1623 + array->dim = dim;
1.1624 + array->size = 0;
1.1625 + array->head = NULL;
1.1626 + array->tail = NULL;
1.1627 + array->tree = NULL;
1.1628 + array->prev = NULL;
1.1629 + array->next = mpl->a_list;
1.1630 + /* include the array in the global array list */
1.1631 + if (array->next != NULL) array->next->prev = array;
1.1632 + mpl->a_list = array;
1.1633 + return array;
1.1634 +}
1.1635 +
1.1636 +/*----------------------------------------------------------------------
1.1637 +-- find_member - find array member with given n-tuple.
1.1638 +--
1.1639 +-- This routine finds an array member, which has given n-tuple. If the
1.1640 +-- array is short, the linear search is used. Otherwise the routine
1.1641 +-- autimatically creates the search tree (i.e. the array index) to find
1.1642 +-- members for logarithmic time. */
1.1643 +
1.1644 +static int compare_member_tuples(void *info, const void *key1,
1.1645 + const void *key2)
1.1646 +{ /* this is an auxiliary routine used to compare keys, which are
1.1647 + n-tuples assigned to array members */
1.1648 + return compare_tuples((MPL *)info, (TUPLE *)key1, (TUPLE *)key2);
1.1649 +}
1.1650 +
1.1651 +MEMBER *find_member
1.1652 +( MPL *mpl,
1.1653 + ARRAY *array, /* not changed */
1.1654 + TUPLE *tuple /* not changed */
1.1655 +)
1.1656 +{ MEMBER *memb;
1.1657 + xassert(array != NULL);
1.1658 + /* the n-tuple must have the same dimension as the array */
1.1659 + xassert(tuple_dimen(mpl, tuple) == array->dim);
1.1660 + /* if the array is large enough, create the search tree and index
1.1661 + all existing members of the array */
1.1662 + if (array->size > 30 && array->tree == NULL)
1.1663 + { array->tree = avl_create_tree(compare_member_tuples, mpl);
1.1664 + for (memb = array->head; memb != NULL; memb = memb->next)
1.1665 +avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
1.1666 + (void *)memb);
1.1667 + }
1.1668 + /* find a member, which has the given tuple */
1.1669 + if (array->tree == NULL)
1.1670 + { /* the search tree doesn't exist; use the linear search */
1.1671 + for (memb = array->head; memb != NULL; memb = memb->next)
1.1672 + if (compare_tuples(mpl, memb->tuple, tuple) == 0) break;
1.1673 + }
1.1674 + else
1.1675 + { /* the search tree exists; use the binary search */
1.1676 + AVLNODE *node;
1.1677 + node = avl_find_node(array->tree, tuple);
1.1678 +memb = (MEMBER *)(node == NULL ? NULL : avl_get_node_link(node));
1.1679 + }
1.1680 + return memb;
1.1681 +}
1.1682 +
1.1683 +/*----------------------------------------------------------------------
1.1684 +-- add_member - add new member to array.
1.1685 +--
1.1686 +-- This routine creates a new member with given n-tuple and adds it to
1.1687 +-- specified array.
1.1688 +--
1.1689 +-- For the sake of efficiency this routine doesn't check whether the
1.1690 +-- array already contains a member with the given n-tuple or not. Thus,
1.1691 +-- if necessary, the calling program should use the routine find_member
1.1692 +-- in order to be sure that the array contains no member with the same
1.1693 +-- n-tuple, because members with duplicate n-tuples are not allowed.
1.1694 +--
1.1695 +-- This routine assigns no generic value to the new member, because the
1.1696 +-- calling program must do that. */
1.1697 +
1.1698 +MEMBER *add_member
1.1699 +( MPL *mpl,
1.1700 + ARRAY *array, /* modified */
1.1701 + TUPLE *tuple /* destroyed */
1.1702 +)
1.1703 +{ MEMBER *memb;
1.1704 + xassert(array != NULL);
1.1705 + /* the n-tuple must have the same dimension as the array */
1.1706 + xassert(tuple_dimen(mpl, tuple) == array->dim);
1.1707 + /* create new member */
1.1708 + memb = dmp_get_atom(mpl->members, sizeof(MEMBER));
1.1709 + memb->tuple = tuple;
1.1710 + memb->next = NULL;
1.1711 + memset(&memb->value, '?', sizeof(VALUE));
1.1712 + /* and append it to the member list */
1.1713 + array->size++;
1.1714 + if (array->head == NULL)
1.1715 + array->head = memb;
1.1716 + else
1.1717 + array->tail->next = memb;
1.1718 + array->tail = memb;
1.1719 + /* if the search tree exists, index the new member */
1.1720 + if (array->tree != NULL)
1.1721 +avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
1.1722 + (void *)memb);
1.1723 + return memb;
1.1724 +}
1.1725 +
1.1726 +/*----------------------------------------------------------------------
1.1727 +-- delete_array - delete array.
1.1728 +--
1.1729 +-- This routine deletes specified array.
1.1730 +--
1.1731 +-- Generic values assigned to the array members are not deleted by this
1.1732 +-- routine. The calling program itself must delete all assigned generic
1.1733 +-- values before deleting the array. */
1.1734 +
1.1735 +void delete_array
1.1736 +( MPL *mpl,
1.1737 + ARRAY *array /* destroyed */
1.1738 +)
1.1739 +{ MEMBER *memb;
1.1740 + xassert(array != NULL);
1.1741 + /* delete all existing array members */
1.1742 + while (array->head != NULL)
1.1743 + { memb = array->head;
1.1744 + array->head = memb->next;
1.1745 + delete_tuple(mpl, memb->tuple);
1.1746 + dmp_free_atom(mpl->members, memb, sizeof(MEMBER));
1.1747 + }
1.1748 + /* if the search tree exists, also delete it */
1.1749 + if (array->tree != NULL) avl_delete_tree(array->tree);
1.1750 + /* remove the array from the global array list */
1.1751 + if (array->prev == NULL)
1.1752 + mpl->a_list = array->next;
1.1753 + else
1.1754 + array->prev->next = array->next;
1.1755 + if (array->next == NULL)
1.1756 + ;
1.1757 + else
1.1758 + array->next->prev = array->prev;
1.1759 + /* delete the array descriptor */
1.1760 + dmp_free_atom(mpl->arrays, array, sizeof(ARRAY));
1.1761 + return;
1.1762 +}
1.1763 +
1.1764 +/**********************************************************************/
1.1765 +/* * * DOMAINS AND DUMMY INDICES * * */
1.1766 +/**********************************************************************/
1.1767 +
1.1768 +/*----------------------------------------------------------------------
1.1769 +-- assign_dummy_index - assign new value to dummy index.
1.1770 +--
1.1771 +-- This routine assigns new value to specified dummy index and, that is
1.1772 +-- important, invalidates all temporary resultant values, which depends
1.1773 +-- on that dummy index. */
1.1774 +
1.1775 +void assign_dummy_index
1.1776 +( MPL *mpl,
1.1777 + DOMAIN_SLOT *slot, /* modified */
1.1778 + SYMBOL *value /* not changed */
1.1779 +)
1.1780 +{ CODE *leaf, *code;
1.1781 + xassert(slot != NULL);
1.1782 + xassert(value != NULL);
1.1783 + /* delete the current value assigned to the dummy index */
1.1784 + if (slot->value != NULL)
1.1785 + { /* if the current value and the new one are identical, actual
1.1786 + assignment is not needed */
1.1787 + if (compare_symbols(mpl, slot->value, value) == 0) goto done;
1.1788 + /* delete a symbol, which is the current value */
1.1789 + delete_symbol(mpl, slot->value), slot->value = NULL;
1.1790 + }
1.1791 + /* now walk through all the pseudo-codes with op = O_INDEX, which
1.1792 + refer to the dummy index to be changed (these pseudo-codes are
1.1793 + leaves in the forest of *all* expressions in the database) */
1.1794 + for (leaf = slot->list; leaf != NULL; leaf = leaf->arg.index.
1.1795 + next)
1.1796 + { xassert(leaf->op == O_INDEX);
1.1797 + /* invalidate all resultant values, which depend on the dummy
1.1798 + index, walking from the current leaf toward the root of the
1.1799 + corresponding expression tree */
1.1800 + for (code = leaf; code != NULL; code = code->up)
1.1801 + { if (code->valid)
1.1802 + { /* invalidate and delete resultant value */
1.1803 + code->valid = 0;
1.1804 + delete_value(mpl, code->type, &code->value);
1.1805 + }
1.1806 + }
1.1807 + }
1.1808 + /* assign new value to the dummy index */
1.1809 + slot->value = copy_symbol(mpl, value);
1.1810 +done: return;
1.1811 +}
1.1812 +
1.1813 +/*----------------------------------------------------------------------
1.1814 +-- update_dummy_indices - update current values of dummy indices.
1.1815 +--
1.1816 +-- This routine assigns components of "backup" n-tuple to dummy indices
1.1817 +-- of specified domain block. If no "backup" n-tuple is defined for the
1.1818 +-- domain block, values of the dummy indices remain untouched. */
1.1819 +
1.1820 +void update_dummy_indices
1.1821 +( MPL *mpl,
1.1822 + DOMAIN_BLOCK *block /* not changed */
1.1823 +)
1.1824 +{ DOMAIN_SLOT *slot;
1.1825 + TUPLE *temp;
1.1826 + if (block->backup != NULL)
1.1827 + { for (slot = block->list, temp = block->backup; slot != NULL;
1.1828 + slot = slot->next, temp = temp->next)
1.1829 + { xassert(temp != NULL);
1.1830 + xassert(temp->sym != NULL);
1.1831 + assign_dummy_index(mpl, slot, temp->sym);
1.1832 + }
1.1833 + }
1.1834 + return;
1.1835 +}
1.1836 +
1.1837 +/*----------------------------------------------------------------------
1.1838 +-- enter_domain_block - enter domain block.
1.1839 +--
1.1840 +-- Let specified domain block have the form:
1.1841 +--
1.1842 +-- { ..., (j1, j2, ..., jn) in J, ... }
1.1843 +--
1.1844 +-- where j1, j2, ..., jn are dummy indices, J is a basic set.
1.1845 +--
1.1846 +-- This routine does the following:
1.1847 +--
1.1848 +-- 1. Checks if the given n-tuple is a member of the basic set J. Note
1.1849 +-- that J being *out of the scope* of the domain block cannot depend
1.1850 +-- on the dummy indices in the same and inner domain blocks, so it
1.1851 +-- can be computed before the dummy indices are assigned new values.
1.1852 +-- If this check fails, the routine returns with non-zero code.
1.1853 +--
1.1854 +-- 2. Saves current values of the dummy indices j1, j2, ..., jn.
1.1855 +--
1.1856 +-- 3. Assigns new values, which are components of the given n-tuple, to
1.1857 +-- the dummy indices j1, j2, ..., jn. If dimension of the n-tuple is
1.1858 +-- larger than n, its extra components n+1, n+2, ... are not used.
1.1859 +--
1.1860 +-- 4. Calls the formal routine func which either enters the next domain
1.1861 +-- block or evaluates some code within the domain scope.
1.1862 +--
1.1863 +-- 5. Restores former values of the dummy indices j1, j2, ..., jn.
1.1864 +--
1.1865 +-- Since current values assigned to the dummy indices on entry to this
1.1866 +-- routine are restored on exit, the formal routine func is allowed to
1.1867 +-- call this routine recursively. */
1.1868 +
1.1869 +int enter_domain_block
1.1870 +( MPL *mpl,
1.1871 + DOMAIN_BLOCK *block, /* not changed */
1.1872 + TUPLE *tuple, /* not changed */
1.1873 + void *info, void (*func)(MPL *mpl, void *info)
1.1874 +)
1.1875 +{ TUPLE *backup;
1.1876 + int ret = 0;
1.1877 + /* check if the given n-tuple is a member of the basic set */
1.1878 + xassert(block->code != NULL);
1.1879 + if (!is_member(mpl, block->code, tuple))
1.1880 + { ret = 1;
1.1881 + goto done;
1.1882 + }
1.1883 + /* save reference to "backup" n-tuple, which was used to assign
1.1884 + current values of the dummy indices (it is sufficient to save
1.1885 + reference, not value, because that n-tuple is defined in some
1.1886 + outer level of recursion and therefore cannot be changed on
1.1887 + this and deeper recursive calls) */
1.1888 + backup = block->backup;
1.1889 + /* set up new "backup" n-tuple, which defines new values of the
1.1890 + dummy indices */
1.1891 + block->backup = tuple;
1.1892 + /* assign new values to the dummy indices */
1.1893 + update_dummy_indices(mpl, block);
1.1894 + /* call the formal routine that does the rest part of the job */
1.1895 + func(mpl, info);
1.1896 + /* restore reference to the former "backup" n-tuple */
1.1897 + block->backup = backup;
1.1898 + /* restore former values of the dummy indices; note that if the
1.1899 + domain block just escaped has no other active instances which
1.1900 + may exist due to recursion (it is indicated by a null pointer
1.1901 + to the former n-tuple), former values of the dummy indices are
1.1902 + undefined; therefore in this case the routine keeps currently
1.1903 + assigned values of the dummy indices that involves keeping all
1.1904 + dependent temporary results and thereby, if this domain block
1.1905 + is not used recursively, allows improving efficiency */
1.1906 + update_dummy_indices(mpl, block);
1.1907 +done: return ret;
1.1908 +}
1.1909 +
1.1910 +/*----------------------------------------------------------------------
1.1911 +-- eval_within_domain - perform evaluation within domain scope.
1.1912 +--
1.1913 +-- This routine assigns new values (symbols) to all dummy indices of
1.1914 +-- specified domain and calls the formal routine func, which is used to
1.1915 +-- evaluate some code in the domain scope. Each free dummy index in the
1.1916 +-- domain is assigned a value specified in the corresponding component
1.1917 +-- of given n-tuple. Non-free dummy indices are assigned values, which
1.1918 +-- are computed by this routine.
1.1919 +--
1.1920 +-- Number of components in the given n-tuple must be the same as number
1.1921 +-- of free indices in the domain.
1.1922 +--
1.1923 +-- If the given n-tuple is not a member of the domain set, the routine
1.1924 +-- func is not called, and non-zero code is returned.
1.1925 +--
1.1926 +-- For the sake of convenience it is allowed to specify domain as NULL
1.1927 +-- (then n-tuple also must be 0-tuple, i.e. empty), in which case this
1.1928 +-- routine just calls the routine func and returns zero.
1.1929 +--
1.1930 +-- This routine allows recursive calls from the routine func providing
1.1931 +-- correct values of dummy indices for each instance.
1.1932 +--
1.1933 +-- NOTE: The n-tuple passed to this routine must not be changed by any
1.1934 +-- other routines called from the formal routine func until this
1.1935 +-- routine has returned. */
1.1936 +
1.1937 +struct eval_domain_info
1.1938 +{ /* working info used by the routine eval_within_domain */
1.1939 + DOMAIN *domain;
1.1940 + /* domain, which has to be entered */
1.1941 + DOMAIN_BLOCK *block;
1.1942 + /* domain block, which is currently processed */
1.1943 + TUPLE *tuple;
1.1944 + /* tail of original n-tuple, whose components have to be assigned
1.1945 + to free dummy indices in the current domain block */
1.1946 + void *info;
1.1947 + /* transit pointer passed to the formal routine func */
1.1948 + void (*func)(MPL *mpl, void *info);
1.1949 + /* routine, which has to be executed in the domain scope */
1.1950 + int failure;
1.1951 + /* this flag indicates that given n-tuple is not a member of the
1.1952 + domain set */
1.1953 +};
1.1954 +
1.1955 +static void eval_domain_func(MPL *mpl, void *_my_info)
1.1956 +{ /* this routine recursively enters into the domain scope and then
1.1957 + calls the routine func */
1.1958 + struct eval_domain_info *my_info = _my_info;
1.1959 + if (my_info->block != NULL)
1.1960 + { /* the current domain block to be entered exists */
1.1961 + DOMAIN_BLOCK *block;
1.1962 + DOMAIN_SLOT *slot;
1.1963 + TUPLE *tuple = NULL, *temp = NULL;
1.1964 + /* save pointer to the current domain block */
1.1965 + block = my_info->block;
1.1966 + /* and get ready to enter the next block (if it exists) */
1.1967 + my_info->block = block->next;
1.1968 + /* construct temporary n-tuple, whose components correspond to
1.1969 + dummy indices (slots) of the current domain; components of
1.1970 + the temporary n-tuple that correspond to free dummy indices
1.1971 + are assigned references (not values!) to symbols specified
1.1972 + in the corresponding components of the given n-tuple, while
1.1973 + other components that correspond to non-free dummy indices
1.1974 + are assigned symbolic values computed here */
1.1975 + for (slot = block->list; slot != NULL; slot = slot->next)
1.1976 + { /* create component that corresponds to the current slot */
1.1977 + if (tuple == NULL)
1.1978 + tuple = temp = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
1.1979 + else
1.1980 +temp = (temp->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));
1.1981 + if (slot->code == NULL)
1.1982 + { /* dummy index is free; take reference to symbol, which
1.1983 + is specified in the corresponding component of given
1.1984 + n-tuple */
1.1985 + xassert(my_info->tuple != NULL);
1.1986 + temp->sym = my_info->tuple->sym;
1.1987 + xassert(temp->sym != NULL);
1.1988 + my_info->tuple = my_info->tuple->next;
1.1989 + }
1.1990 + else
1.1991 + { /* dummy index is non-free; compute symbolic value to be
1.1992 + temporarily assigned to the dummy index */
1.1993 + temp->sym = eval_symbolic(mpl, slot->code);
1.1994 + }
1.1995 + }
1.1996 + temp->next = NULL;
1.1997 + /* enter the current domain block */
1.1998 + if (enter_domain_block(mpl, block, tuple, my_info,
1.1999 + eval_domain_func)) my_info->failure = 1;
1.2000 + /* delete temporary n-tuple as well as symbols that correspond
1.2001 + to non-free dummy indices (they were computed here) */
1.2002 + for (slot = block->list; slot != NULL; slot = slot->next)
1.2003 + { xassert(tuple != NULL);
1.2004 + temp = tuple;
1.2005 + tuple = tuple->next;
1.2006 + if (slot->code != NULL)
1.2007 + { /* dummy index is non-free; delete symbolic value */
1.2008 + delete_symbol(mpl, temp->sym);
1.2009 + }
1.2010 + /* delete component that corresponds to the current slot */
1.2011 + dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));
1.2012 + }
1.2013 + }
1.2014 + else
1.2015 + { /* there are no more domain blocks, i.e. we have reached the
1.2016 + domain scope */
1.2017 + xassert(my_info->tuple == NULL);
1.2018 + /* check optional predicate specified for the domain */
1.2019 + if (my_info->domain->code != NULL && !eval_logical(mpl,
1.2020 + my_info->domain->code))
1.2021 + { /* the predicate is false */
1.2022 + my_info->failure = 2;
1.2023 + }
1.2024 + else
1.2025 + { /* the predicate is true; do the job */
1.2026 + my_info->func(mpl, my_info->info);
1.2027 + }
1.2028 + }
1.2029 + return;
1.2030 +}
1.2031 +
1.2032 +int eval_within_domain
1.2033 +( MPL *mpl,
1.2034 + DOMAIN *domain, /* not changed */
1.2035 + TUPLE *tuple, /* not changed */
1.2036 + void *info, void (*func)(MPL *mpl, void *info)
1.2037 +)
1.2038 +{ /* this routine performs evaluation within domain scope */
1.2039 + struct eval_domain_info _my_info, *my_info = &_my_info;
1.2040 + if (domain == NULL)
1.2041 + { xassert(tuple == NULL);
1.2042 + func(mpl, info);
1.2043 + my_info->failure = 0;
1.2044 + }
1.2045 + else
1.2046 + { xassert(tuple != NULL);
1.2047 + my_info->domain = domain;
1.2048 + my_info->block = domain->list;
1.2049 + my_info->tuple = tuple;
1.2050 + my_info->info = info;
1.2051 + my_info->func = func;
1.2052 + my_info->failure = 0;
1.2053 + /* enter the very first domain block */
1.2054 + eval_domain_func(mpl, my_info);
1.2055 + }
1.2056 + return my_info->failure;
1.2057 +}
1.2058 +
1.2059 +/*----------------------------------------------------------------------
1.2060 +-- loop_within_domain - perform iterations within domain scope.
1.2061 +--
1.2062 +-- This routine iteratively assigns new values (symbols) to the dummy
1.2063 +-- indices of specified domain by enumerating all n-tuples, which are
1.2064 +-- members of the domain set, and for every n-tuple it calls the formal
1.2065 +-- routine func to evaluate some code within the domain scope.
1.2066 +--
1.2067 +-- If the routine func returns non-zero, enumeration within the domain
1.2068 +-- is prematurely terminated.
1.2069 +--
1.2070 +-- For the sake of convenience it is allowed to specify domain as NULL,
1.2071 +-- in which case this routine just calls the routine func only once and
1.2072 +-- returns zero.
1.2073 +--
1.2074 +-- This routine allows recursive calls from the routine func providing
1.2075 +-- correct values of dummy indices for each instance. */
1.2076 +
1.2077 +struct loop_domain_info
1.2078 +{ /* working info used by the routine loop_within_domain */
1.2079 + DOMAIN *domain;
1.2080 + /* domain, which has to be entered */
1.2081 + DOMAIN_BLOCK *block;
1.2082 + /* domain block, which is currently processed */
1.2083 + int looping;
1.2084 + /* clearing this flag leads to terminating enumeration */
1.2085 + void *info;
1.2086 + /* transit pointer passed to the formal routine func */
1.2087 + int (*func)(MPL *mpl, void *info);
1.2088 + /* routine, which needs to be executed in the domain scope */
1.2089 +};
1.2090 +
1.2091 +static void loop_domain_func(MPL *mpl, void *_my_info)
1.2092 +{ /* this routine enumerates all n-tuples in the basic set of the
1.2093 + current domain block, enters recursively into the domain scope
1.2094 + for every n-tuple, and then calls the routine func */
1.2095 + struct loop_domain_info *my_info = _my_info;
1.2096 + if (my_info->block != NULL)
1.2097 + { /* the current domain block to be entered exists */
1.2098 + DOMAIN_BLOCK *block;
1.2099 + DOMAIN_SLOT *slot;
1.2100 + TUPLE *bound;
1.2101 + /* save pointer to the current domain block */
1.2102 + block = my_info->block;
1.2103 + /* and get ready to enter the next block (if it exists) */
1.2104 + my_info->block = block->next;
1.2105 + /* compute symbolic values, at which non-free dummy indices of
1.2106 + the current domain block are bound; since that values don't
1.2107 + depend on free dummy indices of the current block, they can
1.2108 + be computed once out of the enumeration loop */
1.2109 + bound = create_tuple(mpl);
1.2110 + for (slot = block->list; slot != NULL; slot = slot->next)
1.2111 + { if (slot->code != NULL)
1.2112 + bound = expand_tuple(mpl, bound, eval_symbolic(mpl,
1.2113 + slot->code));
1.2114 + }
1.2115 + /* start enumeration */
1.2116 + xassert(block->code != NULL);
1.2117 + if (block->code->op == O_DOTS)
1.2118 + { /* the basic set is "arithmetic", in which case it doesn't
1.2119 + need to be computed explicitly */
1.2120 + TUPLE *tuple;
1.2121 + int n, j;
1.2122 + double t0, tf, dt;
1.2123 + /* compute "parameters" of the basic set */
1.2124 + t0 = eval_numeric(mpl, block->code->arg.arg.x);
1.2125 + tf = eval_numeric(mpl, block->code->arg.arg.y);
1.2126 + if (block->code->arg.arg.z == NULL)
1.2127 + dt = 1.0;
1.2128 + else
1.2129 + dt = eval_numeric(mpl, block->code->arg.arg.z);
1.2130 + /* determine cardinality of the basic set */
1.2131 + n = arelset_size(mpl, t0, tf, dt);
1.2132 + /* create dummy 1-tuple for members of the basic set */
1.2133 + tuple = expand_tuple(mpl, create_tuple(mpl),
1.2134 + create_symbol_num(mpl, 0.0));
1.2135 + /* in case of "arithmetic" set there is exactly one dummy
1.2136 + index, which cannot be non-free */
1.2137 + xassert(bound == NULL);
1.2138 + /* walk through 1-tuples of the basic set */
1.2139 + for (j = 1; j <= n && my_info->looping; j++)
1.2140 + { /* construct dummy 1-tuple for the current member */
1.2141 + tuple->sym->num = arelset_member(mpl, t0, tf, dt, j);
1.2142 + /* enter the current domain block */
1.2143 + enter_domain_block(mpl, block, tuple, my_info,
1.2144 + loop_domain_func);
1.2145 + }
1.2146 + /* delete dummy 1-tuple */
1.2147 + delete_tuple(mpl, tuple);
1.2148 + }
1.2149 + else
1.2150 + { /* the basic set is of general kind, in which case it needs
1.2151 + to be explicitly computed */
1.2152 + ELEMSET *set;
1.2153 + MEMBER *memb;
1.2154 + TUPLE *temp1, *temp2;
1.2155 + /* compute the basic set */
1.2156 + set = eval_elemset(mpl, block->code);
1.2157 + /* walk through all n-tuples of the basic set */
1.2158 + for (memb = set->head; memb != NULL && my_info->looping;
1.2159 + memb = memb->next)
1.2160 + { /* all components of the current n-tuple that correspond
1.2161 + to non-free dummy indices must be feasible; otherwise
1.2162 + the n-tuple is not in the basic set */
1.2163 + temp1 = memb->tuple;
1.2164 + temp2 = bound;
1.2165 + for (slot = block->list; slot != NULL; slot = slot->next)
1.2166 + { xassert(temp1 != NULL);
1.2167 + if (slot->code != NULL)
1.2168 + { /* non-free dummy index */
1.2169 + xassert(temp2 != NULL);
1.2170 + if (compare_symbols(mpl, temp1->sym, temp2->sym)
1.2171 + != 0)
1.2172 + { /* the n-tuple is not in the basic set */
1.2173 + goto skip;
1.2174 + }
1.2175 + temp2 = temp2->next;
1.2176 + }
1.2177 + temp1 = temp1->next;
1.2178 + }
1.2179 + xassert(temp1 == NULL);
1.2180 + xassert(temp2 == NULL);
1.2181 + /* enter the current domain block */
1.2182 + enter_domain_block(mpl, block, memb->tuple, my_info,
1.2183 + loop_domain_func);
1.2184 +skip: ;
1.2185 + }
1.2186 + /* delete the basic set */
1.2187 + delete_elemset(mpl, set);
1.2188 + }
1.2189 + /* delete symbolic values binding non-free dummy indices */
1.2190 + delete_tuple(mpl, bound);
1.2191 + /* restore pointer to the current domain block */
1.2192 + my_info->block = block;
1.2193 + }
1.2194 + else
1.2195 + { /* there are no more domain blocks, i.e. we have reached the
1.2196 + domain scope */
1.2197 + /* check optional predicate specified for the domain */
1.2198 + if (my_info->domain->code != NULL && !eval_logical(mpl,
1.2199 + my_info->domain->code))
1.2200 + { /* the predicate is false */
1.2201 + /* nop */;
1.2202 + }
1.2203 + else
1.2204 + { /* the predicate is true; do the job */
1.2205 + my_info->looping = !my_info->func(mpl, my_info->info);
1.2206 + }
1.2207 + }
1.2208 + return;
1.2209 +}
1.2210 +
1.2211 +void loop_within_domain
1.2212 +( MPL *mpl,
1.2213 + DOMAIN *domain, /* not changed */
1.2214 + void *info, int (*func)(MPL *mpl, void *info)
1.2215 +)
1.2216 +{ /* this routine performs iterations within domain scope */
1.2217 + struct loop_domain_info _my_info, *my_info = &_my_info;
1.2218 + if (domain == NULL)
1.2219 + func(mpl, info);
1.2220 + else
1.2221 + { my_info->domain = domain;
1.2222 + my_info->block = domain->list;
1.2223 + my_info->looping = 1;
1.2224 + my_info->info = info;
1.2225 + my_info->func = func;
1.2226 + /* enter the very first domain block */
1.2227 + loop_domain_func(mpl, my_info);
1.2228 + }
1.2229 + return;
1.2230 +}
1.2231 +
1.2232 +/*----------------------------------------------------------------------
1.2233 +-- out_of_domain - raise domain exception.
1.2234 +--
1.2235 +-- This routine is called when a reference is made to a member of some
1.2236 +-- model object, but its n-tuple is out of the object domain. */
1.2237 +
1.2238 +void out_of_domain
1.2239 +( MPL *mpl,
1.2240 + char *name, /* not changed */
1.2241 + TUPLE *tuple /* not changed */
1.2242 +)
1.2243 +{ xassert(name != NULL);
1.2244 + xassert(tuple != NULL);
1.2245 + error(mpl, "%s%s out of domain", name, format_tuple(mpl, '[',
1.2246 + tuple));
1.2247 + /* no return */
1.2248 +}
1.2249 +
1.2250 +/*----------------------------------------------------------------------
1.2251 +-- get_domain_tuple - obtain current n-tuple from domain.
1.2252 +--
1.2253 +-- This routine constructs n-tuple, whose components are current values
1.2254 +-- assigned to *free* dummy indices of specified domain.
1.2255 +--
1.2256 +-- For the sake of convenience it is allowed to specify domain as NULL,
1.2257 +-- in which case this routine returns 0-tuple.
1.2258 +--
1.2259 +-- NOTE: This routine must not be called out of domain scope. */
1.2260 +
1.2261 +TUPLE *get_domain_tuple
1.2262 +( MPL *mpl,
1.2263 + DOMAIN *domain /* not changed */
1.2264 +)
1.2265 +{ DOMAIN_BLOCK *block;
1.2266 + DOMAIN_SLOT *slot;
1.2267 + TUPLE *tuple;
1.2268 + tuple = create_tuple(mpl);
1.2269 + if (domain != NULL)
1.2270 + { for (block = domain->list; block != NULL; block = block->next)
1.2271 + { for (slot = block->list; slot != NULL; slot = slot->next)
1.2272 + { if (slot->code == NULL)
1.2273 + { xassert(slot->value != NULL);
1.2274 + tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,
1.2275 + slot->value));
1.2276 + }
1.2277 + }
1.2278 + }
1.2279 + }
1.2280 + return tuple;
1.2281 +}
1.2282 +
1.2283 +/*----------------------------------------------------------------------
1.2284 +-- clean_domain - clean domain.
1.2285 +--
1.2286 +-- This routine cleans specified domain that assumes deleting all stuff
1.2287 +-- dynamically allocated during the generation phase. */
1.2288 +
1.2289 +void clean_domain(MPL *mpl, DOMAIN *domain)
1.2290 +{ DOMAIN_BLOCK *block;
1.2291 + DOMAIN_SLOT *slot;
1.2292 + /* if no domain is specified, do nothing */
1.2293 + if (domain == NULL) goto done;
1.2294 + /* clean all domain blocks */
1.2295 + for (block = domain->list; block != NULL; block = block->next)
1.2296 + { /* clean all domain slots */
1.2297 + for (slot = block->list; slot != NULL; slot = slot->next)
1.2298 + { /* clean pseudo-code for computing bound value */
1.2299 + clean_code(mpl, slot->code);
1.2300 + /* delete symbolic value assigned to dummy index */
1.2301 + if (slot->value != NULL)
1.2302 + delete_symbol(mpl, slot->value), slot->value = NULL;
1.2303 + }
1.2304 + /* clean pseudo-code for computing basic set */
1.2305 + clean_code(mpl, block->code);
1.2306 + }
1.2307 + /* clean pseudo-code for computing domain predicate */
1.2308 + clean_code(mpl, domain->code);
1.2309 +done: return;
1.2310 +}
1.2311 +
1.2312 +/**********************************************************************/
1.2313 +/* * * MODEL SETS * * */
1.2314 +/**********************************************************************/
1.2315 +
1.2316 +/*----------------------------------------------------------------------
1.2317 +-- check_elem_set - check elemental set assigned to set member.
1.2318 +--
1.2319 +-- This routine checks if given elemental set being assigned to member
1.2320 +-- of specified model set satisfies to all restrictions.
1.2321 +--
1.2322 +-- NOTE: This routine must not be called out of domain scope. */
1.2323 +
1.2324 +void check_elem_set
1.2325 +( MPL *mpl,
1.2326 + SET *set, /* not changed */
1.2327 + TUPLE *tuple, /* not changed */
1.2328 + ELEMSET *refer /* not changed */
1.2329 +)
1.2330 +{ WITHIN *within;
1.2331 + MEMBER *memb;
1.2332 + int eqno;
1.2333 + /* elemental set must be within all specified supersets */
1.2334 + for (within = set->within, eqno = 1; within != NULL; within =
1.2335 + within->next, eqno++)
1.2336 + { xassert(within->code != NULL);
1.2337 + for (memb = refer->head; memb != NULL; memb = memb->next)
1.2338 + { if (!is_member(mpl, within->code, memb->tuple))
1.2339 + { char buf[255+1];
1.2340 + strcpy(buf, format_tuple(mpl, '(', memb->tuple));
1.2341 + xassert(strlen(buf) < sizeof(buf));
1.2342 + error(mpl, "%s%s contains %s which not within specified "
1.2343 + "set; see (%d)", set->name, format_tuple(mpl, '[',
1.2344 + tuple), buf, eqno);
1.2345 + }
1.2346 + }
1.2347 + }
1.2348 + return;
1.2349 +}
1.2350 +
1.2351 +/*----------------------------------------------------------------------
1.2352 +-- take_member_set - obtain elemental set assigned to set member.
1.2353 +--
1.2354 +-- This routine obtains a reference to elemental set assigned to given
1.2355 +-- member of specified model set and returns it on exit.
1.2356 +--
1.2357 +-- NOTE: This routine must not be called out of domain scope. */
1.2358 +
1.2359 +ELEMSET *take_member_set /* returns reference, not value */
1.2360 +( MPL *mpl,
1.2361 + SET *set, /* not changed */
1.2362 + TUPLE *tuple /* not changed */
1.2363 +)
1.2364 +{ MEMBER *memb;
1.2365 + ELEMSET *refer;
1.2366 + /* find member in the set array */
1.2367 + memb = find_member(mpl, set->array, tuple);
1.2368 + if (memb != NULL)
1.2369 + { /* member exists, so just take the reference */
1.2370 + refer = memb->value.set;
1.2371 + }
1.2372 + else if (set->assign != NULL)
1.2373 + { /* compute value using assignment expression */
1.2374 + refer = eval_elemset(mpl, set->assign);
1.2375 +add: /* check that the elemental set satisfies to all restrictions,
1.2376 + assign it to new member, and add the member to the array */
1.2377 + check_elem_set(mpl, set, tuple, refer);
1.2378 + memb = add_member(mpl, set->array, copy_tuple(mpl, tuple));
1.2379 + memb->value.set = refer;
1.2380 + }
1.2381 + else if (set->option != NULL)
1.2382 + { /* compute default elemental set */
1.2383 + refer = eval_elemset(mpl, set->option);
1.2384 + goto add;
1.2385 + }
1.2386 + else
1.2387 + { /* no value (elemental set) is provided */
1.2388 + error(mpl, "no value for %s%s", set->name, format_tuple(mpl,
1.2389 + '[', tuple));
1.2390 + }
1.2391 + return refer;
1.2392 +}
1.2393 +
1.2394 +/*----------------------------------------------------------------------
1.2395 +-- eval_member_set - evaluate elemental set assigned to set member.
1.2396 +--
1.2397 +-- This routine evaluates a reference to elemental set assigned to given
1.2398 +-- member of specified model set and returns it on exit. */
1.2399 +
1.2400 +struct eval_set_info
1.2401 +{ /* working info used by the routine eval_member_set */
1.2402 + SET *set;
1.2403 + /* model set */
1.2404 + TUPLE *tuple;
1.2405 + /* n-tuple, which defines set member */
1.2406 + MEMBER *memb;
1.2407 + /* normally this pointer is NULL; the routine uses this pointer
1.2408 + to check data provided in the data section, in which case it
1.2409 + points to a member currently checked; this check is performed
1.2410 + automatically only once when a reference to any member occurs
1.2411 + for the first time */
1.2412 + ELEMSET *refer;
1.2413 + /* evaluated reference to elemental set */
1.2414 +};
1.2415 +
1.2416 +static void eval_set_func(MPL *mpl, void *_info)
1.2417 +{ /* this is auxiliary routine to work within domain scope */
1.2418 + struct eval_set_info *info = _info;
1.2419 + if (info->memb != NULL)
1.2420 + { /* checking call; check elemental set being assigned */
1.2421 + check_elem_set(mpl, info->set, info->memb->tuple,
1.2422 + info->memb->value.set);
1.2423 + }
1.2424 + else
1.2425 + { /* normal call; evaluate member, which has given n-tuple */
1.2426 + info->refer = take_member_set(mpl, info->set, info->tuple);
1.2427 + }
1.2428 + return;
1.2429 +}
1.2430 +
1.2431 +#if 1 /* 12/XII-2008 */
1.2432 +static void saturate_set(MPL *mpl, SET *set)
1.2433 +{ GADGET *gadget = set->gadget;
1.2434 + ELEMSET *data;
1.2435 + MEMBER *elem, *memb;
1.2436 + TUPLE *tuple, *work[20];
1.2437 + int i;
1.2438 + xprintf("Generating %s...\n", set->name);
1.2439 + eval_whole_set(mpl, gadget->set);
1.2440 + /* gadget set must have exactly one member */
1.2441 + xassert(gadget->set->array != NULL);
1.2442 + xassert(gadget->set->array->head != NULL);
1.2443 + xassert(gadget->set->array->head == gadget->set->array->tail);
1.2444 + data = gadget->set->array->head->value.set;
1.2445 + xassert(data->type == A_NONE);
1.2446 + xassert(data->dim == gadget->set->dimen);
1.2447 + /* walk thru all elements of the plain set */
1.2448 + for (elem = data->head; elem != NULL; elem = elem->next)
1.2449 + { /* create a copy of n-tuple */
1.2450 + tuple = copy_tuple(mpl, elem->tuple);
1.2451 + /* rearrange component of the n-tuple */
1.2452 + for (i = 0; i < gadget->set->dimen; i++)
1.2453 + work[i] = NULL;
1.2454 + for (i = 0; tuple != NULL; tuple = tuple->next)
1.2455 + work[gadget->ind[i++]-1] = tuple;
1.2456 + xassert(i == gadget->set->dimen);
1.2457 + for (i = 0; i < gadget->set->dimen; i++)
1.2458 + { xassert(work[i] != NULL);
1.2459 + work[i]->next = work[i+1];
1.2460 + }
1.2461 + /* construct subscript list from first set->dim components */
1.2462 + if (set->dim == 0)
1.2463 + tuple = NULL;
1.2464 + else
1.2465 + tuple = work[0], work[set->dim-1]->next = NULL;
1.2466 + /* find corresponding member of the set to be initialized */
1.2467 + memb = find_member(mpl, set->array, tuple);
1.2468 + if (memb == NULL)
1.2469 + { /* not found; add new member to the set and assign it empty
1.2470 + elemental set */
1.2471 + memb = add_member(mpl, set->array, tuple);
1.2472 + memb->value.set = create_elemset(mpl, set->dimen);
1.2473 + }
1.2474 + else
1.2475 + { /* found; free subscript list */
1.2476 + delete_tuple(mpl, tuple);
1.2477 + }
1.2478 + /* construct new n-tuple from rest set->dimen components */
1.2479 + tuple = work[set->dim];
1.2480 + xassert(set->dim + set->dimen == gadget->set->dimen);
1.2481 + work[gadget->set->dimen-1]->next = NULL;
1.2482 + /* and add it to the elemental set assigned to the member
1.2483 + (no check for duplicates is needed) */
1.2484 + add_tuple(mpl, memb->value.set, tuple);
1.2485 + }
1.2486 + /* the set has been saturated with data */
1.2487 + set->data = 1;
1.2488 + return;
1.2489 +}
1.2490 +#endif
1.2491 +
1.2492 +ELEMSET *eval_member_set /* returns reference, not value */
1.2493 +( MPL *mpl,
1.2494 + SET *set, /* not changed */
1.2495 + TUPLE *tuple /* not changed */
1.2496 +)
1.2497 +{ /* this routine evaluates set member */
1.2498 + struct eval_set_info _info, *info = &_info;
1.2499 + xassert(set->dim == tuple_dimen(mpl, tuple));
1.2500 + info->set = set;
1.2501 + info->tuple = tuple;
1.2502 +#if 1 /* 12/XII-2008 */
1.2503 + if (set->gadget != NULL && set->data == 0)
1.2504 + { /* initialize the set with data from a plain set */
1.2505 + saturate_set(mpl, set);
1.2506 + }
1.2507 +#endif
1.2508 + if (set->data == 1)
1.2509 + { /* check data, which are provided in the data section, but not
1.2510 + checked yet */
1.2511 + /* save pointer to the last array member; note that during the
1.2512 + check new members may be added beyond the last member due to
1.2513 + references to the same parameter from default expression as
1.2514 + well as from expressions that define restricting supersets;
1.2515 + however, values assigned to the new members will be checked
1.2516 + by other routine, so we don't need to check them here */
1.2517 + MEMBER *tail = set->array->tail;
1.2518 + /* change the data status to prevent infinite recursive loop
1.2519 + due to references to the same set during the check */
1.2520 + set->data = 2;
1.2521 + /* check elemental sets assigned to array members in the data
1.2522 + section until the marked member has been reached */
1.2523 + for (info->memb = set->array->head; info->memb != NULL;
1.2524 + info->memb = info->memb->next)
1.2525 + { if (eval_within_domain(mpl, set->domain, info->memb->tuple,
1.2526 + info, eval_set_func))
1.2527 + out_of_domain(mpl, set->name, info->memb->tuple);
1.2528 + if (info->memb == tail) break;
1.2529 + }
1.2530 + /* the check has been finished */
1.2531 + }
1.2532 + /* evaluate member, which has given n-tuple */
1.2533 + info->memb = NULL;
1.2534 + if (eval_within_domain(mpl, info->set->domain, info->tuple, info,
1.2535 + eval_set_func))
1.2536 + out_of_domain(mpl, set->name, info->tuple);
1.2537 + /* bring evaluated reference to the calling program */
1.2538 + return info->refer;
1.2539 +}
1.2540 +
1.2541 +/*----------------------------------------------------------------------
1.2542 +-- eval_whole_set - evaluate model set over entire domain.
1.2543 +--
1.2544 +-- This routine evaluates all members of specified model set over entire
1.2545 +-- domain. */
1.2546 +
1.2547 +static int whole_set_func(MPL *mpl, void *info)
1.2548 +{ /* this is auxiliary routine to work within domain scope */
1.2549 + SET *set = (SET *)info;
1.2550 + TUPLE *tuple = get_domain_tuple(mpl, set->domain);
1.2551 + eval_member_set(mpl, set, tuple);
1.2552 + delete_tuple(mpl, tuple);
1.2553 + return 0;
1.2554 +}
1.2555 +
1.2556 +void eval_whole_set(MPL *mpl, SET *set)
1.2557 +{ loop_within_domain(mpl, set->domain, set, whole_set_func);
1.2558 + return;
1.2559 +}
1.2560 +
1.2561 +/*----------------------------------------------------------------------
1.2562 +-- clean set - clean model set.
1.2563 +--
1.2564 +-- This routine cleans specified model set that assumes deleting all
1.2565 +-- stuff dynamically allocated during the generation phase. */
1.2566 +
1.2567 +void clean_set(MPL *mpl, SET *set)
1.2568 +{ WITHIN *within;
1.2569 + MEMBER *memb;
1.2570 + /* clean subscript domain */
1.2571 + clean_domain(mpl, set->domain);
1.2572 + /* clean pseudo-code for computing supersets */
1.2573 + for (within = set->within; within != NULL; within = within->next)
1.2574 + clean_code(mpl, within->code);
1.2575 + /* clean pseudo-code for computing assigned value */
1.2576 + clean_code(mpl, set->assign);
1.2577 + /* clean pseudo-code for computing default value */
1.2578 + clean_code(mpl, set->option);
1.2579 + /* reset data status flag */
1.2580 + set->data = 0;
1.2581 + /* delete content array */
1.2582 + for (memb = set->array->head; memb != NULL; memb = memb->next)
1.2583 + delete_value(mpl, set->array->type, &memb->value);
1.2584 + delete_array(mpl, set->array), set->array = NULL;
1.2585 + return;
1.2586 +}
1.2587 +
1.2588 +/**********************************************************************/
1.2589 +/* * * MODEL PARAMETERS * * */
1.2590 +/**********************************************************************/
1.2591 +
1.2592 +/*----------------------------------------------------------------------
1.2593 +-- check_value_num - check numeric value assigned to parameter member.
1.2594 +--
1.2595 +-- This routine checks if numeric value being assigned to some member
1.2596 +-- of specified numeric model parameter satisfies to all restrictions.
1.2597 +--
1.2598 +-- NOTE: This routine must not be called out of domain scope. */
1.2599 +
1.2600 +void check_value_num
1.2601 +( MPL *mpl,
1.2602 + PARAMETER *par, /* not changed */
1.2603 + TUPLE *tuple, /* not changed */
1.2604 + double value
1.2605 +)
1.2606 +{ CONDITION *cond;
1.2607 + WITHIN *in;
1.2608 + int eqno;
1.2609 + /* the value must satisfy to the parameter type */
1.2610 + switch (par->type)
1.2611 + { case A_NUMERIC:
1.2612 + break;
1.2613 + case A_INTEGER:
1.2614 + if (value != floor(value))
1.2615 + error(mpl, "%s%s = %.*g not integer", par->name,
1.2616 + format_tuple(mpl, '[', tuple), DBL_DIG, value);
1.2617 + break;
1.2618 + case A_BINARY:
1.2619 + if (!(value == 0.0 || value == 1.0))
1.2620 + error(mpl, "%s%s = %.*g not binary", par->name,
1.2621 + format_tuple(mpl, '[', tuple), DBL_DIG, value);
1.2622 + break;
1.2623 + default:
1.2624 + xassert(par != par);
1.2625 + }
1.2626 + /* the value must satisfy to all specified conditions */
1.2627 + for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,
1.2628 + eqno++)
1.2629 + { double bound;
1.2630 + char *rho;
1.2631 + xassert(cond->code != NULL);
1.2632 + bound = eval_numeric(mpl, cond->code);
1.2633 + switch (cond->rho)
1.2634 + { case O_LT:
1.2635 + if (!(value < bound))
1.2636 + { rho = "<";
1.2637 +err: error(mpl, "%s%s = %.*g not %s %.*g; see (%d)",
1.2638 + par->name, format_tuple(mpl, '[', tuple), DBL_DIG,
1.2639 + value, rho, DBL_DIG, bound, eqno);
1.2640 + }
1.2641 + break;
1.2642 + case O_LE:
1.2643 + if (!(value <= bound)) { rho = "<="; goto err; }
1.2644 + break;
1.2645 + case O_EQ:
1.2646 + if (!(value == bound)) { rho = "="; goto err; }
1.2647 + break;
1.2648 + case O_GE:
1.2649 + if (!(value >= bound)) { rho = ">="; goto err; }
1.2650 + break;
1.2651 + case O_GT:
1.2652 + if (!(value > bound)) { rho = ">"; goto err; }
1.2653 + break;
1.2654 + case O_NE:
1.2655 + if (!(value != bound)) { rho = "<>"; goto err; }
1.2656 + break;
1.2657 + default:
1.2658 + xassert(cond != cond);
1.2659 + }
1.2660 + }
1.2661 + /* the value must be in all specified supersets */
1.2662 + for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)
1.2663 + { TUPLE *dummy;
1.2664 + xassert(in->code != NULL);
1.2665 + xassert(in->code->dim == 1);
1.2666 + dummy = expand_tuple(mpl, create_tuple(mpl),
1.2667 + create_symbol_num(mpl, value));
1.2668 + if (!is_member(mpl, in->code, dummy))
1.2669 + error(mpl, "%s%s = %.*g not in specified set; see (%d)",
1.2670 + par->name, format_tuple(mpl, '[', tuple), DBL_DIG,
1.2671 + value, eqno);
1.2672 + delete_tuple(mpl, dummy);
1.2673 + }
1.2674 + return;
1.2675 +}
1.2676 +
1.2677 +/*----------------------------------------------------------------------
1.2678 +-- take_member_num - obtain num. value assigned to parameter member.
1.2679 +--
1.2680 +-- This routine obtains a numeric value assigned to member of specified
1.2681 +-- numeric model parameter and returns it on exit.
1.2682 +--
1.2683 +-- NOTE: This routine must not be called out of domain scope. */
1.2684 +
1.2685 +double take_member_num
1.2686 +( MPL *mpl,
1.2687 + PARAMETER *par, /* not changed */
1.2688 + TUPLE *tuple /* not changed */
1.2689 +)
1.2690 +{ MEMBER *memb;
1.2691 + double value;
1.2692 + /* find member in the parameter array */
1.2693 + memb = find_member(mpl, par->array, tuple);
1.2694 + if (memb != NULL)
1.2695 + { /* member exists, so just take its value */
1.2696 + value = memb->value.num;
1.2697 + }
1.2698 + else if (par->assign != NULL)
1.2699 + { /* compute value using assignment expression */
1.2700 + value = eval_numeric(mpl, par->assign);
1.2701 +add: /* check that the value satisfies to all restrictions, assign
1.2702 + it to new member, and add the member to the array */
1.2703 + check_value_num(mpl, par, tuple, value);
1.2704 + memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));
1.2705 + memb->value.num = value;
1.2706 + }
1.2707 + else if (par->option != NULL)
1.2708 + { /* compute default value */
1.2709 + value = eval_numeric(mpl, par->option);
1.2710 + goto add;
1.2711 + }
1.2712 + else if (par->defval != NULL)
1.2713 + { /* take default value provided in the data section */
1.2714 + if (par->defval->str != NULL)
1.2715 + error(mpl, "cannot convert %s to floating-point number",
1.2716 + format_symbol(mpl, par->defval));
1.2717 + value = par->defval->num;
1.2718 + goto add;
1.2719 + }
1.2720 + else
1.2721 + { /* no value is provided */
1.2722 + error(mpl, "no value for %s%s", par->name, format_tuple(mpl,
1.2723 + '[', tuple));
1.2724 + }
1.2725 + return value;
1.2726 +}
1.2727 +
1.2728 +/*----------------------------------------------------------------------
1.2729 +-- eval_member_num - evaluate num. value assigned to parameter member.
1.2730 +--
1.2731 +-- This routine evaluates a numeric value assigned to given member of
1.2732 +-- specified numeric model parameter and returns it on exit. */
1.2733 +
1.2734 +struct eval_num_info
1.2735 +{ /* working info used by the routine eval_member_num */
1.2736 + PARAMETER *par;
1.2737 + /* model parameter */
1.2738 + TUPLE *tuple;
1.2739 + /* n-tuple, which defines parameter member */
1.2740 + MEMBER *memb;
1.2741 + /* normally this pointer is NULL; the routine uses this pointer
1.2742 + to check data provided in the data section, in which case it
1.2743 + points to a member currently checked; this check is performed
1.2744 + automatically only once when a reference to any member occurs
1.2745 + for the first time */
1.2746 + double value;
1.2747 + /* evaluated numeric value */
1.2748 +};
1.2749 +
1.2750 +static void eval_num_func(MPL *mpl, void *_info)
1.2751 +{ /* this is auxiliary routine to work within domain scope */
1.2752 + struct eval_num_info *info = _info;
1.2753 + if (info->memb != NULL)
1.2754 + { /* checking call; check numeric value being assigned */
1.2755 + check_value_num(mpl, info->par, info->memb->tuple,
1.2756 + info->memb->value.num);
1.2757 + }
1.2758 + else
1.2759 + { /* normal call; evaluate member, which has given n-tuple */
1.2760 + info->value = take_member_num(mpl, info->par, info->tuple);
1.2761 + }
1.2762 + return;
1.2763 +}
1.2764 +
1.2765 +double eval_member_num
1.2766 +( MPL *mpl,
1.2767 + PARAMETER *par, /* not changed */
1.2768 + TUPLE *tuple /* not changed */
1.2769 +)
1.2770 +{ /* this routine evaluates numeric parameter member */
1.2771 + struct eval_num_info _info, *info = &_info;
1.2772 + xassert(par->type == A_NUMERIC || par->type == A_INTEGER ||
1.2773 + par->type == A_BINARY);
1.2774 + xassert(par->dim == tuple_dimen(mpl, tuple));
1.2775 + info->par = par;
1.2776 + info->tuple = tuple;
1.2777 + if (par->data == 1)
1.2778 + { /* check data, which are provided in the data section, but not
1.2779 + checked yet */
1.2780 + /* save pointer to the last array member; note that during the
1.2781 + check new members may be added beyond the last member due to
1.2782 + references to the same parameter from default expression as
1.2783 + well as from expressions that define restricting conditions;
1.2784 + however, values assigned to the new members will be checked
1.2785 + by other routine, so we don't need to check them here */
1.2786 + MEMBER *tail = par->array->tail;
1.2787 + /* change the data status to prevent infinite recursive loop
1.2788 + due to references to the same parameter during the check */
1.2789 + par->data = 2;
1.2790 + /* check values assigned to array members in the data section
1.2791 + until the marked member has been reached */
1.2792 + for (info->memb = par->array->head; info->memb != NULL;
1.2793 + info->memb = info->memb->next)
1.2794 + { if (eval_within_domain(mpl, par->domain, info->memb->tuple,
1.2795 + info, eval_num_func))
1.2796 + out_of_domain(mpl, par->name, info->memb->tuple);
1.2797 + if (info->memb == tail) break;
1.2798 + }
1.2799 + /* the check has been finished */
1.2800 + }
1.2801 + /* evaluate member, which has given n-tuple */
1.2802 + info->memb = NULL;
1.2803 + if (eval_within_domain(mpl, info->par->domain, info->tuple, info,
1.2804 + eval_num_func))
1.2805 + out_of_domain(mpl, par->name, info->tuple);
1.2806 + /* bring evaluated value to the calling program */
1.2807 + return info->value;
1.2808 +}
1.2809 +
1.2810 +/*----------------------------------------------------------------------
1.2811 +-- check_value_sym - check symbolic value assigned to parameter member.
1.2812 +--
1.2813 +-- This routine checks if symbolic value being assigned to some member
1.2814 +-- of specified symbolic model parameter satisfies to all restrictions.
1.2815 +--
1.2816 +-- NOTE: This routine must not be called out of domain scope. */
1.2817 +
1.2818 +void check_value_sym
1.2819 +( MPL *mpl,
1.2820 + PARAMETER *par, /* not changed */
1.2821 + TUPLE *tuple, /* not changed */
1.2822 + SYMBOL *value /* not changed */
1.2823 +)
1.2824 +{ CONDITION *cond;
1.2825 + WITHIN *in;
1.2826 + int eqno;
1.2827 + /* the value must satisfy to all specified conditions */
1.2828 + for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,
1.2829 + eqno++)
1.2830 + { SYMBOL *bound;
1.2831 + char buf[255+1];
1.2832 + xassert(cond->code != NULL);
1.2833 + bound = eval_symbolic(mpl, cond->code);
1.2834 + switch (cond->rho)
1.2835 + {
1.2836 +#if 1 /* 13/VIII-2008 */
1.2837 + case O_LT:
1.2838 + if (!(compare_symbols(mpl, value, bound) < 0))
1.2839 + { strcpy(buf, format_symbol(mpl, bound));
1.2840 + xassert(strlen(buf) < sizeof(buf));
1.2841 + error(mpl, "%s%s = %s not < %s",
1.2842 + par->name, format_tuple(mpl, '[', tuple),
1.2843 + format_symbol(mpl, value), buf, eqno);
1.2844 + }
1.2845 + break;
1.2846 + case O_LE:
1.2847 + if (!(compare_symbols(mpl, value, bound) <= 0))
1.2848 + { strcpy(buf, format_symbol(mpl, bound));
1.2849 + xassert(strlen(buf) < sizeof(buf));
1.2850 + error(mpl, "%s%s = %s not <= %s",
1.2851 + par->name, format_tuple(mpl, '[', tuple),
1.2852 + format_symbol(mpl, value), buf, eqno);
1.2853 + }
1.2854 + break;
1.2855 +#endif
1.2856 + case O_EQ:
1.2857 + if (!(compare_symbols(mpl, value, bound) == 0))
1.2858 + { strcpy(buf, format_symbol(mpl, bound));
1.2859 + xassert(strlen(buf) < sizeof(buf));
1.2860 + error(mpl, "%s%s = %s not = %s",
1.2861 + par->name, format_tuple(mpl, '[', tuple),
1.2862 + format_symbol(mpl, value), buf, eqno);
1.2863 + }
1.2864 + break;
1.2865 +#if 1 /* 13/VIII-2008 */
1.2866 + case O_GE:
1.2867 + if (!(compare_symbols(mpl, value, bound) >= 0))
1.2868 + { strcpy(buf, format_symbol(mpl, bound));
1.2869 + xassert(strlen(buf) < sizeof(buf));
1.2870 + error(mpl, "%s%s = %s not >= %s",
1.2871 + par->name, format_tuple(mpl, '[', tuple),
1.2872 + format_symbol(mpl, value), buf, eqno);
1.2873 + }
1.2874 + break;
1.2875 + case O_GT:
1.2876 + if (!(compare_symbols(mpl, value, bound) > 0))
1.2877 + { strcpy(buf, format_symbol(mpl, bound));
1.2878 + xassert(strlen(buf) < sizeof(buf));
1.2879 + error(mpl, "%s%s = %s not > %s",
1.2880 + par->name, format_tuple(mpl, '[', tuple),
1.2881 + format_symbol(mpl, value), buf, eqno);
1.2882 + }
1.2883 + break;
1.2884 +#endif
1.2885 + case O_NE:
1.2886 + if (!(compare_symbols(mpl, value, bound) != 0))
1.2887 + { strcpy(buf, format_symbol(mpl, bound));
1.2888 + xassert(strlen(buf) < sizeof(buf));
1.2889 + error(mpl, "%s%s = %s not <> %s",
1.2890 + par->name, format_tuple(mpl, '[', tuple),
1.2891 + format_symbol(mpl, value), buf, eqno);
1.2892 + }
1.2893 + break;
1.2894 + default:
1.2895 + xassert(cond != cond);
1.2896 + }
1.2897 + delete_symbol(mpl, bound);
1.2898 + }
1.2899 + /* the value must be in all specified supersets */
1.2900 + for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)
1.2901 + { TUPLE *dummy;
1.2902 + xassert(in->code != NULL);
1.2903 + xassert(in->code->dim == 1);
1.2904 + dummy = expand_tuple(mpl, create_tuple(mpl), copy_symbol(mpl,
1.2905 + value));
1.2906 + if (!is_member(mpl, in->code, dummy))
1.2907 + error(mpl, "%s%s = %s not in specified set; see (%d)",
1.2908 + par->name, format_tuple(mpl, '[', tuple),
1.2909 + format_symbol(mpl, value), eqno);
1.2910 + delete_tuple(mpl, dummy);
1.2911 + }
1.2912 + return;
1.2913 +}
1.2914 +
1.2915 +/*----------------------------------------------------------------------
1.2916 +-- take_member_sym - obtain symb. value assigned to parameter member.
1.2917 +--
1.2918 +-- This routine obtains a symbolic value assigned to member of specified
1.2919 +-- symbolic model parameter and returns it on exit.
1.2920 +--
1.2921 +-- NOTE: This routine must not be called out of domain scope. */
1.2922 +
1.2923 +SYMBOL *take_member_sym /* returns value, not reference */
1.2924 +( MPL *mpl,
1.2925 + PARAMETER *par, /* not changed */
1.2926 + TUPLE *tuple /* not changed */
1.2927 +)
1.2928 +{ MEMBER *memb;
1.2929 + SYMBOL *value;
1.2930 + /* find member in the parameter array */
1.2931 + memb = find_member(mpl, par->array, tuple);
1.2932 + if (memb != NULL)
1.2933 + { /* member exists, so just take its value */
1.2934 + value = copy_symbol(mpl, memb->value.sym);
1.2935 + }
1.2936 + else if (par->assign != NULL)
1.2937 + { /* compute value using assignment expression */
1.2938 + value = eval_symbolic(mpl, par->assign);
1.2939 +add: /* check that the value satisfies to all restrictions, assign
1.2940 + it to new member, and add the member to the array */
1.2941 + check_value_sym(mpl, par, tuple, value);
1.2942 + memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));
1.2943 + memb->value.sym = copy_symbol(mpl, value);
1.2944 + }
1.2945 + else if (par->option != NULL)
1.2946 + { /* compute default value */
1.2947 + value = eval_symbolic(mpl, par->option);
1.2948 + goto add;
1.2949 + }
1.2950 + else if (par->defval != NULL)
1.2951 + { /* take default value provided in the data section */
1.2952 + value = copy_symbol(mpl, par->defval);
1.2953 + goto add;
1.2954 + }
1.2955 + else
1.2956 + { /* no value is provided */
1.2957 + error(mpl, "no value for %s%s", par->name, format_tuple(mpl,
1.2958 + '[', tuple));
1.2959 + }
1.2960 + return value;
1.2961 +}
1.2962 +
1.2963 +/*----------------------------------------------------------------------
1.2964 +-- eval_member_sym - evaluate symb. value assigned to parameter member.
1.2965 +--
1.2966 +-- This routine evaluates a symbolic value assigned to given member of
1.2967 +-- specified symbolic model parameter and returns it on exit. */
1.2968 +
1.2969 +struct eval_sym_info
1.2970 +{ /* working info used by the routine eval_member_sym */
1.2971 + PARAMETER *par;
1.2972 + /* model parameter */
1.2973 + TUPLE *tuple;
1.2974 + /* n-tuple, which defines parameter member */
1.2975 + MEMBER *memb;
1.2976 + /* normally this pointer is NULL; the routine uses this pointer
1.2977 + to check data provided in the data section, in which case it
1.2978 + points to a member currently checked; this check is performed
1.2979 + automatically only once when a reference to any member occurs
1.2980 + for the first time */
1.2981 + SYMBOL *value;
1.2982 + /* evaluated symbolic value */
1.2983 +};
1.2984 +
1.2985 +static void eval_sym_func(MPL *mpl, void *_info)
1.2986 +{ /* this is auxiliary routine to work within domain scope */
1.2987 + struct eval_sym_info *info = _info;
1.2988 + if (info->memb != NULL)
1.2989 + { /* checking call; check symbolic value being assigned */
1.2990 + check_value_sym(mpl, info->par, info->memb->tuple,
1.2991 + info->memb->value.sym);
1.2992 + }
1.2993 + else
1.2994 + { /* normal call; evaluate member, which has given n-tuple */
1.2995 + info->value = take_member_sym(mpl, info->par, info->tuple);
1.2996 + }
1.2997 + return;
1.2998 +}
1.2999 +
1.3000 +SYMBOL *eval_member_sym /* returns value, not reference */
1.3001 +( MPL *mpl,
1.3002 + PARAMETER *par, /* not changed */
1.3003 + TUPLE *tuple /* not changed */
1.3004 +)
1.3005 +{ /* this routine evaluates symbolic parameter member */
1.3006 + struct eval_sym_info _info, *info = &_info;
1.3007 + xassert(par->type == A_SYMBOLIC);
1.3008 + xassert(par->dim == tuple_dimen(mpl, tuple));
1.3009 + info->par = par;
1.3010 + info->tuple = tuple;
1.3011 + if (par->data == 1)
1.3012 + { /* check data, which are provided in the data section, but not
1.3013 + checked yet */
1.3014 + /* save pointer to the last array member; note that during the
1.3015 + check new members may be added beyond the last member due to
1.3016 + references to the same parameter from default expression as
1.3017 + well as from expressions that define restricting conditions;
1.3018 + however, values assigned to the new members will be checked
1.3019 + by other routine, so we don't need to check them here */
1.3020 + MEMBER *tail = par->array->tail;
1.3021 + /* change the data status to prevent infinite recursive loop
1.3022 + due to references to the same parameter during the check */
1.3023 + par->data = 2;
1.3024 + /* check values assigned to array members in the data section
1.3025 + until the marked member has been reached */
1.3026 + for (info->memb = par->array->head; info->memb != NULL;
1.3027 + info->memb = info->memb->next)
1.3028 + { if (eval_within_domain(mpl, par->domain, info->memb->tuple,
1.3029 + info, eval_sym_func))
1.3030 + out_of_domain(mpl, par->name, info->memb->tuple);
1.3031 + if (info->memb == tail) break;
1.3032 + }
1.3033 + /* the check has been finished */
1.3034 + }
1.3035 + /* evaluate member, which has given n-tuple */
1.3036 + info->memb = NULL;
1.3037 + if (eval_within_domain(mpl, info->par->domain, info->tuple, info,
1.3038 + eval_sym_func))
1.3039 + out_of_domain(mpl, par->name, info->tuple);
1.3040 + /* bring evaluated value to the calling program */
1.3041 + return info->value;
1.3042 +}
1.3043 +
1.3044 +/*----------------------------------------------------------------------
1.3045 +-- eval_whole_par - evaluate model parameter over entire domain.
1.3046 +--
1.3047 +-- This routine evaluates all members of specified model parameter over
1.3048 +-- entire domain. */
1.3049 +
1.3050 +static int whole_par_func(MPL *mpl, void *info)
1.3051 +{ /* this is auxiliary routine to work within domain scope */
1.3052 + PARAMETER *par = (PARAMETER *)info;
1.3053 + TUPLE *tuple = get_domain_tuple(mpl, par->domain);
1.3054 + switch (par->type)
1.3055 + { case A_NUMERIC:
1.3056 + case A_INTEGER:
1.3057 + case A_BINARY:
1.3058 + eval_member_num(mpl, par, tuple);
1.3059 + break;
1.3060 + case A_SYMBOLIC:
1.3061 + delete_symbol(mpl, eval_member_sym(mpl, par, tuple));
1.3062 + break;
1.3063 + default:
1.3064 + xassert(par != par);
1.3065 + }
1.3066 + delete_tuple(mpl, tuple);
1.3067 + return 0;
1.3068 +}
1.3069 +
1.3070 +void eval_whole_par(MPL *mpl, PARAMETER *par)
1.3071 +{ loop_within_domain(mpl, par->domain, par, whole_par_func);
1.3072 + return;
1.3073 +}
1.3074 +
1.3075 +/*----------------------------------------------------------------------
1.3076 +-- clean_parameter - clean model parameter.
1.3077 +--
1.3078 +-- This routine cleans specified model parameter that assumes deleting
1.3079 +-- all stuff dynamically allocated during the generation phase. */
1.3080 +
1.3081 +void clean_parameter(MPL *mpl, PARAMETER *par)
1.3082 +{ CONDITION *cond;
1.3083 + WITHIN *in;
1.3084 + MEMBER *memb;
1.3085 + /* clean subscript domain */
1.3086 + clean_domain(mpl, par->domain);
1.3087 + /* clean pseudo-code for computing restricting conditions */
1.3088 + for (cond = par->cond; cond != NULL; cond = cond->next)
1.3089 + clean_code(mpl, cond->code);
1.3090 + /* clean pseudo-code for computing restricting supersets */
1.3091 + for (in = par->in; in != NULL; in = in->next)
1.3092 + clean_code(mpl, in->code);
1.3093 + /* clean pseudo-code for computing assigned value */
1.3094 + clean_code(mpl, par->assign);
1.3095 + /* clean pseudo-code for computing default value */
1.3096 + clean_code(mpl, par->option);
1.3097 + /* reset data status flag */
1.3098 + par->data = 0;
1.3099 + /* delete default symbolic value */
1.3100 + if (par->defval != NULL)
1.3101 + delete_symbol(mpl, par->defval), par->defval = NULL;
1.3102 + /* delete content array */
1.3103 + for (memb = par->array->head; memb != NULL; memb = memb->next)
1.3104 + delete_value(mpl, par->array->type, &memb->value);
1.3105 + delete_array(mpl, par->array), par->array = NULL;
1.3106 + return;
1.3107 +}
1.3108 +
1.3109 +/**********************************************************************/
1.3110 +/* * * MODEL VARIABLES * * */
1.3111 +/**********************************************************************/
1.3112 +
1.3113 +/*----------------------------------------------------------------------
1.3114 +-- take_member_var - obtain reference to elemental variable.
1.3115 +--
1.3116 +-- This routine obtains a reference to elemental variable assigned to
1.3117 +-- given member of specified model variable and returns it on exit. If
1.3118 +-- necessary, new elemental variable is created.
1.3119 +--
1.3120 +-- NOTE: This routine must not be called out of domain scope. */
1.3121 +
1.3122 +ELEMVAR *take_member_var /* returns reference */
1.3123 +( MPL *mpl,
1.3124 + VARIABLE *var, /* not changed */
1.3125 + TUPLE *tuple /* not changed */
1.3126 +)
1.3127 +{ MEMBER *memb;
1.3128 + ELEMVAR *refer;
1.3129 + /* find member in the variable array */
1.3130 + memb = find_member(mpl, var->array, tuple);
1.3131 + if (memb != NULL)
1.3132 + { /* member exists, so just take the reference */
1.3133 + refer = memb->value.var;
1.3134 + }
1.3135 + else
1.3136 + { /* member is referenced for the first time and therefore does
1.3137 + not exist; create new elemental variable, assign it to new
1.3138 + member, and add the member to the variable array */
1.3139 + memb = add_member(mpl, var->array, copy_tuple(mpl, tuple));
1.3140 + refer = (memb->value.var =
1.3141 + dmp_get_atom(mpl->elemvars, sizeof(ELEMVAR)));
1.3142 + refer->j = 0;
1.3143 + refer->var = var;
1.3144 + refer->memb = memb;
1.3145 + /* compute lower bound */
1.3146 + if (var->lbnd == NULL)
1.3147 + refer->lbnd = 0.0;
1.3148 + else
1.3149 + refer->lbnd = eval_numeric(mpl, var->lbnd);
1.3150 + /* compute upper bound */
1.3151 + if (var->ubnd == NULL)
1.3152 + refer->ubnd = 0.0;
1.3153 + else if (var->ubnd == var->lbnd)
1.3154 + refer->ubnd = refer->lbnd;
1.3155 + else
1.3156 + refer->ubnd = eval_numeric(mpl, var->ubnd);
1.3157 + /* nullify working quantity */
1.3158 + refer->temp = 0.0;
1.3159 +#if 1 /* 15/V-2010 */
1.3160 + /* solution has not been obtained by the solver yet */
1.3161 + refer->stat = 0;
1.3162 + refer->prim = refer->dual = 0.0;
1.3163 +#endif
1.3164 + }
1.3165 + return refer;
1.3166 +}
1.3167 +
1.3168 +/*----------------------------------------------------------------------
1.3169 +-- eval_member_var - evaluate reference to elemental variable.
1.3170 +--
1.3171 +-- This routine evaluates a reference to elemental variable assigned to
1.3172 +-- member of specified model variable and returns it on exit. */
1.3173 +
1.3174 +struct eval_var_info
1.3175 +{ /* working info used by the routine eval_member_var */
1.3176 + VARIABLE *var;
1.3177 + /* model variable */
1.3178 + TUPLE *tuple;
1.3179 + /* n-tuple, which defines variable member */
1.3180 + ELEMVAR *refer;
1.3181 + /* evaluated reference to elemental variable */
1.3182 +};
1.3183 +
1.3184 +static void eval_var_func(MPL *mpl, void *_info)
1.3185 +{ /* this is auxiliary routine to work within domain scope */
1.3186 + struct eval_var_info *info = _info;
1.3187 + info->refer = take_member_var(mpl, info->var, info->tuple);
1.3188 + return;
1.3189 +}
1.3190 +
1.3191 +ELEMVAR *eval_member_var /* returns reference */
1.3192 +( MPL *mpl,
1.3193 + VARIABLE *var, /* not changed */
1.3194 + TUPLE *tuple /* not changed */
1.3195 +)
1.3196 +{ /* this routine evaluates variable member */
1.3197 + struct eval_var_info _info, *info = &_info;
1.3198 + xassert(var->dim == tuple_dimen(mpl, tuple));
1.3199 + info->var = var;
1.3200 + info->tuple = tuple;
1.3201 + /* evaluate member, which has given n-tuple */
1.3202 + if (eval_within_domain(mpl, info->var->domain, info->tuple, info,
1.3203 + eval_var_func))
1.3204 + out_of_domain(mpl, var->name, info->tuple);
1.3205 + /* bring evaluated reference to the calling program */
1.3206 + return info->refer;
1.3207 +}
1.3208 +
1.3209 +/*----------------------------------------------------------------------
1.3210 +-- eval_whole_var - evaluate model variable over entire domain.
1.3211 +--
1.3212 +-- This routine evaluates all members of specified model variable over
1.3213 +-- entire domain. */
1.3214 +
1.3215 +static int whole_var_func(MPL *mpl, void *info)
1.3216 +{ /* this is auxiliary routine to work within domain scope */
1.3217 + VARIABLE *var = (VARIABLE *)info;
1.3218 + TUPLE *tuple = get_domain_tuple(mpl, var->domain);
1.3219 + eval_member_var(mpl, var, tuple);
1.3220 + delete_tuple(mpl, tuple);
1.3221 + return 0;
1.3222 +}
1.3223 +
1.3224 +void eval_whole_var(MPL *mpl, VARIABLE *var)
1.3225 +{ loop_within_domain(mpl, var->domain, var, whole_var_func);
1.3226 + return;
1.3227 +}
1.3228 +
1.3229 +/*----------------------------------------------------------------------
1.3230 +-- clean_variable - clean model variable.
1.3231 +--
1.3232 +-- This routine cleans specified model variable that assumes deleting
1.3233 +-- all stuff dynamically allocated during the generation phase. */
1.3234 +
1.3235 +void clean_variable(MPL *mpl, VARIABLE *var)
1.3236 +{ MEMBER *memb;
1.3237 + /* clean subscript domain */
1.3238 + clean_domain(mpl, var->domain);
1.3239 + /* clean code for computing lower bound */
1.3240 + clean_code(mpl, var->lbnd);
1.3241 + /* clean code for computing upper bound */
1.3242 + if (var->ubnd != var->lbnd) clean_code(mpl, var->ubnd);
1.3243 + /* delete content array */
1.3244 + for (memb = var->array->head; memb != NULL; memb = memb->next)
1.3245 + dmp_free_atom(mpl->elemvars, memb->value.var, sizeof(ELEMVAR));
1.3246 + delete_array(mpl, var->array), var->array = NULL;
1.3247 + return;
1.3248 +}
1.3249 +
1.3250 +/**********************************************************************/
1.3251 +/* * * MODEL CONSTRAINTS AND OBJECTIVES * * */
1.3252 +/**********************************************************************/
1.3253 +
1.3254 +/*----------------------------------------------------------------------
1.3255 +-- take_member_con - obtain reference to elemental constraint.
1.3256 +--
1.3257 +-- This routine obtains a reference to elemental constraint assigned
1.3258 +-- to given member of specified model constraint and returns it on exit.
1.3259 +-- If necessary, new elemental constraint is created.
1.3260 +--
1.3261 +-- NOTE: This routine must not be called out of domain scope. */
1.3262 +
1.3263 +ELEMCON *take_member_con /* returns reference */
1.3264 +( MPL *mpl,
1.3265 + CONSTRAINT *con, /* not changed */
1.3266 + TUPLE *tuple /* not changed */
1.3267 +)
1.3268 +{ MEMBER *memb;
1.3269 + ELEMCON *refer;
1.3270 + /* find member in the constraint array */
1.3271 + memb = find_member(mpl, con->array, tuple);
1.3272 + if (memb != NULL)
1.3273 + { /* member exists, so just take the reference */
1.3274 + refer = memb->value.con;
1.3275 + }
1.3276 + else
1.3277 + { /* member is referenced for the first time and therefore does
1.3278 + not exist; create new elemental constraint, assign it to new
1.3279 + member, and add the member to the constraint array */
1.3280 + memb = add_member(mpl, con->array, copy_tuple(mpl, tuple));
1.3281 + refer = (memb->value.con =
1.3282 + dmp_get_atom(mpl->elemcons, sizeof(ELEMCON)));
1.3283 + refer->i = 0;
1.3284 + refer->con = con;
1.3285 + refer->memb = memb;
1.3286 + /* compute linear form */
1.3287 + xassert(con->code != NULL);
1.3288 + refer->form = eval_formula(mpl, con->code);
1.3289 + /* compute lower and upper bounds */
1.3290 + if (con->lbnd == NULL && con->ubnd == NULL)
1.3291 + { /* objective has no bounds */
1.3292 + double temp;
1.3293 + xassert(con->type == A_MINIMIZE || con->type == A_MAXIMIZE);
1.3294 + /* carry the constant term to the right-hand side */
1.3295 + refer->form = remove_constant(mpl, refer->form, &temp);
1.3296 + refer->lbnd = refer->ubnd = - temp;
1.3297 + }
1.3298 + else if (con->lbnd != NULL && con->ubnd == NULL)
1.3299 + { /* constraint a * x + b >= c * y + d is transformed to the
1.3300 + standard form a * x - c * y >= d - b */
1.3301 + double temp;
1.3302 + xassert(con->type == A_CONSTRAINT);
1.3303 + refer->form = linear_comb(mpl,
1.3304 + +1.0, refer->form,
1.3305 + -1.0, eval_formula(mpl, con->lbnd));
1.3306 + refer->form = remove_constant(mpl, refer->form, &temp);
1.3307 + refer->lbnd = - temp;
1.3308 + refer->ubnd = 0.0;
1.3309 + }
1.3310 + else if (con->lbnd == NULL && con->ubnd != NULL)
1.3311 + { /* constraint a * x + b <= c * y + d is transformed to the
1.3312 + standard form a * x - c * y <= d - b */
1.3313 + double temp;
1.3314 + xassert(con->type == A_CONSTRAINT);
1.3315 + refer->form = linear_comb(mpl,
1.3316 + +1.0, refer->form,
1.3317 + -1.0, eval_formula(mpl, con->ubnd));
1.3318 + refer->form = remove_constant(mpl, refer->form, &temp);
1.3319 + refer->lbnd = 0.0;
1.3320 + refer->ubnd = - temp;
1.3321 + }
1.3322 + else if (con->lbnd == con->ubnd)
1.3323 + { /* constraint a * x + b = c * y + d is transformed to the
1.3324 + standard form a * x - c * y = d - b */
1.3325 + double temp;
1.3326 + xassert(con->type == A_CONSTRAINT);
1.3327 + refer->form = linear_comb(mpl,
1.3328 + +1.0, refer->form,
1.3329 + -1.0, eval_formula(mpl, con->lbnd));
1.3330 + refer->form = remove_constant(mpl, refer->form, &temp);
1.3331 + refer->lbnd = refer->ubnd = - temp;
1.3332 + }
1.3333 + else
1.3334 + { /* ranged constraint c <= a * x + b <= d is transformed to
1.3335 + the standard form c - b <= a * x <= d - b */
1.3336 + double temp, temp1, temp2;
1.3337 + xassert(con->type == A_CONSTRAINT);
1.3338 + refer->form = remove_constant(mpl, refer->form, &temp);
1.3339 + xassert(remove_constant(mpl, eval_formula(mpl, con->lbnd),
1.3340 + &temp1) == NULL);
1.3341 + xassert(remove_constant(mpl, eval_formula(mpl, con->ubnd),
1.3342 + &temp2) == NULL);
1.3343 + refer->lbnd = fp_sub(mpl, temp1, temp);
1.3344 + refer->ubnd = fp_sub(mpl, temp2, temp);
1.3345 + }
1.3346 +#if 1 /* 15/V-2010 */
1.3347 + /* solution has not been obtained by the solver yet */
1.3348 + refer->stat = 0;
1.3349 + refer->prim = refer->dual = 0.0;
1.3350 +#endif
1.3351 + }
1.3352 + return refer;
1.3353 +}
1.3354 +
1.3355 +/*----------------------------------------------------------------------
1.3356 +-- eval_member_con - evaluate reference to elemental constraint.
1.3357 +--
1.3358 +-- This routine evaluates a reference to elemental constraint assigned
1.3359 +-- to member of specified model constraint and returns it on exit. */
1.3360 +
1.3361 +struct eval_con_info
1.3362 +{ /* working info used by the routine eval_member_con */
1.3363 + CONSTRAINT *con;
1.3364 + /* model constraint */
1.3365 + TUPLE *tuple;
1.3366 + /* n-tuple, which defines constraint member */
1.3367 + ELEMCON *refer;
1.3368 + /* evaluated reference to elemental constraint */
1.3369 +};
1.3370 +
1.3371 +static void eval_con_func(MPL *mpl, void *_info)
1.3372 +{ /* this is auxiliary routine to work within domain scope */
1.3373 + struct eval_con_info *info = _info;
1.3374 + info->refer = take_member_con(mpl, info->con, info->tuple);
1.3375 + return;
1.3376 +}
1.3377 +
1.3378 +ELEMCON *eval_member_con /* returns reference */
1.3379 +( MPL *mpl,
1.3380 + CONSTRAINT *con, /* not changed */
1.3381 + TUPLE *tuple /* not changed */
1.3382 +)
1.3383 +{ /* this routine evaluates constraint member */
1.3384 + struct eval_con_info _info, *info = &_info;
1.3385 + xassert(con->dim == tuple_dimen(mpl, tuple));
1.3386 + info->con = con;
1.3387 + info->tuple = tuple;
1.3388 + /* evaluate member, which has given n-tuple */
1.3389 + if (eval_within_domain(mpl, info->con->domain, info->tuple, info,
1.3390 + eval_con_func))
1.3391 + out_of_domain(mpl, con->name, info->tuple);
1.3392 + /* bring evaluated reference to the calling program */
1.3393 + return info->refer;
1.3394 +}
1.3395 +
1.3396 +/*----------------------------------------------------------------------
1.3397 +-- eval_whole_con - evaluate model constraint over entire domain.
1.3398 +--
1.3399 +-- This routine evaluates all members of specified model constraint over
1.3400 +-- entire domain. */
1.3401 +
1.3402 +static int whole_con_func(MPL *mpl, void *info)
1.3403 +{ /* this is auxiliary routine to work within domain scope */
1.3404 + CONSTRAINT *con = (CONSTRAINT *)info;
1.3405 + TUPLE *tuple = get_domain_tuple(mpl, con->domain);
1.3406 + eval_member_con(mpl, con, tuple);
1.3407 + delete_tuple(mpl, tuple);
1.3408 + return 0;
1.3409 +}
1.3410 +
1.3411 +void eval_whole_con(MPL *mpl, CONSTRAINT *con)
1.3412 +{ loop_within_domain(mpl, con->domain, con, whole_con_func);
1.3413 + return;
1.3414 +}
1.3415 +
1.3416 +/*----------------------------------------------------------------------
1.3417 +-- clean_constraint - clean model constraint.
1.3418 +--
1.3419 +-- This routine cleans specified model constraint that assumes deleting
1.3420 +-- all stuff dynamically allocated during the generation phase. */
1.3421 +
1.3422 +void clean_constraint(MPL *mpl, CONSTRAINT *con)
1.3423 +{ MEMBER *memb;
1.3424 + /* clean subscript domain */
1.3425 + clean_domain(mpl, con->domain);
1.3426 + /* clean code for computing main linear form */
1.3427 + clean_code(mpl, con->code);
1.3428 + /* clean code for computing lower bound */
1.3429 + clean_code(mpl, con->lbnd);
1.3430 + /* clean code for computing upper bound */
1.3431 + if (con->ubnd != con->lbnd) clean_code(mpl, con->ubnd);
1.3432 + /* delete content array */
1.3433 + for (memb = con->array->head; memb != NULL; memb = memb->next)
1.3434 + { delete_formula(mpl, memb->value.con->form);
1.3435 + dmp_free_atom(mpl->elemcons, memb->value.con, sizeof(ELEMCON));
1.3436 + }
1.3437 + delete_array(mpl, con->array), con->array = NULL;
1.3438 + return;
1.3439 +}
1.3440 +
1.3441 +/**********************************************************************/
1.3442 +/* * * PSEUDO-CODE * * */
1.3443 +/**********************************************************************/
1.3444 +
1.3445 +/*----------------------------------------------------------------------
1.3446 +-- eval_numeric - evaluate pseudo-code to determine numeric value.
1.3447 +--
1.3448 +-- This routine evaluates specified pseudo-code to determine resultant
1.3449 +-- numeric value, which is returned on exit. */
1.3450 +
1.3451 +struct iter_num_info
1.3452 +{ /* working info used by the routine iter_num_func */
1.3453 + CODE *code;
1.3454 + /* pseudo-code for iterated operation to be performed */
1.3455 + double value;
1.3456 + /* resultant value */
1.3457 +};
1.3458 +
1.3459 +static int iter_num_func(MPL *mpl, void *_info)
1.3460 +{ /* this is auxiliary routine used to perform iterated operation
1.3461 + on numeric "integrand" within domain scope */
1.3462 + struct iter_num_info *info = _info;
1.3463 + double temp;
1.3464 + temp = eval_numeric(mpl, info->code->arg.loop.x);
1.3465 + switch (info->code->op)
1.3466 + { case O_SUM:
1.3467 + /* summation over domain */
1.3468 + info->value = fp_add(mpl, info->value, temp);
1.3469 + break;
1.3470 + case O_PROD:
1.3471 + /* multiplication over domain */
1.3472 + info->value = fp_mul(mpl, info->value, temp);
1.3473 + break;
1.3474 + case O_MINIMUM:
1.3475 + /* minimum over domain */
1.3476 + if (info->value > temp) info->value = temp;
1.3477 + break;
1.3478 + case O_MAXIMUM:
1.3479 + /* maximum over domain */
1.3480 + if (info->value < temp) info->value = temp;
1.3481 + break;
1.3482 + default:
1.3483 + xassert(info != info);
1.3484 + }
1.3485 + return 0;
1.3486 +}
1.3487 +
1.3488 +double eval_numeric(MPL *mpl, CODE *code)
1.3489 +{ double value;
1.3490 + xassert(code != NULL);
1.3491 + xassert(code->type == A_NUMERIC);
1.3492 + xassert(code->dim == 0);
1.3493 + /* if the operation has a side effect, invalidate and delete the
1.3494 + resultant value */
1.3495 + if (code->vflag && code->valid)
1.3496 + { code->valid = 0;
1.3497 + delete_value(mpl, code->type, &code->value);
1.3498 + }
1.3499 + /* if resultant value is valid, no evaluation is needed */
1.3500 + if (code->valid)
1.3501 + { value = code->value.num;
1.3502 + goto done;
1.3503 + }
1.3504 + /* evaluate pseudo-code recursively */
1.3505 + switch (code->op)
1.3506 + { case O_NUMBER:
1.3507 + /* take floating-point number */
1.3508 + value = code->arg.num;
1.3509 + break;
1.3510 + case O_MEMNUM:
1.3511 + /* take member of numeric parameter */
1.3512 + { TUPLE *tuple;
1.3513 + ARG_LIST *e;
1.3514 + tuple = create_tuple(mpl);
1.3515 + for (e = code->arg.par.list; e != NULL; e = e->next)
1.3516 + tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
1.3517 + e->x));
1.3518 + value = eval_member_num(mpl, code->arg.par.par, tuple);
1.3519 + delete_tuple(mpl, tuple);
1.3520 + }
1.3521 + break;
1.3522 + case O_MEMVAR:
1.3523 + /* take computed value of elemental variable */
1.3524 + { TUPLE *tuple;
1.3525 + ARG_LIST *e;
1.3526 +#if 1 /* 15/V-2010 */
1.3527 + ELEMVAR *var;
1.3528 +#endif
1.3529 + tuple = create_tuple(mpl);
1.3530 + for (e = code->arg.var.list; e != NULL; e = e->next)
1.3531 + tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
1.3532 + e->x));
1.3533 +#if 0 /* 15/V-2010 */
1.3534 + value = eval_member_var(mpl, code->arg.var.var, tuple)
1.3535 + ->value;
1.3536 +#else
1.3537 + var = eval_member_var(mpl, code->arg.var.var, tuple);
1.3538 + switch (code->arg.var.suff)
1.3539 + { case DOT_LB:
1.3540 + if (var->var->lbnd == NULL)
1.3541 + value = -DBL_MAX;
1.3542 + else
1.3543 + value = var->lbnd;
1.3544 + break;
1.3545 + case DOT_UB:
1.3546 + if (var->var->ubnd == NULL)
1.3547 + value = +DBL_MAX;
1.3548 + else
1.3549 + value = var->ubnd;
1.3550 + break;
1.3551 + case DOT_STATUS:
1.3552 + value = var->stat;
1.3553 + break;
1.3554 + case DOT_VAL:
1.3555 + value = var->prim;
1.3556 + break;
1.3557 + case DOT_DUAL:
1.3558 + value = var->dual;
1.3559 + break;
1.3560 + default:
1.3561 + xassert(code != code);
1.3562 + }
1.3563 +#endif
1.3564 + delete_tuple(mpl, tuple);
1.3565 + }
1.3566 + break;
1.3567 +#if 1 /* 15/V-2010 */
1.3568 + case O_MEMCON:
1.3569 + /* take computed value of elemental constraint */
1.3570 + { TUPLE *tuple;
1.3571 + ARG_LIST *e;
1.3572 + ELEMCON *con;
1.3573 + tuple = create_tuple(mpl);
1.3574 + for (e = code->arg.con.list; e != NULL; e = e->next)
1.3575 + tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
1.3576 + e->x));
1.3577 + con = eval_member_con(mpl, code->arg.con.con, tuple);
1.3578 + switch (code->arg.con.suff)
1.3579 + { case DOT_LB:
1.3580 + if (con->con->lbnd == NULL)
1.3581 + value = -DBL_MAX;
1.3582 + else
1.3583 + value = con->lbnd;
1.3584 + break;
1.3585 + case DOT_UB:
1.3586 + if (con->con->ubnd == NULL)
1.3587 + value = +DBL_MAX;
1.3588 + else
1.3589 + value = con->ubnd;
1.3590 + break;
1.3591 + case DOT_STATUS:
1.3592 + value = con->stat;
1.3593 + break;
1.3594 + case DOT_VAL:
1.3595 + value = con->prim;
1.3596 + break;
1.3597 + case DOT_DUAL:
1.3598 + value = con->dual;
1.3599 + break;
1.3600 + default:
1.3601 + xassert(code != code);
1.3602 + }
1.3603 + delete_tuple(mpl, tuple);
1.3604 + }
1.3605 + break;
1.3606 +#endif
1.3607 + case O_IRAND224:
1.3608 + /* pseudo-random in [0, 2^24-1] */
1.3609 + value = fp_irand224(mpl);
1.3610 + break;
1.3611 + case O_UNIFORM01:
1.3612 + /* pseudo-random in [0, 1) */
1.3613 + value = fp_uniform01(mpl);
1.3614 + break;
1.3615 + case O_NORMAL01:
1.3616 + /* gaussian random, mu = 0, sigma = 1 */
1.3617 + value = fp_normal01(mpl);
1.3618 + break;
1.3619 + case O_GMTIME:
1.3620 + /* current calendar time */
1.3621 + value = fn_gmtime(mpl);
1.3622 + break;
1.3623 + case O_CVTNUM:
1.3624 + /* conversion to numeric */
1.3625 + { SYMBOL *sym;
1.3626 + sym = eval_symbolic(mpl, code->arg.arg.x);
1.3627 +#if 0 /* 23/XI-2008 */
1.3628 + if (sym->str != NULL)
1.3629 + error(mpl, "cannot convert %s to floating-point numbe"
1.3630 + "r", format_symbol(mpl, sym));
1.3631 + value = sym->num;
1.3632 +#else
1.3633 + if (sym->str == NULL)
1.3634 + value = sym->num;
1.3635 + else
1.3636 + { if (str2num(sym->str, &value))
1.3637 + error(mpl, "cannot convert %s to floating-point nu"
1.3638 + "mber", format_symbol(mpl, sym));
1.3639 + }
1.3640 +#endif
1.3641 + delete_symbol(mpl, sym);
1.3642 + }
1.3643 + break;
1.3644 + case O_PLUS:
1.3645 + /* unary plus */
1.3646 + value = + eval_numeric(mpl, code->arg.arg.x);
1.3647 + break;
1.3648 + case O_MINUS:
1.3649 + /* unary minus */
1.3650 + value = - eval_numeric(mpl, code->arg.arg.x);
1.3651 + break;
1.3652 + case O_ABS:
1.3653 + /* absolute value */
1.3654 + value = fabs(eval_numeric(mpl, code->arg.arg.x));
1.3655 + break;
1.3656 + case O_CEIL:
1.3657 + /* round upward ("ceiling of x") */
1.3658 + value = ceil(eval_numeric(mpl, code->arg.arg.x));
1.3659 + break;
1.3660 + case O_FLOOR:
1.3661 + /* round downward ("floor of x") */
1.3662 + value = floor(eval_numeric(mpl, code->arg.arg.x));
1.3663 + break;
1.3664 + case O_EXP:
1.3665 + /* base-e exponential */
1.3666 + value = fp_exp(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3667 + break;
1.3668 + case O_LOG:
1.3669 + /* natural logarithm */
1.3670 + value = fp_log(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3671 + break;
1.3672 + case O_LOG10:
1.3673 + /* common (decimal) logarithm */
1.3674 + value = fp_log10(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3675 + break;
1.3676 + case O_SQRT:
1.3677 + /* square root */
1.3678 + value = fp_sqrt(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3679 + break;
1.3680 + case O_SIN:
1.3681 + /* trigonometric sine */
1.3682 + value = fp_sin(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3683 + break;
1.3684 + case O_COS:
1.3685 + /* trigonometric cosine */
1.3686 + value = fp_cos(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3687 + break;
1.3688 + case O_ATAN:
1.3689 + /* trigonometric arctangent (one argument) */
1.3690 + value = fp_atan(mpl, eval_numeric(mpl, code->arg.arg.x));
1.3691 + break;
1.3692 + case O_ATAN2:
1.3693 + /* trigonometric arctangent (two arguments) */
1.3694 + value = fp_atan2(mpl,
1.3695 + eval_numeric(mpl, code->arg.arg.x),
1.3696 + eval_numeric(mpl, code->arg.arg.y));
1.3697 + break;
1.3698 + case O_ROUND:
1.3699 + /* round to nearest integer */
1.3700 + value = fp_round(mpl,
1.3701 + eval_numeric(mpl, code->arg.arg.x), 0.0);
1.3702 + break;
1.3703 + case O_ROUND2:
1.3704 + /* round to n fractional digits */
1.3705 + value = fp_round(mpl,
1.3706 + eval_numeric(mpl, code->arg.arg.x),
1.3707 + eval_numeric(mpl, code->arg.arg.y));
1.3708 + break;
1.3709 + case O_TRUNC:
1.3710 + /* truncate to nearest integer */
1.3711 + value = fp_trunc(mpl,
1.3712 + eval_numeric(mpl, code->arg.arg.x), 0.0);
1.3713 + break;
1.3714 + case O_TRUNC2:
1.3715 + /* truncate to n fractional digits */
1.3716 + value = fp_trunc(mpl,
1.3717 + eval_numeric(mpl, code->arg.arg.x),
1.3718 + eval_numeric(mpl, code->arg.arg.y));
1.3719 + break;
1.3720 + case O_ADD:
1.3721 + /* addition */
1.3722 + value = fp_add(mpl,
1.3723 + eval_numeric(mpl, code->arg.arg.x),
1.3724 + eval_numeric(mpl, code->arg.arg.y));
1.3725 + break;
1.3726 + case O_SUB:
1.3727 + /* subtraction */
1.3728 + value = fp_sub(mpl,
1.3729 + eval_numeric(mpl, code->arg.arg.x),
1.3730 + eval_numeric(mpl, code->arg.arg.y));
1.3731 + break;
1.3732 + case O_LESS:
1.3733 + /* non-negative subtraction */
1.3734 + value = fp_less(mpl,
1.3735 + eval_numeric(mpl, code->arg.arg.x),
1.3736 + eval_numeric(mpl, code->arg.arg.y));
1.3737 + break;
1.3738 + case O_MUL:
1.3739 + /* multiplication */
1.3740 + value = fp_mul(mpl,
1.3741 + eval_numeric(mpl, code->arg.arg.x),
1.3742 + eval_numeric(mpl, code->arg.arg.y));
1.3743 + break;
1.3744 + case O_DIV:
1.3745 + /* division */
1.3746 + value = fp_div(mpl,
1.3747 + eval_numeric(mpl, code->arg.arg.x),
1.3748 + eval_numeric(mpl, code->arg.arg.y));
1.3749 + break;
1.3750 + case O_IDIV:
1.3751 + /* quotient of exact division */
1.3752 + value = fp_idiv(mpl,
1.3753 + eval_numeric(mpl, code->arg.arg.x),
1.3754 + eval_numeric(mpl, code->arg.arg.y));
1.3755 + break;
1.3756 + case O_MOD:
1.3757 + /* remainder of exact division */
1.3758 + value = fp_mod(mpl,
1.3759 + eval_numeric(mpl, code->arg.arg.x),
1.3760 + eval_numeric(mpl, code->arg.arg.y));
1.3761 + break;
1.3762 + case O_POWER:
1.3763 + /* exponentiation (raise to power) */
1.3764 + value = fp_power(mpl,
1.3765 + eval_numeric(mpl, code->arg.arg.x),
1.3766 + eval_numeric(mpl, code->arg.arg.y));
1.3767 + break;
1.3768 + case O_UNIFORM:
1.3769 + /* pseudo-random in [a, b) */
1.3770 + value = fp_uniform(mpl,
1.3771 + eval_numeric(mpl, code->arg.arg.x),
1.3772 + eval_numeric(mpl, code->arg.arg.y));
1.3773 + break;
1.3774 + case O_NORMAL:
1.3775 + /* gaussian random, given mu and sigma */
1.3776 + value = fp_normal(mpl,
1.3777 + eval_numeric(mpl, code->arg.arg.x),
1.3778 + eval_numeric(mpl, code->arg.arg.y));
1.3779 + break;
1.3780 + case O_CARD:
1.3781 + { ELEMSET *set;
1.3782 + set = eval_elemset(mpl, code->arg.arg.x);
1.3783 + value = set->size;
1.3784 + delete_array(mpl, set);
1.3785 + }
1.3786 + break;
1.3787 + case O_LENGTH:
1.3788 + { SYMBOL *sym;
1.3789 + char str[MAX_LENGTH+1];
1.3790 + sym = eval_symbolic(mpl, code->arg.arg.x);
1.3791 + if (sym->str == NULL)
1.3792 + sprintf(str, "%.*g", DBL_DIG, sym->num);
1.3793 + else
1.3794 + fetch_string(mpl, sym->str, str);
1.3795 + delete_symbol(mpl, sym);
1.3796 + value = strlen(str);
1.3797 + }
1.3798 + break;
1.3799 + case O_STR2TIME:
1.3800 + { SYMBOL *sym;
1.3801 + char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];
1.3802 + sym = eval_symbolic(mpl, code->arg.arg.x);
1.3803 + if (sym->str == NULL)
1.3804 + sprintf(str, "%.*g", DBL_DIG, sym->num);
1.3805 + else
1.3806 + fetch_string(mpl, sym->str, str);
1.3807 + delete_symbol(mpl, sym);
1.3808 + sym = eval_symbolic(mpl, code->arg.arg.y);
1.3809 + if (sym->str == NULL)
1.3810 + sprintf(fmt, "%.*g", DBL_DIG, sym->num);
1.3811 + else
1.3812 + fetch_string(mpl, sym->str, fmt);
1.3813 + delete_symbol(mpl, sym);
1.3814 + value = fn_str2time(mpl, str, fmt);
1.3815 + }
1.3816 + break;
1.3817 + case O_FORK:
1.3818 + /* if-then-else */
1.3819 + if (eval_logical(mpl, code->arg.arg.x))
1.3820 + value = eval_numeric(mpl, code->arg.arg.y);
1.3821 + else if (code->arg.arg.z == NULL)
1.3822 + value = 0.0;
1.3823 + else
1.3824 + value = eval_numeric(mpl, code->arg.arg.z);
1.3825 + break;
1.3826 + case O_MIN:
1.3827 + /* minimal value (n-ary) */
1.3828 + { ARG_LIST *e;
1.3829 + double temp;
1.3830 + value = +DBL_MAX;
1.3831 + for (e = code->arg.list; e != NULL; e = e->next)
1.3832 + { temp = eval_numeric(mpl, e->x);
1.3833 + if (value > temp) value = temp;
1.3834 + }
1.3835 + }
1.3836 + break;
1.3837 + case O_MAX:
1.3838 + /* maximal value (n-ary) */
1.3839 + { ARG_LIST *e;
1.3840 + double temp;
1.3841 + value = -DBL_MAX;
1.3842 + for (e = code->arg.list; e != NULL; e = e->next)
1.3843 + { temp = eval_numeric(mpl, e->x);
1.3844 + if (value < temp) value = temp;
1.3845 + }
1.3846 + }
1.3847 + break;
1.3848 + case O_SUM:
1.3849 + /* summation over domain */
1.3850 + { struct iter_num_info _info, *info = &_info;
1.3851 + info->code = code;
1.3852 + info->value = 0.0;
1.3853 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.3854 + iter_num_func);
1.3855 + value = info->value;
1.3856 + }
1.3857 + break;
1.3858 + case O_PROD:
1.3859 + /* multiplication over domain */
1.3860 + { struct iter_num_info _info, *info = &_info;
1.3861 + info->code = code;
1.3862 + info->value = 1.0;
1.3863 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.3864 + iter_num_func);
1.3865 + value = info->value;
1.3866 + }
1.3867 + break;
1.3868 + case O_MINIMUM:
1.3869 + /* minimum over domain */
1.3870 + { struct iter_num_info _info, *info = &_info;
1.3871 + info->code = code;
1.3872 + info->value = +DBL_MAX;
1.3873 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.3874 + iter_num_func);
1.3875 + if (info->value == +DBL_MAX)
1.3876 + error(mpl, "min{} over empty set; result undefined");
1.3877 + value = info->value;
1.3878 + }
1.3879 + break;
1.3880 + case O_MAXIMUM:
1.3881 + /* maximum over domain */
1.3882 + { struct iter_num_info _info, *info = &_info;
1.3883 + info->code = code;
1.3884 + info->value = -DBL_MAX;
1.3885 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.3886 + iter_num_func);
1.3887 + if (info->value == -DBL_MAX)
1.3888 + error(mpl, "max{} over empty set; result undefined");
1.3889 + value = info->value;
1.3890 + }
1.3891 + break;
1.3892 + default:
1.3893 + xassert(code != code);
1.3894 + }
1.3895 + /* save resultant value */
1.3896 + xassert(!code->valid);
1.3897 + code->valid = 1;
1.3898 + code->value.num = value;
1.3899 +done: return value;
1.3900 +}
1.3901 +
1.3902 +/*----------------------------------------------------------------------
1.3903 +-- eval_symbolic - evaluate pseudo-code to determine symbolic value.
1.3904 +--
1.3905 +-- This routine evaluates specified pseudo-code to determine resultant
1.3906 +-- symbolic value, which is returned on exit. */
1.3907 +
1.3908 +SYMBOL *eval_symbolic(MPL *mpl, CODE *code)
1.3909 +{ SYMBOL *value;
1.3910 + xassert(code != NULL);
1.3911 + xassert(code->type == A_SYMBOLIC);
1.3912 + xassert(code->dim == 0);
1.3913 + /* if the operation has a side effect, invalidate and delete the
1.3914 + resultant value */
1.3915 + if (code->vflag && code->valid)
1.3916 + { code->valid = 0;
1.3917 + delete_value(mpl, code->type, &code->value);
1.3918 + }
1.3919 + /* if resultant value is valid, no evaluation is needed */
1.3920 + if (code->valid)
1.3921 + { value = copy_symbol(mpl, code->value.sym);
1.3922 + goto done;
1.3923 + }
1.3924 + /* evaluate pseudo-code recursively */
1.3925 + switch (code->op)
1.3926 + { case O_STRING:
1.3927 + /* take character string */
1.3928 + value = create_symbol_str(mpl, create_string(mpl,
1.3929 + code->arg.str));
1.3930 + break;
1.3931 + case O_INDEX:
1.3932 + /* take dummy index */
1.3933 + xassert(code->arg.index.slot->value != NULL);
1.3934 + value = copy_symbol(mpl, code->arg.index.slot->value);
1.3935 + break;
1.3936 + case O_MEMSYM:
1.3937 + /* take member of symbolic parameter */
1.3938 + { TUPLE *tuple;
1.3939 + ARG_LIST *e;
1.3940 + tuple = create_tuple(mpl);
1.3941 + for (e = code->arg.par.list; e != NULL; e = e->next)
1.3942 + tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
1.3943 + e->x));
1.3944 + value = eval_member_sym(mpl, code->arg.par.par, tuple);
1.3945 + delete_tuple(mpl, tuple);
1.3946 + }
1.3947 + break;
1.3948 + case O_CVTSYM:
1.3949 + /* conversion to symbolic */
1.3950 + value = create_symbol_num(mpl, eval_numeric(mpl,
1.3951 + code->arg.arg.x));
1.3952 + break;
1.3953 + case O_CONCAT:
1.3954 + /* concatenation */
1.3955 + value = concat_symbols(mpl,
1.3956 + eval_symbolic(mpl, code->arg.arg.x),
1.3957 + eval_symbolic(mpl, code->arg.arg.y));
1.3958 + break;
1.3959 + case O_FORK:
1.3960 + /* if-then-else */
1.3961 + if (eval_logical(mpl, code->arg.arg.x))
1.3962 + value = eval_symbolic(mpl, code->arg.arg.y);
1.3963 + else if (code->arg.arg.z == NULL)
1.3964 + value = create_symbol_num(mpl, 0.0);
1.3965 + else
1.3966 + value = eval_symbolic(mpl, code->arg.arg.z);
1.3967 + break;
1.3968 + case O_SUBSTR:
1.3969 + case O_SUBSTR3:
1.3970 + { double pos, len;
1.3971 + char str[MAX_LENGTH+1];
1.3972 + value = eval_symbolic(mpl, code->arg.arg.x);
1.3973 + if (value->str == NULL)
1.3974 + sprintf(str, "%.*g", DBL_DIG, value->num);
1.3975 + else
1.3976 + fetch_string(mpl, value->str, str);
1.3977 + delete_symbol(mpl, value);
1.3978 + if (code->op == O_SUBSTR)
1.3979 + { pos = eval_numeric(mpl, code->arg.arg.y);
1.3980 + if (pos != floor(pos))
1.3981 + error(mpl, "substr('...', %.*g); non-integer secon"
1.3982 + "d argument", DBL_DIG, pos);
1.3983 + if (pos < 1 || pos > strlen(str) + 1)
1.3984 + error(mpl, "substr('...', %.*g); substring out of "
1.3985 + "range", DBL_DIG, pos);
1.3986 + }
1.3987 + else
1.3988 + { pos = eval_numeric(mpl, code->arg.arg.y);
1.3989 + len = eval_numeric(mpl, code->arg.arg.z);
1.3990 + if (pos != floor(pos) || len != floor(len))
1.3991 + error(mpl, "substr('...', %.*g, %.*g); non-integer"
1.3992 + " second and/or third argument", DBL_DIG, pos,
1.3993 + DBL_DIG, len);
1.3994 + if (pos < 1 || len < 0 || pos + len > strlen(str) + 1)
1.3995 + error(mpl, "substr('...', %.*g, %.*g); substring o"
1.3996 + "ut of range", DBL_DIG, pos, DBL_DIG, len);
1.3997 + str[(int)pos + (int)len - 1] = '\0';
1.3998 + }
1.3999 + value = create_symbol_str(mpl, create_string(mpl, str +
1.4000 + (int)pos - 1));
1.4001 + }
1.4002 + break;
1.4003 + case O_TIME2STR:
1.4004 + { double num;
1.4005 + SYMBOL *sym;
1.4006 + char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];
1.4007 + num = eval_numeric(mpl, code->arg.arg.x);
1.4008 + sym = eval_symbolic(mpl, code->arg.arg.y);
1.4009 + if (sym->str == NULL)
1.4010 + sprintf(fmt, "%.*g", DBL_DIG, sym->num);
1.4011 + else
1.4012 + fetch_string(mpl, sym->str, fmt);
1.4013 + delete_symbol(mpl, sym);
1.4014 + fn_time2str(mpl, str, num, fmt);
1.4015 + value = create_symbol_str(mpl, create_string(mpl, str));
1.4016 + }
1.4017 + break;
1.4018 + default:
1.4019 + xassert(code != code);
1.4020 + }
1.4021 + /* save resultant value */
1.4022 + xassert(!code->valid);
1.4023 + code->valid = 1;
1.4024 + code->value.sym = copy_symbol(mpl, value);
1.4025 +done: return value;
1.4026 +}
1.4027 +
1.4028 +/*----------------------------------------------------------------------
1.4029 +-- eval_logical - evaluate pseudo-code to determine logical value.
1.4030 +--
1.4031 +-- This routine evaluates specified pseudo-code to determine resultant
1.4032 +-- logical value, which is returned on exit. */
1.4033 +
1.4034 +struct iter_log_info
1.4035 +{ /* working info used by the routine iter_log_func */
1.4036 + CODE *code;
1.4037 + /* pseudo-code for iterated operation to be performed */
1.4038 + int value;
1.4039 + /* resultant value */
1.4040 +};
1.4041 +
1.4042 +static int iter_log_func(MPL *mpl, void *_info)
1.4043 +{ /* this is auxiliary routine used to perform iterated operation
1.4044 + on logical "integrand" within domain scope */
1.4045 + struct iter_log_info *info = _info;
1.4046 + int ret = 0;
1.4047 + switch (info->code->op)
1.4048 + { case O_FORALL:
1.4049 + /* conjunction over domain */
1.4050 + info->value &= eval_logical(mpl, info->code->arg.loop.x);
1.4051 + if (!info->value) ret = 1;
1.4052 + break;
1.4053 + case O_EXISTS:
1.4054 + /* disjunction over domain */
1.4055 + info->value |= eval_logical(mpl, info->code->arg.loop.x);
1.4056 + if (info->value) ret = 1;
1.4057 + break;
1.4058 + default:
1.4059 + xassert(info != info);
1.4060 + }
1.4061 + return ret;
1.4062 +}
1.4063 +
1.4064 +int eval_logical(MPL *mpl, CODE *code)
1.4065 +{ int value;
1.4066 + xassert(code->type == A_LOGICAL);
1.4067 + xassert(code->dim == 0);
1.4068 + /* if the operation has a side effect, invalidate and delete the
1.4069 + resultant value */
1.4070 + if (code->vflag && code->valid)
1.4071 + { code->valid = 0;
1.4072 + delete_value(mpl, code->type, &code->value);
1.4073 + }
1.4074 + /* if resultant value is valid, no evaluation is needed */
1.4075 + if (code->valid)
1.4076 + { value = code->value.bit;
1.4077 + goto done;
1.4078 + }
1.4079 + /* evaluate pseudo-code recursively */
1.4080 + switch (code->op)
1.4081 + { case O_CVTLOG:
1.4082 + /* conversion to logical */
1.4083 + value = (eval_numeric(mpl, code->arg.arg.x) != 0.0);
1.4084 + break;
1.4085 + case O_NOT:
1.4086 + /* negation (logical "not") */
1.4087 + value = !eval_logical(mpl, code->arg.arg.x);
1.4088 + break;
1.4089 + case O_LT:
1.4090 + /* comparison on 'less than' */
1.4091 +#if 0 /* 02/VIII-2008 */
1.4092 + value = (eval_numeric(mpl, code->arg.arg.x) <
1.4093 + eval_numeric(mpl, code->arg.arg.y));
1.4094 +#else
1.4095 + xassert(code->arg.arg.x != NULL);
1.4096 + if (code->arg.arg.x->type == A_NUMERIC)
1.4097 + value = (eval_numeric(mpl, code->arg.arg.x) <
1.4098 + eval_numeric(mpl, code->arg.arg.y));
1.4099 + else
1.4100 + { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
1.4101 + SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
1.4102 + value = (compare_symbols(mpl, sym1, sym2) < 0);
1.4103 + delete_symbol(mpl, sym1);
1.4104 + delete_symbol(mpl, sym2);
1.4105 + }
1.4106 +#endif
1.4107 + break;
1.4108 + case O_LE:
1.4109 + /* comparison on 'not greater than' */
1.4110 +#if 0 /* 02/VIII-2008 */
1.4111 + value = (eval_numeric(mpl, code->arg.arg.x) <=
1.4112 + eval_numeric(mpl, code->arg.arg.y));
1.4113 +#else
1.4114 + xassert(code->arg.arg.x != NULL);
1.4115 + if (code->arg.arg.x->type == A_NUMERIC)
1.4116 + value = (eval_numeric(mpl, code->arg.arg.x) <=
1.4117 + eval_numeric(mpl, code->arg.arg.y));
1.4118 + else
1.4119 + { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
1.4120 + SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
1.4121 + value = (compare_symbols(mpl, sym1, sym2) <= 0);
1.4122 + delete_symbol(mpl, sym1);
1.4123 + delete_symbol(mpl, sym2);
1.4124 + }
1.4125 +#endif
1.4126 + break;
1.4127 + case O_EQ:
1.4128 + /* comparison on 'equal to' */
1.4129 + xassert(code->arg.arg.x != NULL);
1.4130 + if (code->arg.arg.x->type == A_NUMERIC)
1.4131 + value = (eval_numeric(mpl, code->arg.arg.x) ==
1.4132 + eval_numeric(mpl, code->arg.arg.y));
1.4133 + else
1.4134 + { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
1.4135 + SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
1.4136 + value = (compare_symbols(mpl, sym1, sym2) == 0);
1.4137 + delete_symbol(mpl, sym1);
1.4138 + delete_symbol(mpl, sym2);
1.4139 + }
1.4140 + break;
1.4141 + case O_GE:
1.4142 + /* comparison on 'not less than' */
1.4143 +#if 0 /* 02/VIII-2008 */
1.4144 + value = (eval_numeric(mpl, code->arg.arg.x) >=
1.4145 + eval_numeric(mpl, code->arg.arg.y));
1.4146 +#else
1.4147 + xassert(code->arg.arg.x != NULL);
1.4148 + if (code->arg.arg.x->type == A_NUMERIC)
1.4149 + value = (eval_numeric(mpl, code->arg.arg.x) >=
1.4150 + eval_numeric(mpl, code->arg.arg.y));
1.4151 + else
1.4152 + { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
1.4153 + SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
1.4154 + value = (compare_symbols(mpl, sym1, sym2) >= 0);
1.4155 + delete_symbol(mpl, sym1);
1.4156 + delete_symbol(mpl, sym2);
1.4157 + }
1.4158 +#endif
1.4159 + break;
1.4160 + case O_GT:
1.4161 + /* comparison on 'greater than' */
1.4162 +#if 0 /* 02/VIII-2008 */
1.4163 + value = (eval_numeric(mpl, code->arg.arg.x) >
1.4164 + eval_numeric(mpl, code->arg.arg.y));
1.4165 +#else
1.4166 + xassert(code->arg.arg.x != NULL);
1.4167 + if (code->arg.arg.x->type == A_NUMERIC)
1.4168 + value = (eval_numeric(mpl, code->arg.arg.x) >
1.4169 + eval_numeric(mpl, code->arg.arg.y));
1.4170 + else
1.4171 + { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
1.4172 + SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
1.4173 + value = (compare_symbols(mpl, sym1, sym2) > 0);
1.4174 + delete_symbol(mpl, sym1);
1.4175 + delete_symbol(mpl, sym2);
1.4176 + }
1.4177 +#endif
1.4178 + break;
1.4179 + case O_NE:
1.4180 + /* comparison on 'not equal to' */
1.4181 + xassert(code->arg.arg.x != NULL);
1.4182 + if (code->arg.arg.x->type == A_NUMERIC)
1.4183 + value = (eval_numeric(mpl, code->arg.arg.x) !=
1.4184 + eval_numeric(mpl, code->arg.arg.y));
1.4185 + else
1.4186 + { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
1.4187 + SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
1.4188 + value = (compare_symbols(mpl, sym1, sym2) != 0);
1.4189 + delete_symbol(mpl, sym1);
1.4190 + delete_symbol(mpl, sym2);
1.4191 + }
1.4192 + break;
1.4193 + case O_AND:
1.4194 + /* conjunction (logical "and") */
1.4195 + value = eval_logical(mpl, code->arg.arg.x) &&
1.4196 + eval_logical(mpl, code->arg.arg.y);
1.4197 + break;
1.4198 + case O_OR:
1.4199 + /* disjunction (logical "or") */
1.4200 + value = eval_logical(mpl, code->arg.arg.x) ||
1.4201 + eval_logical(mpl, code->arg.arg.y);
1.4202 + break;
1.4203 + case O_IN:
1.4204 + /* test on 'x in Y' */
1.4205 + { TUPLE *tuple;
1.4206 + tuple = eval_tuple(mpl, code->arg.arg.x);
1.4207 + value = is_member(mpl, code->arg.arg.y, tuple);
1.4208 + delete_tuple(mpl, tuple);
1.4209 + }
1.4210 + break;
1.4211 + case O_NOTIN:
1.4212 + /* test on 'x not in Y' */
1.4213 + { TUPLE *tuple;
1.4214 + tuple = eval_tuple(mpl, code->arg.arg.x);
1.4215 + value = !is_member(mpl, code->arg.arg.y, tuple);
1.4216 + delete_tuple(mpl, tuple);
1.4217 + }
1.4218 + break;
1.4219 + case O_WITHIN:
1.4220 + /* test on 'X within Y' */
1.4221 + { ELEMSET *set;
1.4222 + MEMBER *memb;
1.4223 + set = eval_elemset(mpl, code->arg.arg.x);
1.4224 + value = 1;
1.4225 + for (memb = set->head; memb != NULL; memb = memb->next)
1.4226 + { if (!is_member(mpl, code->arg.arg.y, memb->tuple))
1.4227 + { value = 0;
1.4228 + break;
1.4229 + }
1.4230 + }
1.4231 + delete_elemset(mpl, set);
1.4232 + }
1.4233 + break;
1.4234 + case O_NOTWITHIN:
1.4235 + /* test on 'X not within Y' */
1.4236 + { ELEMSET *set;
1.4237 + MEMBER *memb;
1.4238 + set = eval_elemset(mpl, code->arg.arg.x);
1.4239 + value = 1;
1.4240 + for (memb = set->head; memb != NULL; memb = memb->next)
1.4241 + { if (is_member(mpl, code->arg.arg.y, memb->tuple))
1.4242 + { value = 0;
1.4243 + break;
1.4244 + }
1.4245 + }
1.4246 + delete_elemset(mpl, set);
1.4247 + }
1.4248 + break;
1.4249 + case O_FORALL:
1.4250 + /* conjunction (A-quantification) */
1.4251 + { struct iter_log_info _info, *info = &_info;
1.4252 + info->code = code;
1.4253 + info->value = 1;
1.4254 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.4255 + iter_log_func);
1.4256 + value = info->value;
1.4257 + }
1.4258 + break;
1.4259 + case O_EXISTS:
1.4260 + /* disjunction (E-quantification) */
1.4261 + { struct iter_log_info _info, *info = &_info;
1.4262 + info->code = code;
1.4263 + info->value = 0;
1.4264 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.4265 + iter_log_func);
1.4266 + value = info->value;
1.4267 + }
1.4268 + break;
1.4269 + default:
1.4270 + xassert(code != code);
1.4271 + }
1.4272 + /* save resultant value */
1.4273 + xassert(!code->valid);
1.4274 + code->valid = 1;
1.4275 + code->value.bit = value;
1.4276 +done: return value;
1.4277 +}
1.4278 +
1.4279 +/*----------------------------------------------------------------------
1.4280 +-- eval_tuple - evaluate pseudo-code to construct n-tuple.
1.4281 +--
1.4282 +-- This routine evaluates specified pseudo-code to construct resultant
1.4283 +-- n-tuple, which is returned on exit. */
1.4284 +
1.4285 +TUPLE *eval_tuple(MPL *mpl, CODE *code)
1.4286 +{ TUPLE *value;
1.4287 + xassert(code != NULL);
1.4288 + xassert(code->type == A_TUPLE);
1.4289 + xassert(code->dim > 0);
1.4290 + /* if the operation has a side effect, invalidate and delete the
1.4291 + resultant value */
1.4292 + if (code->vflag && code->valid)
1.4293 + { code->valid = 0;
1.4294 + delete_value(mpl, code->type, &code->value);
1.4295 + }
1.4296 + /* if resultant value is valid, no evaluation is needed */
1.4297 + if (code->valid)
1.4298 + { value = copy_tuple(mpl, code->value.tuple);
1.4299 + goto done;
1.4300 + }
1.4301 + /* evaluate pseudo-code recursively */
1.4302 + switch (code->op)
1.4303 + { case O_TUPLE:
1.4304 + /* make n-tuple */
1.4305 + { ARG_LIST *e;
1.4306 + value = create_tuple(mpl);
1.4307 + for (e = code->arg.list; e != NULL; e = e->next)
1.4308 + value = expand_tuple(mpl, value, eval_symbolic(mpl,
1.4309 + e->x));
1.4310 + }
1.4311 + break;
1.4312 + case O_CVTTUP:
1.4313 + /* convert to 1-tuple */
1.4314 + value = expand_tuple(mpl, create_tuple(mpl),
1.4315 + eval_symbolic(mpl, code->arg.arg.x));
1.4316 + break;
1.4317 + default:
1.4318 + xassert(code != code);
1.4319 + }
1.4320 + /* save resultant value */
1.4321 + xassert(!code->valid);
1.4322 + code->valid = 1;
1.4323 + code->value.tuple = copy_tuple(mpl, value);
1.4324 +done: return value;
1.4325 +}
1.4326 +
1.4327 +/*----------------------------------------------------------------------
1.4328 +-- eval_elemset - evaluate pseudo-code to construct elemental set.
1.4329 +--
1.4330 +-- This routine evaluates specified pseudo-code to construct resultant
1.4331 +-- elemental set, which is returned on exit. */
1.4332 +
1.4333 +struct iter_set_info
1.4334 +{ /* working info used by the routine iter_set_func */
1.4335 + CODE *code;
1.4336 + /* pseudo-code for iterated operation to be performed */
1.4337 + ELEMSET *value;
1.4338 + /* resultant value */
1.4339 +};
1.4340 +
1.4341 +static int iter_set_func(MPL *mpl, void *_info)
1.4342 +{ /* this is auxiliary routine used to perform iterated operation
1.4343 + on n-tuple "integrand" within domain scope */
1.4344 + struct iter_set_info *info = _info;
1.4345 + TUPLE *tuple;
1.4346 + switch (info->code->op)
1.4347 + { case O_SETOF:
1.4348 + /* compute next n-tuple and add it to the set; in this case
1.4349 + duplicate n-tuples are silently ignored */
1.4350 + tuple = eval_tuple(mpl, info->code->arg.loop.x);
1.4351 + if (find_tuple(mpl, info->value, tuple) == NULL)
1.4352 + add_tuple(mpl, info->value, tuple);
1.4353 + else
1.4354 + delete_tuple(mpl, tuple);
1.4355 + break;
1.4356 + case O_BUILD:
1.4357 + /* construct next n-tuple using current values assigned to
1.4358 + *free* dummy indices as its components and add it to the
1.4359 + set; in this case duplicate n-tuples cannot appear */
1.4360 + add_tuple(mpl, info->value, get_domain_tuple(mpl,
1.4361 + info->code->arg.loop.domain));
1.4362 + break;
1.4363 + default:
1.4364 + xassert(info != info);
1.4365 + }
1.4366 + return 0;
1.4367 +}
1.4368 +
1.4369 +ELEMSET *eval_elemset(MPL *mpl, CODE *code)
1.4370 +{ ELEMSET *value;
1.4371 + xassert(code != NULL);
1.4372 + xassert(code->type == A_ELEMSET);
1.4373 + xassert(code->dim > 0);
1.4374 + /* if the operation has a side effect, invalidate and delete the
1.4375 + resultant value */
1.4376 + if (code->vflag && code->valid)
1.4377 + { code->valid = 0;
1.4378 + delete_value(mpl, code->type, &code->value);
1.4379 + }
1.4380 + /* if resultant value is valid, no evaluation is needed */
1.4381 + if (code->valid)
1.4382 + { value = copy_elemset(mpl, code->value.set);
1.4383 + goto done;
1.4384 + }
1.4385 + /* evaluate pseudo-code recursively */
1.4386 + switch (code->op)
1.4387 + { case O_MEMSET:
1.4388 + /* take member of set */
1.4389 + { TUPLE *tuple;
1.4390 + ARG_LIST *e;
1.4391 + tuple = create_tuple(mpl);
1.4392 + for (e = code->arg.set.list; e != NULL; e = e->next)
1.4393 + tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
1.4394 + e->x));
1.4395 + value = copy_elemset(mpl,
1.4396 + eval_member_set(mpl, code->arg.set.set, tuple));
1.4397 + delete_tuple(mpl, tuple);
1.4398 + }
1.4399 + break;
1.4400 + case O_MAKE:
1.4401 + /* make elemental set of n-tuples */
1.4402 + { ARG_LIST *e;
1.4403 + value = create_elemset(mpl, code->dim);
1.4404 + for (e = code->arg.list; e != NULL; e = e->next)
1.4405 + check_then_add(mpl, value, eval_tuple(mpl, e->x));
1.4406 + }
1.4407 + break;
1.4408 + case O_UNION:
1.4409 + /* union of two elemental sets */
1.4410 + value = set_union(mpl,
1.4411 + eval_elemset(mpl, code->arg.arg.x),
1.4412 + eval_elemset(mpl, code->arg.arg.y));
1.4413 + break;
1.4414 + case O_DIFF:
1.4415 + /* difference between two elemental sets */
1.4416 + value = set_diff(mpl,
1.4417 + eval_elemset(mpl, code->arg.arg.x),
1.4418 + eval_elemset(mpl, code->arg.arg.y));
1.4419 + break;
1.4420 + case O_SYMDIFF:
1.4421 + /* symmetric difference between two elemental sets */
1.4422 + value = set_symdiff(mpl,
1.4423 + eval_elemset(mpl, code->arg.arg.x),
1.4424 + eval_elemset(mpl, code->arg.arg.y));
1.4425 + break;
1.4426 + case O_INTER:
1.4427 + /* intersection of two elemental sets */
1.4428 + value = set_inter(mpl,
1.4429 + eval_elemset(mpl, code->arg.arg.x),
1.4430 + eval_elemset(mpl, code->arg.arg.y));
1.4431 + break;
1.4432 + case O_CROSS:
1.4433 + /* cross (Cartesian) product of two elemental sets */
1.4434 + value = set_cross(mpl,
1.4435 + eval_elemset(mpl, code->arg.arg.x),
1.4436 + eval_elemset(mpl, code->arg.arg.y));
1.4437 + break;
1.4438 + case O_DOTS:
1.4439 + /* build "arithmetic" elemental set */
1.4440 + value = create_arelset(mpl,
1.4441 + eval_numeric(mpl, code->arg.arg.x),
1.4442 + eval_numeric(mpl, code->arg.arg.y),
1.4443 + code->arg.arg.z == NULL ? 1.0 : eval_numeric(mpl,
1.4444 + code->arg.arg.z));
1.4445 + break;
1.4446 + case O_FORK:
1.4447 + /* if-then-else */
1.4448 + if (eval_logical(mpl, code->arg.arg.x))
1.4449 + value = eval_elemset(mpl, code->arg.arg.y);
1.4450 + else
1.4451 + value = eval_elemset(mpl, code->arg.arg.z);
1.4452 + break;
1.4453 + case O_SETOF:
1.4454 + /* compute elemental set */
1.4455 + { struct iter_set_info _info, *info = &_info;
1.4456 + info->code = code;
1.4457 + info->value = create_elemset(mpl, code->dim);
1.4458 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.4459 + iter_set_func);
1.4460 + value = info->value;
1.4461 + }
1.4462 + break;
1.4463 + case O_BUILD:
1.4464 + /* build elemental set identical to domain set */
1.4465 + { struct iter_set_info _info, *info = &_info;
1.4466 + info->code = code;
1.4467 + info->value = create_elemset(mpl, code->dim);
1.4468 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.4469 + iter_set_func);
1.4470 + value = info->value;
1.4471 + }
1.4472 + break;
1.4473 + default:
1.4474 + xassert(code != code);
1.4475 + }
1.4476 + /* save resultant value */
1.4477 + xassert(!code->valid);
1.4478 + code->valid = 1;
1.4479 + code->value.set = copy_elemset(mpl, value);
1.4480 +done: return value;
1.4481 +}
1.4482 +
1.4483 +/*----------------------------------------------------------------------
1.4484 +-- is_member - check if n-tuple is in set specified by pseudo-code.
1.4485 +--
1.4486 +-- This routine checks if given n-tuple is a member of elemental set
1.4487 +-- specified in the form of pseudo-code (i.e. by expression).
1.4488 +--
1.4489 +-- The n-tuple may have more components that dimension of the elemental
1.4490 +-- set, in which case the extra components are ignored. */
1.4491 +
1.4492 +static void null_func(MPL *mpl, void *info)
1.4493 +{ /* this is dummy routine used to enter the domain scope */
1.4494 + xassert(mpl == mpl);
1.4495 + xassert(info == NULL);
1.4496 + return;
1.4497 +}
1.4498 +
1.4499 +int is_member(MPL *mpl, CODE *code, TUPLE *tuple)
1.4500 +{ int value;
1.4501 + xassert(code != NULL);
1.4502 + xassert(code->type == A_ELEMSET);
1.4503 + xassert(code->dim > 0);
1.4504 + xassert(tuple != NULL);
1.4505 + switch (code->op)
1.4506 + { case O_MEMSET:
1.4507 + /* check if given n-tuple is member of elemental set, which
1.4508 + is assigned to member of model set */
1.4509 + { ARG_LIST *e;
1.4510 + TUPLE *temp;
1.4511 + ELEMSET *set;
1.4512 + /* evaluate reference to elemental set */
1.4513 + temp = create_tuple(mpl);
1.4514 + for (e = code->arg.set.list; e != NULL; e = e->next)
1.4515 + temp = expand_tuple(mpl, temp, eval_symbolic(mpl,
1.4516 + e->x));
1.4517 + set = eval_member_set(mpl, code->arg.set.set, temp);
1.4518 + delete_tuple(mpl, temp);
1.4519 + /* check if the n-tuple is contained in the set array */
1.4520 + temp = build_subtuple(mpl, tuple, set->dim);
1.4521 + value = (find_tuple(mpl, set, temp) != NULL);
1.4522 + delete_tuple(mpl, temp);
1.4523 + }
1.4524 + break;
1.4525 + case O_MAKE:
1.4526 + /* check if given n-tuple is member of literal set */
1.4527 + { ARG_LIST *e;
1.4528 + TUPLE *temp, *that;
1.4529 + value = 0;
1.4530 + temp = build_subtuple(mpl, tuple, code->dim);
1.4531 + for (e = code->arg.list; e != NULL; e = e->next)
1.4532 + { that = eval_tuple(mpl, e->x);
1.4533 + value = (compare_tuples(mpl, temp, that) == 0);
1.4534 + delete_tuple(mpl, that);
1.4535 + if (value) break;
1.4536 + }
1.4537 + delete_tuple(mpl, temp);
1.4538 + }
1.4539 + break;
1.4540 + case O_UNION:
1.4541 + value = is_member(mpl, code->arg.arg.x, tuple) ||
1.4542 + is_member(mpl, code->arg.arg.y, tuple);
1.4543 + break;
1.4544 + case O_DIFF:
1.4545 + value = is_member(mpl, code->arg.arg.x, tuple) &&
1.4546 + !is_member(mpl, code->arg.arg.y, tuple);
1.4547 + break;
1.4548 + case O_SYMDIFF:
1.4549 + { int in1 = is_member(mpl, code->arg.arg.x, tuple);
1.4550 + int in2 = is_member(mpl, code->arg.arg.y, tuple);
1.4551 + value = (in1 && !in2) || (!in1 && in2);
1.4552 + }
1.4553 + break;
1.4554 + case O_INTER:
1.4555 + value = is_member(mpl, code->arg.arg.x, tuple) &&
1.4556 + is_member(mpl, code->arg.arg.y, tuple);
1.4557 + break;
1.4558 + case O_CROSS:
1.4559 + { int j;
1.4560 + value = is_member(mpl, code->arg.arg.x, tuple);
1.4561 + if (value)
1.4562 + { for (j = 1; j <= code->arg.arg.x->dim; j++)
1.4563 + { xassert(tuple != NULL);
1.4564 + tuple = tuple->next;
1.4565 + }
1.4566 + value = is_member(mpl, code->arg.arg.y, tuple);
1.4567 + }
1.4568 + }
1.4569 + break;
1.4570 + case O_DOTS:
1.4571 + /* check if given 1-tuple is member of "arithmetic" set */
1.4572 + { int j;
1.4573 + double x, t0, tf, dt;
1.4574 + xassert(code->dim == 1);
1.4575 + /* compute "parameters" of the "arithmetic" set */
1.4576 + t0 = eval_numeric(mpl, code->arg.arg.x);
1.4577 + tf = eval_numeric(mpl, code->arg.arg.y);
1.4578 + if (code->arg.arg.z == NULL)
1.4579 + dt = 1.0;
1.4580 + else
1.4581 + dt = eval_numeric(mpl, code->arg.arg.z);
1.4582 + /* make sure the parameters are correct */
1.4583 + arelset_size(mpl, t0, tf, dt);
1.4584 + /* if component of 1-tuple is symbolic, not numeric, the
1.4585 + 1-tuple cannot be member of "arithmetic" set */
1.4586 + xassert(tuple->sym != NULL);
1.4587 + if (tuple->sym->str != NULL)
1.4588 + { value = 0;
1.4589 + break;
1.4590 + }
1.4591 + /* determine numeric value of the component */
1.4592 + x = tuple->sym->num;
1.4593 + /* if the component value is out of the set range, the
1.4594 + 1-tuple is not in the set */
1.4595 + if (dt > 0.0 && !(t0 <= x && x <= tf) ||
1.4596 + dt < 0.0 && !(tf <= x && x <= t0))
1.4597 + { value = 0;
1.4598 + break;
1.4599 + }
1.4600 + /* estimate ordinal number of the 1-tuple in the set */
1.4601 + j = (int)(((x - t0) / dt) + 0.5) + 1;
1.4602 + /* perform the main check */
1.4603 + value = (arelset_member(mpl, t0, tf, dt, j) == x);
1.4604 + }
1.4605 + break;
1.4606 + case O_FORK:
1.4607 + /* check if given n-tuple is member of conditional set */
1.4608 + if (eval_logical(mpl, code->arg.arg.x))
1.4609 + value = is_member(mpl, code->arg.arg.y, tuple);
1.4610 + else
1.4611 + value = is_member(mpl, code->arg.arg.z, tuple);
1.4612 + break;
1.4613 + case O_SETOF:
1.4614 + /* check if given n-tuple is member of computed set */
1.4615 + /* it is not clear how to efficiently perform the check not
1.4616 + computing the entire elemental set :+( */
1.4617 + error(mpl, "implementation restriction; in/within setof{} n"
1.4618 + "ot allowed");
1.4619 + break;
1.4620 + case O_BUILD:
1.4621 + /* check if given n-tuple is member of domain set */
1.4622 + { TUPLE *temp;
1.4623 + temp = build_subtuple(mpl, tuple, code->dim);
1.4624 + /* try to enter the domain scope; if it is successful,
1.4625 + the n-tuple is in the domain set */
1.4626 + value = (eval_within_domain(mpl, code->arg.loop.domain,
1.4627 + temp, NULL, null_func) == 0);
1.4628 + delete_tuple(mpl, temp);
1.4629 + }
1.4630 + break;
1.4631 + default:
1.4632 + xassert(code != code);
1.4633 + }
1.4634 + return value;
1.4635 +}
1.4636 +
1.4637 +/*----------------------------------------------------------------------
1.4638 +-- eval_formula - evaluate pseudo-code to construct linear form.
1.4639 +--
1.4640 +-- This routine evaluates specified pseudo-code to construct resultant
1.4641 +-- linear form, which is returned on exit. */
1.4642 +
1.4643 +struct iter_form_info
1.4644 +{ /* working info used by the routine iter_form_func */
1.4645 + CODE *code;
1.4646 + /* pseudo-code for iterated operation to be performed */
1.4647 + FORMULA *value;
1.4648 + /* resultant value */
1.4649 + FORMULA *tail;
1.4650 + /* pointer to the last term */
1.4651 +};
1.4652 +
1.4653 +static int iter_form_func(MPL *mpl, void *_info)
1.4654 +{ /* this is auxiliary routine used to perform iterated operation
1.4655 + on linear form "integrand" within domain scope */
1.4656 + struct iter_form_info *info = _info;
1.4657 + switch (info->code->op)
1.4658 + { case O_SUM:
1.4659 + /* summation over domain */
1.4660 +#if 0
1.4661 + info->value =
1.4662 + linear_comb(mpl,
1.4663 + +1.0, info->value,
1.4664 + +1.0, eval_formula(mpl, info->code->arg.loop.x));
1.4665 +#else
1.4666 + /* the routine linear_comb needs to look through all terms
1.4667 + of both linear forms to reduce identical terms, so using
1.4668 + it here is not a good idea (for example, evaluation of
1.4669 + sum{i in 1..n} x[i] required quadratic time); the better
1.4670 + idea is to gather all terms of the integrand in one list
1.4671 + and reduce identical terms only once after all terms of
1.4672 + the resultant linear form have been evaluated */
1.4673 + { FORMULA *form, *term;
1.4674 + form = eval_formula(mpl, info->code->arg.loop.x);
1.4675 + if (info->value == NULL)
1.4676 + { xassert(info->tail == NULL);
1.4677 + info->value = form;
1.4678 + }
1.4679 + else
1.4680 + { xassert(info->tail != NULL);
1.4681 + info->tail->next = form;
1.4682 + }
1.4683 + for (term = form; term != NULL; term = term->next)
1.4684 + info->tail = term;
1.4685 + }
1.4686 +#endif
1.4687 + break;
1.4688 + default:
1.4689 + xassert(info != info);
1.4690 + }
1.4691 + return 0;
1.4692 +}
1.4693 +
1.4694 +FORMULA *eval_formula(MPL *mpl, CODE *code)
1.4695 +{ FORMULA *value;
1.4696 + xassert(code != NULL);
1.4697 + xassert(code->type == A_FORMULA);
1.4698 + xassert(code->dim == 0);
1.4699 + /* if the operation has a side effect, invalidate and delete the
1.4700 + resultant value */
1.4701 + if (code->vflag && code->valid)
1.4702 + { code->valid = 0;
1.4703 + delete_value(mpl, code->type, &code->value);
1.4704 + }
1.4705 + /* if resultant value is valid, no evaluation is needed */
1.4706 + if (code->valid)
1.4707 + { value = copy_formula(mpl, code->value.form);
1.4708 + goto done;
1.4709 + }
1.4710 + /* evaluate pseudo-code recursively */
1.4711 + switch (code->op)
1.4712 + { case O_MEMVAR:
1.4713 + /* take member of variable */
1.4714 + { TUPLE *tuple;
1.4715 + ARG_LIST *e;
1.4716 + tuple = create_tuple(mpl);
1.4717 + for (e = code->arg.var.list; e != NULL; e = e->next)
1.4718 + tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
1.4719 + e->x));
1.4720 +#if 1 /* 15/V-2010 */
1.4721 + xassert(code->arg.var.suff == DOT_NONE);
1.4722 +#endif
1.4723 + value = single_variable(mpl,
1.4724 + eval_member_var(mpl, code->arg.var.var, tuple));
1.4725 + delete_tuple(mpl, tuple);
1.4726 + }
1.4727 + break;
1.4728 + case O_CVTLFM:
1.4729 + /* convert to linear form */
1.4730 + value = constant_term(mpl, eval_numeric(mpl,
1.4731 + code->arg.arg.x));
1.4732 + break;
1.4733 + case O_PLUS:
1.4734 + /* unary plus */
1.4735 + value = linear_comb(mpl,
1.4736 + 0.0, constant_term(mpl, 0.0),
1.4737 + +1.0, eval_formula(mpl, code->arg.arg.x));
1.4738 + break;
1.4739 + case O_MINUS:
1.4740 + /* unary minus */
1.4741 + value = linear_comb(mpl,
1.4742 + 0.0, constant_term(mpl, 0.0),
1.4743 + -1.0, eval_formula(mpl, code->arg.arg.x));
1.4744 + break;
1.4745 + case O_ADD:
1.4746 + /* addition */
1.4747 + value = linear_comb(mpl,
1.4748 + +1.0, eval_formula(mpl, code->arg.arg.x),
1.4749 + +1.0, eval_formula(mpl, code->arg.arg.y));
1.4750 + break;
1.4751 + case O_SUB:
1.4752 + /* subtraction */
1.4753 + value = linear_comb(mpl,
1.4754 + +1.0, eval_formula(mpl, code->arg.arg.x),
1.4755 + -1.0, eval_formula(mpl, code->arg.arg.y));
1.4756 + break;
1.4757 + case O_MUL:
1.4758 + /* multiplication */
1.4759 + xassert(code->arg.arg.x != NULL);
1.4760 + xassert(code->arg.arg.y != NULL);
1.4761 + if (code->arg.arg.x->type == A_NUMERIC)
1.4762 + { xassert(code->arg.arg.y->type == A_FORMULA);
1.4763 + value = linear_comb(mpl,
1.4764 + eval_numeric(mpl, code->arg.arg.x),
1.4765 + eval_formula(mpl, code->arg.arg.y),
1.4766 + 0.0, constant_term(mpl, 0.0));
1.4767 + }
1.4768 + else
1.4769 + { xassert(code->arg.arg.x->type == A_FORMULA);
1.4770 + xassert(code->arg.arg.y->type == A_NUMERIC);
1.4771 + value = linear_comb(mpl,
1.4772 + eval_numeric(mpl, code->arg.arg.y),
1.4773 + eval_formula(mpl, code->arg.arg.x),
1.4774 + 0.0, constant_term(mpl, 0.0));
1.4775 + }
1.4776 + break;
1.4777 + case O_DIV:
1.4778 + /* division */
1.4779 + value = linear_comb(mpl,
1.4780 + fp_div(mpl, 1.0, eval_numeric(mpl, code->arg.arg.y)),
1.4781 + eval_formula(mpl, code->arg.arg.x),
1.4782 + 0.0, constant_term(mpl, 0.0));
1.4783 + break;
1.4784 + case O_FORK:
1.4785 + /* if-then-else */
1.4786 + if (eval_logical(mpl, code->arg.arg.x))
1.4787 + value = eval_formula(mpl, code->arg.arg.y);
1.4788 + else if (code->arg.arg.z == NULL)
1.4789 + value = constant_term(mpl, 0.0);
1.4790 + else
1.4791 + value = eval_formula(mpl, code->arg.arg.z);
1.4792 + break;
1.4793 + case O_SUM:
1.4794 + /* summation over domain */
1.4795 + { struct iter_form_info _info, *info = &_info;
1.4796 + info->code = code;
1.4797 + info->value = constant_term(mpl, 0.0);
1.4798 + info->tail = NULL;
1.4799 + loop_within_domain(mpl, code->arg.loop.domain, info,
1.4800 + iter_form_func);
1.4801 + value = reduce_terms(mpl, info->value);
1.4802 + }
1.4803 + break;
1.4804 + default:
1.4805 + xassert(code != code);
1.4806 + }
1.4807 + /* save resultant value */
1.4808 + xassert(!code->valid);
1.4809 + code->valid = 1;
1.4810 + code->value.form = copy_formula(mpl, value);
1.4811 +done: return value;
1.4812 +}
1.4813 +
1.4814 +/*----------------------------------------------------------------------
1.4815 +-- clean_code - clean pseudo-code.
1.4816 +--
1.4817 +-- This routine recursively cleans specified pseudo-code that assumes
1.4818 +-- deleting all temporary resultant values. */
1.4819 +
1.4820 +void clean_code(MPL *mpl, CODE *code)
1.4821 +{ ARG_LIST *e;
1.4822 + /* if no pseudo-code is specified, do nothing */
1.4823 + if (code == NULL) goto done;
1.4824 + /* if resultant value is valid (exists), delete it */
1.4825 + if (code->valid)
1.4826 + { code->valid = 0;
1.4827 + delete_value(mpl, code->type, &code->value);
1.4828 + }
1.4829 + /* recursively clean pseudo-code for operands */
1.4830 + switch (code->op)
1.4831 + { case O_NUMBER:
1.4832 + case O_STRING:
1.4833 + case O_INDEX:
1.4834 + break;
1.4835 + case O_MEMNUM:
1.4836 + case O_MEMSYM:
1.4837 + for (e = code->arg.par.list; e != NULL; e = e->next)
1.4838 + clean_code(mpl, e->x);
1.4839 + break;
1.4840 + case O_MEMSET:
1.4841 + for (e = code->arg.set.list; e != NULL; e = e->next)
1.4842 + clean_code(mpl, e->x);
1.4843 + break;
1.4844 + case O_MEMVAR:
1.4845 + for (e = code->arg.var.list; e != NULL; e = e->next)
1.4846 + clean_code(mpl, e->x);
1.4847 + break;
1.4848 +#if 1 /* 15/V-2010 */
1.4849 + case O_MEMCON:
1.4850 + for (e = code->arg.con.list; e != NULL; e = e->next)
1.4851 + clean_code(mpl, e->x);
1.4852 + break;
1.4853 +#endif
1.4854 + case O_TUPLE:
1.4855 + case O_MAKE:
1.4856 + for (e = code->arg.list; e != NULL; e = e->next)
1.4857 + clean_code(mpl, e->x);
1.4858 + break;
1.4859 + case O_SLICE:
1.4860 + xassert(code != code);
1.4861 + case O_IRAND224:
1.4862 + case O_UNIFORM01:
1.4863 + case O_NORMAL01:
1.4864 + case O_GMTIME:
1.4865 + break;
1.4866 + case O_CVTNUM:
1.4867 + case O_CVTSYM:
1.4868 + case O_CVTLOG:
1.4869 + case O_CVTTUP:
1.4870 + case O_CVTLFM:
1.4871 + case O_PLUS:
1.4872 + case O_MINUS:
1.4873 + case O_NOT:
1.4874 + case O_ABS:
1.4875 + case O_CEIL:
1.4876 + case O_FLOOR:
1.4877 + case O_EXP:
1.4878 + case O_LOG:
1.4879 + case O_LOG10:
1.4880 + case O_SQRT:
1.4881 + case O_SIN:
1.4882 + case O_COS:
1.4883 + case O_ATAN:
1.4884 + case O_ROUND:
1.4885 + case O_TRUNC:
1.4886 + case O_CARD:
1.4887 + case O_LENGTH:
1.4888 + /* unary operation */
1.4889 + clean_code(mpl, code->arg.arg.x);
1.4890 + break;
1.4891 + case O_ADD:
1.4892 + case O_SUB:
1.4893 + case O_LESS:
1.4894 + case O_MUL:
1.4895 + case O_DIV:
1.4896 + case O_IDIV:
1.4897 + case O_MOD:
1.4898 + case O_POWER:
1.4899 + case O_ATAN2:
1.4900 + case O_ROUND2:
1.4901 + case O_TRUNC2:
1.4902 + case O_UNIFORM:
1.4903 + case O_NORMAL:
1.4904 + case O_CONCAT:
1.4905 + case O_LT:
1.4906 + case O_LE:
1.4907 + case O_EQ:
1.4908 + case O_GE:
1.4909 + case O_GT:
1.4910 + case O_NE:
1.4911 + case O_AND:
1.4912 + case O_OR:
1.4913 + case O_UNION:
1.4914 + case O_DIFF:
1.4915 + case O_SYMDIFF:
1.4916 + case O_INTER:
1.4917 + case O_CROSS:
1.4918 + case O_IN:
1.4919 + case O_NOTIN:
1.4920 + case O_WITHIN:
1.4921 + case O_NOTWITHIN:
1.4922 + case O_SUBSTR:
1.4923 + case O_STR2TIME:
1.4924 + case O_TIME2STR:
1.4925 + /* binary operation */
1.4926 + clean_code(mpl, code->arg.arg.x);
1.4927 + clean_code(mpl, code->arg.arg.y);
1.4928 + break;
1.4929 + case O_DOTS:
1.4930 + case O_FORK:
1.4931 + case O_SUBSTR3:
1.4932 + /* ternary operation */
1.4933 + clean_code(mpl, code->arg.arg.x);
1.4934 + clean_code(mpl, code->arg.arg.y);
1.4935 + clean_code(mpl, code->arg.arg.z);
1.4936 + break;
1.4937 + case O_MIN:
1.4938 + case O_MAX:
1.4939 + /* n-ary operation */
1.4940 + for (e = code->arg.list; e != NULL; e = e->next)
1.4941 + clean_code(mpl, e->x);
1.4942 + break;
1.4943 + case O_SUM:
1.4944 + case O_PROD:
1.4945 + case O_MINIMUM:
1.4946 + case O_MAXIMUM:
1.4947 + case O_FORALL:
1.4948 + case O_EXISTS:
1.4949 + case O_SETOF:
1.4950 + case O_BUILD:
1.4951 + /* iterated operation */
1.4952 + clean_domain(mpl, code->arg.loop.domain);
1.4953 + clean_code(mpl, code->arg.loop.x);
1.4954 + break;
1.4955 + default:
1.4956 + xassert(code->op != code->op);
1.4957 + }
1.4958 +done: return;
1.4959 +}
1.4960 +
1.4961 +#if 1 /* 11/II-2008 */
1.4962 +/**********************************************************************/
1.4963 +/* * * DATA TABLES * * */
1.4964 +/**********************************************************************/
1.4965 +
1.4966 +int mpl_tab_num_args(TABDCA *dca)
1.4967 +{ /* returns the number of arguments */
1.4968 + return dca->na;
1.4969 +}
1.4970 +
1.4971 +const char *mpl_tab_get_arg(TABDCA *dca, int k)
1.4972 +{ /* returns pointer to k-th argument */
1.4973 + xassert(1 <= k && k <= dca->na);
1.4974 + return dca->arg[k];
1.4975 +}
1.4976 +
1.4977 +int mpl_tab_num_flds(TABDCA *dca)
1.4978 +{ /* returns the number of fields */
1.4979 + return dca->nf;
1.4980 +}
1.4981 +
1.4982 +const char *mpl_tab_get_name(TABDCA *dca, int k)
1.4983 +{ /* returns pointer to name of k-th field */
1.4984 + xassert(1 <= k && k <= dca->nf);
1.4985 + return dca->name[k];
1.4986 +}
1.4987 +
1.4988 +int mpl_tab_get_type(TABDCA *dca, int k)
1.4989 +{ /* returns type of k-th field */
1.4990 + xassert(1 <= k && k <= dca->nf);
1.4991 + return dca->type[k];
1.4992 +}
1.4993 +
1.4994 +double mpl_tab_get_num(TABDCA *dca, int k)
1.4995 +{ /* returns numeric value of k-th field */
1.4996 + xassert(1 <= k && k <= dca->nf);
1.4997 + xassert(dca->type[k] == 'N');
1.4998 + return dca->num[k];
1.4999 +}
1.5000 +
1.5001 +const char *mpl_tab_get_str(TABDCA *dca, int k)
1.5002 +{ /* returns pointer to string value of k-th field */
1.5003 + xassert(1 <= k && k <= dca->nf);
1.5004 + xassert(dca->type[k] == 'S');
1.5005 + xassert(dca->str[k] != NULL);
1.5006 + return dca->str[k];
1.5007 +}
1.5008 +
1.5009 +void mpl_tab_set_num(TABDCA *dca, int k, double num)
1.5010 +{ /* assign numeric value to k-th field */
1.5011 + xassert(1 <= k && k <= dca->nf);
1.5012 + xassert(dca->type[k] == '?');
1.5013 + dca->type[k] = 'N';
1.5014 + dca->num[k] = num;
1.5015 + return;
1.5016 +}
1.5017 +
1.5018 +void mpl_tab_set_str(TABDCA *dca, int k, const char *str)
1.5019 +{ /* assign string value to k-th field */
1.5020 + xassert(1 <= k && k <= dca->nf);
1.5021 + xassert(dca->type[k] == '?');
1.5022 + xassert(strlen(str) <= MAX_LENGTH);
1.5023 + xassert(dca->str[k] != NULL);
1.5024 + dca->type[k] = 'S';
1.5025 + strcpy(dca->str[k], str);
1.5026 + return;
1.5027 +}
1.5028 +
1.5029 +static int write_func(MPL *mpl, void *info)
1.5030 +{ /* this is auxiliary routine to work within domain scope */
1.5031 + TABLE *tab = info;
1.5032 + TABDCA *dca = mpl->dca;
1.5033 + TABOUT *out;
1.5034 + SYMBOL *sym;
1.5035 + int k;
1.5036 + char buf[MAX_LENGTH+1];
1.5037 + /* evaluate field values */
1.5038 + k = 0;
1.5039 + for (out = tab->u.out.list; out != NULL; out = out->next)
1.5040 + { k++;
1.5041 + switch (out->code->type)
1.5042 + { case A_NUMERIC:
1.5043 + dca->type[k] = 'N';
1.5044 + dca->num[k] = eval_numeric(mpl, out->code);
1.5045 + dca->str[k][0] = '\0';
1.5046 + break;
1.5047 + case A_SYMBOLIC:
1.5048 + sym = eval_symbolic(mpl, out->code);
1.5049 + if (sym->str == NULL)
1.5050 + { dca->type[k] = 'N';
1.5051 + dca->num[k] = sym->num;
1.5052 + dca->str[k][0] = '\0';
1.5053 + }
1.5054 + else
1.5055 + { dca->type[k] = 'S';
1.5056 + dca->num[k] = 0.0;
1.5057 + fetch_string(mpl, sym->str, buf);
1.5058 + strcpy(dca->str[k], buf);
1.5059 + }
1.5060 + delete_symbol(mpl, sym);
1.5061 + break;
1.5062 + default:
1.5063 + xassert(out != out);
1.5064 + }
1.5065 + }
1.5066 + /* write record to output table */
1.5067 + mpl_tab_drv_write(mpl);
1.5068 + return 0;
1.5069 +}
1.5070 +
1.5071 +void execute_table(MPL *mpl, TABLE *tab)
1.5072 +{ /* execute table statement */
1.5073 + TABARG *arg;
1.5074 + TABFLD *fld;
1.5075 + TABIN *in;
1.5076 + TABOUT *out;
1.5077 + TABDCA *dca;
1.5078 + SET *set;
1.5079 + int k;
1.5080 + char buf[MAX_LENGTH+1];
1.5081 + /* allocate table driver communication area */
1.5082 + xassert(mpl->dca == NULL);
1.5083 + mpl->dca = dca = xmalloc(sizeof(TABDCA));
1.5084 + dca->id = 0;
1.5085 + dca->link = NULL;
1.5086 + dca->na = 0;
1.5087 + dca->arg = NULL;
1.5088 + dca->nf = 0;
1.5089 + dca->name = NULL;
1.5090 + dca->type = NULL;
1.5091 + dca->num = NULL;
1.5092 + dca->str = NULL;
1.5093 + /* allocate arguments */
1.5094 + xassert(dca->na == 0);
1.5095 + for (arg = tab->arg; arg != NULL; arg = arg->next)
1.5096 + dca->na++;
1.5097 + dca->arg = xcalloc(1+dca->na, sizeof(char *));
1.5098 +#if 1 /* 28/IX-2008 */
1.5099 + for (k = 1; k <= dca->na; k++) dca->arg[k] = NULL;
1.5100 +#endif
1.5101 + /* evaluate argument values */
1.5102 + k = 0;
1.5103 + for (arg = tab->arg; arg != NULL; arg = arg->next)
1.5104 + { SYMBOL *sym;
1.5105 + k++;
1.5106 + xassert(arg->code->type == A_SYMBOLIC);
1.5107 + sym = eval_symbolic(mpl, arg->code);
1.5108 + if (sym->str == NULL)
1.5109 + sprintf(buf, "%.*g", DBL_DIG, sym->num);
1.5110 + else
1.5111 + fetch_string(mpl, sym->str, buf);
1.5112 + delete_symbol(mpl, sym);
1.5113 + dca->arg[k] = xmalloc(strlen(buf)+1);
1.5114 + strcpy(dca->arg[k], buf);
1.5115 + }
1.5116 + /* perform table input/output */
1.5117 + switch (tab->type)
1.5118 + { case A_INPUT: goto read_table;
1.5119 + case A_OUTPUT: goto write_table;
1.5120 + default: xassert(tab != tab);
1.5121 + }
1.5122 +read_table:
1.5123 + /* read data from input table */
1.5124 + /* add the only member to the control set and assign it empty
1.5125 + elemental set */
1.5126 + set = tab->u.in.set;
1.5127 + if (set != NULL)
1.5128 + { if (set->data)
1.5129 + error(mpl, "%s already provided with data", set->name);
1.5130 + xassert(set->array->head == NULL);
1.5131 + add_member(mpl, set->array, NULL)->value.set =
1.5132 + create_elemset(mpl, set->dimen);
1.5133 + set->data = 1;
1.5134 + }
1.5135 + /* check parameters specified in the input list */
1.5136 + for (in = tab->u.in.list; in != NULL; in = in->next)
1.5137 + { if (in->par->data)
1.5138 + error(mpl, "%s already provided with data", in->par->name);
1.5139 + in->par->data = 1;
1.5140 + }
1.5141 + /* allocate and initialize fields */
1.5142 + xassert(dca->nf == 0);
1.5143 + for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
1.5144 + dca->nf++;
1.5145 + for (in = tab->u.in.list; in != NULL; in = in->next)
1.5146 + dca->nf++;
1.5147 + dca->name = xcalloc(1+dca->nf, sizeof(char *));
1.5148 + dca->type = xcalloc(1+dca->nf, sizeof(int));
1.5149 + dca->num = xcalloc(1+dca->nf, sizeof(double));
1.5150 + dca->str = xcalloc(1+dca->nf, sizeof(char *));
1.5151 + k = 0;
1.5152 + for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
1.5153 + { k++;
1.5154 + dca->name[k] = fld->name;
1.5155 + dca->type[k] = '?';
1.5156 + dca->num[k] = 0.0;
1.5157 + dca->str[k] = xmalloc(MAX_LENGTH+1);
1.5158 + dca->str[k][0] = '\0';
1.5159 + }
1.5160 + for (in = tab->u.in.list; in != NULL; in = in->next)
1.5161 + { k++;
1.5162 + dca->name[k] = in->name;
1.5163 + dca->type[k] = '?';
1.5164 + dca->num[k] = 0.0;
1.5165 + dca->str[k] = xmalloc(MAX_LENGTH+1);
1.5166 + dca->str[k][0] = '\0';
1.5167 + }
1.5168 + /* open input table */
1.5169 + mpl_tab_drv_open(mpl, 'R');
1.5170 + /* read and process records */
1.5171 + for (;;)
1.5172 + { TUPLE *tup;
1.5173 + /* reset field types */
1.5174 + for (k = 1; k <= dca->nf; k++)
1.5175 + dca->type[k] = '?';
1.5176 + /* read next record */
1.5177 + if (mpl_tab_drv_read(mpl)) break;
1.5178 + /* all fields must be set by the driver */
1.5179 + for (k = 1; k <= dca->nf; k++)
1.5180 + { if (dca->type[k] == '?')
1.5181 + error(mpl, "field %s missing in input table",
1.5182 + dca->name[k]);
1.5183 + }
1.5184 + /* construct n-tuple */
1.5185 + tup = create_tuple(mpl);
1.5186 + k = 0;
1.5187 + for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
1.5188 + { k++;
1.5189 + xassert(k <= dca->nf);
1.5190 + switch (dca->type[k])
1.5191 + { case 'N':
1.5192 + tup = expand_tuple(mpl, tup, create_symbol_num(mpl,
1.5193 + dca->num[k]));
1.5194 + break;
1.5195 + case 'S':
1.5196 + xassert(strlen(dca->str[k]) <= MAX_LENGTH);
1.5197 + tup = expand_tuple(mpl, tup, create_symbol_str(mpl,
1.5198 + create_string(mpl, dca->str[k])));
1.5199 + break;
1.5200 + default:
1.5201 + xassert(dca != dca);
1.5202 + }
1.5203 + }
1.5204 + /* add n-tuple just read to the control set */
1.5205 + if (tab->u.in.set != NULL)
1.5206 + check_then_add(mpl, tab->u.in.set->array->head->value.set,
1.5207 + copy_tuple(mpl, tup));
1.5208 + /* assign values to the parameters in the input list */
1.5209 + for (in = tab->u.in.list; in != NULL; in = in->next)
1.5210 + { MEMBER *memb;
1.5211 + k++;
1.5212 + xassert(k <= dca->nf);
1.5213 + /* there must be no member with the same n-tuple */
1.5214 + if (find_member(mpl, in->par->array, tup) != NULL)
1.5215 + error(mpl, "%s%s already defined", in->par->name,
1.5216 + format_tuple(mpl, '[', tup));
1.5217 + /* create new parameter member with given n-tuple */
1.5218 + memb = add_member(mpl, in->par->array, copy_tuple(mpl, tup))
1.5219 + ;
1.5220 + /* assign value to the parameter member */
1.5221 + switch (in->par->type)
1.5222 + { case A_NUMERIC:
1.5223 + case A_INTEGER:
1.5224 + case A_BINARY:
1.5225 + if (dca->type[k] != 'N')
1.5226 + error(mpl, "%s requires numeric data",
1.5227 + in->par->name);
1.5228 + memb->value.num = dca->num[k];
1.5229 + break;
1.5230 + case A_SYMBOLIC:
1.5231 + switch (dca->type[k])
1.5232 + { case 'N':
1.5233 + memb->value.sym = create_symbol_num(mpl,
1.5234 + dca->num[k]);
1.5235 + break;
1.5236 + case 'S':
1.5237 + xassert(strlen(dca->str[k]) <= MAX_LENGTH);
1.5238 + memb->value.sym = create_symbol_str(mpl,
1.5239 + create_string(mpl,dca->str[k]));
1.5240 + break;
1.5241 + default:
1.5242 + xassert(dca != dca);
1.5243 + }
1.5244 + break;
1.5245 + default:
1.5246 + xassert(in != in);
1.5247 + }
1.5248 + }
1.5249 + /* n-tuple is no more needed */
1.5250 + delete_tuple(mpl, tup);
1.5251 + }
1.5252 + /* close input table */
1.5253 + mpl_tab_drv_close(mpl);
1.5254 + goto done;
1.5255 +write_table:
1.5256 + /* write data to output table */
1.5257 + /* allocate and initialize fields */
1.5258 + xassert(dca->nf == 0);
1.5259 + for (out = tab->u.out.list; out != NULL; out = out->next)
1.5260 + dca->nf++;
1.5261 + dca->name = xcalloc(1+dca->nf, sizeof(char *));
1.5262 + dca->type = xcalloc(1+dca->nf, sizeof(int));
1.5263 + dca->num = xcalloc(1+dca->nf, sizeof(double));
1.5264 + dca->str = xcalloc(1+dca->nf, sizeof(char *));
1.5265 + k = 0;
1.5266 + for (out = tab->u.out.list; out != NULL; out = out->next)
1.5267 + { k++;
1.5268 + dca->name[k] = out->name;
1.5269 + dca->type[k] = '?';
1.5270 + dca->num[k] = 0.0;
1.5271 + dca->str[k] = xmalloc(MAX_LENGTH+1);
1.5272 + dca->str[k][0] = '\0';
1.5273 + }
1.5274 + /* open output table */
1.5275 + mpl_tab_drv_open(mpl, 'W');
1.5276 + /* evaluate fields and write records */
1.5277 + loop_within_domain(mpl, tab->u.out.domain, tab, write_func);
1.5278 + /* close output table */
1.5279 + mpl_tab_drv_close(mpl);
1.5280 +done: /* free table driver communication area */
1.5281 + free_dca(mpl);
1.5282 + return;
1.5283 +}
1.5284 +
1.5285 +void free_dca(MPL *mpl)
1.5286 +{ /* free table driver communucation area */
1.5287 + TABDCA *dca = mpl->dca;
1.5288 + int k;
1.5289 + if (dca != NULL)
1.5290 + { if (dca->link != NULL)
1.5291 + mpl_tab_drv_close(mpl);
1.5292 + if (dca->arg != NULL)
1.5293 + { for (k = 1; k <= dca->na; k++)
1.5294 +#if 1 /* 28/IX-2008 */
1.5295 + if (dca->arg[k] != NULL)
1.5296 +#endif
1.5297 + xfree(dca->arg[k]);
1.5298 + xfree(dca->arg);
1.5299 + }
1.5300 + if (dca->name != NULL) xfree(dca->name);
1.5301 + if (dca->type != NULL) xfree(dca->type);
1.5302 + if (dca->num != NULL) xfree(dca->num);
1.5303 + if (dca->str != NULL)
1.5304 + { for (k = 1; k <= dca->nf; k++)
1.5305 + xfree(dca->str[k]);
1.5306 + xfree(dca->str);
1.5307 + }
1.5308 + xfree(dca), mpl->dca = NULL;
1.5309 + }
1.5310 + return;
1.5311 +}
1.5312 +
1.5313 +void clean_table(MPL *mpl, TABLE *tab)
1.5314 +{ /* clean table statement */
1.5315 + TABARG *arg;
1.5316 + TABOUT *out;
1.5317 + /* clean string list */
1.5318 + for (arg = tab->arg; arg != NULL; arg = arg->next)
1.5319 + clean_code(mpl, arg->code);
1.5320 + switch (tab->type)
1.5321 + { case A_INPUT:
1.5322 + break;
1.5323 + case A_OUTPUT:
1.5324 + /* clean subscript domain */
1.5325 + clean_domain(mpl, tab->u.out.domain);
1.5326 + /* clean output list */
1.5327 + for (out = tab->u.out.list; out != NULL; out = out->next)
1.5328 + clean_code(mpl, out->code);
1.5329 + break;
1.5330 + default:
1.5331 + xassert(tab != tab);
1.5332 + }
1.5333 + return;
1.5334 +}
1.5335 +#endif
1.5336 +
1.5337 +/**********************************************************************/
1.5338 +/* * * MODEL STATEMENTS * * */
1.5339 +/**********************************************************************/
1.5340 +
1.5341 +/*----------------------------------------------------------------------
1.5342 +-- execute_check - execute check statement.
1.5343 +--
1.5344 +-- This routine executes specified check statement. */
1.5345 +
1.5346 +static int check_func(MPL *mpl, void *info)
1.5347 +{ /* this is auxiliary routine to work within domain scope */
1.5348 + CHECK *chk = (CHECK *)info;
1.5349 + if (!eval_logical(mpl, chk->code))
1.5350 + error(mpl, "check%s failed", format_tuple(mpl, '[',
1.5351 + get_domain_tuple(mpl, chk->domain)));
1.5352 + return 0;
1.5353 +}
1.5354 +
1.5355 +void execute_check(MPL *mpl, CHECK *chk)
1.5356 +{ loop_within_domain(mpl, chk->domain, chk, check_func);
1.5357 + return;
1.5358 +}
1.5359 +
1.5360 +/*----------------------------------------------------------------------
1.5361 +-- clean_check - clean check statement.
1.5362 +--
1.5363 +-- This routine cleans specified check statement that assumes deleting
1.5364 +-- all stuff dynamically allocated on generating/postsolving phase. */
1.5365 +
1.5366 +void clean_check(MPL *mpl, CHECK *chk)
1.5367 +{ /* clean subscript domain */
1.5368 + clean_domain(mpl, chk->domain);
1.5369 + /* clean pseudo-code for computing predicate */
1.5370 + clean_code(mpl, chk->code);
1.5371 + return;
1.5372 +}
1.5373 +
1.5374 +/*----------------------------------------------------------------------
1.5375 +-- execute_display - execute display statement.
1.5376 +--
1.5377 +-- This routine executes specified display statement. */
1.5378 +
1.5379 +static void display_set(MPL *mpl, SET *set, MEMBER *memb)
1.5380 +{ /* display member of model set */
1.5381 + ELEMSET *s = memb->value.set;
1.5382 + MEMBER *m;
1.5383 + write_text(mpl, "%s%s%s\n", set->name,
1.5384 + format_tuple(mpl, '[', memb->tuple),
1.5385 + s->head == NULL ? " is empty" : ":");
1.5386 + for (m = s->head; m != NULL; m = m->next)
1.5387 + write_text(mpl, " %s\n", format_tuple(mpl, '(', m->tuple));
1.5388 + return;
1.5389 +}
1.5390 +
1.5391 +static void display_par(MPL *mpl, PARAMETER *par, MEMBER *memb)
1.5392 +{ /* display member of model parameter */
1.5393 + switch (par->type)
1.5394 + { case A_NUMERIC:
1.5395 + case A_INTEGER:
1.5396 + case A_BINARY:
1.5397 + write_text(mpl, "%s%s = %.*g\n", par->name,
1.5398 + format_tuple(mpl, '[', memb->tuple),
1.5399 + DBL_DIG, memb->value.num);
1.5400 + break;
1.5401 + case A_SYMBOLIC:
1.5402 + write_text(mpl, "%s%s = %s\n", par->name,
1.5403 + format_tuple(mpl, '[', memb->tuple),
1.5404 + format_symbol(mpl, memb->value.sym));
1.5405 + break;
1.5406 + default:
1.5407 + xassert(par != par);
1.5408 + }
1.5409 + return;
1.5410 +}
1.5411 +
1.5412 +#if 1 /* 15/V-2010 */
1.5413 +static void display_var(MPL *mpl, VARIABLE *var, MEMBER *memb,
1.5414 + int suff)
1.5415 +{ /* display member of model variable */
1.5416 + if (suff == DOT_NONE || suff == DOT_VAL)
1.5417 + write_text(mpl, "%s%s.val = %.*g\n", var->name,
1.5418 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5419 + memb->value.var->prim);
1.5420 + else if (suff == DOT_LB)
1.5421 + write_text(mpl, "%s%s.lb = %.*g\n", var->name,
1.5422 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5423 + memb->value.var->var->lbnd == NULL ? -DBL_MAX :
1.5424 + memb->value.var->lbnd);
1.5425 + else if (suff == DOT_UB)
1.5426 + write_text(mpl, "%s%s.ub = %.*g\n", var->name,
1.5427 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5428 + memb->value.var->var->ubnd == NULL ? +DBL_MAX :
1.5429 + memb->value.var->ubnd);
1.5430 + else if (suff == DOT_STATUS)
1.5431 + write_text(mpl, "%s%s.status = %d\n", var->name, format_tuple
1.5432 + (mpl, '[', memb->tuple), memb->value.var->stat);
1.5433 + else if (suff == DOT_DUAL)
1.5434 + write_text(mpl, "%s%s.dual = %.*g\n", var->name,
1.5435 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5436 + memb->value.var->dual);
1.5437 + else
1.5438 + xassert(suff != suff);
1.5439 + return;
1.5440 +}
1.5441 +#endif
1.5442 +
1.5443 +#if 1 /* 15/V-2010 */
1.5444 +static void display_con(MPL *mpl, CONSTRAINT *con, MEMBER *memb,
1.5445 + int suff)
1.5446 +{ /* display member of model constraint */
1.5447 + if (suff == DOT_NONE || suff == DOT_VAL)
1.5448 + write_text(mpl, "%s%s.val = %.*g\n", con->name,
1.5449 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5450 + memb->value.con->prim);
1.5451 + else if (suff == DOT_LB)
1.5452 + write_text(mpl, "%s%s.lb = %.*g\n", con->name,
1.5453 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5454 + memb->value.con->con->lbnd == NULL ? -DBL_MAX :
1.5455 + memb->value.con->lbnd);
1.5456 + else if (suff == DOT_UB)
1.5457 + write_text(mpl, "%s%s.ub = %.*g\n", con->name,
1.5458 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5459 + memb->value.con->con->ubnd == NULL ? +DBL_MAX :
1.5460 + memb->value.con->ubnd);
1.5461 + else if (suff == DOT_STATUS)
1.5462 + write_text(mpl, "%s%s.status = %d\n", con->name, format_tuple
1.5463 + (mpl, '[', memb->tuple), memb->value.con->stat);
1.5464 + else if (suff == DOT_DUAL)
1.5465 + write_text(mpl, "%s%s.dual = %.*g\n", con->name,
1.5466 + format_tuple(mpl, '[', memb->tuple), DBL_DIG,
1.5467 + memb->value.con->dual);
1.5468 + else
1.5469 + xassert(suff != suff);
1.5470 + return;
1.5471 +}
1.5472 +#endif
1.5473 +
1.5474 +static void display_memb(MPL *mpl, CODE *code)
1.5475 +{ /* display member specified by pseudo-code */
1.5476 + MEMBER memb;
1.5477 + ARG_LIST *e;
1.5478 + xassert(code->op == O_MEMNUM || code->op == O_MEMSYM
1.5479 + || code->op == O_MEMSET || code->op == O_MEMVAR
1.5480 + || code->op == O_MEMCON);
1.5481 + memb.tuple = create_tuple(mpl);
1.5482 + for (e = code->arg.par.list; e != NULL; e = e->next)
1.5483 + memb.tuple = expand_tuple(mpl, memb.tuple, eval_symbolic(mpl,
1.5484 + e->x));
1.5485 + switch (code->op)
1.5486 + { case O_MEMNUM:
1.5487 + memb.value.num = eval_member_num(mpl, code->arg.par.par,
1.5488 + memb.tuple);
1.5489 + display_par(mpl, code->arg.par.par, &memb);
1.5490 + break;
1.5491 + case O_MEMSYM:
1.5492 + memb.value.sym = eval_member_sym(mpl, code->arg.par.par,
1.5493 + memb.tuple);
1.5494 + display_par(mpl, code->arg.par.par, &memb);
1.5495 + delete_symbol(mpl, memb.value.sym);
1.5496 + break;
1.5497 + case O_MEMSET:
1.5498 + memb.value.set = eval_member_set(mpl, code->arg.set.set,
1.5499 + memb.tuple);
1.5500 + display_set(mpl, code->arg.set.set, &memb);
1.5501 + break;
1.5502 + case O_MEMVAR:
1.5503 + memb.value.var = eval_member_var(mpl, code->arg.var.var,
1.5504 + memb.tuple);
1.5505 + display_var
1.5506 + (mpl, code->arg.var.var, &memb, code->arg.var.suff);
1.5507 + break;
1.5508 + case O_MEMCON:
1.5509 + memb.value.con = eval_member_con(mpl, code->arg.con.con,
1.5510 + memb.tuple);
1.5511 + display_con
1.5512 + (mpl, code->arg.con.con, &memb, code->arg.con.suff);
1.5513 + break;
1.5514 + default:
1.5515 + xassert(code != code);
1.5516 + }
1.5517 + delete_tuple(mpl, memb.tuple);
1.5518 + return;
1.5519 +}
1.5520 +
1.5521 +static void display_code(MPL *mpl, CODE *code)
1.5522 +{ /* display value of expression */
1.5523 + switch (code->type)
1.5524 + { case A_NUMERIC:
1.5525 + /* numeric value */
1.5526 + { double num;
1.5527 + num = eval_numeric(mpl, code);
1.5528 + write_text(mpl, "%.*g\n", DBL_DIG, num);
1.5529 + }
1.5530 + break;
1.5531 + case A_SYMBOLIC:
1.5532 + /* symbolic value */
1.5533 + { SYMBOL *sym;
1.5534 + sym = eval_symbolic(mpl, code);
1.5535 + write_text(mpl, "%s\n", format_symbol(mpl, sym));
1.5536 + delete_symbol(mpl, sym);
1.5537 + }
1.5538 + break;
1.5539 + case A_LOGICAL:
1.5540 + /* logical value */
1.5541 + { int bit;
1.5542 + bit = eval_logical(mpl, code);
1.5543 + write_text(mpl, "%s\n", bit ? "true" : "false");
1.5544 + }
1.5545 + break;
1.5546 + case A_TUPLE:
1.5547 + /* n-tuple */
1.5548 + { TUPLE *tuple;
1.5549 + tuple = eval_tuple(mpl, code);
1.5550 + write_text(mpl, "%s\n", format_tuple(mpl, '(', tuple));
1.5551 + delete_tuple(mpl, tuple);
1.5552 + }
1.5553 + break;
1.5554 + case A_ELEMSET:
1.5555 + /* elemental set */
1.5556 + { ELEMSET *set;
1.5557 + MEMBER *memb;
1.5558 + set = eval_elemset(mpl, code);
1.5559 + if (set->head == 0)
1.5560 + write_text(mpl, "set is empty\n");
1.5561 + for (memb = set->head; memb != NULL; memb = memb->next)
1.5562 + write_text(mpl, " %s\n", format_tuple(mpl, '(',
1.5563 + memb->tuple));
1.5564 + delete_elemset(mpl, set);
1.5565 + }
1.5566 + break;
1.5567 + case A_FORMULA:
1.5568 + /* linear form */
1.5569 + { FORMULA *form, *term;
1.5570 + form = eval_formula(mpl, code);
1.5571 + if (form == NULL)
1.5572 + write_text(mpl, "linear form is empty\n");
1.5573 + for (term = form; term != NULL; term = term->next)
1.5574 + { if (term->var == NULL)
1.5575 + write_text(mpl, " %.*g\n", term->coef);
1.5576 + else
1.5577 + write_text(mpl, " %.*g %s%s\n", DBL_DIG,
1.5578 + term->coef, term->var->var->name,
1.5579 + format_tuple(mpl, '[', term->var->memb->tuple));
1.5580 + }
1.5581 + delete_formula(mpl, form);
1.5582 + }
1.5583 + break;
1.5584 + default:
1.5585 + xassert(code != code);
1.5586 + }
1.5587 + return;
1.5588 +}
1.5589 +
1.5590 +static int display_func(MPL *mpl, void *info)
1.5591 +{ /* this is auxiliary routine to work within domain scope */
1.5592 + DISPLAY *dpy = (DISPLAY *)info;
1.5593 + DISPLAY1 *entry;
1.5594 + for (entry = dpy->list; entry != NULL; entry = entry->next)
1.5595 + { if (entry->type == A_INDEX)
1.5596 + { /* dummy index */
1.5597 + DOMAIN_SLOT *slot = entry->u.slot;
1.5598 + write_text(mpl, "%s = %s\n", slot->name,
1.5599 + format_symbol(mpl, slot->value));
1.5600 + }
1.5601 + else if (entry->type == A_SET)
1.5602 + { /* model set */
1.5603 + SET *set = entry->u.set;
1.5604 + MEMBER *memb;
1.5605 + if (set->assign != NULL)
1.5606 + { /* the set has assignment expression; evaluate all its
1.5607 + members over entire domain */
1.5608 + eval_whole_set(mpl, set);
1.5609 + }
1.5610 + else
1.5611 + { /* the set has no assignment expression; refer to its
1.5612 + any existing member ignoring resultant value to check
1.5613 + the data provided the data section */
1.5614 +#if 1 /* 12/XII-2008 */
1.5615 + if (set->gadget != NULL && set->data == 0)
1.5616 + { /* initialize the set with data from a plain set */
1.5617 + saturate_set(mpl, set);
1.5618 + }
1.5619 +#endif
1.5620 + if (set->array->head != NULL)
1.5621 + eval_member_set(mpl, set, set->array->head->tuple);
1.5622 + }
1.5623 + /* display all members of the set array */
1.5624 + if (set->array->head == NULL)
1.5625 + write_text(mpl, "%s has empty content\n", set->name);
1.5626 + for (memb = set->array->head; memb != NULL; memb =
1.5627 + memb->next) display_set(mpl, set, memb);
1.5628 + }
1.5629 + else if (entry->type == A_PARAMETER)
1.5630 + { /* model parameter */
1.5631 + PARAMETER *par = entry->u.par;
1.5632 + MEMBER *memb;
1.5633 + if (par->assign != NULL)
1.5634 + { /* the parameter has an assignment expression; evaluate
1.5635 + all its member over entire domain */
1.5636 + eval_whole_par(mpl, par);
1.5637 + }
1.5638 + else
1.5639 + { /* the parameter has no assignment expression; refer to
1.5640 + its any existing member ignoring resultant value to
1.5641 + check the data provided in the data section */
1.5642 + if (par->array->head != NULL)
1.5643 + { if (par->type != A_SYMBOLIC)
1.5644 + eval_member_num(mpl, par, par->array->head->tuple);
1.5645 + else
1.5646 + delete_symbol(mpl, eval_member_sym(mpl, par,
1.5647 + par->array->head->tuple));
1.5648 + }
1.5649 + }
1.5650 + /* display all members of the parameter array */
1.5651 + if (par->array->head == NULL)
1.5652 + write_text(mpl, "%s has empty content\n", par->name);
1.5653 + for (memb = par->array->head; memb != NULL; memb =
1.5654 + memb->next) display_par(mpl, par, memb);
1.5655 + }
1.5656 + else if (entry->type == A_VARIABLE)
1.5657 + { /* model variable */
1.5658 + VARIABLE *var = entry->u.var;
1.5659 + MEMBER *memb;
1.5660 + xassert(mpl->flag_p);
1.5661 + /* display all members of the variable array */
1.5662 + if (var->array->head == NULL)
1.5663 + write_text(mpl, "%s has empty content\n", var->name);
1.5664 + for (memb = var->array->head; memb != NULL; memb =
1.5665 + memb->next) display_var(mpl, var, memb, DOT_NONE);
1.5666 + }
1.5667 + else if (entry->type == A_CONSTRAINT)
1.5668 + { /* model constraint */
1.5669 + CONSTRAINT *con = entry->u.con;
1.5670 + MEMBER *memb;
1.5671 + xassert(mpl->flag_p);
1.5672 + /* display all members of the constraint array */
1.5673 + if (con->array->head == NULL)
1.5674 + write_text(mpl, "%s has empty content\n", con->name);
1.5675 + for (memb = con->array->head; memb != NULL; memb =
1.5676 + memb->next) display_con(mpl, con, memb, DOT_NONE);
1.5677 + }
1.5678 + else if (entry->type == A_EXPRESSION)
1.5679 + { /* expression */
1.5680 + CODE *code = entry->u.code;
1.5681 + if (code->op == O_MEMNUM || code->op == O_MEMSYM ||
1.5682 + code->op == O_MEMSET || code->op == O_MEMVAR ||
1.5683 + code->op == O_MEMCON)
1.5684 + display_memb(mpl, code);
1.5685 + else
1.5686 + display_code(mpl, code);
1.5687 + }
1.5688 + else
1.5689 + xassert(entry != entry);
1.5690 + }
1.5691 + return 0;
1.5692 +}
1.5693 +
1.5694 +void execute_display(MPL *mpl, DISPLAY *dpy)
1.5695 +{ loop_within_domain(mpl, dpy->domain, dpy, display_func);
1.5696 + return;
1.5697 +}
1.5698 +
1.5699 +/*----------------------------------------------------------------------
1.5700 +-- clean_display - clean display statement.
1.5701 +--
1.5702 +-- This routine cleans specified display statement that assumes deleting
1.5703 +-- all stuff dynamically allocated on generating/postsolving phase. */
1.5704 +
1.5705 +void clean_display(MPL *mpl, DISPLAY *dpy)
1.5706 +{ DISPLAY1 *d;
1.5707 +#if 0 /* 15/V-2010 */
1.5708 + ARG_LIST *e;
1.5709 +#endif
1.5710 + /* clean subscript domain */
1.5711 + clean_domain(mpl, dpy->domain);
1.5712 + /* clean display list */
1.5713 + for (d = dpy->list; d != NULL; d = d->next)
1.5714 + { /* clean pseudo-code for computing expression */
1.5715 + if (d->type == A_EXPRESSION)
1.5716 + clean_code(mpl, d->u.code);
1.5717 +#if 0 /* 15/V-2010 */
1.5718 + /* clean pseudo-code for computing subscripts */
1.5719 + for (e = d->list; e != NULL; e = e->next)
1.5720 + clean_code(mpl, e->x);
1.5721 +#endif
1.5722 + }
1.5723 + return;
1.5724 +}
1.5725 +
1.5726 +/*----------------------------------------------------------------------
1.5727 +-- execute_printf - execute printf statement.
1.5728 +--
1.5729 +-- This routine executes specified printf statement. */
1.5730 +
1.5731 +#if 1 /* 14/VII-2006 */
1.5732 +static void print_char(MPL *mpl, int c)
1.5733 +{ if (mpl->prt_fp == NULL)
1.5734 + write_char(mpl, c);
1.5735 + else
1.5736 + xfputc(c, mpl->prt_fp);
1.5737 + return;
1.5738 +}
1.5739 +
1.5740 +static void print_text(MPL *mpl, char *fmt, ...)
1.5741 +{ va_list arg;
1.5742 + char buf[OUTBUF_SIZE], *c;
1.5743 + va_start(arg, fmt);
1.5744 + vsprintf(buf, fmt, arg);
1.5745 + xassert(strlen(buf) < sizeof(buf));
1.5746 + va_end(arg);
1.5747 + for (c = buf; *c != '\0'; c++) print_char(mpl, *c);
1.5748 + return;
1.5749 +}
1.5750 +#endif
1.5751 +
1.5752 +static int printf_func(MPL *mpl, void *info)
1.5753 +{ /* this is auxiliary routine to work within domain scope */
1.5754 + PRINTF *prt = (PRINTF *)info;
1.5755 + PRINTF1 *entry;
1.5756 + SYMBOL *sym;
1.5757 + char fmt[MAX_LENGTH+1], *c, *from, save;
1.5758 + /* evaluate format control string */
1.5759 + sym = eval_symbolic(mpl, prt->fmt);
1.5760 + if (sym->str == NULL)
1.5761 + sprintf(fmt, "%.*g", DBL_DIG, sym->num);
1.5762 + else
1.5763 + fetch_string(mpl, sym->str, fmt);
1.5764 + delete_symbol(mpl, sym);
1.5765 + /* scan format control string and perform formatting output */
1.5766 + entry = prt->list;
1.5767 + for (c = fmt; *c != '\0'; c++)
1.5768 + { if (*c == '%')
1.5769 + { /* scan format specifier */
1.5770 + from = c++;
1.5771 + if (*c == '%')
1.5772 + { print_char(mpl, '%');
1.5773 + continue;
1.5774 + }
1.5775 + if (entry == NULL) break;
1.5776 + /* scan optional flags */
1.5777 + while (*c == '-' || *c == '+' || *c == ' ' || *c == '#' ||
1.5778 + *c == '0') c++;
1.5779 + /* scan optional minimum field width */
1.5780 + while (isdigit((unsigned char)*c)) c++;
1.5781 + /* scan optional precision */
1.5782 + if (*c == '.')
1.5783 + { c++;
1.5784 + while (isdigit((unsigned char)*c)) c++;
1.5785 + }
1.5786 + /* scan conversion specifier and perform formatting */
1.5787 + save = *(c+1), *(c+1) = '\0';
1.5788 + if (*c == 'd' || *c == 'i' || *c == 'e' || *c == 'E' ||
1.5789 + *c == 'f' || *c == 'F' || *c == 'g' || *c == 'G')
1.5790 + { /* the specifier requires numeric value */
1.5791 + double value;
1.5792 + xassert(entry != NULL);
1.5793 + switch (entry->code->type)
1.5794 + { case A_NUMERIC:
1.5795 + value = eval_numeric(mpl, entry->code);
1.5796 + break;
1.5797 + case A_SYMBOLIC:
1.5798 + sym = eval_symbolic(mpl, entry->code);
1.5799 + if (sym->str != NULL)
1.5800 + error(mpl, "cannot convert %s to floating-point"
1.5801 + " number", format_symbol(mpl, sym));
1.5802 + value = sym->num;
1.5803 + delete_symbol(mpl, sym);
1.5804 + break;
1.5805 + case A_LOGICAL:
1.5806 + if (eval_logical(mpl, entry->code))
1.5807 + value = 1.0;
1.5808 + else
1.5809 + value = 0.0;
1.5810 + break;
1.5811 + default:
1.5812 + xassert(entry != entry);
1.5813 + }
1.5814 + if (*c == 'd' || *c == 'i')
1.5815 + { double int_max = (double)INT_MAX;
1.5816 + if (!(-int_max <= value && value <= +int_max))
1.5817 + error(mpl, "cannot convert %.*g to integer",
1.5818 + DBL_DIG, value);
1.5819 + print_text(mpl, from, (int)floor(value + 0.5));
1.5820 + }
1.5821 + else
1.5822 + print_text(mpl, from, value);
1.5823 + }
1.5824 + else if (*c == 's')
1.5825 + { /* the specifier requires symbolic value */
1.5826 + char value[MAX_LENGTH+1];
1.5827 + switch (entry->code->type)
1.5828 + { case A_NUMERIC:
1.5829 + sprintf(value, "%.*g", DBL_DIG, eval_numeric(mpl,
1.5830 + entry->code));
1.5831 + break;
1.5832 + case A_LOGICAL:
1.5833 + if (eval_logical(mpl, entry->code))
1.5834 + strcpy(value, "T");
1.5835 + else
1.5836 + strcpy(value, "F");
1.5837 + break;
1.5838 + case A_SYMBOLIC:
1.5839 + sym = eval_symbolic(mpl, entry->code);
1.5840 + if (sym->str == NULL)
1.5841 + sprintf(value, "%.*g", DBL_DIG, sym->num);
1.5842 + else
1.5843 + fetch_string(mpl, sym->str, value);
1.5844 + delete_symbol(mpl, sym);
1.5845 + break;
1.5846 + default:
1.5847 + xassert(entry != entry);
1.5848 + }
1.5849 + print_text(mpl, from, value);
1.5850 + }
1.5851 + else
1.5852 + error(mpl, "format specifier missing or invalid");
1.5853 + *(c+1) = save;
1.5854 + entry = entry->next;
1.5855 + }
1.5856 + else if (*c == '\\')
1.5857 + { /* write some control character */
1.5858 + c++;
1.5859 + if (*c == 't')
1.5860 + print_char(mpl, '\t');
1.5861 + else if (*c == 'n')
1.5862 + print_char(mpl, '\n');
1.5863 +#if 1 /* 28/X-2010 */
1.5864 + else if (*c == '\0')
1.5865 + { /* format string ends with backslash */
1.5866 + error(mpl, "invalid use of escape character \\ in format"
1.5867 + " control string");
1.5868 + }
1.5869 +#endif
1.5870 + else
1.5871 + print_char(mpl, *c);
1.5872 + }
1.5873 + else
1.5874 + { /* write character without formatting */
1.5875 + print_char(mpl, *c);
1.5876 + }
1.5877 + }
1.5878 + return 0;
1.5879 +}
1.5880 +
1.5881 +#if 0 /* 14/VII-2006 */
1.5882 +void execute_printf(MPL *mpl, PRINTF *prt)
1.5883 +{ loop_within_domain(mpl, prt->domain, prt, printf_func);
1.5884 + return;
1.5885 +}
1.5886 +#else
1.5887 +void execute_printf(MPL *mpl, PRINTF *prt)
1.5888 +{ if (prt->fname == NULL)
1.5889 + { /* switch to the standard output */
1.5890 + if (mpl->prt_fp != NULL)
1.5891 + { xfclose(mpl->prt_fp), mpl->prt_fp = NULL;
1.5892 + xfree(mpl->prt_file), mpl->prt_file = NULL;
1.5893 + }
1.5894 + }
1.5895 + else
1.5896 + { /* evaluate file name string */
1.5897 + SYMBOL *sym;
1.5898 + char fname[MAX_LENGTH+1];
1.5899 + sym = eval_symbolic(mpl, prt->fname);
1.5900 + if (sym->str == NULL)
1.5901 + sprintf(fname, "%.*g", DBL_DIG, sym->num);
1.5902 + else
1.5903 + fetch_string(mpl, sym->str, fname);
1.5904 + delete_symbol(mpl, sym);
1.5905 + /* close the current print file, if necessary */
1.5906 + if (mpl->prt_fp != NULL &&
1.5907 + (!prt->app || strcmp(mpl->prt_file, fname) != 0))
1.5908 + { xfclose(mpl->prt_fp), mpl->prt_fp = NULL;
1.5909 + xfree(mpl->prt_file), mpl->prt_file = NULL;
1.5910 + }
1.5911 + /* open the specified print file, if necessary */
1.5912 + if (mpl->prt_fp == NULL)
1.5913 + { mpl->prt_fp = xfopen(fname, prt->app ? "a" : "w");
1.5914 + if (mpl->prt_fp == NULL)
1.5915 + error(mpl, "unable to open `%s' for writing - %s",
1.5916 + fname, xerrmsg());
1.5917 + mpl->prt_file = xmalloc(strlen(fname)+1);
1.5918 + strcpy(mpl->prt_file, fname);
1.5919 + }
1.5920 + }
1.5921 + loop_within_domain(mpl, prt->domain, prt, printf_func);
1.5922 + if (mpl->prt_fp != NULL)
1.5923 + { xfflush(mpl->prt_fp);
1.5924 + if (xferror(mpl->prt_fp))
1.5925 + error(mpl, "writing error to `%s' - %s", mpl->prt_file,
1.5926 + xerrmsg());
1.5927 + }
1.5928 + return;
1.5929 +}
1.5930 +#endif
1.5931 +
1.5932 +/*----------------------------------------------------------------------
1.5933 +-- clean_printf - clean printf statement.
1.5934 +--
1.5935 +-- This routine cleans specified printf statement that assumes deleting
1.5936 +-- all stuff dynamically allocated on generating/postsolving phase. */
1.5937 +
1.5938 +void clean_printf(MPL *mpl, PRINTF *prt)
1.5939 +{ PRINTF1 *p;
1.5940 + /* clean subscript domain */
1.5941 + clean_domain(mpl, prt->domain);
1.5942 + /* clean pseudo-code for computing format string */
1.5943 + clean_code(mpl, prt->fmt);
1.5944 + /* clean printf list */
1.5945 + for (p = prt->list; p != NULL; p = p->next)
1.5946 + { /* clean pseudo-code for computing value to be printed */
1.5947 + clean_code(mpl, p->code);
1.5948 + }
1.5949 +#if 1 /* 14/VII-2006 */
1.5950 + /* clean pseudo-code for computing file name string */
1.5951 + clean_code(mpl, prt->fname);
1.5952 +#endif
1.5953 + return;
1.5954 +}
1.5955 +
1.5956 +/*----------------------------------------------------------------------
1.5957 +-- execute_for - execute for statement.
1.5958 +--
1.5959 +-- This routine executes specified for statement. */
1.5960 +
1.5961 +static int for_func(MPL *mpl, void *info)
1.5962 +{ /* this is auxiliary routine to work within domain scope */
1.5963 + FOR *fur = (FOR *)info;
1.5964 + STATEMENT *stmt, *save;
1.5965 + save = mpl->stmt;
1.5966 + for (stmt = fur->list; stmt != NULL; stmt = stmt->next)
1.5967 + execute_statement(mpl, stmt);
1.5968 + mpl->stmt = save;
1.5969 + return 0;
1.5970 +}
1.5971 +
1.5972 +void execute_for(MPL *mpl, FOR *fur)
1.5973 +{ loop_within_domain(mpl, fur->domain, fur, for_func);
1.5974 + return;
1.5975 +}
1.5976 +
1.5977 +/*----------------------------------------------------------------------
1.5978 +-- clean_for - clean for statement.
1.5979 +--
1.5980 +-- This routine cleans specified for statement that assumes deleting all
1.5981 +-- stuff dynamically allocated on generating/postsolving phase. */
1.5982 +
1.5983 +void clean_for(MPL *mpl, FOR *fur)
1.5984 +{ STATEMENT *stmt;
1.5985 + /* clean subscript domain */
1.5986 + clean_domain(mpl, fur->domain);
1.5987 + /* clean all sub-statements */
1.5988 + for (stmt = fur->list; stmt != NULL; stmt = stmt->next)
1.5989 + clean_statement(mpl, stmt);
1.5990 + return;
1.5991 +}
1.5992 +
1.5993 +/*----------------------------------------------------------------------
1.5994 +-- execute_statement - execute specified model statement.
1.5995 +--
1.5996 +-- This routine executes specified model statement. */
1.5997 +
1.5998 +void execute_statement(MPL *mpl, STATEMENT *stmt)
1.5999 +{ mpl->stmt = stmt;
1.6000 + switch (stmt->type)
1.6001 + { case A_SET:
1.6002 + case A_PARAMETER:
1.6003 + case A_VARIABLE:
1.6004 + break;
1.6005 + case A_CONSTRAINT:
1.6006 + xprintf("Generating %s...\n", stmt->u.con->name);
1.6007 + eval_whole_con(mpl, stmt->u.con);
1.6008 + break;
1.6009 + case A_TABLE:
1.6010 + switch (stmt->u.tab->type)
1.6011 + { case A_INPUT:
1.6012 + xprintf("Reading %s...\n", stmt->u.tab->name);
1.6013 + break;
1.6014 + case A_OUTPUT:
1.6015 + xprintf("Writing %s...\n", stmt->u.tab->name);
1.6016 + break;
1.6017 + default:
1.6018 + xassert(stmt != stmt);
1.6019 + }
1.6020 + execute_table(mpl, stmt->u.tab);
1.6021 + break;
1.6022 + case A_SOLVE:
1.6023 + break;
1.6024 + case A_CHECK:
1.6025 + xprintf("Checking (line %d)...\n", stmt->line);
1.6026 + execute_check(mpl, stmt->u.chk);
1.6027 + break;
1.6028 + case A_DISPLAY:
1.6029 + write_text(mpl, "Display statement at line %d\n",
1.6030 + stmt->line);
1.6031 + execute_display(mpl, stmt->u.dpy);
1.6032 + break;
1.6033 + case A_PRINTF:
1.6034 + execute_printf(mpl, stmt->u.prt);
1.6035 + break;
1.6036 + case A_FOR:
1.6037 + execute_for(mpl, stmt->u.fur);
1.6038 + break;
1.6039 + default:
1.6040 + xassert(stmt != stmt);
1.6041 + }
1.6042 + return;
1.6043 +}
1.6044 +
1.6045 +/*----------------------------------------------------------------------
1.6046 +-- clean_statement - clean specified model statement.
1.6047 +--
1.6048 +-- This routine cleans specified model statement that assumes deleting
1.6049 +-- all stuff dynamically allocated on generating/postsolving phase. */
1.6050 +
1.6051 +void clean_statement(MPL *mpl, STATEMENT *stmt)
1.6052 +{ switch(stmt->type)
1.6053 + { case A_SET:
1.6054 + clean_set(mpl, stmt->u.set); break;
1.6055 + case A_PARAMETER:
1.6056 + clean_parameter(mpl, stmt->u.par); break;
1.6057 + case A_VARIABLE:
1.6058 + clean_variable(mpl, stmt->u.var); break;
1.6059 + case A_CONSTRAINT:
1.6060 + clean_constraint(mpl, stmt->u.con); break;
1.6061 +#if 1 /* 11/II-2008 */
1.6062 + case A_TABLE:
1.6063 + clean_table(mpl, stmt->u.tab); break;
1.6064 +#endif
1.6065 + case A_SOLVE:
1.6066 + break;
1.6067 + case A_CHECK:
1.6068 + clean_check(mpl, stmt->u.chk); break;
1.6069 + case A_DISPLAY:
1.6070 + clean_display(mpl, stmt->u.dpy); break;
1.6071 + case A_PRINTF:
1.6072 + clean_printf(mpl, stmt->u.prt); break;
1.6073 + case A_FOR:
1.6074 + clean_for(mpl, stmt->u.fur); break;
1.6075 + default:
1.6076 + xassert(stmt != stmt);
1.6077 + }
1.6078 + return;
1.6079 +}
1.6080 +
1.6081 +/* eof */