lemon-project-template-glpk
diff deps/glpk/src/glpmpl03.c @ 9:33de93886c88
Import GLPK 4.47
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 20:59:10 +0100 |
parents | |
children |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/deps/glpk/src/glpmpl03.c Sun Nov 06 20:59:10 2011 +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, 2011 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 */