lemon-project-template-glpk
view 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 source
1 /* glpmpl03.c */
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 * E-mail: <mao@gnu.org>.
10 *
11 * GLPK is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
25 #define _GLPSTD_ERRNO
26 #define _GLPSTD_STDIO
27 #include "glpenv.h"
28 #include "glpmpl.h"
30 /**********************************************************************/
31 /* * * FLOATING-POINT NUMBERS * * */
32 /**********************************************************************/
34 /*----------------------------------------------------------------------
35 -- fp_add - floating-point addition.
36 --
37 -- This routine computes the sum x + y. */
39 double fp_add(MPL *mpl, double x, double y)
40 { if (x > 0.0 && y > 0.0 && x > + 0.999 * DBL_MAX - y ||
41 x < 0.0 && y < 0.0 && x < - 0.999 * DBL_MAX - y)
42 error(mpl, "%.*g + %.*g; floating-point overflow",
43 DBL_DIG, x, DBL_DIG, y);
44 return x + y;
45 }
47 /*----------------------------------------------------------------------
48 -- fp_sub - floating-point subtraction.
49 --
50 -- This routine computes the difference x - y. */
52 double fp_sub(MPL *mpl, double x, double y)
53 { if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y ||
54 x < 0.0 && y > 0.0 && x < - 0.999 * DBL_MAX + y)
55 error(mpl, "%.*g - %.*g; floating-point overflow",
56 DBL_DIG, x, DBL_DIG, y);
57 return x - y;
58 }
60 /*----------------------------------------------------------------------
61 -- fp_less - floating-point non-negative subtraction.
62 --
63 -- This routine computes the non-negative difference max(0, x - y). */
65 double fp_less(MPL *mpl, double x, double y)
66 { if (x < y) return 0.0;
67 if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y)
68 error(mpl, "%.*g less %.*g; floating-point overflow",
69 DBL_DIG, x, DBL_DIG, y);
70 return x - y;
71 }
73 /*----------------------------------------------------------------------
74 -- fp_mul - floating-point multiplication.
75 --
76 -- This routine computes the product x * y. */
78 double fp_mul(MPL *mpl, double x, double y)
79 { if (fabs(y) > 1.0 && fabs(x) > (0.999 * DBL_MAX) / fabs(y))
80 error(mpl, "%.*g * %.*g; floating-point overflow",
81 DBL_DIG, x, DBL_DIG, y);
82 return x * y;
83 }
85 /*----------------------------------------------------------------------
86 -- fp_div - floating-point division.
87 --
88 -- This routine computes the quotient x / y. */
90 double fp_div(MPL *mpl, double x, double y)
91 { if (fabs(y) < DBL_MIN)
92 error(mpl, "%.*g / %.*g; floating-point zero divide",
93 DBL_DIG, x, DBL_DIG, y);
94 if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))
95 error(mpl, "%.*g / %.*g; floating-point overflow",
96 DBL_DIG, x, DBL_DIG, y);
97 return x / y;
98 }
100 /*----------------------------------------------------------------------
101 -- fp_idiv - floating-point quotient of exact division.
102 --
103 -- This routine computes the quotient of exact division x div y. */
105 double fp_idiv(MPL *mpl, double x, double y)
106 { if (fabs(y) < DBL_MIN)
107 error(mpl, "%.*g div %.*g; floating-point zero divide",
108 DBL_DIG, x, DBL_DIG, y);
109 if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))
110 error(mpl, "%.*g div %.*g; floating-point overflow",
111 DBL_DIG, x, DBL_DIG, y);
112 x /= y;
113 return x > 0.0 ? floor(x) : x < 0.0 ? ceil(x) : 0.0;
114 }
116 /*----------------------------------------------------------------------
117 -- fp_mod - floating-point remainder of exact division.
118 --
119 -- This routine computes the remainder of exact division x mod y.
120 --
121 -- NOTE: By definition x mod y = x - y * floor(x / y). */
123 double fp_mod(MPL *mpl, double x, double y)
124 { double r;
125 xassert(mpl == mpl);
126 if (x == 0.0)
127 r = 0.0;
128 else if (y == 0.0)
129 r = x;
130 else
131 { r = fmod(fabs(x), fabs(y));
132 if (r != 0.0)
133 { if (x < 0.0) r = - r;
134 if (x > 0.0 && y < 0.0 || x < 0.0 && y > 0.0) r += y;
135 }
136 }
137 return r;
138 }
140 /*----------------------------------------------------------------------
141 -- fp_power - floating-point exponentiation (raise to power).
142 --
143 -- This routine computes the exponentiation x ** y. */
145 double fp_power(MPL *mpl, double x, double y)
146 { double r;
147 if (x == 0.0 && y <= 0.0 || x < 0.0 && y != floor(y))
148 error(mpl, "%.*g ** %.*g; result undefined",
149 DBL_DIG, x, DBL_DIG, y);
150 if (x == 0.0) goto eval;
151 if (fabs(x) > 1.0 && y > +1.0 &&
152 +log(fabs(x)) > (0.999 * log(DBL_MAX)) / y ||
153 fabs(x) < 1.0 && y < -1.0 &&
154 +log(fabs(x)) < (0.999 * log(DBL_MAX)) / y)
155 error(mpl, "%.*g ** %.*g; floating-point overflow",
156 DBL_DIG, x, DBL_DIG, y);
157 if (fabs(x) > 1.0 && y < -1.0 &&
158 -log(fabs(x)) < (0.999 * log(DBL_MAX)) / y ||
159 fabs(x) < 1.0 && y > +1.0 &&
160 -log(fabs(x)) > (0.999 * log(DBL_MAX)) / y)
161 r = 0.0;
162 else
163 eval: r = pow(x, y);
164 return r;
165 }
167 /*----------------------------------------------------------------------
168 -- fp_exp - floating-point base-e exponential.
169 --
170 -- This routine computes the base-e exponential e ** x. */
172 double fp_exp(MPL *mpl, double x)
173 { if (x > 0.999 * log(DBL_MAX))
174 error(mpl, "exp(%.*g); floating-point overflow", DBL_DIG, x);
175 return exp(x);
176 }
178 /*----------------------------------------------------------------------
179 -- fp_log - floating-point natural logarithm.
180 --
181 -- This routine computes the natural logarithm log x. */
183 double fp_log(MPL *mpl, double x)
184 { if (x <= 0.0)
185 error(mpl, "log(%.*g); non-positive argument", DBL_DIG, x);
186 return log(x);
187 }
189 /*----------------------------------------------------------------------
190 -- fp_log10 - floating-point common (decimal) logarithm.
191 --
192 -- This routine computes the common (decimal) logarithm lg x. */
194 double fp_log10(MPL *mpl, double x)
195 { if (x <= 0.0)
196 error(mpl, "log10(%.*g); non-positive argument", DBL_DIG, x);
197 return log10(x);
198 }
200 /*----------------------------------------------------------------------
201 -- fp_sqrt - floating-point square root.
202 --
203 -- This routine computes the square root x ** 0.5. */
205 double fp_sqrt(MPL *mpl, double x)
206 { if (x < 0.0)
207 error(mpl, "sqrt(%.*g); negative argument", DBL_DIG, x);
208 return sqrt(x);
209 }
211 /*----------------------------------------------------------------------
212 -- fp_sin - floating-point trigonometric sine.
213 --
214 -- This routine computes the trigonometric sine sin(x). */
216 double fp_sin(MPL *mpl, double x)
217 { if (!(-1e6 <= x && x <= +1e6))
218 error(mpl, "sin(%.*g); argument too large", DBL_DIG, x);
219 return sin(x);
220 }
222 /*----------------------------------------------------------------------
223 -- fp_cos - floating-point trigonometric cosine.
224 --
225 -- This routine computes the trigonometric cosine cos(x). */
227 double fp_cos(MPL *mpl, double x)
228 { if (!(-1e6 <= x && x <= +1e6))
229 error(mpl, "cos(%.*g); argument too large", DBL_DIG, x);
230 return cos(x);
231 }
233 /*----------------------------------------------------------------------
234 -- fp_atan - floating-point trigonometric arctangent.
235 --
236 -- This routine computes the trigonometric arctangent atan(x). */
238 double fp_atan(MPL *mpl, double x)
239 { xassert(mpl == mpl);
240 return atan(x);
241 }
243 /*----------------------------------------------------------------------
244 -- fp_atan2 - floating-point trigonometric arctangent.
245 --
246 -- This routine computes the trigonometric arctangent atan(y / x). */
248 double fp_atan2(MPL *mpl, double y, double x)
249 { xassert(mpl == mpl);
250 return atan2(y, x);
251 }
253 /*----------------------------------------------------------------------
254 -- fp_round - round floating-point value to n fractional digits.
255 --
256 -- This routine rounds given floating-point value x to n fractional
257 -- digits with the formula:
258 --
259 -- round(x, n) = floor(x * 10^n + 0.5) / 10^n.
260 --
261 -- The parameter n is assumed to be integer. */
263 double fp_round(MPL *mpl, double x, double n)
264 { double ten_to_n;
265 if (n != floor(n))
266 error(mpl, "round(%.*g, %.*g); non-integer second argument",
267 DBL_DIG, x, DBL_DIG, n);
268 if (n <= DBL_DIG + 2)
269 { ten_to_n = pow(10.0, n);
270 if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)
271 { x = floor(x * ten_to_n + 0.5);
272 if (x != 0.0) x /= ten_to_n;
273 }
274 }
275 return x;
276 }
278 /*----------------------------------------------------------------------
279 -- fp_trunc - truncate floating-point value to n fractional digits.
280 --
281 -- This routine truncates given floating-point value x to n fractional
282 -- digits with the formula:
283 --
284 -- ( floor(x * 10^n) / 10^n, if x >= 0
285 -- trunc(x, n) = <
286 -- ( ceil(x * 10^n) / 10^n, if x < 0
287 --
288 -- The parameter n is assumed to be integer. */
290 double fp_trunc(MPL *mpl, double x, double n)
291 { double ten_to_n;
292 if (n != floor(n))
293 error(mpl, "trunc(%.*g, %.*g); non-integer second argument",
294 DBL_DIG, x, DBL_DIG, n);
295 if (n <= DBL_DIG + 2)
296 { ten_to_n = pow(10.0, n);
297 if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)
298 { x = (x >= 0.0 ? floor(x * ten_to_n) : ceil(x * ten_to_n));
299 if (x != 0.0) x /= ten_to_n;
300 }
301 }
302 return x;
303 }
305 /**********************************************************************/
306 /* * * PSEUDO-RANDOM NUMBER GENERATORS * * */
307 /**********************************************************************/
309 /*----------------------------------------------------------------------
310 -- fp_irand224 - pseudo-random integer in the range [0, 2^24).
311 --
312 -- This routine returns a next pseudo-random integer (converted to
313 -- floating-point) which is uniformly distributed between 0 and 2^24-1,
314 -- inclusive. */
316 #define two_to_the_24 0x1000000
318 double fp_irand224(MPL *mpl)
319 { return
320 (double)rng_unif_rand(mpl->rand, two_to_the_24);
321 }
323 /*----------------------------------------------------------------------
324 -- fp_uniform01 - pseudo-random number in the range [0, 1).
325 --
326 -- This routine returns a next pseudo-random number which is uniformly
327 -- distributed in the range [0, 1). */
329 #define two_to_the_31 ((unsigned int)0x80000000)
331 double fp_uniform01(MPL *mpl)
332 { return
333 (double)rng_next_rand(mpl->rand) / (double)two_to_the_31;
334 }
336 /*----------------------------------------------------------------------
337 -- fp_uniform - pseudo-random number in the range [a, b).
338 --
339 -- This routine returns a next pseudo-random number which is uniformly
340 -- distributed in the range [a, b). */
342 double fp_uniform(MPL *mpl, double a, double b)
343 { double x;
344 if (a >= b)
345 error(mpl, "Uniform(%.*g, %.*g); invalid range",
346 DBL_DIG, a, DBL_DIG, b);
347 x = fp_uniform01(mpl);
348 #if 0
349 x = a * (1.0 - x) + b * x;
350 #else
351 x = fp_add(mpl, a * (1.0 - x), b * x);
352 #endif
353 return x;
354 }
356 /*----------------------------------------------------------------------
357 -- fp_normal01 - Gaussian random variate with mu = 0 and sigma = 1.
358 --
359 -- This routine returns a Gaussian random variate with zero mean and
360 -- unit standard deviation. The polar (Box-Mueller) method is used.
361 --
362 -- This code is a modified version of the routine gsl_ran_gaussian from
363 -- the GNU Scientific Library Version 1.0. */
365 double fp_normal01(MPL *mpl)
366 { double x, y, r2;
367 do
368 { /* choose x, y in uniform square (-1,-1) to (+1,+1) */
369 x = -1.0 + 2.0 * fp_uniform01(mpl);
370 y = -1.0 + 2.0 * fp_uniform01(mpl);
371 /* see if it is in the unit circle */
372 r2 = x * x + y * y;
373 } while (r2 > 1.0 || r2 == 0.0);
374 /* Box-Muller transform */
375 return y * sqrt(-2.0 * log (r2) / r2);
376 }
378 /*----------------------------------------------------------------------
379 -- fp_normal - Gaussian random variate with specified mu and sigma.
380 --
381 -- This routine returns a Gaussian random variate with mean mu and
382 -- standard deviation sigma. */
384 double fp_normal(MPL *mpl, double mu, double sigma)
385 { double x;
386 #if 0
387 x = mu + sigma * fp_normal01(mpl);
388 #else
389 x = fp_add(mpl, mu, fp_mul(mpl, sigma, fp_normal01(mpl)));
390 #endif
391 return x;
392 }
394 /**********************************************************************/
395 /* * * SEGMENTED CHARACTER STRINGS * * */
396 /**********************************************************************/
398 /*----------------------------------------------------------------------
399 -- create_string - create character string.
400 --
401 -- This routine creates a segmented character string, which is exactly
402 -- equivalent to specified character string. */
404 STRING *create_string
405 ( MPL *mpl,
406 char buf[MAX_LENGTH+1] /* not changed */
407 )
408 #if 0
409 { STRING *head, *tail;
410 int i, j;
411 xassert(buf != NULL);
412 xassert(strlen(buf) <= MAX_LENGTH);
413 head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));
414 for (i = j = 0; ; i++)
415 { if ((tail->seg[j++] = buf[i]) == '\0') break;
416 if (j == STRSEG_SIZE)
417 tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING))), j = 0;
418 }
419 tail->next = NULL;
420 return head;
421 }
422 #else
423 { STRING *str;
424 xassert(strlen(buf) <= MAX_LENGTH);
425 str = dmp_get_atom(mpl->strings, strlen(buf)+1);
426 strcpy(str, buf);
427 return str;
428 }
429 #endif
431 /*----------------------------------------------------------------------
432 -- copy_string - make copy of character string.
433 --
434 -- This routine returns an exact copy of segmented character string. */
436 STRING *copy_string
437 ( MPL *mpl,
438 STRING *str /* not changed */
439 )
440 #if 0
441 { STRING *head, *tail;
442 xassert(str != NULL);
443 head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));
444 for (; str != NULL; str = str->next)
445 { memcpy(tail->seg, str->seg, STRSEG_SIZE);
446 if (str->next != NULL)
447 tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING)));
448 }
449 tail->next = NULL;
450 return head;
451 }
452 #else
453 { xassert(mpl == mpl);
454 return create_string(mpl, str);
455 }
456 #endif
458 /*----------------------------------------------------------------------
459 -- compare_strings - compare one character string with another.
460 --
461 -- This routine compares one segmented character strings with another
462 -- and returns the result of comparison as follows:
463 --
464 -- = 0 - both strings are identical;
465 -- < 0 - the first string precedes the second one;
466 -- > 0 - the first string follows the second one. */
468 int compare_strings
469 ( MPL *mpl,
470 STRING *str1, /* not changed */
471 STRING *str2 /* not changed */
472 )
473 #if 0
474 { int j, c1, c2;
475 xassert(mpl == mpl);
476 for (;; str1 = str1->next, str2 = str2->next)
477 { xassert(str1 != NULL);
478 xassert(str2 != NULL);
479 for (j = 0; j < STRSEG_SIZE; j++)
480 { c1 = (unsigned char)str1->seg[j];
481 c2 = (unsigned char)str2->seg[j];
482 if (c1 < c2) return -1;
483 if (c1 > c2) return +1;
484 if (c1 == '\0') goto done;
485 }
486 }
487 done: return 0;
488 }
489 #else
490 { xassert(mpl == mpl);
491 return strcmp(str1, str2);
492 }
493 #endif
495 /*----------------------------------------------------------------------
496 -- fetch_string - extract content of character string.
497 --
498 -- This routine returns a character string, which is exactly equivalent
499 -- to specified segmented character string. */
501 char *fetch_string
502 ( MPL *mpl,
503 STRING *str, /* not changed */
504 char buf[MAX_LENGTH+1] /* modified */
505 )
506 #if 0
507 { int i, j;
508 xassert(mpl == mpl);
509 xassert(buf != NULL);
510 for (i = 0; ; str = str->next)
511 { xassert(str != NULL);
512 for (j = 0; j < STRSEG_SIZE; j++)
513 if ((buf[i++] = str->seg[j]) == '\0') goto done;
514 }
515 done: xassert(strlen(buf) <= MAX_LENGTH);
516 return buf;
517 }
518 #else
519 { xassert(mpl == mpl);
520 return strcpy(buf, str);
521 }
522 #endif
524 /*----------------------------------------------------------------------
525 -- delete_string - delete character string.
526 --
527 -- This routine deletes specified segmented character string. */
529 void delete_string
530 ( MPL *mpl,
531 STRING *str /* destroyed */
532 )
533 #if 0
534 { STRING *temp;
535 xassert(str != NULL);
536 while (str != NULL)
537 { temp = str;
538 str = str->next;
539 dmp_free_atom(mpl->strings, temp, sizeof(STRING));
540 }
541 return;
542 }
543 #else
544 { dmp_free_atom(mpl->strings, str, strlen(str)+1);
545 return;
546 }
547 #endif
549 /**********************************************************************/
550 /* * * SYMBOLS * * */
551 /**********************************************************************/
553 /*----------------------------------------------------------------------
554 -- create_symbol_num - create symbol of numeric type.
555 --
556 -- This routine creates a symbol, which has a numeric value specified
557 -- as floating-point number. */
559 SYMBOL *create_symbol_num(MPL *mpl, double num)
560 { SYMBOL *sym;
561 sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));
562 sym->num = num;
563 sym->str = NULL;
564 return sym;
565 }
567 /*----------------------------------------------------------------------
568 -- create_symbol_str - create symbol of abstract type.
569 --
570 -- This routine creates a symbol, which has an abstract value specified
571 -- as segmented character string. */
573 SYMBOL *create_symbol_str
574 ( MPL *mpl,
575 STRING *str /* destroyed */
576 )
577 { SYMBOL *sym;
578 xassert(str != NULL);
579 sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));
580 sym->num = 0.0;
581 sym->str = str;
582 return sym;
583 }
585 /*----------------------------------------------------------------------
586 -- copy_symbol - make copy of symbol.
587 --
588 -- This routine returns an exact copy of symbol. */
590 SYMBOL *copy_symbol
591 ( MPL *mpl,
592 SYMBOL *sym /* not changed */
593 )
594 { SYMBOL *copy;
595 xassert(sym != NULL);
596 copy = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));
597 if (sym->str == NULL)
598 { copy->num = sym->num;
599 copy->str = NULL;
600 }
601 else
602 { copy->num = 0.0;
603 copy->str = copy_string(mpl, sym->str);
604 }
605 return copy;
606 }
608 /*----------------------------------------------------------------------
609 -- compare_symbols - compare one symbol with another.
610 --
611 -- This routine compares one symbol with another and returns the result
612 -- of comparison as follows:
613 --
614 -- = 0 - both symbols are identical;
615 -- < 0 - the first symbol precedes the second one;
616 -- > 0 - the first symbol follows the second one.
617 --
618 -- Note that the linear order, in which symbols follow each other, is
619 -- implementation-dependent. It may be not an alphabetical order. */
621 int compare_symbols
622 ( MPL *mpl,
623 SYMBOL *sym1, /* not changed */
624 SYMBOL *sym2 /* not changed */
625 )
626 { xassert(sym1 != NULL);
627 xassert(sym2 != NULL);
628 /* let all numeric quantities precede all symbolic quantities */
629 if (sym1->str == NULL && sym2->str == NULL)
630 { if (sym1->num < sym2->num) return -1;
631 if (sym1->num > sym2->num) return +1;
632 return 0;
633 }
634 if (sym1->str == NULL) return -1;
635 if (sym2->str == NULL) return +1;
636 return compare_strings(mpl, sym1->str, sym2->str);
637 }
639 /*----------------------------------------------------------------------
640 -- delete_symbol - delete symbol.
641 --
642 -- This routine deletes specified symbol. */
644 void delete_symbol
645 ( MPL *mpl,
646 SYMBOL *sym /* destroyed */
647 )
648 { xassert(sym != NULL);
649 if (sym->str != NULL) delete_string(mpl, sym->str);
650 dmp_free_atom(mpl->symbols, sym, sizeof(SYMBOL));
651 return;
652 }
654 /*----------------------------------------------------------------------
655 -- format_symbol - format symbol for displaying or printing.
656 --
657 -- This routine converts specified symbol to a charater string, which
658 -- is suitable for displaying or printing.
659 --
660 -- The resultant string is never longer than 255 characters. If it gets
661 -- longer, it is truncated from the right and appended by dots. */
663 char *format_symbol
664 ( MPL *mpl,
665 SYMBOL *sym /* not changed */
666 )
667 { char *buf = mpl->sym_buf;
668 xassert(sym != NULL);
669 if (sym->str == NULL)
670 sprintf(buf, "%.*g", DBL_DIG, sym->num);
671 else
672 { char str[MAX_LENGTH+1];
673 int quoted, j, len;
674 fetch_string(mpl, sym->str, str);
675 if (!(isalpha((unsigned char)str[0]) || str[0] == '_'))
676 quoted = 1;
677 else
678 { quoted = 0;
679 for (j = 1; str[j] != '\0'; j++)
680 { if (!(isalnum((unsigned char)str[j]) ||
681 strchr("+-._", (unsigned char)str[j]) != NULL))
682 { quoted = 1;
683 break;
684 }
685 }
686 }
687 # define safe_append(c) \
688 (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)
689 buf[0] = '\0', len = 0;
690 if (quoted) safe_append('\'');
691 for (j = 0; str[j] != '\0'; j++)
692 { if (quoted && str[j] == '\'') safe_append('\'');
693 safe_append(str[j]);
694 }
695 if (quoted) safe_append('\'');
696 # undef safe_append
697 buf[len] = '\0';
698 if (len == 255) strcpy(buf+252, "...");
699 }
700 xassert(strlen(buf) <= 255);
701 return buf;
702 }
704 /*----------------------------------------------------------------------
705 -- concat_symbols - concatenate one symbol with another.
706 --
707 -- This routine concatenates values of two given symbols and assigns
708 -- the resultant character string to a new symbol, which is returned on
709 -- exit. Both original symbols are destroyed. */
711 SYMBOL *concat_symbols
712 ( MPL *mpl,
713 SYMBOL *sym1, /* destroyed */
714 SYMBOL *sym2 /* destroyed */
715 )
716 { char str1[MAX_LENGTH+1], str2[MAX_LENGTH+1];
717 xassert(MAX_LENGTH >= DBL_DIG + DBL_DIG);
718 if (sym1->str == NULL)
719 sprintf(str1, "%.*g", DBL_DIG, sym1->num);
720 else
721 fetch_string(mpl, sym1->str, str1);
722 if (sym2->str == NULL)
723 sprintf(str2, "%.*g", DBL_DIG, sym2->num);
724 else
725 fetch_string(mpl, sym2->str, str2);
726 if (strlen(str1) + strlen(str2) > MAX_LENGTH)
727 { char buf[255+1];
728 strcpy(buf, format_symbol(mpl, sym1));
729 xassert(strlen(buf) < sizeof(buf));
730 error(mpl, "%s & %s; resultant symbol exceeds %d characters",
731 buf, format_symbol(mpl, sym2), MAX_LENGTH);
732 }
733 delete_symbol(mpl, sym1);
734 delete_symbol(mpl, sym2);
735 return create_symbol_str(mpl, create_string(mpl, strcat(str1,
736 str2)));
737 }
739 /**********************************************************************/
740 /* * * N-TUPLES * * */
741 /**********************************************************************/
743 /*----------------------------------------------------------------------
744 -- create_tuple - create n-tuple.
745 --
746 -- This routine creates a n-tuple, which initially has no components,
747 -- i.e. which is 0-tuple. */
749 TUPLE *create_tuple(MPL *mpl)
750 { TUPLE *tuple;
751 xassert(mpl == mpl);
752 tuple = NULL;
753 return tuple;
754 }
756 /*----------------------------------------------------------------------
757 -- expand_tuple - append symbol to n-tuple.
758 --
759 -- This routine expands n-tuple appending to it a given symbol, which
760 -- becomes its new last component. */
762 TUPLE *expand_tuple
763 ( MPL *mpl,
764 TUPLE *tuple, /* destroyed */
765 SYMBOL *sym /* destroyed */
766 )
767 { TUPLE *tail, *temp;
768 xassert(sym != NULL);
769 /* create a new component */
770 tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
771 tail->sym = sym;
772 tail->next = NULL;
773 /* and append it to the component list */
774 if (tuple == NULL)
775 tuple = tail;
776 else
777 { for (temp = tuple; temp->next != NULL; temp = temp->next);
778 temp->next = tail;
779 }
780 return tuple;
781 }
783 /*----------------------------------------------------------------------
784 -- tuple_dimen - determine dimension of n-tuple.
785 --
786 -- This routine returns dimension of n-tuple, i.e. number of components
787 -- in the n-tuple. */
789 int tuple_dimen
790 ( MPL *mpl,
791 TUPLE *tuple /* not changed */
792 )
793 { TUPLE *temp;
794 int dim = 0;
795 xassert(mpl == mpl);
796 for (temp = tuple; temp != NULL; temp = temp->next) dim++;
797 return dim;
798 }
800 /*----------------------------------------------------------------------
801 -- copy_tuple - make copy of n-tuple.
802 --
803 -- This routine returns an exact copy of n-tuple. */
805 TUPLE *copy_tuple
806 ( MPL *mpl,
807 TUPLE *tuple /* not changed */
808 )
809 { TUPLE *head, *tail;
810 if (tuple == NULL)
811 head = NULL;
812 else
813 { head = tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
814 for (; tuple != NULL; tuple = tuple->next)
815 { xassert(tuple->sym != NULL);
816 tail->sym = copy_symbol(mpl, tuple->sym);
817 if (tuple->next != NULL)
818 tail = (tail->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));
819 }
820 tail->next = NULL;
821 }
822 return head;
823 }
825 /*----------------------------------------------------------------------
826 -- compare_tuples - compare one n-tuple with another.
827 --
828 -- This routine compares two given n-tuples, which must have the same
829 -- dimension (not checked for the sake of efficiency), and returns one
830 -- of the following codes:
831 --
832 -- = 0 - both n-tuples are identical;
833 -- < 0 - the first n-tuple precedes the second one;
834 -- > 0 - the first n-tuple follows the second one.
835 --
836 -- Note that the linear order, in which n-tuples follow each other, is
837 -- implementation-dependent. It may be not an alphabetical order. */
839 int compare_tuples
840 ( MPL *mpl,
841 TUPLE *tuple1, /* not changed */
842 TUPLE *tuple2 /* not changed */
843 )
844 { TUPLE *item1, *item2;
845 int ret;
846 xassert(mpl == mpl);
847 for (item1 = tuple1, item2 = tuple2; item1 != NULL;
848 item1 = item1->next, item2 = item2->next)
849 { xassert(item2 != NULL);
850 xassert(item1->sym != NULL);
851 xassert(item2->sym != NULL);
852 ret = compare_symbols(mpl, item1->sym, item2->sym);
853 if (ret != 0) return ret;
854 }
855 xassert(item2 == NULL);
856 return 0;
857 }
859 /*----------------------------------------------------------------------
860 -- build_subtuple - build subtuple of given n-tuple.
861 --
862 -- This routine builds subtuple, which consists of first dim components
863 -- of given n-tuple. */
865 TUPLE *build_subtuple
866 ( MPL *mpl,
867 TUPLE *tuple, /* not changed */
868 int dim
869 )
870 { TUPLE *head, *temp;
871 int j;
872 head = create_tuple(mpl);
873 for (j = 1, temp = tuple; j <= dim; j++, temp = temp->next)
874 { xassert(temp != NULL);
875 head = expand_tuple(mpl, head, copy_symbol(mpl, temp->sym));
876 }
877 return head;
878 }
880 /*----------------------------------------------------------------------
881 -- delete_tuple - delete n-tuple.
882 --
883 -- This routine deletes specified n-tuple. */
885 void delete_tuple
886 ( MPL *mpl,
887 TUPLE *tuple /* destroyed */
888 )
889 { TUPLE *temp;
890 while (tuple != NULL)
891 { temp = tuple;
892 tuple = temp->next;
893 xassert(temp->sym != NULL);
894 delete_symbol(mpl, temp->sym);
895 dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));
896 }
897 return;
898 }
900 /*----------------------------------------------------------------------
901 -- format_tuple - format n-tuple for displaying or printing.
902 --
903 -- This routine converts specified n-tuple to a character string, which
904 -- is suitable for displaying or printing.
905 --
906 -- The resultant string is never longer than 255 characters. If it gets
907 -- longer, it is truncated from the right and appended by dots. */
909 char *format_tuple
910 ( MPL *mpl,
911 int c,
912 TUPLE *tuple /* not changed */
913 )
914 { TUPLE *temp;
915 int dim, j, len;
916 char *buf = mpl->tup_buf, str[255+1], *save;
917 # define safe_append(c) \
918 (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)
919 buf[0] = '\0', len = 0;
920 dim = tuple_dimen(mpl, tuple);
921 if (c == '[' && dim > 0) safe_append('[');
922 if (c == '(' && dim > 1) safe_append('(');
923 for (temp = tuple; temp != NULL; temp = temp->next)
924 { if (temp != tuple) safe_append(',');
925 xassert(temp->sym != NULL);
926 save = mpl->sym_buf;
927 mpl->sym_buf = str;
928 format_symbol(mpl, temp->sym);
929 mpl->sym_buf = save;
930 xassert(strlen(str) < sizeof(str));
931 for (j = 0; str[j] != '\0'; j++) safe_append(str[j]);
932 }
933 if (c == '[' && dim > 0) safe_append(']');
934 if (c == '(' && dim > 1) safe_append(')');
935 # undef safe_append
936 buf[len] = '\0';
937 if (len == 255) strcpy(buf+252, "...");
938 xassert(strlen(buf) <= 255);
939 return buf;
940 }
942 /**********************************************************************/
943 /* * * ELEMENTAL SETS * * */
944 /**********************************************************************/
946 /*----------------------------------------------------------------------
947 -- create_elemset - create elemental set.
948 --
949 -- This routine creates an elemental set, whose members are n-tuples of
950 -- specified dimension. Being created the set is initially empty. */
952 ELEMSET *create_elemset(MPL *mpl, int dim)
953 { ELEMSET *set;
954 xassert(dim > 0);
955 set = create_array(mpl, A_NONE, dim);
956 return set;
957 }
959 /*----------------------------------------------------------------------
960 -- find_tuple - check if elemental set contains given n-tuple.
961 --
962 -- This routine finds given n-tuple in specified elemental set in order
963 -- to check if the set contains that n-tuple. If the n-tuple is found,
964 -- the routine returns pointer to corresponding array member. Otherwise
965 -- null pointer is returned. */
967 MEMBER *find_tuple
968 ( MPL *mpl,
969 ELEMSET *set, /* not changed */
970 TUPLE *tuple /* not changed */
971 )
972 { xassert(set != NULL);
973 xassert(set->type == A_NONE);
974 xassert(set->dim == tuple_dimen(mpl, tuple));
975 return find_member(mpl, set, tuple);
976 }
978 /*----------------------------------------------------------------------
979 -- add_tuple - add new n-tuple to elemental set.
980 --
981 -- This routine adds given n-tuple to specified elemental set.
982 --
983 -- For the sake of efficiency this routine doesn't check whether the
984 -- set already contains the same n-tuple or not. Therefore the calling
985 -- program should use the routine find_tuple (if necessary) in order to
986 -- make sure that the given n-tuple is not contained in the set, since
987 -- duplicate n-tuples within the same set are not allowed. */
989 MEMBER *add_tuple
990 ( MPL *mpl,
991 ELEMSET *set, /* modified */
992 TUPLE *tuple /* destroyed */
993 )
994 { MEMBER *memb;
995 xassert(set != NULL);
996 xassert(set->type == A_NONE);
997 xassert(set->dim == tuple_dimen(mpl, tuple));
998 memb = add_member(mpl, set, tuple);
999 memb->value.none = NULL;
1000 return memb;
1001 }
1003 /*----------------------------------------------------------------------
1004 -- check_then_add - check and add new n-tuple to elemental set.
1005 --
1006 -- This routine is equivalent to the routine add_tuple except that it
1007 -- does check for duplicate n-tuples. */
1009 MEMBER *check_then_add
1010 ( MPL *mpl,
1011 ELEMSET *set, /* modified */
1012 TUPLE *tuple /* destroyed */
1013 )
1014 { if (find_tuple(mpl, set, tuple) != NULL)
1015 error(mpl, "duplicate tuple %s detected", format_tuple(mpl,
1016 '(', tuple));
1017 return add_tuple(mpl, set, tuple);
1018 }
1020 /*----------------------------------------------------------------------
1021 -- copy_elemset - make copy of elemental set.
1022 --
1023 -- This routine makes an exact copy of elemental set. */
1025 ELEMSET *copy_elemset
1026 ( MPL *mpl,
1027 ELEMSET *set /* not changed */
1028 )
1029 { ELEMSET *copy;
1030 MEMBER *memb;
1031 xassert(set != NULL);
1032 xassert(set->type == A_NONE);
1033 xassert(set->dim > 0);
1034 copy = create_elemset(mpl, set->dim);
1035 for (memb = set->head; memb != NULL; memb = memb->next)
1036 add_tuple(mpl, copy, copy_tuple(mpl, memb->tuple));
1037 return copy;
1038 }
1040 /*----------------------------------------------------------------------
1041 -- delete_elemset - delete elemental set.
1042 --
1043 -- This routine deletes specified elemental set. */
1045 void delete_elemset
1046 ( MPL *mpl,
1047 ELEMSET *set /* destroyed */
1048 )
1049 { xassert(set != NULL);
1050 xassert(set->type == A_NONE);
1051 delete_array(mpl, set);
1052 return;
1053 }
1055 /*----------------------------------------------------------------------
1056 -- arelset_size - compute size of "arithmetic" elemental set.
1057 --
1058 -- This routine computes the size of "arithmetic" elemental set, which
1059 -- is specified in the form of arithmetic progression:
1060 --
1061 -- { t0 .. tf by dt }.
1062 --
1063 -- The size is computed using the formula:
1064 --
1065 -- n = max(0, floor((tf - t0) / dt) + 1). */
1067 int arelset_size(MPL *mpl, double t0, double tf, double dt)
1068 { double temp;
1069 if (dt == 0.0)
1070 error(mpl, "%.*g .. %.*g by %.*g; zero stride not allowed",
1071 DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);
1072 if (tf > 0.0 && t0 < 0.0 && tf > + 0.999 * DBL_MAX + t0)
1073 temp = +DBL_MAX;
1074 else if (tf < 0.0 && t0 > 0.0 && tf < - 0.999 * DBL_MAX + t0)
1075 temp = -DBL_MAX;
1076 else
1077 temp = tf - t0;
1078 if (fabs(dt) < 1.0 && fabs(temp) > (0.999 * DBL_MAX) * fabs(dt))
1079 { if (temp > 0.0 && dt > 0.0 || temp < 0.0 && dt < 0.0)
1080 temp = +DBL_MAX;
1081 else
1082 temp = 0.0;
1083 }
1084 else
1085 { temp = floor(temp / dt) + 1.0;
1086 if (temp < 0.0) temp = 0.0;
1087 }
1088 xassert(temp >= 0.0);
1089 if (temp > (double)(INT_MAX - 1))
1090 error(mpl, "%.*g .. %.*g by %.*g; set too large",
1091 DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);
1092 return (int)(temp + 0.5);
1093 }
1095 /*----------------------------------------------------------------------
1096 -- arelset_member - compute member of "arithmetic" elemental set.
1097 --
1098 -- This routine returns a numeric value of symbol, which is equivalent
1099 -- to j-th member of given "arithmetic" elemental set specified in the
1100 -- form of arithmetic progression:
1101 --
1102 -- { t0 .. tf by dt }.
1103 --
1104 -- The symbol value is computed with the formula:
1105 --
1106 -- j-th member = t0 + (j - 1) * dt,
1107 --
1108 -- The number j must satisfy to the restriction 1 <= j <= n, where n is
1109 -- the set size computed by the routine arelset_size. */
1111 double arelset_member(MPL *mpl, double t0, double tf, double dt, int j)
1112 { xassert(1 <= j && j <= arelset_size(mpl, t0, tf, dt));
1113 return t0 + (double)(j - 1) * dt;
1114 }
1116 /*----------------------------------------------------------------------
1117 -- create_arelset - create "arithmetic" elemental set.
1118 --
1119 -- This routine creates "arithmetic" elemental set, which is specified
1120 -- in the form of arithmetic progression:
1121 --
1122 -- { t0 .. tf by dt }.
1123 --
1124 -- Components of this set are 1-tuples. */
1126 ELEMSET *create_arelset(MPL *mpl, double t0, double tf, double dt)
1127 { ELEMSET *set;
1128 int j, n;
1129 set = create_elemset(mpl, 1);
1130 n = arelset_size(mpl, t0, tf, dt);
1131 for (j = 1; j <= n; j++)
1132 { add_tuple
1133 ( mpl,
1134 set,
1135 expand_tuple
1136 ( mpl,
1137 create_tuple(mpl),
1138 create_symbol_num
1139 ( mpl,
1140 arelset_member(mpl, t0, tf, dt, j)
1141 )
1142 )
1143 );
1144 }
1145 return set;
1146 }
1148 /*----------------------------------------------------------------------
1149 -- set_union - union of two elemental sets.
1150 --
1151 -- This routine computes the union:
1152 --
1153 -- X U Y = { j | (j in X) or (j in Y) },
1154 --
1155 -- where X and Y are given elemental sets (destroyed on exit). */
1157 ELEMSET *set_union
1158 ( MPL *mpl,
1159 ELEMSET *X, /* destroyed */
1160 ELEMSET *Y /* destroyed */
1161 )
1162 { MEMBER *memb;
1163 xassert(X != NULL);
1164 xassert(X->type == A_NONE);
1165 xassert(X->dim > 0);
1166 xassert(Y != NULL);
1167 xassert(Y->type == A_NONE);
1168 xassert(Y->dim > 0);
1169 xassert(X->dim == Y->dim);
1170 for (memb = Y->head; memb != NULL; memb = memb->next)
1171 { if (find_tuple(mpl, X, memb->tuple) == NULL)
1172 add_tuple(mpl, X, copy_tuple(mpl, memb->tuple));
1173 }
1174 delete_elemset(mpl, Y);
1175 return X;
1176 }
1178 /*----------------------------------------------------------------------
1179 -- set_diff - difference between two elemental sets.
1180 --
1181 -- This routine computes the difference:
1182 --
1183 -- X \ Y = { j | (j in X) and (j not in Y) },
1184 --
1185 -- where X and Y are given elemental sets (destroyed on exit). */
1187 ELEMSET *set_diff
1188 ( MPL *mpl,
1189 ELEMSET *X, /* destroyed */
1190 ELEMSET *Y /* destroyed */
1191 )
1192 { ELEMSET *Z;
1193 MEMBER *memb;
1194 xassert(X != NULL);
1195 xassert(X->type == A_NONE);
1196 xassert(X->dim > 0);
1197 xassert(Y != NULL);
1198 xassert(Y->type == A_NONE);
1199 xassert(Y->dim > 0);
1200 xassert(X->dim == Y->dim);
1201 Z = create_elemset(mpl, X->dim);
1202 for (memb = X->head; memb != NULL; memb = memb->next)
1203 { if (find_tuple(mpl, Y, memb->tuple) == NULL)
1204 add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1205 }
1206 delete_elemset(mpl, X);
1207 delete_elemset(mpl, Y);
1208 return Z;
1209 }
1211 /*----------------------------------------------------------------------
1212 -- set_symdiff - symmetric difference between two elemental sets.
1213 --
1214 -- This routine computes the symmetric difference:
1215 --
1216 -- X (+) Y = (X \ Y) U (Y \ X),
1217 --
1218 -- where X and Y are given elemental sets (destroyed on exit). */
1220 ELEMSET *set_symdiff
1221 ( MPL *mpl,
1222 ELEMSET *X, /* destroyed */
1223 ELEMSET *Y /* destroyed */
1224 )
1225 { ELEMSET *Z;
1226 MEMBER *memb;
1227 xassert(X != NULL);
1228 xassert(X->type == A_NONE);
1229 xassert(X->dim > 0);
1230 xassert(Y != NULL);
1231 xassert(Y->type == A_NONE);
1232 xassert(Y->dim > 0);
1233 xassert(X->dim == Y->dim);
1234 /* Z := X \ Y */
1235 Z = create_elemset(mpl, X->dim);
1236 for (memb = X->head; memb != NULL; memb = memb->next)
1237 { if (find_tuple(mpl, Y, memb->tuple) == NULL)
1238 add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1239 }
1240 /* Z := Z U (Y \ X) */
1241 for (memb = Y->head; memb != NULL; memb = memb->next)
1242 { if (find_tuple(mpl, X, memb->tuple) == NULL)
1243 add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1244 }
1245 delete_elemset(mpl, X);
1246 delete_elemset(mpl, Y);
1247 return Z;
1248 }
1250 /*----------------------------------------------------------------------
1251 -- set_inter - intersection of two elemental sets.
1252 --
1253 -- This routine computes the intersection:
1254 --
1255 -- X ^ Y = { j | (j in X) and (j in Y) },
1256 --
1257 -- where X and Y are given elemental sets (destroyed on exit). */
1259 ELEMSET *set_inter
1260 ( MPL *mpl,
1261 ELEMSET *X, /* destroyed */
1262 ELEMSET *Y /* destroyed */
1263 )
1264 { ELEMSET *Z;
1265 MEMBER *memb;
1266 xassert(X != NULL);
1267 xassert(X->type == A_NONE);
1268 xassert(X->dim > 0);
1269 xassert(Y != NULL);
1270 xassert(Y->type == A_NONE);
1271 xassert(Y->dim > 0);
1272 xassert(X->dim == Y->dim);
1273 Z = create_elemset(mpl, X->dim);
1274 for (memb = X->head; memb != NULL; memb = memb->next)
1275 { if (find_tuple(mpl, Y, memb->tuple) != NULL)
1276 add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));
1277 }
1278 delete_elemset(mpl, X);
1279 delete_elemset(mpl, Y);
1280 return Z;
1281 }
1283 /*----------------------------------------------------------------------
1284 -- set_cross - cross (Cartesian) product of two elemental sets.
1285 --
1286 -- This routine computes the cross (Cartesian) product:
1287 --
1288 -- X x Y = { (i,j) | (i in X) and (j in Y) },
1289 --
1290 -- where X and Y are given elemental sets (destroyed on exit). */
1292 ELEMSET *set_cross
1293 ( MPL *mpl,
1294 ELEMSET *X, /* destroyed */
1295 ELEMSET *Y /* destroyed */
1296 )
1297 { ELEMSET *Z;
1298 MEMBER *memx, *memy;
1299 TUPLE *tuple, *temp;
1300 xassert(X != NULL);
1301 xassert(X->type == A_NONE);
1302 xassert(X->dim > 0);
1303 xassert(Y != NULL);
1304 xassert(Y->type == A_NONE);
1305 xassert(Y->dim > 0);
1306 Z = create_elemset(mpl, X->dim + Y->dim);
1307 for (memx = X->head; memx != NULL; memx = memx->next)
1308 { for (memy = Y->head; memy != NULL; memy = memy->next)
1309 { tuple = copy_tuple(mpl, memx->tuple);
1310 for (temp = memy->tuple; temp != NULL; temp = temp->next)
1311 tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,
1312 temp->sym));
1313 add_tuple(mpl, Z, tuple);
1314 }
1315 }
1316 delete_elemset(mpl, X);
1317 delete_elemset(mpl, Y);
1318 return Z;
1319 }
1321 /**********************************************************************/
1322 /* * * ELEMENTAL VARIABLES * * */
1323 /**********************************************************************/
1325 /* (there are no specific routines for elemental variables) */
1327 /**********************************************************************/
1328 /* * * LINEAR FORMS * * */
1329 /**********************************************************************/
1331 /*----------------------------------------------------------------------
1332 -- constant_term - create constant term.
1333 --
1334 -- This routine creates the linear form, which is a constant term. */
1336 FORMULA *constant_term(MPL *mpl, double coef)
1337 { FORMULA *form;
1338 if (coef == 0.0)
1339 form = NULL;
1340 else
1341 { form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1342 form->coef = coef;
1343 form->var = NULL;
1344 form->next = NULL;
1345 }
1346 return form;
1347 }
1349 /*----------------------------------------------------------------------
1350 -- single_variable - create single variable.
1351 --
1352 -- This routine creates the linear form, which is a single elemental
1353 -- variable. */
1355 FORMULA *single_variable
1356 ( MPL *mpl,
1357 ELEMVAR *var /* referenced */
1358 )
1359 { FORMULA *form;
1360 xassert(var != NULL);
1361 form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1362 form->coef = 1.0;
1363 form->var = var;
1364 form->next = NULL;
1365 return form;
1366 }
1368 /*----------------------------------------------------------------------
1369 -- copy_formula - make copy of linear form.
1370 --
1371 -- This routine returns an exact copy of linear form. */
1373 FORMULA *copy_formula
1374 ( MPL *mpl,
1375 FORMULA *form /* not changed */
1376 )
1377 { FORMULA *head, *tail;
1378 if (form == NULL)
1379 head = NULL;
1380 else
1381 { head = tail = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1382 for (; form != NULL; form = form->next)
1383 { tail->coef = form->coef;
1384 tail->var = form->var;
1385 if (form->next != NULL)
1386 tail = (tail->next = dmp_get_atom(mpl->formulae, sizeof(FORMULA)));
1387 }
1388 tail->next = NULL;
1389 }
1390 return head;
1391 }
1393 /*----------------------------------------------------------------------
1394 -- delete_formula - delete linear form.
1395 --
1396 -- This routine deletes specified linear form. */
1398 void delete_formula
1399 ( MPL *mpl,
1400 FORMULA *form /* destroyed */
1401 )
1402 { FORMULA *temp;
1403 while (form != NULL)
1404 { temp = form;
1405 form = form->next;
1406 dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));
1407 }
1408 return;
1409 }
1411 /*----------------------------------------------------------------------
1412 -- linear_comb - linear combination of two linear forms.
1413 --
1414 -- This routine computes the linear combination:
1415 --
1416 -- a * fx + b * fy,
1417 --
1418 -- where a and b are numeric coefficients, fx and fy are linear forms
1419 -- (destroyed on exit). */
1421 FORMULA *linear_comb
1422 ( MPL *mpl,
1423 double a, FORMULA *fx, /* destroyed */
1424 double b, FORMULA *fy /* destroyed */
1425 )
1426 { FORMULA *form = NULL, *term, *temp;
1427 double c0 = 0.0;
1428 for (term = fx; term != NULL; term = term->next)
1429 { if (term->var == NULL)
1430 c0 = fp_add(mpl, c0, fp_mul(mpl, a, term->coef));
1431 else
1432 term->var->temp =
1433 fp_add(mpl, term->var->temp, fp_mul(mpl, a, term->coef));
1434 }
1435 for (term = fy; term != NULL; term = term->next)
1436 { if (term->var == NULL)
1437 c0 = fp_add(mpl, c0, fp_mul(mpl, b, term->coef));
1438 else
1439 term->var->temp =
1440 fp_add(mpl, term->var->temp, fp_mul(mpl, b, term->coef));
1441 }
1442 for (term = fx; term != NULL; term = term->next)
1443 { if (term->var != NULL && term->var->temp != 0.0)
1444 { temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1445 temp->coef = term->var->temp, temp->var = term->var;
1446 temp->next = form, form = temp;
1447 term->var->temp = 0.0;
1448 }
1449 }
1450 for (term = fy; term != NULL; term = term->next)
1451 { if (term->var != NULL && term->var->temp != 0.0)
1452 { temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1453 temp->coef = term->var->temp, temp->var = term->var;
1454 temp->next = form, form = temp;
1455 term->var->temp = 0.0;
1456 }
1457 }
1458 if (c0 != 0.0)
1459 { temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));
1460 temp->coef = c0, temp->var = NULL;
1461 temp->next = form, form = temp;
1462 }
1463 delete_formula(mpl, fx);
1464 delete_formula(mpl, fy);
1465 return form;
1466 }
1468 /*----------------------------------------------------------------------
1469 -- remove_constant - remove constant term from linear form.
1470 --
1471 -- This routine removes constant term from linear form and stores its
1472 -- value to given location. */
1474 FORMULA *remove_constant
1475 ( MPL *mpl,
1476 FORMULA *form, /* destroyed */
1477 double *coef /* modified */
1478 )
1479 { FORMULA *head = NULL, *temp;
1480 *coef = 0.0;
1481 while (form != NULL)
1482 { temp = form;
1483 form = form->next;
1484 if (temp->var == NULL)
1485 { /* constant term */
1486 *coef = fp_add(mpl, *coef, temp->coef);
1487 dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));
1488 }
1489 else
1490 { /* linear term */
1491 temp->next = head;
1492 head = temp;
1493 }
1494 }
1495 return head;
1496 }
1498 /*----------------------------------------------------------------------
1499 -- reduce_terms - reduce identical terms in linear form.
1500 --
1501 -- This routine reduces identical terms in specified linear form. */
1503 FORMULA *reduce_terms
1504 ( MPL *mpl,
1505 FORMULA *form /* destroyed */
1506 )
1507 { FORMULA *term, *next_term;
1508 double c0 = 0.0;
1509 for (term = form; term != NULL; term = term->next)
1510 { if (term->var == NULL)
1511 c0 = fp_add(mpl, c0, term->coef);
1512 else
1513 term->var->temp = fp_add(mpl, term->var->temp, term->coef);
1514 }
1515 next_term = form, form = NULL;
1516 for (term = next_term; term != NULL; term = next_term)
1517 { next_term = term->next;
1518 if (term->var == NULL && c0 != 0.0)
1519 { term->coef = c0, c0 = 0.0;
1520 term->next = form, form = term;
1521 }
1522 else if (term->var != NULL && term->var->temp != 0.0)
1523 { term->coef = term->var->temp, term->var->temp = 0.0;
1524 term->next = form, form = term;
1525 }
1526 else
1527 dmp_free_atom(mpl->formulae, term, sizeof(FORMULA));
1528 }
1529 return form;
1530 }
1532 /**********************************************************************/
1533 /* * * ELEMENTAL CONSTRAINTS * * */
1534 /**********************************************************************/
1536 /* (there are no specific routines for elemental constraints) */
1538 /**********************************************************************/
1539 /* * * GENERIC VALUES * * */
1540 /**********************************************************************/
1542 /*----------------------------------------------------------------------
1543 -- delete_value - delete generic value.
1544 --
1545 -- This routine deletes specified generic value.
1546 --
1547 -- NOTE: The generic value to be deleted must be valid. */
1549 void delete_value
1550 ( MPL *mpl,
1551 int type,
1552 VALUE *value /* content destroyed */
1553 )
1554 { xassert(value != NULL);
1555 switch (type)
1556 { case A_NONE:
1557 value->none = NULL;
1558 break;
1559 case A_NUMERIC:
1560 value->num = 0.0;
1561 break;
1562 case A_SYMBOLIC:
1563 delete_symbol(mpl, value->sym), value->sym = NULL;
1564 break;
1565 case A_LOGICAL:
1566 value->bit = 0;
1567 break;
1568 case A_TUPLE:
1569 delete_tuple(mpl, value->tuple), value->tuple = NULL;
1570 break;
1571 case A_ELEMSET:
1572 delete_elemset(mpl, value->set), value->set = NULL;
1573 break;
1574 case A_ELEMVAR:
1575 value->var = NULL;
1576 break;
1577 case A_FORMULA:
1578 delete_formula(mpl, value->form), value->form = NULL;
1579 break;
1580 case A_ELEMCON:
1581 value->con = NULL;
1582 break;
1583 default:
1584 xassert(type != type);
1585 }
1586 return;
1587 }
1589 /**********************************************************************/
1590 /* * * SYMBOLICALLY INDEXED ARRAYS * * */
1591 /**********************************************************************/
1593 /*----------------------------------------------------------------------
1594 -- create_array - create array.
1595 --
1596 -- This routine creates an array of specified type and dimension. Being
1597 -- created the array is initially empty.
1598 --
1599 -- The type indicator determines generic values, which can be assigned
1600 -- to the array members:
1601 --
1602 -- A_NONE - none (members have no assigned values)
1603 -- A_NUMERIC - floating-point numbers
1604 -- A_SYMBOLIC - symbols
1605 -- A_ELEMSET - elemental sets
1606 -- A_ELEMVAR - elemental variables
1607 -- A_ELEMCON - elemental constraints
1608 --
1609 -- The dimension may be 0, in which case the array consists of the only
1610 -- member (such arrays represent 0-dimensional objects). */
1612 ARRAY *create_array(MPL *mpl, int type, int dim)
1613 { ARRAY *array;
1614 xassert(type == A_NONE || type == A_NUMERIC ||
1615 type == A_SYMBOLIC || type == A_ELEMSET ||
1616 type == A_ELEMVAR || type == A_ELEMCON);
1617 xassert(dim >= 0);
1618 array = dmp_get_atom(mpl->arrays, sizeof(ARRAY));
1619 array->type = type;
1620 array->dim = dim;
1621 array->size = 0;
1622 array->head = NULL;
1623 array->tail = NULL;
1624 array->tree = NULL;
1625 array->prev = NULL;
1626 array->next = mpl->a_list;
1627 /* include the array in the global array list */
1628 if (array->next != NULL) array->next->prev = array;
1629 mpl->a_list = array;
1630 return array;
1631 }
1633 /*----------------------------------------------------------------------
1634 -- find_member - find array member with given n-tuple.
1635 --
1636 -- This routine finds an array member, which has given n-tuple. If the
1637 -- array is short, the linear search is used. Otherwise the routine
1638 -- autimatically creates the search tree (i.e. the array index) to find
1639 -- members for logarithmic time. */
1641 static int compare_member_tuples(void *info, const void *key1,
1642 const void *key2)
1643 { /* this is an auxiliary routine used to compare keys, which are
1644 n-tuples assigned to array members */
1645 return compare_tuples((MPL *)info, (TUPLE *)key1, (TUPLE *)key2);
1646 }
1648 MEMBER *find_member
1649 ( MPL *mpl,
1650 ARRAY *array, /* not changed */
1651 TUPLE *tuple /* not changed */
1652 )
1653 { MEMBER *memb;
1654 xassert(array != NULL);
1655 /* the n-tuple must have the same dimension as the array */
1656 xassert(tuple_dimen(mpl, tuple) == array->dim);
1657 /* if the array is large enough, create the search tree and index
1658 all existing members of the array */
1659 if (array->size > 30 && array->tree == NULL)
1660 { array->tree = avl_create_tree(compare_member_tuples, mpl);
1661 for (memb = array->head; memb != NULL; memb = memb->next)
1662 avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
1663 (void *)memb);
1664 }
1665 /* find a member, which has the given tuple */
1666 if (array->tree == NULL)
1667 { /* the search tree doesn't exist; use the linear search */
1668 for (memb = array->head; memb != NULL; memb = memb->next)
1669 if (compare_tuples(mpl, memb->tuple, tuple) == 0) break;
1670 }
1671 else
1672 { /* the search tree exists; use the binary search */
1673 AVLNODE *node;
1674 node = avl_find_node(array->tree, tuple);
1675 memb = (MEMBER *)(node == NULL ? NULL : avl_get_node_link(node));
1676 }
1677 return memb;
1678 }
1680 /*----------------------------------------------------------------------
1681 -- add_member - add new member to array.
1682 --
1683 -- This routine creates a new member with given n-tuple and adds it to
1684 -- specified array.
1685 --
1686 -- For the sake of efficiency this routine doesn't check whether the
1687 -- array already contains a member with the given n-tuple or not. Thus,
1688 -- if necessary, the calling program should use the routine find_member
1689 -- in order to be sure that the array contains no member with the same
1690 -- n-tuple, because members with duplicate n-tuples are not allowed.
1691 --
1692 -- This routine assigns no generic value to the new member, because the
1693 -- calling program must do that. */
1695 MEMBER *add_member
1696 ( MPL *mpl,
1697 ARRAY *array, /* modified */
1698 TUPLE *tuple /* destroyed */
1699 )
1700 { MEMBER *memb;
1701 xassert(array != NULL);
1702 /* the n-tuple must have the same dimension as the array */
1703 xassert(tuple_dimen(mpl, tuple) == array->dim);
1704 /* create new member */
1705 memb = dmp_get_atom(mpl->members, sizeof(MEMBER));
1706 memb->tuple = tuple;
1707 memb->next = NULL;
1708 memset(&memb->value, '?', sizeof(VALUE));
1709 /* and append it to the member list */
1710 array->size++;
1711 if (array->head == NULL)
1712 array->head = memb;
1713 else
1714 array->tail->next = memb;
1715 array->tail = memb;
1716 /* if the search tree exists, index the new member */
1717 if (array->tree != NULL)
1718 avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
1719 (void *)memb);
1720 return memb;
1721 }
1723 /*----------------------------------------------------------------------
1724 -- delete_array - delete array.
1725 --
1726 -- This routine deletes specified array.
1727 --
1728 -- Generic values assigned to the array members are not deleted by this
1729 -- routine. The calling program itself must delete all assigned generic
1730 -- values before deleting the array. */
1732 void delete_array
1733 ( MPL *mpl,
1734 ARRAY *array /* destroyed */
1735 )
1736 { MEMBER *memb;
1737 xassert(array != NULL);
1738 /* delete all existing array members */
1739 while (array->head != NULL)
1740 { memb = array->head;
1741 array->head = memb->next;
1742 delete_tuple(mpl, memb->tuple);
1743 dmp_free_atom(mpl->members, memb, sizeof(MEMBER));
1744 }
1745 /* if the search tree exists, also delete it */
1746 if (array->tree != NULL) avl_delete_tree(array->tree);
1747 /* remove the array from the global array list */
1748 if (array->prev == NULL)
1749 mpl->a_list = array->next;
1750 else
1751 array->prev->next = array->next;
1752 if (array->next == NULL)
1753 ;
1754 else
1755 array->next->prev = array->prev;
1756 /* delete the array descriptor */
1757 dmp_free_atom(mpl->arrays, array, sizeof(ARRAY));
1758 return;
1759 }
1761 /**********************************************************************/
1762 /* * * DOMAINS AND DUMMY INDICES * * */
1763 /**********************************************************************/
1765 /*----------------------------------------------------------------------
1766 -- assign_dummy_index - assign new value to dummy index.
1767 --
1768 -- This routine assigns new value to specified dummy index and, that is
1769 -- important, invalidates all temporary resultant values, which depends
1770 -- on that dummy index. */
1772 void assign_dummy_index
1773 ( MPL *mpl,
1774 DOMAIN_SLOT *slot, /* modified */
1775 SYMBOL *value /* not changed */
1776 )
1777 { CODE *leaf, *code;
1778 xassert(slot != NULL);
1779 xassert(value != NULL);
1780 /* delete the current value assigned to the dummy index */
1781 if (slot->value != NULL)
1782 { /* if the current value and the new one are identical, actual
1783 assignment is not needed */
1784 if (compare_symbols(mpl, slot->value, value) == 0) goto done;
1785 /* delete a symbol, which is the current value */
1786 delete_symbol(mpl, slot->value), slot->value = NULL;
1787 }
1788 /* now walk through all the pseudo-codes with op = O_INDEX, which
1789 refer to the dummy index to be changed (these pseudo-codes are
1790 leaves in the forest of *all* expressions in the database) */
1791 for (leaf = slot->list; leaf != NULL; leaf = leaf->arg.index.
1792 next)
1793 { xassert(leaf->op == O_INDEX);
1794 /* invalidate all resultant values, which depend on the dummy
1795 index, walking from the current leaf toward the root of the
1796 corresponding expression tree */
1797 for (code = leaf; code != NULL; code = code->up)
1798 { if (code->valid)
1799 { /* invalidate and delete resultant value */
1800 code->valid = 0;
1801 delete_value(mpl, code->type, &code->value);
1802 }
1803 }
1804 }
1805 /* assign new value to the dummy index */
1806 slot->value = copy_symbol(mpl, value);
1807 done: return;
1808 }
1810 /*----------------------------------------------------------------------
1811 -- update_dummy_indices - update current values of dummy indices.
1812 --
1813 -- This routine assigns components of "backup" n-tuple to dummy indices
1814 -- of specified domain block. If no "backup" n-tuple is defined for the
1815 -- domain block, values of the dummy indices remain untouched. */
1817 void update_dummy_indices
1818 ( MPL *mpl,
1819 DOMAIN_BLOCK *block /* not changed */
1820 )
1821 { DOMAIN_SLOT *slot;
1822 TUPLE *temp;
1823 if (block->backup != NULL)
1824 { for (slot = block->list, temp = block->backup; slot != NULL;
1825 slot = slot->next, temp = temp->next)
1826 { xassert(temp != NULL);
1827 xassert(temp->sym != NULL);
1828 assign_dummy_index(mpl, slot, temp->sym);
1829 }
1830 }
1831 return;
1832 }
1834 /*----------------------------------------------------------------------
1835 -- enter_domain_block - enter domain block.
1836 --
1837 -- Let specified domain block have the form:
1838 --
1839 -- { ..., (j1, j2, ..., jn) in J, ... }
1840 --
1841 -- where j1, j2, ..., jn are dummy indices, J is a basic set.
1842 --
1843 -- This routine does the following:
1844 --
1845 -- 1. Checks if the given n-tuple is a member of the basic set J. Note
1846 -- that J being *out of the scope* of the domain block cannot depend
1847 -- on the dummy indices in the same and inner domain blocks, so it
1848 -- can be computed before the dummy indices are assigned new values.
1849 -- If this check fails, the routine returns with non-zero code.
1850 --
1851 -- 2. Saves current values of the dummy indices j1, j2, ..., jn.
1852 --
1853 -- 3. Assigns new values, which are components of the given n-tuple, to
1854 -- the dummy indices j1, j2, ..., jn. If dimension of the n-tuple is
1855 -- larger than n, its extra components n+1, n+2, ... are not used.
1856 --
1857 -- 4. Calls the formal routine func which either enters the next domain
1858 -- block or evaluates some code within the domain scope.
1859 --
1860 -- 5. Restores former values of the dummy indices j1, j2, ..., jn.
1861 --
1862 -- Since current values assigned to the dummy indices on entry to this
1863 -- routine are restored on exit, the formal routine func is allowed to
1864 -- call this routine recursively. */
1866 int enter_domain_block
1867 ( MPL *mpl,
1868 DOMAIN_BLOCK *block, /* not changed */
1869 TUPLE *tuple, /* not changed */
1870 void *info, void (*func)(MPL *mpl, void *info)
1871 )
1872 { TUPLE *backup;
1873 int ret = 0;
1874 /* check if the given n-tuple is a member of the basic set */
1875 xassert(block->code != NULL);
1876 if (!is_member(mpl, block->code, tuple))
1877 { ret = 1;
1878 goto done;
1879 }
1880 /* save reference to "backup" n-tuple, which was used to assign
1881 current values of the dummy indices (it is sufficient to save
1882 reference, not value, because that n-tuple is defined in some
1883 outer level of recursion and therefore cannot be changed on
1884 this and deeper recursive calls) */
1885 backup = block->backup;
1886 /* set up new "backup" n-tuple, which defines new values of the
1887 dummy indices */
1888 block->backup = tuple;
1889 /* assign new values to the dummy indices */
1890 update_dummy_indices(mpl, block);
1891 /* call the formal routine that does the rest part of the job */
1892 func(mpl, info);
1893 /* restore reference to the former "backup" n-tuple */
1894 block->backup = backup;
1895 /* restore former values of the dummy indices; note that if the
1896 domain block just escaped has no other active instances which
1897 may exist due to recursion (it is indicated by a null pointer
1898 to the former n-tuple), former values of the dummy indices are
1899 undefined; therefore in this case the routine keeps currently
1900 assigned values of the dummy indices that involves keeping all
1901 dependent temporary results and thereby, if this domain block
1902 is not used recursively, allows improving efficiency */
1903 update_dummy_indices(mpl, block);
1904 done: return ret;
1905 }
1907 /*----------------------------------------------------------------------
1908 -- eval_within_domain - perform evaluation within domain scope.
1909 --
1910 -- This routine assigns new values (symbols) to all dummy indices of
1911 -- specified domain and calls the formal routine func, which is used to
1912 -- evaluate some code in the domain scope. Each free dummy index in the
1913 -- domain is assigned a value specified in the corresponding component
1914 -- of given n-tuple. Non-free dummy indices are assigned values, which
1915 -- are computed by this routine.
1916 --
1917 -- Number of components in the given n-tuple must be the same as number
1918 -- of free indices in the domain.
1919 --
1920 -- If the given n-tuple is not a member of the domain set, the routine
1921 -- func is not called, and non-zero code is returned.
1922 --
1923 -- For the sake of convenience it is allowed to specify domain as NULL
1924 -- (then n-tuple also must be 0-tuple, i.e. empty), in which case this
1925 -- routine just calls the routine func and returns zero.
1926 --
1927 -- This routine allows recursive calls from the routine func providing
1928 -- correct values of dummy indices for each instance.
1929 --
1930 -- NOTE: The n-tuple passed to this routine must not be changed by any
1931 -- other routines called from the formal routine func until this
1932 -- routine has returned. */
1934 struct eval_domain_info
1935 { /* working info used by the routine eval_within_domain */
1936 DOMAIN *domain;
1937 /* domain, which has to be entered */
1938 DOMAIN_BLOCK *block;
1939 /* domain block, which is currently processed */
1940 TUPLE *tuple;
1941 /* tail of original n-tuple, whose components have to be assigned
1942 to free dummy indices in the current domain block */
1943 void *info;
1944 /* transit pointer passed to the formal routine func */
1945 void (*func)(MPL *mpl, void *info);
1946 /* routine, which has to be executed in the domain scope */
1947 int failure;
1948 /* this flag indicates that given n-tuple is not a member of the
1949 domain set */
1950 };
1952 static void eval_domain_func(MPL *mpl, void *_my_info)
1953 { /* this routine recursively enters into the domain scope and then
1954 calls the routine func */
1955 struct eval_domain_info *my_info = _my_info;
1956 if (my_info->block != NULL)
1957 { /* the current domain block to be entered exists */
1958 DOMAIN_BLOCK *block;
1959 DOMAIN_SLOT *slot;
1960 TUPLE *tuple = NULL, *temp = NULL;
1961 /* save pointer to the current domain block */
1962 block = my_info->block;
1963 /* and get ready to enter the next block (if it exists) */
1964 my_info->block = block->next;
1965 /* construct temporary n-tuple, whose components correspond to
1966 dummy indices (slots) of the current domain; components of
1967 the temporary n-tuple that correspond to free dummy indices
1968 are assigned references (not values!) to symbols specified
1969 in the corresponding components of the given n-tuple, while
1970 other components that correspond to non-free dummy indices
1971 are assigned symbolic values computed here */
1972 for (slot = block->list; slot != NULL; slot = slot->next)
1973 { /* create component that corresponds to the current slot */
1974 if (tuple == NULL)
1975 tuple = temp = dmp_get_atom(mpl->tuples, sizeof(TUPLE));
1976 else
1977 temp = (temp->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));
1978 if (slot->code == NULL)
1979 { /* dummy index is free; take reference to symbol, which
1980 is specified in the corresponding component of given
1981 n-tuple */
1982 xassert(my_info->tuple != NULL);
1983 temp->sym = my_info->tuple->sym;
1984 xassert(temp->sym != NULL);
1985 my_info->tuple = my_info->tuple->next;
1986 }
1987 else
1988 { /* dummy index is non-free; compute symbolic value to be
1989 temporarily assigned to the dummy index */
1990 temp->sym = eval_symbolic(mpl, slot->code);
1991 }
1992 }
1993 temp->next = NULL;
1994 /* enter the current domain block */
1995 if (enter_domain_block(mpl, block, tuple, my_info,
1996 eval_domain_func)) my_info->failure = 1;
1997 /* delete temporary n-tuple as well as symbols that correspond
1998 to non-free dummy indices (they were computed here) */
1999 for (slot = block->list; slot != NULL; slot = slot->next)
2000 { xassert(tuple != NULL);
2001 temp = tuple;
2002 tuple = tuple->next;
2003 if (slot->code != NULL)
2004 { /* dummy index is non-free; delete symbolic value */
2005 delete_symbol(mpl, temp->sym);
2006 }
2007 /* delete component that corresponds to the current slot */
2008 dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));
2009 }
2010 }
2011 else
2012 { /* there are no more domain blocks, i.e. we have reached the
2013 domain scope */
2014 xassert(my_info->tuple == NULL);
2015 /* check optional predicate specified for the domain */
2016 if (my_info->domain->code != NULL && !eval_logical(mpl,
2017 my_info->domain->code))
2018 { /* the predicate is false */
2019 my_info->failure = 2;
2020 }
2021 else
2022 { /* the predicate is true; do the job */
2023 my_info->func(mpl, my_info->info);
2024 }
2025 }
2026 return;
2027 }
2029 int eval_within_domain
2030 ( MPL *mpl,
2031 DOMAIN *domain, /* not changed */
2032 TUPLE *tuple, /* not changed */
2033 void *info, void (*func)(MPL *mpl, void *info)
2034 )
2035 { /* this routine performs evaluation within domain scope */
2036 struct eval_domain_info _my_info, *my_info = &_my_info;
2037 if (domain == NULL)
2038 { xassert(tuple == NULL);
2039 func(mpl, info);
2040 my_info->failure = 0;
2041 }
2042 else
2043 { xassert(tuple != NULL);
2044 my_info->domain = domain;
2045 my_info->block = domain->list;
2046 my_info->tuple = tuple;
2047 my_info->info = info;
2048 my_info->func = func;
2049 my_info->failure = 0;
2050 /* enter the very first domain block */
2051 eval_domain_func(mpl, my_info);
2052 }
2053 return my_info->failure;
2054 }
2056 /*----------------------------------------------------------------------
2057 -- loop_within_domain - perform iterations within domain scope.
2058 --
2059 -- This routine iteratively assigns new values (symbols) to the dummy
2060 -- indices of specified domain by enumerating all n-tuples, which are
2061 -- members of the domain set, and for every n-tuple it calls the formal
2062 -- routine func to evaluate some code within the domain scope.
2063 --
2064 -- If the routine func returns non-zero, enumeration within the domain
2065 -- is prematurely terminated.
2066 --
2067 -- For the sake of convenience it is allowed to specify domain as NULL,
2068 -- in which case this routine just calls the routine func only once and
2069 -- returns zero.
2070 --
2071 -- This routine allows recursive calls from the routine func providing
2072 -- correct values of dummy indices for each instance. */
2074 struct loop_domain_info
2075 { /* working info used by the routine loop_within_domain */
2076 DOMAIN *domain;
2077 /* domain, which has to be entered */
2078 DOMAIN_BLOCK *block;
2079 /* domain block, which is currently processed */
2080 int looping;
2081 /* clearing this flag leads to terminating enumeration */
2082 void *info;
2083 /* transit pointer passed to the formal routine func */
2084 int (*func)(MPL *mpl, void *info);
2085 /* routine, which needs to be executed in the domain scope */
2086 };
2088 static void loop_domain_func(MPL *mpl, void *_my_info)
2089 { /* this routine enumerates all n-tuples in the basic set of the
2090 current domain block, enters recursively into the domain scope
2091 for every n-tuple, and then calls the routine func */
2092 struct loop_domain_info *my_info = _my_info;
2093 if (my_info->block != NULL)
2094 { /* the current domain block to be entered exists */
2095 DOMAIN_BLOCK *block;
2096 DOMAIN_SLOT *slot;
2097 TUPLE *bound;
2098 /* save pointer to the current domain block */
2099 block = my_info->block;
2100 /* and get ready to enter the next block (if it exists) */
2101 my_info->block = block->next;
2102 /* compute symbolic values, at which non-free dummy indices of
2103 the current domain block are bound; since that values don't
2104 depend on free dummy indices of the current block, they can
2105 be computed once out of the enumeration loop */
2106 bound = create_tuple(mpl);
2107 for (slot = block->list; slot != NULL; slot = slot->next)
2108 { if (slot->code != NULL)
2109 bound = expand_tuple(mpl, bound, eval_symbolic(mpl,
2110 slot->code));
2111 }
2112 /* start enumeration */
2113 xassert(block->code != NULL);
2114 if (block->code->op == O_DOTS)
2115 { /* the basic set is "arithmetic", in which case it doesn't
2116 need to be computed explicitly */
2117 TUPLE *tuple;
2118 int n, j;
2119 double t0, tf, dt;
2120 /* compute "parameters" of the basic set */
2121 t0 = eval_numeric(mpl, block->code->arg.arg.x);
2122 tf = eval_numeric(mpl, block->code->arg.arg.y);
2123 if (block->code->arg.arg.z == NULL)
2124 dt = 1.0;
2125 else
2126 dt = eval_numeric(mpl, block->code->arg.arg.z);
2127 /* determine cardinality of the basic set */
2128 n = arelset_size(mpl, t0, tf, dt);
2129 /* create dummy 1-tuple for members of the basic set */
2130 tuple = expand_tuple(mpl, create_tuple(mpl),
2131 create_symbol_num(mpl, 0.0));
2132 /* in case of "arithmetic" set there is exactly one dummy
2133 index, which cannot be non-free */
2134 xassert(bound == NULL);
2135 /* walk through 1-tuples of the basic set */
2136 for (j = 1; j <= n && my_info->looping; j++)
2137 { /* construct dummy 1-tuple for the current member */
2138 tuple->sym->num = arelset_member(mpl, t0, tf, dt, j);
2139 /* enter the current domain block */
2140 enter_domain_block(mpl, block, tuple, my_info,
2141 loop_domain_func);
2142 }
2143 /* delete dummy 1-tuple */
2144 delete_tuple(mpl, tuple);
2145 }
2146 else
2147 { /* the basic set is of general kind, in which case it needs
2148 to be explicitly computed */
2149 ELEMSET *set;
2150 MEMBER *memb;
2151 TUPLE *temp1, *temp2;
2152 /* compute the basic set */
2153 set = eval_elemset(mpl, block->code);
2154 /* walk through all n-tuples of the basic set */
2155 for (memb = set->head; memb != NULL && my_info->looping;
2156 memb = memb->next)
2157 { /* all components of the current n-tuple that correspond
2158 to non-free dummy indices must be feasible; otherwise
2159 the n-tuple is not in the basic set */
2160 temp1 = memb->tuple;
2161 temp2 = bound;
2162 for (slot = block->list; slot != NULL; slot = slot->next)
2163 { xassert(temp1 != NULL);
2164 if (slot->code != NULL)
2165 { /* non-free dummy index */
2166 xassert(temp2 != NULL);
2167 if (compare_symbols(mpl, temp1->sym, temp2->sym)
2168 != 0)
2169 { /* the n-tuple is not in the basic set */
2170 goto skip;
2171 }
2172 temp2 = temp2->next;
2173 }
2174 temp1 = temp1->next;
2175 }
2176 xassert(temp1 == NULL);
2177 xassert(temp2 == NULL);
2178 /* enter the current domain block */
2179 enter_domain_block(mpl, block, memb->tuple, my_info,
2180 loop_domain_func);
2181 skip: ;
2182 }
2183 /* delete the basic set */
2184 delete_elemset(mpl, set);
2185 }
2186 /* delete symbolic values binding non-free dummy indices */
2187 delete_tuple(mpl, bound);
2188 /* restore pointer to the current domain block */
2189 my_info->block = block;
2190 }
2191 else
2192 { /* there are no more domain blocks, i.e. we have reached the
2193 domain scope */
2194 /* check optional predicate specified for the domain */
2195 if (my_info->domain->code != NULL && !eval_logical(mpl,
2196 my_info->domain->code))
2197 { /* the predicate is false */
2198 /* nop */;
2199 }
2200 else
2201 { /* the predicate is true; do the job */
2202 my_info->looping = !my_info->func(mpl, my_info->info);
2203 }
2204 }
2205 return;
2206 }
2208 void loop_within_domain
2209 ( MPL *mpl,
2210 DOMAIN *domain, /* not changed */
2211 void *info, int (*func)(MPL *mpl, void *info)
2212 )
2213 { /* this routine performs iterations within domain scope */
2214 struct loop_domain_info _my_info, *my_info = &_my_info;
2215 if (domain == NULL)
2216 func(mpl, info);
2217 else
2218 { my_info->domain = domain;
2219 my_info->block = domain->list;
2220 my_info->looping = 1;
2221 my_info->info = info;
2222 my_info->func = func;
2223 /* enter the very first domain block */
2224 loop_domain_func(mpl, my_info);
2225 }
2226 return;
2227 }
2229 /*----------------------------------------------------------------------
2230 -- out_of_domain - raise domain exception.
2231 --
2232 -- This routine is called when a reference is made to a member of some
2233 -- model object, but its n-tuple is out of the object domain. */
2235 void out_of_domain
2236 ( MPL *mpl,
2237 char *name, /* not changed */
2238 TUPLE *tuple /* not changed */
2239 )
2240 { xassert(name != NULL);
2241 xassert(tuple != NULL);
2242 error(mpl, "%s%s out of domain", name, format_tuple(mpl, '[',
2243 tuple));
2244 /* no return */
2245 }
2247 /*----------------------------------------------------------------------
2248 -- get_domain_tuple - obtain current n-tuple from domain.
2249 --
2250 -- This routine constructs n-tuple, whose components are current values
2251 -- assigned to *free* dummy indices of specified domain.
2252 --
2253 -- For the sake of convenience it is allowed to specify domain as NULL,
2254 -- in which case this routine returns 0-tuple.
2255 --
2256 -- NOTE: This routine must not be called out of domain scope. */
2258 TUPLE *get_domain_tuple
2259 ( MPL *mpl,
2260 DOMAIN *domain /* not changed */
2261 )
2262 { DOMAIN_BLOCK *block;
2263 DOMAIN_SLOT *slot;
2264 TUPLE *tuple;
2265 tuple = create_tuple(mpl);
2266 if (domain != NULL)
2267 { for (block = domain->list; block != NULL; block = block->next)
2268 { for (slot = block->list; slot != NULL; slot = slot->next)
2269 { if (slot->code == NULL)
2270 { xassert(slot->value != NULL);
2271 tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,
2272 slot->value));
2273 }
2274 }
2275 }
2276 }
2277 return tuple;
2278 }
2280 /*----------------------------------------------------------------------
2281 -- clean_domain - clean domain.
2282 --
2283 -- This routine cleans specified domain that assumes deleting all stuff
2284 -- dynamically allocated during the generation phase. */
2286 void clean_domain(MPL *mpl, DOMAIN *domain)
2287 { DOMAIN_BLOCK *block;
2288 DOMAIN_SLOT *slot;
2289 /* if no domain is specified, do nothing */
2290 if (domain == NULL) goto done;
2291 /* clean all domain blocks */
2292 for (block = domain->list; block != NULL; block = block->next)
2293 { /* clean all domain slots */
2294 for (slot = block->list; slot != NULL; slot = slot->next)
2295 { /* clean pseudo-code for computing bound value */
2296 clean_code(mpl, slot->code);
2297 /* delete symbolic value assigned to dummy index */
2298 if (slot->value != NULL)
2299 delete_symbol(mpl, slot->value), slot->value = NULL;
2300 }
2301 /* clean pseudo-code for computing basic set */
2302 clean_code(mpl, block->code);
2303 }
2304 /* clean pseudo-code for computing domain predicate */
2305 clean_code(mpl, domain->code);
2306 done: return;
2307 }
2309 /**********************************************************************/
2310 /* * * MODEL SETS * * */
2311 /**********************************************************************/
2313 /*----------------------------------------------------------------------
2314 -- check_elem_set - check elemental set assigned to set member.
2315 --
2316 -- This routine checks if given elemental set being assigned to member
2317 -- of specified model set satisfies to all restrictions.
2318 --
2319 -- NOTE: This routine must not be called out of domain scope. */
2321 void check_elem_set
2322 ( MPL *mpl,
2323 SET *set, /* not changed */
2324 TUPLE *tuple, /* not changed */
2325 ELEMSET *refer /* not changed */
2326 )
2327 { WITHIN *within;
2328 MEMBER *memb;
2329 int eqno;
2330 /* elemental set must be within all specified supersets */
2331 for (within = set->within, eqno = 1; within != NULL; within =
2332 within->next, eqno++)
2333 { xassert(within->code != NULL);
2334 for (memb = refer->head; memb != NULL; memb = memb->next)
2335 { if (!is_member(mpl, within->code, memb->tuple))
2336 { char buf[255+1];
2337 strcpy(buf, format_tuple(mpl, '(', memb->tuple));
2338 xassert(strlen(buf) < sizeof(buf));
2339 error(mpl, "%s%s contains %s which not within specified "
2340 "set; see (%d)", set->name, format_tuple(mpl, '[',
2341 tuple), buf, eqno);
2342 }
2343 }
2344 }
2345 return;
2346 }
2348 /*----------------------------------------------------------------------
2349 -- take_member_set - obtain elemental set assigned to set member.
2350 --
2351 -- This routine obtains a reference to elemental set assigned to given
2352 -- member of specified model set and returns it on exit.
2353 --
2354 -- NOTE: This routine must not be called out of domain scope. */
2356 ELEMSET *take_member_set /* returns reference, not value */
2357 ( MPL *mpl,
2358 SET *set, /* not changed */
2359 TUPLE *tuple /* not changed */
2360 )
2361 { MEMBER *memb;
2362 ELEMSET *refer;
2363 /* find member in the set array */
2364 memb = find_member(mpl, set->array, tuple);
2365 if (memb != NULL)
2366 { /* member exists, so just take the reference */
2367 refer = memb->value.set;
2368 }
2369 else if (set->assign != NULL)
2370 { /* compute value using assignment expression */
2371 refer = eval_elemset(mpl, set->assign);
2372 add: /* check that the elemental set satisfies to all restrictions,
2373 assign it to new member, and add the member to the array */
2374 check_elem_set(mpl, set, tuple, refer);
2375 memb = add_member(mpl, set->array, copy_tuple(mpl, tuple));
2376 memb->value.set = refer;
2377 }
2378 else if (set->option != NULL)
2379 { /* compute default elemental set */
2380 refer = eval_elemset(mpl, set->option);
2381 goto add;
2382 }
2383 else
2384 { /* no value (elemental set) is provided */
2385 error(mpl, "no value for %s%s", set->name, format_tuple(mpl,
2386 '[', tuple));
2387 }
2388 return refer;
2389 }
2391 /*----------------------------------------------------------------------
2392 -- eval_member_set - evaluate elemental set assigned to set member.
2393 --
2394 -- This routine evaluates a reference to elemental set assigned to given
2395 -- member of specified model set and returns it on exit. */
2397 struct eval_set_info
2398 { /* working info used by the routine eval_member_set */
2399 SET *set;
2400 /* model set */
2401 TUPLE *tuple;
2402 /* n-tuple, which defines set member */
2403 MEMBER *memb;
2404 /* normally this pointer is NULL; the routine uses this pointer
2405 to check data provided in the data section, in which case it
2406 points to a member currently checked; this check is performed
2407 automatically only once when a reference to any member occurs
2408 for the first time */
2409 ELEMSET *refer;
2410 /* evaluated reference to elemental set */
2411 };
2413 static void eval_set_func(MPL *mpl, void *_info)
2414 { /* this is auxiliary routine to work within domain scope */
2415 struct eval_set_info *info = _info;
2416 if (info->memb != NULL)
2417 { /* checking call; check elemental set being assigned */
2418 check_elem_set(mpl, info->set, info->memb->tuple,
2419 info->memb->value.set);
2420 }
2421 else
2422 { /* normal call; evaluate member, which has given n-tuple */
2423 info->refer = take_member_set(mpl, info->set, info->tuple);
2424 }
2425 return;
2426 }
2428 #if 1 /* 12/XII-2008 */
2429 static void saturate_set(MPL *mpl, SET *set)
2430 { GADGET *gadget = set->gadget;
2431 ELEMSET *data;
2432 MEMBER *elem, *memb;
2433 TUPLE *tuple, *work[20];
2434 int i;
2435 xprintf("Generating %s...\n", set->name);
2436 eval_whole_set(mpl, gadget->set);
2437 /* gadget set must have exactly one member */
2438 xassert(gadget->set->array != NULL);
2439 xassert(gadget->set->array->head != NULL);
2440 xassert(gadget->set->array->head == gadget->set->array->tail);
2441 data = gadget->set->array->head->value.set;
2442 xassert(data->type == A_NONE);
2443 xassert(data->dim == gadget->set->dimen);
2444 /* walk thru all elements of the plain set */
2445 for (elem = data->head; elem != NULL; elem = elem->next)
2446 { /* create a copy of n-tuple */
2447 tuple = copy_tuple(mpl, elem->tuple);
2448 /* rearrange component of the n-tuple */
2449 for (i = 0; i < gadget->set->dimen; i++)
2450 work[i] = NULL;
2451 for (i = 0; tuple != NULL; tuple = tuple->next)
2452 work[gadget->ind[i++]-1] = tuple;
2453 xassert(i == gadget->set->dimen);
2454 for (i = 0; i < gadget->set->dimen; i++)
2455 { xassert(work[i] != NULL);
2456 work[i]->next = work[i+1];
2457 }
2458 /* construct subscript list from first set->dim components */
2459 if (set->dim == 0)
2460 tuple = NULL;
2461 else
2462 tuple = work[0], work[set->dim-1]->next = NULL;
2463 /* find corresponding member of the set to be initialized */
2464 memb = find_member(mpl, set->array, tuple);
2465 if (memb == NULL)
2466 { /* not found; add new member to the set and assign it empty
2467 elemental set */
2468 memb = add_member(mpl, set->array, tuple);
2469 memb->value.set = create_elemset(mpl, set->dimen);
2470 }
2471 else
2472 { /* found; free subscript list */
2473 delete_tuple(mpl, tuple);
2474 }
2475 /* construct new n-tuple from rest set->dimen components */
2476 tuple = work[set->dim];
2477 xassert(set->dim + set->dimen == gadget->set->dimen);
2478 work[gadget->set->dimen-1]->next = NULL;
2479 /* and add it to the elemental set assigned to the member
2480 (no check for duplicates is needed) */
2481 add_tuple(mpl, memb->value.set, tuple);
2482 }
2483 /* the set has been saturated with data */
2484 set->data = 1;
2485 return;
2486 }
2487 #endif
2489 ELEMSET *eval_member_set /* returns reference, not value */
2490 ( MPL *mpl,
2491 SET *set, /* not changed */
2492 TUPLE *tuple /* not changed */
2493 )
2494 { /* this routine evaluates set member */
2495 struct eval_set_info _info, *info = &_info;
2496 xassert(set->dim == tuple_dimen(mpl, tuple));
2497 info->set = set;
2498 info->tuple = tuple;
2499 #if 1 /* 12/XII-2008 */
2500 if (set->gadget != NULL && set->data == 0)
2501 { /* initialize the set with data from a plain set */
2502 saturate_set(mpl, set);
2503 }
2504 #endif
2505 if (set->data == 1)
2506 { /* check data, which are provided in the data section, but not
2507 checked yet */
2508 /* save pointer to the last array member; note that during the
2509 check new members may be added beyond the last member due to
2510 references to the same parameter from default expression as
2511 well as from expressions that define restricting supersets;
2512 however, values assigned to the new members will be checked
2513 by other routine, so we don't need to check them here */
2514 MEMBER *tail = set->array->tail;
2515 /* change the data status to prevent infinite recursive loop
2516 due to references to the same set during the check */
2517 set->data = 2;
2518 /* check elemental sets assigned to array members in the data
2519 section until the marked member has been reached */
2520 for (info->memb = set->array->head; info->memb != NULL;
2521 info->memb = info->memb->next)
2522 { if (eval_within_domain(mpl, set->domain, info->memb->tuple,
2523 info, eval_set_func))
2524 out_of_domain(mpl, set->name, info->memb->tuple);
2525 if (info->memb == tail) break;
2526 }
2527 /* the check has been finished */
2528 }
2529 /* evaluate member, which has given n-tuple */
2530 info->memb = NULL;
2531 if (eval_within_domain(mpl, info->set->domain, info->tuple, info,
2532 eval_set_func))
2533 out_of_domain(mpl, set->name, info->tuple);
2534 /* bring evaluated reference to the calling program */
2535 return info->refer;
2536 }
2538 /*----------------------------------------------------------------------
2539 -- eval_whole_set - evaluate model set over entire domain.
2540 --
2541 -- This routine evaluates all members of specified model set over entire
2542 -- domain. */
2544 static int whole_set_func(MPL *mpl, void *info)
2545 { /* this is auxiliary routine to work within domain scope */
2546 SET *set = (SET *)info;
2547 TUPLE *tuple = get_domain_tuple(mpl, set->domain);
2548 eval_member_set(mpl, set, tuple);
2549 delete_tuple(mpl, tuple);
2550 return 0;
2551 }
2553 void eval_whole_set(MPL *mpl, SET *set)
2554 { loop_within_domain(mpl, set->domain, set, whole_set_func);
2555 return;
2556 }
2558 /*----------------------------------------------------------------------
2559 -- clean set - clean model set.
2560 --
2561 -- This routine cleans specified model set that assumes deleting all
2562 -- stuff dynamically allocated during the generation phase. */
2564 void clean_set(MPL *mpl, SET *set)
2565 { WITHIN *within;
2566 MEMBER *memb;
2567 /* clean subscript domain */
2568 clean_domain(mpl, set->domain);
2569 /* clean pseudo-code for computing supersets */
2570 for (within = set->within; within != NULL; within = within->next)
2571 clean_code(mpl, within->code);
2572 /* clean pseudo-code for computing assigned value */
2573 clean_code(mpl, set->assign);
2574 /* clean pseudo-code for computing default value */
2575 clean_code(mpl, set->option);
2576 /* reset data status flag */
2577 set->data = 0;
2578 /* delete content array */
2579 for (memb = set->array->head; memb != NULL; memb = memb->next)
2580 delete_value(mpl, set->array->type, &memb->value);
2581 delete_array(mpl, set->array), set->array = NULL;
2582 return;
2583 }
2585 /**********************************************************************/
2586 /* * * MODEL PARAMETERS * * */
2587 /**********************************************************************/
2589 /*----------------------------------------------------------------------
2590 -- check_value_num - check numeric value assigned to parameter member.
2591 --
2592 -- This routine checks if numeric value being assigned to some member
2593 -- of specified numeric model parameter satisfies to all restrictions.
2594 --
2595 -- NOTE: This routine must not be called out of domain scope. */
2597 void check_value_num
2598 ( MPL *mpl,
2599 PARAMETER *par, /* not changed */
2600 TUPLE *tuple, /* not changed */
2601 double value
2602 )
2603 { CONDITION *cond;
2604 WITHIN *in;
2605 int eqno;
2606 /* the value must satisfy to the parameter type */
2607 switch (par->type)
2608 { case A_NUMERIC:
2609 break;
2610 case A_INTEGER:
2611 if (value != floor(value))
2612 error(mpl, "%s%s = %.*g not integer", par->name,
2613 format_tuple(mpl, '[', tuple), DBL_DIG, value);
2614 break;
2615 case A_BINARY:
2616 if (!(value == 0.0 || value == 1.0))
2617 error(mpl, "%s%s = %.*g not binary", par->name,
2618 format_tuple(mpl, '[', tuple), DBL_DIG, value);
2619 break;
2620 default:
2621 xassert(par != par);
2622 }
2623 /* the value must satisfy to all specified conditions */
2624 for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,
2625 eqno++)
2626 { double bound;
2627 char *rho;
2628 xassert(cond->code != NULL);
2629 bound = eval_numeric(mpl, cond->code);
2630 switch (cond->rho)
2631 { case O_LT:
2632 if (!(value < bound))
2633 { rho = "<";
2634 err: error(mpl, "%s%s = %.*g not %s %.*g; see (%d)",
2635 par->name, format_tuple(mpl, '[', tuple), DBL_DIG,
2636 value, rho, DBL_DIG, bound, eqno);
2637 }
2638 break;
2639 case O_LE:
2640 if (!(value <= bound)) { rho = "<="; goto err; }
2641 break;
2642 case O_EQ:
2643 if (!(value == bound)) { rho = "="; goto err; }
2644 break;
2645 case O_GE:
2646 if (!(value >= bound)) { rho = ">="; goto err; }
2647 break;
2648 case O_GT:
2649 if (!(value > bound)) { rho = ">"; goto err; }
2650 break;
2651 case O_NE:
2652 if (!(value != bound)) { rho = "<>"; goto err; }
2653 break;
2654 default:
2655 xassert(cond != cond);
2656 }
2657 }
2658 /* the value must be in all specified supersets */
2659 for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)
2660 { TUPLE *dummy;
2661 xassert(in->code != NULL);
2662 xassert(in->code->dim == 1);
2663 dummy = expand_tuple(mpl, create_tuple(mpl),
2664 create_symbol_num(mpl, value));
2665 if (!is_member(mpl, in->code, dummy))
2666 error(mpl, "%s%s = %.*g not in specified set; see (%d)",
2667 par->name, format_tuple(mpl, '[', tuple), DBL_DIG,
2668 value, eqno);
2669 delete_tuple(mpl, dummy);
2670 }
2671 return;
2672 }
2674 /*----------------------------------------------------------------------
2675 -- take_member_num - obtain num. value assigned to parameter member.
2676 --
2677 -- This routine obtains a numeric value assigned to member of specified
2678 -- numeric model parameter and returns it on exit.
2679 --
2680 -- NOTE: This routine must not be called out of domain scope. */
2682 double take_member_num
2683 ( MPL *mpl,
2684 PARAMETER *par, /* not changed */
2685 TUPLE *tuple /* not changed */
2686 )
2687 { MEMBER *memb;
2688 double value;
2689 /* find member in the parameter array */
2690 memb = find_member(mpl, par->array, tuple);
2691 if (memb != NULL)
2692 { /* member exists, so just take its value */
2693 value = memb->value.num;
2694 }
2695 else if (par->assign != NULL)
2696 { /* compute value using assignment expression */
2697 value = eval_numeric(mpl, par->assign);
2698 add: /* check that the value satisfies to all restrictions, assign
2699 it to new member, and add the member to the array */
2700 check_value_num(mpl, par, tuple, value);
2701 memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));
2702 memb->value.num = value;
2703 }
2704 else if (par->option != NULL)
2705 { /* compute default value */
2706 value = eval_numeric(mpl, par->option);
2707 goto add;
2708 }
2709 else if (par->defval != NULL)
2710 { /* take default value provided in the data section */
2711 if (par->defval->str != NULL)
2712 error(mpl, "cannot convert %s to floating-point number",
2713 format_symbol(mpl, par->defval));
2714 value = par->defval->num;
2715 goto add;
2716 }
2717 else
2718 { /* no value is provided */
2719 error(mpl, "no value for %s%s", par->name, format_tuple(mpl,
2720 '[', tuple));
2721 }
2722 return value;
2723 }
2725 /*----------------------------------------------------------------------
2726 -- eval_member_num - evaluate num. value assigned to parameter member.
2727 --
2728 -- This routine evaluates a numeric value assigned to given member of
2729 -- specified numeric model parameter and returns it on exit. */
2731 struct eval_num_info
2732 { /* working info used by the routine eval_member_num */
2733 PARAMETER *par;
2734 /* model parameter */
2735 TUPLE *tuple;
2736 /* n-tuple, which defines parameter member */
2737 MEMBER *memb;
2738 /* normally this pointer is NULL; the routine uses this pointer
2739 to check data provided in the data section, in which case it
2740 points to a member currently checked; this check is performed
2741 automatically only once when a reference to any member occurs
2742 for the first time */
2743 double value;
2744 /* evaluated numeric value */
2745 };
2747 static void eval_num_func(MPL *mpl, void *_info)
2748 { /* this is auxiliary routine to work within domain scope */
2749 struct eval_num_info *info = _info;
2750 if (info->memb != NULL)
2751 { /* checking call; check numeric value being assigned */
2752 check_value_num(mpl, info->par, info->memb->tuple,
2753 info->memb->value.num);
2754 }
2755 else
2756 { /* normal call; evaluate member, which has given n-tuple */
2757 info->value = take_member_num(mpl, info->par, info->tuple);
2758 }
2759 return;
2760 }
2762 double eval_member_num
2763 ( MPL *mpl,
2764 PARAMETER *par, /* not changed */
2765 TUPLE *tuple /* not changed */
2766 )
2767 { /* this routine evaluates numeric parameter member */
2768 struct eval_num_info _info, *info = &_info;
2769 xassert(par->type == A_NUMERIC || par->type == A_INTEGER ||
2770 par->type == A_BINARY);
2771 xassert(par->dim == tuple_dimen(mpl, tuple));
2772 info->par = par;
2773 info->tuple = tuple;
2774 if (par->data == 1)
2775 { /* check data, which are provided in the data section, but not
2776 checked yet */
2777 /* save pointer to the last array member; note that during the
2778 check new members may be added beyond the last member due to
2779 references to the same parameter from default expression as
2780 well as from expressions that define restricting conditions;
2781 however, values assigned to the new members will be checked
2782 by other routine, so we don't need to check them here */
2783 MEMBER *tail = par->array->tail;
2784 /* change the data status to prevent infinite recursive loop
2785 due to references to the same parameter during the check */
2786 par->data = 2;
2787 /* check values assigned to array members in the data section
2788 until the marked member has been reached */
2789 for (info->memb = par->array->head; info->memb != NULL;
2790 info->memb = info->memb->next)
2791 { if (eval_within_domain(mpl, par->domain, info->memb->tuple,
2792 info, eval_num_func))
2793 out_of_domain(mpl, par->name, info->memb->tuple);
2794 if (info->memb == tail) break;
2795 }
2796 /* the check has been finished */
2797 }
2798 /* evaluate member, which has given n-tuple */
2799 info->memb = NULL;
2800 if (eval_within_domain(mpl, info->par->domain, info->tuple, info,
2801 eval_num_func))
2802 out_of_domain(mpl, par->name, info->tuple);
2803 /* bring evaluated value to the calling program */
2804 return info->value;
2805 }
2807 /*----------------------------------------------------------------------
2808 -- check_value_sym - check symbolic value assigned to parameter member.
2809 --
2810 -- This routine checks if symbolic value being assigned to some member
2811 -- of specified symbolic model parameter satisfies to all restrictions.
2812 --
2813 -- NOTE: This routine must not be called out of domain scope. */
2815 void check_value_sym
2816 ( MPL *mpl,
2817 PARAMETER *par, /* not changed */
2818 TUPLE *tuple, /* not changed */
2819 SYMBOL *value /* not changed */
2820 )
2821 { CONDITION *cond;
2822 WITHIN *in;
2823 int eqno;
2824 /* the value must satisfy to all specified conditions */
2825 for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,
2826 eqno++)
2827 { SYMBOL *bound;
2828 char buf[255+1];
2829 xassert(cond->code != NULL);
2830 bound = eval_symbolic(mpl, cond->code);
2831 switch (cond->rho)
2832 {
2833 #if 1 /* 13/VIII-2008 */
2834 case O_LT:
2835 if (!(compare_symbols(mpl, value, bound) < 0))
2836 { strcpy(buf, format_symbol(mpl, bound));
2837 xassert(strlen(buf) < sizeof(buf));
2838 error(mpl, "%s%s = %s not < %s",
2839 par->name, format_tuple(mpl, '[', tuple),
2840 format_symbol(mpl, value), buf, eqno);
2841 }
2842 break;
2843 case O_LE:
2844 if (!(compare_symbols(mpl, value, bound) <= 0))
2845 { strcpy(buf, format_symbol(mpl, bound));
2846 xassert(strlen(buf) < sizeof(buf));
2847 error(mpl, "%s%s = %s not <= %s",
2848 par->name, format_tuple(mpl, '[', tuple),
2849 format_symbol(mpl, value), buf, eqno);
2850 }
2851 break;
2852 #endif
2853 case O_EQ:
2854 if (!(compare_symbols(mpl, value, bound) == 0))
2855 { strcpy(buf, format_symbol(mpl, bound));
2856 xassert(strlen(buf) < sizeof(buf));
2857 error(mpl, "%s%s = %s not = %s",
2858 par->name, format_tuple(mpl, '[', tuple),
2859 format_symbol(mpl, value), buf, eqno);
2860 }
2861 break;
2862 #if 1 /* 13/VIII-2008 */
2863 case O_GE:
2864 if (!(compare_symbols(mpl, value, bound) >= 0))
2865 { strcpy(buf, format_symbol(mpl, bound));
2866 xassert(strlen(buf) < sizeof(buf));
2867 error(mpl, "%s%s = %s not >= %s",
2868 par->name, format_tuple(mpl, '[', tuple),
2869 format_symbol(mpl, value), buf, eqno);
2870 }
2871 break;
2872 case O_GT:
2873 if (!(compare_symbols(mpl, value, bound) > 0))
2874 { strcpy(buf, format_symbol(mpl, bound));
2875 xassert(strlen(buf) < sizeof(buf));
2876 error(mpl, "%s%s = %s not > %s",
2877 par->name, format_tuple(mpl, '[', tuple),
2878 format_symbol(mpl, value), buf, eqno);
2879 }
2880 break;
2881 #endif
2882 case O_NE:
2883 if (!(compare_symbols(mpl, value, bound) != 0))
2884 { strcpy(buf, format_symbol(mpl, bound));
2885 xassert(strlen(buf) < sizeof(buf));
2886 error(mpl, "%s%s = %s not <> %s",
2887 par->name, format_tuple(mpl, '[', tuple),
2888 format_symbol(mpl, value), buf, eqno);
2889 }
2890 break;
2891 default:
2892 xassert(cond != cond);
2893 }
2894 delete_symbol(mpl, bound);
2895 }
2896 /* the value must be in all specified supersets */
2897 for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)
2898 { TUPLE *dummy;
2899 xassert(in->code != NULL);
2900 xassert(in->code->dim == 1);
2901 dummy = expand_tuple(mpl, create_tuple(mpl), copy_symbol(mpl,
2902 value));
2903 if (!is_member(mpl, in->code, dummy))
2904 error(mpl, "%s%s = %s not in specified set; see (%d)",
2905 par->name, format_tuple(mpl, '[', tuple),
2906 format_symbol(mpl, value), eqno);
2907 delete_tuple(mpl, dummy);
2908 }
2909 return;
2910 }
2912 /*----------------------------------------------------------------------
2913 -- take_member_sym - obtain symb. value assigned to parameter member.
2914 --
2915 -- This routine obtains a symbolic value assigned to member of specified
2916 -- symbolic model parameter and returns it on exit.
2917 --
2918 -- NOTE: This routine must not be called out of domain scope. */
2920 SYMBOL *take_member_sym /* returns value, not reference */
2921 ( MPL *mpl,
2922 PARAMETER *par, /* not changed */
2923 TUPLE *tuple /* not changed */
2924 )
2925 { MEMBER *memb;
2926 SYMBOL *value;
2927 /* find member in the parameter array */
2928 memb = find_member(mpl, par->array, tuple);
2929 if (memb != NULL)
2930 { /* member exists, so just take its value */
2931 value = copy_symbol(mpl, memb->value.sym);
2932 }
2933 else if (par->assign != NULL)
2934 { /* compute value using assignment expression */
2935 value = eval_symbolic(mpl, par->assign);
2936 add: /* check that the value satisfies to all restrictions, assign
2937 it to new member, and add the member to the array */
2938 check_value_sym(mpl, par, tuple, value);
2939 memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));
2940 memb->value.sym = copy_symbol(mpl, value);
2941 }
2942 else if (par->option != NULL)
2943 { /* compute default value */
2944 value = eval_symbolic(mpl, par->option);
2945 goto add;
2946 }
2947 else if (par->defval != NULL)
2948 { /* take default value provided in the data section */
2949 value = copy_symbol(mpl, par->defval);
2950 goto add;
2951 }
2952 else
2953 { /* no value is provided */
2954 error(mpl, "no value for %s%s", par->name, format_tuple(mpl,
2955 '[', tuple));
2956 }
2957 return value;
2958 }
2960 /*----------------------------------------------------------------------
2961 -- eval_member_sym - evaluate symb. value assigned to parameter member.
2962 --
2963 -- This routine evaluates a symbolic value assigned to given member of
2964 -- specified symbolic model parameter and returns it on exit. */
2966 struct eval_sym_info
2967 { /* working info used by the routine eval_member_sym */
2968 PARAMETER *par;
2969 /* model parameter */
2970 TUPLE *tuple;
2971 /* n-tuple, which defines parameter member */
2972 MEMBER *memb;
2973 /* normally this pointer is NULL; the routine uses this pointer
2974 to check data provided in the data section, in which case it
2975 points to a member currently checked; this check is performed
2976 automatically only once when a reference to any member occurs
2977 for the first time */
2978 SYMBOL *value;
2979 /* evaluated symbolic value */
2980 };
2982 static void eval_sym_func(MPL *mpl, void *_info)
2983 { /* this is auxiliary routine to work within domain scope */
2984 struct eval_sym_info *info = _info;
2985 if (info->memb != NULL)
2986 { /* checking call; check symbolic value being assigned */
2987 check_value_sym(mpl, info->par, info->memb->tuple,
2988 info->memb->value.sym);
2989 }
2990 else
2991 { /* normal call; evaluate member, which has given n-tuple */
2992 info->value = take_member_sym(mpl, info->par, info->tuple);
2993 }
2994 return;
2995 }
2997 SYMBOL *eval_member_sym /* returns value, not reference */
2998 ( MPL *mpl,
2999 PARAMETER *par, /* not changed */
3000 TUPLE *tuple /* not changed */
3001 )
3002 { /* this routine evaluates symbolic parameter member */
3003 struct eval_sym_info _info, *info = &_info;
3004 xassert(par->type == A_SYMBOLIC);
3005 xassert(par->dim == tuple_dimen(mpl, tuple));
3006 info->par = par;
3007 info->tuple = tuple;
3008 if (par->data == 1)
3009 { /* check data, which are provided in the data section, but not
3010 checked yet */
3011 /* save pointer to the last array member; note that during the
3012 check new members may be added beyond the last member due to
3013 references to the same parameter from default expression as
3014 well as from expressions that define restricting conditions;
3015 however, values assigned to the new members will be checked
3016 by other routine, so we don't need to check them here */
3017 MEMBER *tail = par->array->tail;
3018 /* change the data status to prevent infinite recursive loop
3019 due to references to the same parameter during the check */
3020 par->data = 2;
3021 /* check values assigned to array members in the data section
3022 until the marked member has been reached */
3023 for (info->memb = par->array->head; info->memb != NULL;
3024 info->memb = info->memb->next)
3025 { if (eval_within_domain(mpl, par->domain, info->memb->tuple,
3026 info, eval_sym_func))
3027 out_of_domain(mpl, par->name, info->memb->tuple);
3028 if (info->memb == tail) break;
3029 }
3030 /* the check has been finished */
3031 }
3032 /* evaluate member, which has given n-tuple */
3033 info->memb = NULL;
3034 if (eval_within_domain(mpl, info->par->domain, info->tuple, info,
3035 eval_sym_func))
3036 out_of_domain(mpl, par->name, info->tuple);
3037 /* bring evaluated value to the calling program */
3038 return info->value;
3039 }
3041 /*----------------------------------------------------------------------
3042 -- eval_whole_par - evaluate model parameter over entire domain.
3043 --
3044 -- This routine evaluates all members of specified model parameter over
3045 -- entire domain. */
3047 static int whole_par_func(MPL *mpl, void *info)
3048 { /* this is auxiliary routine to work within domain scope */
3049 PARAMETER *par = (PARAMETER *)info;
3050 TUPLE *tuple = get_domain_tuple(mpl, par->domain);
3051 switch (par->type)
3052 { case A_NUMERIC:
3053 case A_INTEGER:
3054 case A_BINARY:
3055 eval_member_num(mpl, par, tuple);
3056 break;
3057 case A_SYMBOLIC:
3058 delete_symbol(mpl, eval_member_sym(mpl, par, tuple));
3059 break;
3060 default:
3061 xassert(par != par);
3062 }
3063 delete_tuple(mpl, tuple);
3064 return 0;
3065 }
3067 void eval_whole_par(MPL *mpl, PARAMETER *par)
3068 { loop_within_domain(mpl, par->domain, par, whole_par_func);
3069 return;
3070 }
3072 /*----------------------------------------------------------------------
3073 -- clean_parameter - clean model parameter.
3074 --
3075 -- This routine cleans specified model parameter that assumes deleting
3076 -- all stuff dynamically allocated during the generation phase. */
3078 void clean_parameter(MPL *mpl, PARAMETER *par)
3079 { CONDITION *cond;
3080 WITHIN *in;
3081 MEMBER *memb;
3082 /* clean subscript domain */
3083 clean_domain(mpl, par->domain);
3084 /* clean pseudo-code for computing restricting conditions */
3085 for (cond = par->cond; cond != NULL; cond = cond->next)
3086 clean_code(mpl, cond->code);
3087 /* clean pseudo-code for computing restricting supersets */
3088 for (in = par->in; in != NULL; in = in->next)
3089 clean_code(mpl, in->code);
3090 /* clean pseudo-code for computing assigned value */
3091 clean_code(mpl, par->assign);
3092 /* clean pseudo-code for computing default value */
3093 clean_code(mpl, par->option);
3094 /* reset data status flag */
3095 par->data = 0;
3096 /* delete default symbolic value */
3097 if (par->defval != NULL)
3098 delete_symbol(mpl, par->defval), par->defval = NULL;
3099 /* delete content array */
3100 for (memb = par->array->head; memb != NULL; memb = memb->next)
3101 delete_value(mpl, par->array->type, &memb->value);
3102 delete_array(mpl, par->array), par->array = NULL;
3103 return;
3104 }
3106 /**********************************************************************/
3107 /* * * MODEL VARIABLES * * */
3108 /**********************************************************************/
3110 /*----------------------------------------------------------------------
3111 -- take_member_var - obtain reference to elemental variable.
3112 --
3113 -- This routine obtains a reference to elemental variable assigned to
3114 -- given member of specified model variable and returns it on exit. If
3115 -- necessary, new elemental variable is created.
3116 --
3117 -- NOTE: This routine must not be called out of domain scope. */
3119 ELEMVAR *take_member_var /* returns reference */
3120 ( MPL *mpl,
3121 VARIABLE *var, /* not changed */
3122 TUPLE *tuple /* not changed */
3123 )
3124 { MEMBER *memb;
3125 ELEMVAR *refer;
3126 /* find member in the variable array */
3127 memb = find_member(mpl, var->array, tuple);
3128 if (memb != NULL)
3129 { /* member exists, so just take the reference */
3130 refer = memb->value.var;
3131 }
3132 else
3133 { /* member is referenced for the first time and therefore does
3134 not exist; create new elemental variable, assign it to new
3135 member, and add the member to the variable array */
3136 memb = add_member(mpl, var->array, copy_tuple(mpl, tuple));
3137 refer = (memb->value.var =
3138 dmp_get_atom(mpl->elemvars, sizeof(ELEMVAR)));
3139 refer->j = 0;
3140 refer->var = var;
3141 refer->memb = memb;
3142 /* compute lower bound */
3143 if (var->lbnd == NULL)
3144 refer->lbnd = 0.0;
3145 else
3146 refer->lbnd = eval_numeric(mpl, var->lbnd);
3147 /* compute upper bound */
3148 if (var->ubnd == NULL)
3149 refer->ubnd = 0.0;
3150 else if (var->ubnd == var->lbnd)
3151 refer->ubnd = refer->lbnd;
3152 else
3153 refer->ubnd = eval_numeric(mpl, var->ubnd);
3154 /* nullify working quantity */
3155 refer->temp = 0.0;
3156 #if 1 /* 15/V-2010 */
3157 /* solution has not been obtained by the solver yet */
3158 refer->stat = 0;
3159 refer->prim = refer->dual = 0.0;
3160 #endif
3161 }
3162 return refer;
3163 }
3165 /*----------------------------------------------------------------------
3166 -- eval_member_var - evaluate reference to elemental variable.
3167 --
3168 -- This routine evaluates a reference to elemental variable assigned to
3169 -- member of specified model variable and returns it on exit. */
3171 struct eval_var_info
3172 { /* working info used by the routine eval_member_var */
3173 VARIABLE *var;
3174 /* model variable */
3175 TUPLE *tuple;
3176 /* n-tuple, which defines variable member */
3177 ELEMVAR *refer;
3178 /* evaluated reference to elemental variable */
3179 };
3181 static void eval_var_func(MPL *mpl, void *_info)
3182 { /* this is auxiliary routine to work within domain scope */
3183 struct eval_var_info *info = _info;
3184 info->refer = take_member_var(mpl, info->var, info->tuple);
3185 return;
3186 }
3188 ELEMVAR *eval_member_var /* returns reference */
3189 ( MPL *mpl,
3190 VARIABLE *var, /* not changed */
3191 TUPLE *tuple /* not changed */
3192 )
3193 { /* this routine evaluates variable member */
3194 struct eval_var_info _info, *info = &_info;
3195 xassert(var->dim == tuple_dimen(mpl, tuple));
3196 info->var = var;
3197 info->tuple = tuple;
3198 /* evaluate member, which has given n-tuple */
3199 if (eval_within_domain(mpl, info->var->domain, info->tuple, info,
3200 eval_var_func))
3201 out_of_domain(mpl, var->name, info->tuple);
3202 /* bring evaluated reference to the calling program */
3203 return info->refer;
3204 }
3206 /*----------------------------------------------------------------------
3207 -- eval_whole_var - evaluate model variable over entire domain.
3208 --
3209 -- This routine evaluates all members of specified model variable over
3210 -- entire domain. */
3212 static int whole_var_func(MPL *mpl, void *info)
3213 { /* this is auxiliary routine to work within domain scope */
3214 VARIABLE *var = (VARIABLE *)info;
3215 TUPLE *tuple = get_domain_tuple(mpl, var->domain);
3216 eval_member_var(mpl, var, tuple);
3217 delete_tuple(mpl, tuple);
3218 return 0;
3219 }
3221 void eval_whole_var(MPL *mpl, VARIABLE *var)
3222 { loop_within_domain(mpl, var->domain, var, whole_var_func);
3223 return;
3224 }
3226 /*----------------------------------------------------------------------
3227 -- clean_variable - clean model variable.
3228 --
3229 -- This routine cleans specified model variable that assumes deleting
3230 -- all stuff dynamically allocated during the generation phase. */
3232 void clean_variable(MPL *mpl, VARIABLE *var)
3233 { MEMBER *memb;
3234 /* clean subscript domain */
3235 clean_domain(mpl, var->domain);
3236 /* clean code for computing lower bound */
3237 clean_code(mpl, var->lbnd);
3238 /* clean code for computing upper bound */
3239 if (var->ubnd != var->lbnd) clean_code(mpl, var->ubnd);
3240 /* delete content array */
3241 for (memb = var->array->head; memb != NULL; memb = memb->next)
3242 dmp_free_atom(mpl->elemvars, memb->value.var, sizeof(ELEMVAR));
3243 delete_array(mpl, var->array), var->array = NULL;
3244 return;
3245 }
3247 /**********************************************************************/
3248 /* * * MODEL CONSTRAINTS AND OBJECTIVES * * */
3249 /**********************************************************************/
3251 /*----------------------------------------------------------------------
3252 -- take_member_con - obtain reference to elemental constraint.
3253 --
3254 -- This routine obtains a reference to elemental constraint assigned
3255 -- to given member of specified model constraint and returns it on exit.
3256 -- If necessary, new elemental constraint is created.
3257 --
3258 -- NOTE: This routine must not be called out of domain scope. */
3260 ELEMCON *take_member_con /* returns reference */
3261 ( MPL *mpl,
3262 CONSTRAINT *con, /* not changed */
3263 TUPLE *tuple /* not changed */
3264 )
3265 { MEMBER *memb;
3266 ELEMCON *refer;
3267 /* find member in the constraint array */
3268 memb = find_member(mpl, con->array, tuple);
3269 if (memb != NULL)
3270 { /* member exists, so just take the reference */
3271 refer = memb->value.con;
3272 }
3273 else
3274 { /* member is referenced for the first time and therefore does
3275 not exist; create new elemental constraint, assign it to new
3276 member, and add the member to the constraint array */
3277 memb = add_member(mpl, con->array, copy_tuple(mpl, tuple));
3278 refer = (memb->value.con =
3279 dmp_get_atom(mpl->elemcons, sizeof(ELEMCON)));
3280 refer->i = 0;
3281 refer->con = con;
3282 refer->memb = memb;
3283 /* compute linear form */
3284 xassert(con->code != NULL);
3285 refer->form = eval_formula(mpl, con->code);
3286 /* compute lower and upper bounds */
3287 if (con->lbnd == NULL && con->ubnd == NULL)
3288 { /* objective has no bounds */
3289 double temp;
3290 xassert(con->type == A_MINIMIZE || con->type == A_MAXIMIZE);
3291 /* carry the constant term to the right-hand side */
3292 refer->form = remove_constant(mpl, refer->form, &temp);
3293 refer->lbnd = refer->ubnd = - temp;
3294 }
3295 else if (con->lbnd != NULL && con->ubnd == NULL)
3296 { /* constraint a * x + b >= c * y + d is transformed to the
3297 standard form a * x - c * y >= d - b */
3298 double temp;
3299 xassert(con->type == A_CONSTRAINT);
3300 refer->form = linear_comb(mpl,
3301 +1.0, refer->form,
3302 -1.0, eval_formula(mpl, con->lbnd));
3303 refer->form = remove_constant(mpl, refer->form, &temp);
3304 refer->lbnd = - temp;
3305 refer->ubnd = 0.0;
3306 }
3307 else if (con->lbnd == NULL && con->ubnd != NULL)
3308 { /* constraint a * x + b <= c * y + d is transformed to the
3309 standard form a * x - c * y <= d - b */
3310 double temp;
3311 xassert(con->type == A_CONSTRAINT);
3312 refer->form = linear_comb(mpl,
3313 +1.0, refer->form,
3314 -1.0, eval_formula(mpl, con->ubnd));
3315 refer->form = remove_constant(mpl, refer->form, &temp);
3316 refer->lbnd = 0.0;
3317 refer->ubnd = - temp;
3318 }
3319 else if (con->lbnd == con->ubnd)
3320 { /* constraint a * x + b = c * y + d is transformed to the
3321 standard form a * x - c * y = d - b */
3322 double temp;
3323 xassert(con->type == A_CONSTRAINT);
3324 refer->form = linear_comb(mpl,
3325 +1.0, refer->form,
3326 -1.0, eval_formula(mpl, con->lbnd));
3327 refer->form = remove_constant(mpl, refer->form, &temp);
3328 refer->lbnd = refer->ubnd = - temp;
3329 }
3330 else
3331 { /* ranged constraint c <= a * x + b <= d is transformed to
3332 the standard form c - b <= a * x <= d - b */
3333 double temp, temp1, temp2;
3334 xassert(con->type == A_CONSTRAINT);
3335 refer->form = remove_constant(mpl, refer->form, &temp);
3336 xassert(remove_constant(mpl, eval_formula(mpl, con->lbnd),
3337 &temp1) == NULL);
3338 xassert(remove_constant(mpl, eval_formula(mpl, con->ubnd),
3339 &temp2) == NULL);
3340 refer->lbnd = fp_sub(mpl, temp1, temp);
3341 refer->ubnd = fp_sub(mpl, temp2, temp);
3342 }
3343 #if 1 /* 15/V-2010 */
3344 /* solution has not been obtained by the solver yet */
3345 refer->stat = 0;
3346 refer->prim = refer->dual = 0.0;
3347 #endif
3348 }
3349 return refer;
3350 }
3352 /*----------------------------------------------------------------------
3353 -- eval_member_con - evaluate reference to elemental constraint.
3354 --
3355 -- This routine evaluates a reference to elemental constraint assigned
3356 -- to member of specified model constraint and returns it on exit. */
3358 struct eval_con_info
3359 { /* working info used by the routine eval_member_con */
3360 CONSTRAINT *con;
3361 /* model constraint */
3362 TUPLE *tuple;
3363 /* n-tuple, which defines constraint member */
3364 ELEMCON *refer;
3365 /* evaluated reference to elemental constraint */
3366 };
3368 static void eval_con_func(MPL *mpl, void *_info)
3369 { /* this is auxiliary routine to work within domain scope */
3370 struct eval_con_info *info = _info;
3371 info->refer = take_member_con(mpl, info->con, info->tuple);
3372 return;
3373 }
3375 ELEMCON *eval_member_con /* returns reference */
3376 ( MPL *mpl,
3377 CONSTRAINT *con, /* not changed */
3378 TUPLE *tuple /* not changed */
3379 )
3380 { /* this routine evaluates constraint member */
3381 struct eval_con_info _info, *info = &_info;
3382 xassert(con->dim == tuple_dimen(mpl, tuple));
3383 info->con = con;
3384 info->tuple = tuple;
3385 /* evaluate member, which has given n-tuple */
3386 if (eval_within_domain(mpl, info->con->domain, info->tuple, info,
3387 eval_con_func))
3388 out_of_domain(mpl, con->name, info->tuple);
3389 /* bring evaluated reference to the calling program */
3390 return info->refer;
3391 }
3393 /*----------------------------------------------------------------------
3394 -- eval_whole_con - evaluate model constraint over entire domain.
3395 --
3396 -- This routine evaluates all members of specified model constraint over
3397 -- entire domain. */
3399 static int whole_con_func(MPL *mpl, void *info)
3400 { /* this is auxiliary routine to work within domain scope */
3401 CONSTRAINT *con = (CONSTRAINT *)info;
3402 TUPLE *tuple = get_domain_tuple(mpl, con->domain);
3403 eval_member_con(mpl, con, tuple);
3404 delete_tuple(mpl, tuple);
3405 return 0;
3406 }
3408 void eval_whole_con(MPL *mpl, CONSTRAINT *con)
3409 { loop_within_domain(mpl, con->domain, con, whole_con_func);
3410 return;
3411 }
3413 /*----------------------------------------------------------------------
3414 -- clean_constraint - clean model constraint.
3415 --
3416 -- This routine cleans specified model constraint that assumes deleting
3417 -- all stuff dynamically allocated during the generation phase. */
3419 void clean_constraint(MPL *mpl, CONSTRAINT *con)
3420 { MEMBER *memb;
3421 /* clean subscript domain */
3422 clean_domain(mpl, con->domain);
3423 /* clean code for computing main linear form */
3424 clean_code(mpl, con->code);
3425 /* clean code for computing lower bound */
3426 clean_code(mpl, con->lbnd);
3427 /* clean code for computing upper bound */
3428 if (con->ubnd != con->lbnd) clean_code(mpl, con->ubnd);
3429 /* delete content array */
3430 for (memb = con->array->head; memb != NULL; memb = memb->next)
3431 { delete_formula(mpl, memb->value.con->form);
3432 dmp_free_atom(mpl->elemcons, memb->value.con, sizeof(ELEMCON));
3433 }
3434 delete_array(mpl, con->array), con->array = NULL;
3435 return;
3436 }
3438 /**********************************************************************/
3439 /* * * PSEUDO-CODE * * */
3440 /**********************************************************************/
3442 /*----------------------------------------------------------------------
3443 -- eval_numeric - evaluate pseudo-code to determine numeric value.
3444 --
3445 -- This routine evaluates specified pseudo-code to determine resultant
3446 -- numeric value, which is returned on exit. */
3448 struct iter_num_info
3449 { /* working info used by the routine iter_num_func */
3450 CODE *code;
3451 /* pseudo-code for iterated operation to be performed */
3452 double value;
3453 /* resultant value */
3454 };
3456 static int iter_num_func(MPL *mpl, void *_info)
3457 { /* this is auxiliary routine used to perform iterated operation
3458 on numeric "integrand" within domain scope */
3459 struct iter_num_info *info = _info;
3460 double temp;
3461 temp = eval_numeric(mpl, info->code->arg.loop.x);
3462 switch (info->code->op)
3463 { case O_SUM:
3464 /* summation over domain */
3465 info->value = fp_add(mpl, info->value, temp);
3466 break;
3467 case O_PROD:
3468 /* multiplication over domain */
3469 info->value = fp_mul(mpl, info->value, temp);
3470 break;
3471 case O_MINIMUM:
3472 /* minimum over domain */
3473 if (info->value > temp) info->value = temp;
3474 break;
3475 case O_MAXIMUM:
3476 /* maximum over domain */
3477 if (info->value < temp) info->value = temp;
3478 break;
3479 default:
3480 xassert(info != info);
3481 }
3482 return 0;
3483 }
3485 double eval_numeric(MPL *mpl, CODE *code)
3486 { double value;
3487 xassert(code != NULL);
3488 xassert(code->type == A_NUMERIC);
3489 xassert(code->dim == 0);
3490 /* if the operation has a side effect, invalidate and delete the
3491 resultant value */
3492 if (code->vflag && code->valid)
3493 { code->valid = 0;
3494 delete_value(mpl, code->type, &code->value);
3495 }
3496 /* if resultant value is valid, no evaluation is needed */
3497 if (code->valid)
3498 { value = code->value.num;
3499 goto done;
3500 }
3501 /* evaluate pseudo-code recursively */
3502 switch (code->op)
3503 { case O_NUMBER:
3504 /* take floating-point number */
3505 value = code->arg.num;
3506 break;
3507 case O_MEMNUM:
3508 /* take member of numeric parameter */
3509 { TUPLE *tuple;
3510 ARG_LIST *e;
3511 tuple = create_tuple(mpl);
3512 for (e = code->arg.par.list; e != NULL; e = e->next)
3513 tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
3514 e->x));
3515 value = eval_member_num(mpl, code->arg.par.par, tuple);
3516 delete_tuple(mpl, tuple);
3517 }
3518 break;
3519 case O_MEMVAR:
3520 /* take computed value of elemental variable */
3521 { TUPLE *tuple;
3522 ARG_LIST *e;
3523 #if 1 /* 15/V-2010 */
3524 ELEMVAR *var;
3525 #endif
3526 tuple = create_tuple(mpl);
3527 for (e = code->arg.var.list; e != NULL; e = e->next)
3528 tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
3529 e->x));
3530 #if 0 /* 15/V-2010 */
3531 value = eval_member_var(mpl, code->arg.var.var, tuple)
3532 ->value;
3533 #else
3534 var = eval_member_var(mpl, code->arg.var.var, tuple);
3535 switch (code->arg.var.suff)
3536 { case DOT_LB:
3537 if (var->var->lbnd == NULL)
3538 value = -DBL_MAX;
3539 else
3540 value = var->lbnd;
3541 break;
3542 case DOT_UB:
3543 if (var->var->ubnd == NULL)
3544 value = +DBL_MAX;
3545 else
3546 value = var->ubnd;
3547 break;
3548 case DOT_STATUS:
3549 value = var->stat;
3550 break;
3551 case DOT_VAL:
3552 value = var->prim;
3553 break;
3554 case DOT_DUAL:
3555 value = var->dual;
3556 break;
3557 default:
3558 xassert(code != code);
3559 }
3560 #endif
3561 delete_tuple(mpl, tuple);
3562 }
3563 break;
3564 #if 1 /* 15/V-2010 */
3565 case O_MEMCON:
3566 /* take computed value of elemental constraint */
3567 { TUPLE *tuple;
3568 ARG_LIST *e;
3569 ELEMCON *con;
3570 tuple = create_tuple(mpl);
3571 for (e = code->arg.con.list; e != NULL; e = e->next)
3572 tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
3573 e->x));
3574 con = eval_member_con(mpl, code->arg.con.con, tuple);
3575 switch (code->arg.con.suff)
3576 { case DOT_LB:
3577 if (con->con->lbnd == NULL)
3578 value = -DBL_MAX;
3579 else
3580 value = con->lbnd;
3581 break;
3582 case DOT_UB:
3583 if (con->con->ubnd == NULL)
3584 value = +DBL_MAX;
3585 else
3586 value = con->ubnd;
3587 break;
3588 case DOT_STATUS:
3589 value = con->stat;
3590 break;
3591 case DOT_VAL:
3592 value = con->prim;
3593 break;
3594 case DOT_DUAL:
3595 value = con->dual;
3596 break;
3597 default:
3598 xassert(code != code);
3599 }
3600 delete_tuple(mpl, tuple);
3601 }
3602 break;
3603 #endif
3604 case O_IRAND224:
3605 /* pseudo-random in [0, 2^24-1] */
3606 value = fp_irand224(mpl);
3607 break;
3608 case O_UNIFORM01:
3609 /* pseudo-random in [0, 1) */
3610 value = fp_uniform01(mpl);
3611 break;
3612 case O_NORMAL01:
3613 /* gaussian random, mu = 0, sigma = 1 */
3614 value = fp_normal01(mpl);
3615 break;
3616 case O_GMTIME:
3617 /* current calendar time */
3618 value = fn_gmtime(mpl);
3619 break;
3620 case O_CVTNUM:
3621 /* conversion to numeric */
3622 { SYMBOL *sym;
3623 sym = eval_symbolic(mpl, code->arg.arg.x);
3624 #if 0 /* 23/XI-2008 */
3625 if (sym->str != NULL)
3626 error(mpl, "cannot convert %s to floating-point numbe"
3627 "r", format_symbol(mpl, sym));
3628 value = sym->num;
3629 #else
3630 if (sym->str == NULL)
3631 value = sym->num;
3632 else
3633 { if (str2num(sym->str, &value))
3634 error(mpl, "cannot convert %s to floating-point nu"
3635 "mber", format_symbol(mpl, sym));
3636 }
3637 #endif
3638 delete_symbol(mpl, sym);
3639 }
3640 break;
3641 case O_PLUS:
3642 /* unary plus */
3643 value = + eval_numeric(mpl, code->arg.arg.x);
3644 break;
3645 case O_MINUS:
3646 /* unary minus */
3647 value = - eval_numeric(mpl, code->arg.arg.x);
3648 break;
3649 case O_ABS:
3650 /* absolute value */
3651 value = fabs(eval_numeric(mpl, code->arg.arg.x));
3652 break;
3653 case O_CEIL:
3654 /* round upward ("ceiling of x") */
3655 value = ceil(eval_numeric(mpl, code->arg.arg.x));
3656 break;
3657 case O_FLOOR:
3658 /* round downward ("floor of x") */
3659 value = floor(eval_numeric(mpl, code->arg.arg.x));
3660 break;
3661 case O_EXP:
3662 /* base-e exponential */
3663 value = fp_exp(mpl, eval_numeric(mpl, code->arg.arg.x));
3664 break;
3665 case O_LOG:
3666 /* natural logarithm */
3667 value = fp_log(mpl, eval_numeric(mpl, code->arg.arg.x));
3668 break;
3669 case O_LOG10:
3670 /* common (decimal) logarithm */
3671 value = fp_log10(mpl, eval_numeric(mpl, code->arg.arg.x));
3672 break;
3673 case O_SQRT:
3674 /* square root */
3675 value = fp_sqrt(mpl, eval_numeric(mpl, code->arg.arg.x));
3676 break;
3677 case O_SIN:
3678 /* trigonometric sine */
3679 value = fp_sin(mpl, eval_numeric(mpl, code->arg.arg.x));
3680 break;
3681 case O_COS:
3682 /* trigonometric cosine */
3683 value = fp_cos(mpl, eval_numeric(mpl, code->arg.arg.x));
3684 break;
3685 case O_ATAN:
3686 /* trigonometric arctangent (one argument) */
3687 value = fp_atan(mpl, eval_numeric(mpl, code->arg.arg.x));
3688 break;
3689 case O_ATAN2:
3690 /* trigonometric arctangent (two arguments) */
3691 value = fp_atan2(mpl,
3692 eval_numeric(mpl, code->arg.arg.x),
3693 eval_numeric(mpl, code->arg.arg.y));
3694 break;
3695 case O_ROUND:
3696 /* round to nearest integer */
3697 value = fp_round(mpl,
3698 eval_numeric(mpl, code->arg.arg.x), 0.0);
3699 break;
3700 case O_ROUND2:
3701 /* round to n fractional digits */
3702 value = fp_round(mpl,
3703 eval_numeric(mpl, code->arg.arg.x),
3704 eval_numeric(mpl, code->arg.arg.y));
3705 break;
3706 case O_TRUNC:
3707 /* truncate to nearest integer */
3708 value = fp_trunc(mpl,
3709 eval_numeric(mpl, code->arg.arg.x), 0.0);
3710 break;
3711 case O_TRUNC2:
3712 /* truncate to n fractional digits */
3713 value = fp_trunc(mpl,
3714 eval_numeric(mpl, code->arg.arg.x),
3715 eval_numeric(mpl, code->arg.arg.y));
3716 break;
3717 case O_ADD:
3718 /* addition */
3719 value = fp_add(mpl,
3720 eval_numeric(mpl, code->arg.arg.x),
3721 eval_numeric(mpl, code->arg.arg.y));
3722 break;
3723 case O_SUB:
3724 /* subtraction */
3725 value = fp_sub(mpl,
3726 eval_numeric(mpl, code->arg.arg.x),
3727 eval_numeric(mpl, code->arg.arg.y));
3728 break;
3729 case O_LESS:
3730 /* non-negative subtraction */
3731 value = fp_less(mpl,
3732 eval_numeric(mpl, code->arg.arg.x),
3733 eval_numeric(mpl, code->arg.arg.y));
3734 break;
3735 case O_MUL:
3736 /* multiplication */
3737 value = fp_mul(mpl,
3738 eval_numeric(mpl, code->arg.arg.x),
3739 eval_numeric(mpl, code->arg.arg.y));
3740 break;
3741 case O_DIV:
3742 /* division */
3743 value = fp_div(mpl,
3744 eval_numeric(mpl, code->arg.arg.x),
3745 eval_numeric(mpl, code->arg.arg.y));
3746 break;
3747 case O_IDIV:
3748 /* quotient of exact division */
3749 value = fp_idiv(mpl,
3750 eval_numeric(mpl, code->arg.arg.x),
3751 eval_numeric(mpl, code->arg.arg.y));
3752 break;
3753 case O_MOD:
3754 /* remainder of exact division */
3755 value = fp_mod(mpl,
3756 eval_numeric(mpl, code->arg.arg.x),
3757 eval_numeric(mpl, code->arg.arg.y));
3758 break;
3759 case O_POWER:
3760 /* exponentiation (raise to power) */
3761 value = fp_power(mpl,
3762 eval_numeric(mpl, code->arg.arg.x),
3763 eval_numeric(mpl, code->arg.arg.y));
3764 break;
3765 case O_UNIFORM:
3766 /* pseudo-random in [a, b) */
3767 value = fp_uniform(mpl,
3768 eval_numeric(mpl, code->arg.arg.x),
3769 eval_numeric(mpl, code->arg.arg.y));
3770 break;
3771 case O_NORMAL:
3772 /* gaussian random, given mu and sigma */
3773 value = fp_normal(mpl,
3774 eval_numeric(mpl, code->arg.arg.x),
3775 eval_numeric(mpl, code->arg.arg.y));
3776 break;
3777 case O_CARD:
3778 { ELEMSET *set;
3779 set = eval_elemset(mpl, code->arg.arg.x);
3780 value = set->size;
3781 delete_array(mpl, set);
3782 }
3783 break;
3784 case O_LENGTH:
3785 { SYMBOL *sym;
3786 char str[MAX_LENGTH+1];
3787 sym = eval_symbolic(mpl, code->arg.arg.x);
3788 if (sym->str == NULL)
3789 sprintf(str, "%.*g", DBL_DIG, sym->num);
3790 else
3791 fetch_string(mpl, sym->str, str);
3792 delete_symbol(mpl, sym);
3793 value = strlen(str);
3794 }
3795 break;
3796 case O_STR2TIME:
3797 { SYMBOL *sym;
3798 char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];
3799 sym = eval_symbolic(mpl, code->arg.arg.x);
3800 if (sym->str == NULL)
3801 sprintf(str, "%.*g", DBL_DIG, sym->num);
3802 else
3803 fetch_string(mpl, sym->str, str);
3804 delete_symbol(mpl, sym);
3805 sym = eval_symbolic(mpl, code->arg.arg.y);
3806 if (sym->str == NULL)
3807 sprintf(fmt, "%.*g", DBL_DIG, sym->num);
3808 else
3809 fetch_string(mpl, sym->str, fmt);
3810 delete_symbol(mpl, sym);
3811 value = fn_str2time(mpl, str, fmt);
3812 }
3813 break;
3814 case O_FORK:
3815 /* if-then-else */
3816 if (eval_logical(mpl, code->arg.arg.x))
3817 value = eval_numeric(mpl, code->arg.arg.y);
3818 else if (code->arg.arg.z == NULL)
3819 value = 0.0;
3820 else
3821 value = eval_numeric(mpl, code->arg.arg.z);
3822 break;
3823 case O_MIN:
3824 /* minimal value (n-ary) */
3825 { ARG_LIST *e;
3826 double temp;
3827 value = +DBL_MAX;
3828 for (e = code->arg.list; e != NULL; e = e->next)
3829 { temp = eval_numeric(mpl, e->x);
3830 if (value > temp) value = temp;
3831 }
3832 }
3833 break;
3834 case O_MAX:
3835 /* maximal value (n-ary) */
3836 { ARG_LIST *e;
3837 double temp;
3838 value = -DBL_MAX;
3839 for (e = code->arg.list; e != NULL; e = e->next)
3840 { temp = eval_numeric(mpl, e->x);
3841 if (value < temp) value = temp;
3842 }
3843 }
3844 break;
3845 case O_SUM:
3846 /* summation over domain */
3847 { struct iter_num_info _info, *info = &_info;
3848 info->code = code;
3849 info->value = 0.0;
3850 loop_within_domain(mpl, code->arg.loop.domain, info,
3851 iter_num_func);
3852 value = info->value;
3853 }
3854 break;
3855 case O_PROD:
3856 /* multiplication over domain */
3857 { struct iter_num_info _info, *info = &_info;
3858 info->code = code;
3859 info->value = 1.0;
3860 loop_within_domain(mpl, code->arg.loop.domain, info,
3861 iter_num_func);
3862 value = info->value;
3863 }
3864 break;
3865 case O_MINIMUM:
3866 /* minimum over domain */
3867 { struct iter_num_info _info, *info = &_info;
3868 info->code = code;
3869 info->value = +DBL_MAX;
3870 loop_within_domain(mpl, code->arg.loop.domain, info,
3871 iter_num_func);
3872 if (info->value == +DBL_MAX)
3873 error(mpl, "min{} over empty set; result undefined");
3874 value = info->value;
3875 }
3876 break;
3877 case O_MAXIMUM:
3878 /* maximum over domain */
3879 { struct iter_num_info _info, *info = &_info;
3880 info->code = code;
3881 info->value = -DBL_MAX;
3882 loop_within_domain(mpl, code->arg.loop.domain, info,
3883 iter_num_func);
3884 if (info->value == -DBL_MAX)
3885 error(mpl, "max{} over empty set; result undefined");
3886 value = info->value;
3887 }
3888 break;
3889 default:
3890 xassert(code != code);
3891 }
3892 /* save resultant value */
3893 xassert(!code->valid);
3894 code->valid = 1;
3895 code->value.num = value;
3896 done: return value;
3897 }
3899 /*----------------------------------------------------------------------
3900 -- eval_symbolic - evaluate pseudo-code to determine symbolic value.
3901 --
3902 -- This routine evaluates specified pseudo-code to determine resultant
3903 -- symbolic value, which is returned on exit. */
3905 SYMBOL *eval_symbolic(MPL *mpl, CODE *code)
3906 { SYMBOL *value;
3907 xassert(code != NULL);
3908 xassert(code->type == A_SYMBOLIC);
3909 xassert(code->dim == 0);
3910 /* if the operation has a side effect, invalidate and delete the
3911 resultant value */
3912 if (code->vflag && code->valid)
3913 { code->valid = 0;
3914 delete_value(mpl, code->type, &code->value);
3915 }
3916 /* if resultant value is valid, no evaluation is needed */
3917 if (code->valid)
3918 { value = copy_symbol(mpl, code->value.sym);
3919 goto done;
3920 }
3921 /* evaluate pseudo-code recursively */
3922 switch (code->op)
3923 { case O_STRING:
3924 /* take character string */
3925 value = create_symbol_str(mpl, create_string(mpl,
3926 code->arg.str));
3927 break;
3928 case O_INDEX:
3929 /* take dummy index */
3930 xassert(code->arg.index.slot->value != NULL);
3931 value = copy_symbol(mpl, code->arg.index.slot->value);
3932 break;
3933 case O_MEMSYM:
3934 /* take member of symbolic parameter */
3935 { TUPLE *tuple;
3936 ARG_LIST *e;
3937 tuple = create_tuple(mpl);
3938 for (e = code->arg.par.list; e != NULL; e = e->next)
3939 tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
3940 e->x));
3941 value = eval_member_sym(mpl, code->arg.par.par, tuple);
3942 delete_tuple(mpl, tuple);
3943 }
3944 break;
3945 case O_CVTSYM:
3946 /* conversion to symbolic */
3947 value = create_symbol_num(mpl, eval_numeric(mpl,
3948 code->arg.arg.x));
3949 break;
3950 case O_CONCAT:
3951 /* concatenation */
3952 value = concat_symbols(mpl,
3953 eval_symbolic(mpl, code->arg.arg.x),
3954 eval_symbolic(mpl, code->arg.arg.y));
3955 break;
3956 case O_FORK:
3957 /* if-then-else */
3958 if (eval_logical(mpl, code->arg.arg.x))
3959 value = eval_symbolic(mpl, code->arg.arg.y);
3960 else if (code->arg.arg.z == NULL)
3961 value = create_symbol_num(mpl, 0.0);
3962 else
3963 value = eval_symbolic(mpl, code->arg.arg.z);
3964 break;
3965 case O_SUBSTR:
3966 case O_SUBSTR3:
3967 { double pos, len;
3968 char str[MAX_LENGTH+1];
3969 value = eval_symbolic(mpl, code->arg.arg.x);
3970 if (value->str == NULL)
3971 sprintf(str, "%.*g", DBL_DIG, value->num);
3972 else
3973 fetch_string(mpl, value->str, str);
3974 delete_symbol(mpl, value);
3975 if (code->op == O_SUBSTR)
3976 { pos = eval_numeric(mpl, code->arg.arg.y);
3977 if (pos != floor(pos))
3978 error(mpl, "substr('...', %.*g); non-integer secon"
3979 "d argument", DBL_DIG, pos);
3980 if (pos < 1 || pos > strlen(str) + 1)
3981 error(mpl, "substr('...', %.*g); substring out of "
3982 "range", DBL_DIG, pos);
3983 }
3984 else
3985 { pos = eval_numeric(mpl, code->arg.arg.y);
3986 len = eval_numeric(mpl, code->arg.arg.z);
3987 if (pos != floor(pos) || len != floor(len))
3988 error(mpl, "substr('...', %.*g, %.*g); non-integer"
3989 " second and/or third argument", DBL_DIG, pos,
3990 DBL_DIG, len);
3991 if (pos < 1 || len < 0 || pos + len > strlen(str) + 1)
3992 error(mpl, "substr('...', %.*g, %.*g); substring o"
3993 "ut of range", DBL_DIG, pos, DBL_DIG, len);
3994 str[(int)pos + (int)len - 1] = '\0';
3995 }
3996 value = create_symbol_str(mpl, create_string(mpl, str +
3997 (int)pos - 1));
3998 }
3999 break;
4000 case O_TIME2STR:
4001 { double num;
4002 SYMBOL *sym;
4003 char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];
4004 num = eval_numeric(mpl, code->arg.arg.x);
4005 sym = eval_symbolic(mpl, code->arg.arg.y);
4006 if (sym->str == NULL)
4007 sprintf(fmt, "%.*g", DBL_DIG, sym->num);
4008 else
4009 fetch_string(mpl, sym->str, fmt);
4010 delete_symbol(mpl, sym);
4011 fn_time2str(mpl, str, num, fmt);
4012 value = create_symbol_str(mpl, create_string(mpl, str));
4013 }
4014 break;
4015 default:
4016 xassert(code != code);
4017 }
4018 /* save resultant value */
4019 xassert(!code->valid);
4020 code->valid = 1;
4021 code->value.sym = copy_symbol(mpl, value);
4022 done: return value;
4023 }
4025 /*----------------------------------------------------------------------
4026 -- eval_logical - evaluate pseudo-code to determine logical value.
4027 --
4028 -- This routine evaluates specified pseudo-code to determine resultant
4029 -- logical value, which is returned on exit. */
4031 struct iter_log_info
4032 { /* working info used by the routine iter_log_func */
4033 CODE *code;
4034 /* pseudo-code for iterated operation to be performed */
4035 int value;
4036 /* resultant value */
4037 };
4039 static int iter_log_func(MPL *mpl, void *_info)
4040 { /* this is auxiliary routine used to perform iterated operation
4041 on logical "integrand" within domain scope */
4042 struct iter_log_info *info = _info;
4043 int ret = 0;
4044 switch (info->code->op)
4045 { case O_FORALL:
4046 /* conjunction over domain */
4047 info->value &= eval_logical(mpl, info->code->arg.loop.x);
4048 if (!info->value) ret = 1;
4049 break;
4050 case O_EXISTS:
4051 /* disjunction over domain */
4052 info->value |= eval_logical(mpl, info->code->arg.loop.x);
4053 if (info->value) ret = 1;
4054 break;
4055 default:
4056 xassert(info != info);
4057 }
4058 return ret;
4059 }
4061 int eval_logical(MPL *mpl, CODE *code)
4062 { int value;
4063 xassert(code->type == A_LOGICAL);
4064 xassert(code->dim == 0);
4065 /* if the operation has a side effect, invalidate and delete the
4066 resultant value */
4067 if (code->vflag && code->valid)
4068 { code->valid = 0;
4069 delete_value(mpl, code->type, &code->value);
4070 }
4071 /* if resultant value is valid, no evaluation is needed */
4072 if (code->valid)
4073 { value = code->value.bit;
4074 goto done;
4075 }
4076 /* evaluate pseudo-code recursively */
4077 switch (code->op)
4078 { case O_CVTLOG:
4079 /* conversion to logical */
4080 value = (eval_numeric(mpl, code->arg.arg.x) != 0.0);
4081 break;
4082 case O_NOT:
4083 /* negation (logical "not") */
4084 value = !eval_logical(mpl, code->arg.arg.x);
4085 break;
4086 case O_LT:
4087 /* comparison on 'less than' */
4088 #if 0 /* 02/VIII-2008 */
4089 value = (eval_numeric(mpl, code->arg.arg.x) <
4090 eval_numeric(mpl, code->arg.arg.y));
4091 #else
4092 xassert(code->arg.arg.x != NULL);
4093 if (code->arg.arg.x->type == A_NUMERIC)
4094 value = (eval_numeric(mpl, code->arg.arg.x) <
4095 eval_numeric(mpl, code->arg.arg.y));
4096 else
4097 { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
4098 SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
4099 value = (compare_symbols(mpl, sym1, sym2) < 0);
4100 delete_symbol(mpl, sym1);
4101 delete_symbol(mpl, sym2);
4102 }
4103 #endif
4104 break;
4105 case O_LE:
4106 /* comparison on 'not greater than' */
4107 #if 0 /* 02/VIII-2008 */
4108 value = (eval_numeric(mpl, code->arg.arg.x) <=
4109 eval_numeric(mpl, code->arg.arg.y));
4110 #else
4111 xassert(code->arg.arg.x != NULL);
4112 if (code->arg.arg.x->type == A_NUMERIC)
4113 value = (eval_numeric(mpl, code->arg.arg.x) <=
4114 eval_numeric(mpl, code->arg.arg.y));
4115 else
4116 { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
4117 SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
4118 value = (compare_symbols(mpl, sym1, sym2) <= 0);
4119 delete_symbol(mpl, sym1);
4120 delete_symbol(mpl, sym2);
4121 }
4122 #endif
4123 break;
4124 case O_EQ:
4125 /* comparison on 'equal to' */
4126 xassert(code->arg.arg.x != NULL);
4127 if (code->arg.arg.x->type == A_NUMERIC)
4128 value = (eval_numeric(mpl, code->arg.arg.x) ==
4129 eval_numeric(mpl, code->arg.arg.y));
4130 else
4131 { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
4132 SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
4133 value = (compare_symbols(mpl, sym1, sym2) == 0);
4134 delete_symbol(mpl, sym1);
4135 delete_symbol(mpl, sym2);
4136 }
4137 break;
4138 case O_GE:
4139 /* comparison on 'not less than' */
4140 #if 0 /* 02/VIII-2008 */
4141 value = (eval_numeric(mpl, code->arg.arg.x) >=
4142 eval_numeric(mpl, code->arg.arg.y));
4143 #else
4144 xassert(code->arg.arg.x != NULL);
4145 if (code->arg.arg.x->type == A_NUMERIC)
4146 value = (eval_numeric(mpl, code->arg.arg.x) >=
4147 eval_numeric(mpl, code->arg.arg.y));
4148 else
4149 { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
4150 SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
4151 value = (compare_symbols(mpl, sym1, sym2) >= 0);
4152 delete_symbol(mpl, sym1);
4153 delete_symbol(mpl, sym2);
4154 }
4155 #endif
4156 break;
4157 case O_GT:
4158 /* comparison on 'greater than' */
4159 #if 0 /* 02/VIII-2008 */
4160 value = (eval_numeric(mpl, code->arg.arg.x) >
4161 eval_numeric(mpl, code->arg.arg.y));
4162 #else
4163 xassert(code->arg.arg.x != NULL);
4164 if (code->arg.arg.x->type == A_NUMERIC)
4165 value = (eval_numeric(mpl, code->arg.arg.x) >
4166 eval_numeric(mpl, code->arg.arg.y));
4167 else
4168 { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
4169 SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
4170 value = (compare_symbols(mpl, sym1, sym2) > 0);
4171 delete_symbol(mpl, sym1);
4172 delete_symbol(mpl, sym2);
4173 }
4174 #endif
4175 break;
4176 case O_NE:
4177 /* comparison on 'not equal to' */
4178 xassert(code->arg.arg.x != NULL);
4179 if (code->arg.arg.x->type == A_NUMERIC)
4180 value = (eval_numeric(mpl, code->arg.arg.x) !=
4181 eval_numeric(mpl, code->arg.arg.y));
4182 else
4183 { SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);
4184 SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);
4185 value = (compare_symbols(mpl, sym1, sym2) != 0);
4186 delete_symbol(mpl, sym1);
4187 delete_symbol(mpl, sym2);
4188 }
4189 break;
4190 case O_AND:
4191 /* conjunction (logical "and") */
4192 value = eval_logical(mpl, code->arg.arg.x) &&
4193 eval_logical(mpl, code->arg.arg.y);
4194 break;
4195 case O_OR:
4196 /* disjunction (logical "or") */
4197 value = eval_logical(mpl, code->arg.arg.x) ||
4198 eval_logical(mpl, code->arg.arg.y);
4199 break;
4200 case O_IN:
4201 /* test on 'x in Y' */
4202 { TUPLE *tuple;
4203 tuple = eval_tuple(mpl, code->arg.arg.x);
4204 value = is_member(mpl, code->arg.arg.y, tuple);
4205 delete_tuple(mpl, tuple);
4206 }
4207 break;
4208 case O_NOTIN:
4209 /* test on 'x not in Y' */
4210 { TUPLE *tuple;
4211 tuple = eval_tuple(mpl, code->arg.arg.x);
4212 value = !is_member(mpl, code->arg.arg.y, tuple);
4213 delete_tuple(mpl, tuple);
4214 }
4215 break;
4216 case O_WITHIN:
4217 /* test on 'X within Y' */
4218 { ELEMSET *set;
4219 MEMBER *memb;
4220 set = eval_elemset(mpl, code->arg.arg.x);
4221 value = 1;
4222 for (memb = set->head; memb != NULL; memb = memb->next)
4223 { if (!is_member(mpl, code->arg.arg.y, memb->tuple))
4224 { value = 0;
4225 break;
4226 }
4227 }
4228 delete_elemset(mpl, set);
4229 }
4230 break;
4231 case O_NOTWITHIN:
4232 /* test on 'X not within Y' */
4233 { ELEMSET *set;
4234 MEMBER *memb;
4235 set = eval_elemset(mpl, code->arg.arg.x);
4236 value = 1;
4237 for (memb = set->head; memb != NULL; memb = memb->next)
4238 { if (is_member(mpl, code->arg.arg.y, memb->tuple))
4239 { value = 0;
4240 break;
4241 }
4242 }
4243 delete_elemset(mpl, set);
4244 }
4245 break;
4246 case O_FORALL:
4247 /* conjunction (A-quantification) */
4248 { struct iter_log_info _info, *info = &_info;
4249 info->code = code;
4250 info->value = 1;
4251 loop_within_domain(mpl, code->arg.loop.domain, info,
4252 iter_log_func);
4253 value = info->value;
4254 }
4255 break;
4256 case O_EXISTS:
4257 /* disjunction (E-quantification) */
4258 { struct iter_log_info _info, *info = &_info;
4259 info->code = code;
4260 info->value = 0;
4261 loop_within_domain(mpl, code->arg.loop.domain, info,
4262 iter_log_func);
4263 value = info->value;
4264 }
4265 break;
4266 default:
4267 xassert(code != code);
4268 }
4269 /* save resultant value */
4270 xassert(!code->valid);
4271 code->valid = 1;
4272 code->value.bit = value;
4273 done: return value;
4274 }
4276 /*----------------------------------------------------------------------
4277 -- eval_tuple - evaluate pseudo-code to construct n-tuple.
4278 --
4279 -- This routine evaluates specified pseudo-code to construct resultant
4280 -- n-tuple, which is returned on exit. */
4282 TUPLE *eval_tuple(MPL *mpl, CODE *code)
4283 { TUPLE *value;
4284 xassert(code != NULL);
4285 xassert(code->type == A_TUPLE);
4286 xassert(code->dim > 0);
4287 /* if the operation has a side effect, invalidate and delete the
4288 resultant value */
4289 if (code->vflag && code->valid)
4290 { code->valid = 0;
4291 delete_value(mpl, code->type, &code->value);
4292 }
4293 /* if resultant value is valid, no evaluation is needed */
4294 if (code->valid)
4295 { value = copy_tuple(mpl, code->value.tuple);
4296 goto done;
4297 }
4298 /* evaluate pseudo-code recursively */
4299 switch (code->op)
4300 { case O_TUPLE:
4301 /* make n-tuple */
4302 { ARG_LIST *e;
4303 value = create_tuple(mpl);
4304 for (e = code->arg.list; e != NULL; e = e->next)
4305 value = expand_tuple(mpl, value, eval_symbolic(mpl,
4306 e->x));
4307 }
4308 break;
4309 case O_CVTTUP:
4310 /* convert to 1-tuple */
4311 value = expand_tuple(mpl, create_tuple(mpl),
4312 eval_symbolic(mpl, code->arg.arg.x));
4313 break;
4314 default:
4315 xassert(code != code);
4316 }
4317 /* save resultant value */
4318 xassert(!code->valid);
4319 code->valid = 1;
4320 code->value.tuple = copy_tuple(mpl, value);
4321 done: return value;
4322 }
4324 /*----------------------------------------------------------------------
4325 -- eval_elemset - evaluate pseudo-code to construct elemental set.
4326 --
4327 -- This routine evaluates specified pseudo-code to construct resultant
4328 -- elemental set, which is returned on exit. */
4330 struct iter_set_info
4331 { /* working info used by the routine iter_set_func */
4332 CODE *code;
4333 /* pseudo-code for iterated operation to be performed */
4334 ELEMSET *value;
4335 /* resultant value */
4336 };
4338 static int iter_set_func(MPL *mpl, void *_info)
4339 { /* this is auxiliary routine used to perform iterated operation
4340 on n-tuple "integrand" within domain scope */
4341 struct iter_set_info *info = _info;
4342 TUPLE *tuple;
4343 switch (info->code->op)
4344 { case O_SETOF:
4345 /* compute next n-tuple and add it to the set; in this case
4346 duplicate n-tuples are silently ignored */
4347 tuple = eval_tuple(mpl, info->code->arg.loop.x);
4348 if (find_tuple(mpl, info->value, tuple) == NULL)
4349 add_tuple(mpl, info->value, tuple);
4350 else
4351 delete_tuple(mpl, tuple);
4352 break;
4353 case O_BUILD:
4354 /* construct next n-tuple using current values assigned to
4355 *free* dummy indices as its components and add it to the
4356 set; in this case duplicate n-tuples cannot appear */
4357 add_tuple(mpl, info->value, get_domain_tuple(mpl,
4358 info->code->arg.loop.domain));
4359 break;
4360 default:
4361 xassert(info != info);
4362 }
4363 return 0;
4364 }
4366 ELEMSET *eval_elemset(MPL *mpl, CODE *code)
4367 { ELEMSET *value;
4368 xassert(code != NULL);
4369 xassert(code->type == A_ELEMSET);
4370 xassert(code->dim > 0);
4371 /* if the operation has a side effect, invalidate and delete the
4372 resultant value */
4373 if (code->vflag && code->valid)
4374 { code->valid = 0;
4375 delete_value(mpl, code->type, &code->value);
4376 }
4377 /* if resultant value is valid, no evaluation is needed */
4378 if (code->valid)
4379 { value = copy_elemset(mpl, code->value.set);
4380 goto done;
4381 }
4382 /* evaluate pseudo-code recursively */
4383 switch (code->op)
4384 { case O_MEMSET:
4385 /* take member of set */
4386 { TUPLE *tuple;
4387 ARG_LIST *e;
4388 tuple = create_tuple(mpl);
4389 for (e = code->arg.set.list; e != NULL; e = e->next)
4390 tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
4391 e->x));
4392 value = copy_elemset(mpl,
4393 eval_member_set(mpl, code->arg.set.set, tuple));
4394 delete_tuple(mpl, tuple);
4395 }
4396 break;
4397 case O_MAKE:
4398 /* make elemental set of n-tuples */
4399 { ARG_LIST *e;
4400 value = create_elemset(mpl, code->dim);
4401 for (e = code->arg.list; e != NULL; e = e->next)
4402 check_then_add(mpl, value, eval_tuple(mpl, e->x));
4403 }
4404 break;
4405 case O_UNION:
4406 /* union of two elemental sets */
4407 value = set_union(mpl,
4408 eval_elemset(mpl, code->arg.arg.x),
4409 eval_elemset(mpl, code->arg.arg.y));
4410 break;
4411 case O_DIFF:
4412 /* difference between two elemental sets */
4413 value = set_diff(mpl,
4414 eval_elemset(mpl, code->arg.arg.x),
4415 eval_elemset(mpl, code->arg.arg.y));
4416 break;
4417 case O_SYMDIFF:
4418 /* symmetric difference between two elemental sets */
4419 value = set_symdiff(mpl,
4420 eval_elemset(mpl, code->arg.arg.x),
4421 eval_elemset(mpl, code->arg.arg.y));
4422 break;
4423 case O_INTER:
4424 /* intersection of two elemental sets */
4425 value = set_inter(mpl,
4426 eval_elemset(mpl, code->arg.arg.x),
4427 eval_elemset(mpl, code->arg.arg.y));
4428 break;
4429 case O_CROSS:
4430 /* cross (Cartesian) product of two elemental sets */
4431 value = set_cross(mpl,
4432 eval_elemset(mpl, code->arg.arg.x),
4433 eval_elemset(mpl, code->arg.arg.y));
4434 break;
4435 case O_DOTS:
4436 /* build "arithmetic" elemental set */
4437 value = create_arelset(mpl,
4438 eval_numeric(mpl, code->arg.arg.x),
4439 eval_numeric(mpl, code->arg.arg.y),
4440 code->arg.arg.z == NULL ? 1.0 : eval_numeric(mpl,
4441 code->arg.arg.z));
4442 break;
4443 case O_FORK:
4444 /* if-then-else */
4445 if (eval_logical(mpl, code->arg.arg.x))
4446 value = eval_elemset(mpl, code->arg.arg.y);
4447 else
4448 value = eval_elemset(mpl, code->arg.arg.z);
4449 break;
4450 case O_SETOF:
4451 /* compute elemental set */
4452 { struct iter_set_info _info, *info = &_info;
4453 info->code = code;
4454 info->value = create_elemset(mpl, code->dim);
4455 loop_within_domain(mpl, code->arg.loop.domain, info,
4456 iter_set_func);
4457 value = info->value;
4458 }
4459 break;
4460 case O_BUILD:
4461 /* build elemental set identical to domain set */
4462 { struct iter_set_info _info, *info = &_info;
4463 info->code = code;
4464 info->value = create_elemset(mpl, code->dim);
4465 loop_within_domain(mpl, code->arg.loop.domain, info,
4466 iter_set_func);
4467 value = info->value;
4468 }
4469 break;
4470 default:
4471 xassert(code != code);
4472 }
4473 /* save resultant value */
4474 xassert(!code->valid);
4475 code->valid = 1;
4476 code->value.set = copy_elemset(mpl, value);
4477 done: return value;
4478 }
4480 /*----------------------------------------------------------------------
4481 -- is_member - check if n-tuple is in set specified by pseudo-code.
4482 --
4483 -- This routine checks if given n-tuple is a member of elemental set
4484 -- specified in the form of pseudo-code (i.e. by expression).
4485 --
4486 -- The n-tuple may have more components that dimension of the elemental
4487 -- set, in which case the extra components are ignored. */
4489 static void null_func(MPL *mpl, void *info)
4490 { /* this is dummy routine used to enter the domain scope */
4491 xassert(mpl == mpl);
4492 xassert(info == NULL);
4493 return;
4494 }
4496 int is_member(MPL *mpl, CODE *code, TUPLE *tuple)
4497 { int value;
4498 xassert(code != NULL);
4499 xassert(code->type == A_ELEMSET);
4500 xassert(code->dim > 0);
4501 xassert(tuple != NULL);
4502 switch (code->op)
4503 { case O_MEMSET:
4504 /* check if given n-tuple is member of elemental set, which
4505 is assigned to member of model set */
4506 { ARG_LIST *e;
4507 TUPLE *temp;
4508 ELEMSET *set;
4509 /* evaluate reference to elemental set */
4510 temp = create_tuple(mpl);
4511 for (e = code->arg.set.list; e != NULL; e = e->next)
4512 temp = expand_tuple(mpl, temp, eval_symbolic(mpl,
4513 e->x));
4514 set = eval_member_set(mpl, code->arg.set.set, temp);
4515 delete_tuple(mpl, temp);
4516 /* check if the n-tuple is contained in the set array */
4517 temp = build_subtuple(mpl, tuple, set->dim);
4518 value = (find_tuple(mpl, set, temp) != NULL);
4519 delete_tuple(mpl, temp);
4520 }
4521 break;
4522 case O_MAKE:
4523 /* check if given n-tuple is member of literal set */
4524 { ARG_LIST *e;
4525 TUPLE *temp, *that;
4526 value = 0;
4527 temp = build_subtuple(mpl, tuple, code->dim);
4528 for (e = code->arg.list; e != NULL; e = e->next)
4529 { that = eval_tuple(mpl, e->x);
4530 value = (compare_tuples(mpl, temp, that) == 0);
4531 delete_tuple(mpl, that);
4532 if (value) break;
4533 }
4534 delete_tuple(mpl, temp);
4535 }
4536 break;
4537 case O_UNION:
4538 value = is_member(mpl, code->arg.arg.x, tuple) ||
4539 is_member(mpl, code->arg.arg.y, tuple);
4540 break;
4541 case O_DIFF:
4542 value = is_member(mpl, code->arg.arg.x, tuple) &&
4543 !is_member(mpl, code->arg.arg.y, tuple);
4544 break;
4545 case O_SYMDIFF:
4546 { int in1 = is_member(mpl, code->arg.arg.x, tuple);
4547 int in2 = is_member(mpl, code->arg.arg.y, tuple);
4548 value = (in1 && !in2) || (!in1 && in2);
4549 }
4550 break;
4551 case O_INTER:
4552 value = is_member(mpl, code->arg.arg.x, tuple) &&
4553 is_member(mpl, code->arg.arg.y, tuple);
4554 break;
4555 case O_CROSS:
4556 { int j;
4557 value = is_member(mpl, code->arg.arg.x, tuple);
4558 if (value)
4559 { for (j = 1; j <= code->arg.arg.x->dim; j++)
4560 { xassert(tuple != NULL);
4561 tuple = tuple->next;
4562 }
4563 value = is_member(mpl, code->arg.arg.y, tuple);
4564 }
4565 }
4566 break;
4567 case O_DOTS:
4568 /* check if given 1-tuple is member of "arithmetic" set */
4569 { int j;
4570 double x, t0, tf, dt;
4571 xassert(code->dim == 1);
4572 /* compute "parameters" of the "arithmetic" set */
4573 t0 = eval_numeric(mpl, code->arg.arg.x);
4574 tf = eval_numeric(mpl, code->arg.arg.y);
4575 if (code->arg.arg.z == NULL)
4576 dt = 1.0;
4577 else
4578 dt = eval_numeric(mpl, code->arg.arg.z);
4579 /* make sure the parameters are correct */
4580 arelset_size(mpl, t0, tf, dt);
4581 /* if component of 1-tuple is symbolic, not numeric, the
4582 1-tuple cannot be member of "arithmetic" set */
4583 xassert(tuple->sym != NULL);
4584 if (tuple->sym->str != NULL)
4585 { value = 0;
4586 break;
4587 }
4588 /* determine numeric value of the component */
4589 x = tuple->sym->num;
4590 /* if the component value is out of the set range, the
4591 1-tuple is not in the set */
4592 if (dt > 0.0 && !(t0 <= x && x <= tf) ||
4593 dt < 0.0 && !(tf <= x && x <= t0))
4594 { value = 0;
4595 break;
4596 }
4597 /* estimate ordinal number of the 1-tuple in the set */
4598 j = (int)(((x - t0) / dt) + 0.5) + 1;
4599 /* perform the main check */
4600 value = (arelset_member(mpl, t0, tf, dt, j) == x);
4601 }
4602 break;
4603 case O_FORK:
4604 /* check if given n-tuple is member of conditional set */
4605 if (eval_logical(mpl, code->arg.arg.x))
4606 value = is_member(mpl, code->arg.arg.y, tuple);
4607 else
4608 value = is_member(mpl, code->arg.arg.z, tuple);
4609 break;
4610 case O_SETOF:
4611 /* check if given n-tuple is member of computed set */
4612 /* it is not clear how to efficiently perform the check not
4613 computing the entire elemental set :+( */
4614 error(mpl, "implementation restriction; in/within setof{} n"
4615 "ot allowed");
4616 break;
4617 case O_BUILD:
4618 /* check if given n-tuple is member of domain set */
4619 { TUPLE *temp;
4620 temp = build_subtuple(mpl, tuple, code->dim);
4621 /* try to enter the domain scope; if it is successful,
4622 the n-tuple is in the domain set */
4623 value = (eval_within_domain(mpl, code->arg.loop.domain,
4624 temp, NULL, null_func) == 0);
4625 delete_tuple(mpl, temp);
4626 }
4627 break;
4628 default:
4629 xassert(code != code);
4630 }
4631 return value;
4632 }
4634 /*----------------------------------------------------------------------
4635 -- eval_formula - evaluate pseudo-code to construct linear form.
4636 --
4637 -- This routine evaluates specified pseudo-code to construct resultant
4638 -- linear form, which is returned on exit. */
4640 struct iter_form_info
4641 { /* working info used by the routine iter_form_func */
4642 CODE *code;
4643 /* pseudo-code for iterated operation to be performed */
4644 FORMULA *value;
4645 /* resultant value */
4646 FORMULA *tail;
4647 /* pointer to the last term */
4648 };
4650 static int iter_form_func(MPL *mpl, void *_info)
4651 { /* this is auxiliary routine used to perform iterated operation
4652 on linear form "integrand" within domain scope */
4653 struct iter_form_info *info = _info;
4654 switch (info->code->op)
4655 { case O_SUM:
4656 /* summation over domain */
4657 #if 0
4658 info->value =
4659 linear_comb(mpl,
4660 +1.0, info->value,
4661 +1.0, eval_formula(mpl, info->code->arg.loop.x));
4662 #else
4663 /* the routine linear_comb needs to look through all terms
4664 of both linear forms to reduce identical terms, so using
4665 it here is not a good idea (for example, evaluation of
4666 sum{i in 1..n} x[i] required quadratic time); the better
4667 idea is to gather all terms of the integrand in one list
4668 and reduce identical terms only once after all terms of
4669 the resultant linear form have been evaluated */
4670 { FORMULA *form, *term;
4671 form = eval_formula(mpl, info->code->arg.loop.x);
4672 if (info->value == NULL)
4673 { xassert(info->tail == NULL);
4674 info->value = form;
4675 }
4676 else
4677 { xassert(info->tail != NULL);
4678 info->tail->next = form;
4679 }
4680 for (term = form; term != NULL; term = term->next)
4681 info->tail = term;
4682 }
4683 #endif
4684 break;
4685 default:
4686 xassert(info != info);
4687 }
4688 return 0;
4689 }
4691 FORMULA *eval_formula(MPL *mpl, CODE *code)
4692 { FORMULA *value;
4693 xassert(code != NULL);
4694 xassert(code->type == A_FORMULA);
4695 xassert(code->dim == 0);
4696 /* if the operation has a side effect, invalidate and delete the
4697 resultant value */
4698 if (code->vflag && code->valid)
4699 { code->valid = 0;
4700 delete_value(mpl, code->type, &code->value);
4701 }
4702 /* if resultant value is valid, no evaluation is needed */
4703 if (code->valid)
4704 { value = copy_formula(mpl, code->value.form);
4705 goto done;
4706 }
4707 /* evaluate pseudo-code recursively */
4708 switch (code->op)
4709 { case O_MEMVAR:
4710 /* take member of variable */
4711 { TUPLE *tuple;
4712 ARG_LIST *e;
4713 tuple = create_tuple(mpl);
4714 for (e = code->arg.var.list; e != NULL; e = e->next)
4715 tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,
4716 e->x));
4717 #if 1 /* 15/V-2010 */
4718 xassert(code->arg.var.suff == DOT_NONE);
4719 #endif
4720 value = single_variable(mpl,
4721 eval_member_var(mpl, code->arg.var.var, tuple));
4722 delete_tuple(mpl, tuple);
4723 }
4724 break;
4725 case O_CVTLFM:
4726 /* convert to linear form */
4727 value = constant_term(mpl, eval_numeric(mpl,
4728 code->arg.arg.x));
4729 break;
4730 case O_PLUS:
4731 /* unary plus */
4732 value = linear_comb(mpl,
4733 0.0, constant_term(mpl, 0.0),
4734 +1.0, eval_formula(mpl, code->arg.arg.x));
4735 break;
4736 case O_MINUS:
4737 /* unary minus */
4738 value = linear_comb(mpl,
4739 0.0, constant_term(mpl, 0.0),
4740 -1.0, eval_formula(mpl, code->arg.arg.x));
4741 break;
4742 case O_ADD:
4743 /* addition */
4744 value = linear_comb(mpl,
4745 +1.0, eval_formula(mpl, code->arg.arg.x),
4746 +1.0, eval_formula(mpl, code->arg.arg.y));
4747 break;
4748 case O_SUB:
4749 /* subtraction */
4750 value = linear_comb(mpl,
4751 +1.0, eval_formula(mpl, code->arg.arg.x),
4752 -1.0, eval_formula(mpl, code->arg.arg.y));
4753 break;
4754 case O_MUL:
4755 /* multiplication */
4756 xassert(code->arg.arg.x != NULL);
4757 xassert(code->arg.arg.y != NULL);
4758 if (code->arg.arg.x->type == A_NUMERIC)
4759 { xassert(code->arg.arg.y->type == A_FORMULA);
4760 value = linear_comb(mpl,
4761 eval_numeric(mpl, code->arg.arg.x),
4762 eval_formula(mpl, code->arg.arg.y),
4763 0.0, constant_term(mpl, 0.0));
4764 }
4765 else
4766 { xassert(code->arg.arg.x->type == A_FORMULA);
4767 xassert(code->arg.arg.y->type == A_NUMERIC);
4768 value = linear_comb(mpl,
4769 eval_numeric(mpl, code->arg.arg.y),
4770 eval_formula(mpl, code->arg.arg.x),
4771 0.0, constant_term(mpl, 0.0));
4772 }
4773 break;
4774 case O_DIV:
4775 /* division */
4776 value = linear_comb(mpl,
4777 fp_div(mpl, 1.0, eval_numeric(mpl, code->arg.arg.y)),
4778 eval_formula(mpl, code->arg.arg.x),
4779 0.0, constant_term(mpl, 0.0));
4780 break;
4781 case O_FORK:
4782 /* if-then-else */
4783 if (eval_logical(mpl, code->arg.arg.x))
4784 value = eval_formula(mpl, code->arg.arg.y);
4785 else if (code->arg.arg.z == NULL)
4786 value = constant_term(mpl, 0.0);
4787 else
4788 value = eval_formula(mpl, code->arg.arg.z);
4789 break;
4790 case O_SUM:
4791 /* summation over domain */
4792 { struct iter_form_info _info, *info = &_info;
4793 info->code = code;
4794 info->value = constant_term(mpl, 0.0);
4795 info->tail = NULL;
4796 loop_within_domain(mpl, code->arg.loop.domain, info,
4797 iter_form_func);
4798 value = reduce_terms(mpl, info->value);
4799 }
4800 break;
4801 default:
4802 xassert(code != code);
4803 }
4804 /* save resultant value */
4805 xassert(!code->valid);
4806 code->valid = 1;
4807 code->value.form = copy_formula(mpl, value);
4808 done: return value;
4809 }
4811 /*----------------------------------------------------------------------
4812 -- clean_code - clean pseudo-code.
4813 --
4814 -- This routine recursively cleans specified pseudo-code that assumes
4815 -- deleting all temporary resultant values. */
4817 void clean_code(MPL *mpl, CODE *code)
4818 { ARG_LIST *e;
4819 /* if no pseudo-code is specified, do nothing */
4820 if (code == NULL) goto done;
4821 /* if resultant value is valid (exists), delete it */
4822 if (code->valid)
4823 { code->valid = 0;
4824 delete_value(mpl, code->type, &code->value);
4825 }
4826 /* recursively clean pseudo-code for operands */
4827 switch (code->op)
4828 { case O_NUMBER:
4829 case O_STRING:
4830 case O_INDEX:
4831 break;
4832 case O_MEMNUM:
4833 case O_MEMSYM:
4834 for (e = code->arg.par.list; e != NULL; e = e->next)
4835 clean_code(mpl, e->x);
4836 break;
4837 case O_MEMSET:
4838 for (e = code->arg.set.list; e != NULL; e = e->next)
4839 clean_code(mpl, e->x);
4840 break;
4841 case O_MEMVAR:
4842 for (e = code->arg.var.list; e != NULL; e = e->next)
4843 clean_code(mpl, e->x);
4844 break;
4845 #if 1 /* 15/V-2010 */
4846 case O_MEMCON:
4847 for (e = code->arg.con.list; e != NULL; e = e->next)
4848 clean_code(mpl, e->x);
4849 break;
4850 #endif
4851 case O_TUPLE:
4852 case O_MAKE:
4853 for (e = code->arg.list; e != NULL; e = e->next)
4854 clean_code(mpl, e->x);
4855 break;
4856 case O_SLICE:
4857 xassert(code != code);
4858 case O_IRAND224:
4859 case O_UNIFORM01:
4860 case O_NORMAL01:
4861 case O_GMTIME:
4862 break;
4863 case O_CVTNUM:
4864 case O_CVTSYM:
4865 case O_CVTLOG:
4866 case O_CVTTUP:
4867 case O_CVTLFM:
4868 case O_PLUS:
4869 case O_MINUS:
4870 case O_NOT:
4871 case O_ABS:
4872 case O_CEIL:
4873 case O_FLOOR:
4874 case O_EXP:
4875 case O_LOG:
4876 case O_LOG10:
4877 case O_SQRT:
4878 case O_SIN:
4879 case O_COS:
4880 case O_ATAN:
4881 case O_ROUND:
4882 case O_TRUNC:
4883 case O_CARD:
4884 case O_LENGTH:
4885 /* unary operation */
4886 clean_code(mpl, code->arg.arg.x);
4887 break;
4888 case O_ADD:
4889 case O_SUB:
4890 case O_LESS:
4891 case O_MUL:
4892 case O_DIV:
4893 case O_IDIV:
4894 case O_MOD:
4895 case O_POWER:
4896 case O_ATAN2:
4897 case O_ROUND2:
4898 case O_TRUNC2:
4899 case O_UNIFORM:
4900 case O_NORMAL:
4901 case O_CONCAT:
4902 case O_LT:
4903 case O_LE:
4904 case O_EQ:
4905 case O_GE:
4906 case O_GT:
4907 case O_NE:
4908 case O_AND:
4909 case O_OR:
4910 case O_UNION:
4911 case O_DIFF:
4912 case O_SYMDIFF:
4913 case O_INTER:
4914 case O_CROSS:
4915 case O_IN:
4916 case O_NOTIN:
4917 case O_WITHIN:
4918 case O_NOTWITHIN:
4919 case O_SUBSTR:
4920 case O_STR2TIME:
4921 case O_TIME2STR:
4922 /* binary operation */
4923 clean_code(mpl, code->arg.arg.x);
4924 clean_code(mpl, code->arg.arg.y);
4925 break;
4926 case O_DOTS:
4927 case O_FORK:
4928 case O_SUBSTR3:
4929 /* ternary operation */
4930 clean_code(mpl, code->arg.arg.x);
4931 clean_code(mpl, code->arg.arg.y);
4932 clean_code(mpl, code->arg.arg.z);
4933 break;
4934 case O_MIN:
4935 case O_MAX:
4936 /* n-ary operation */
4937 for (e = code->arg.list; e != NULL; e = e->next)
4938 clean_code(mpl, e->x);
4939 break;
4940 case O_SUM:
4941 case O_PROD:
4942 case O_MINIMUM:
4943 case O_MAXIMUM:
4944 case O_FORALL:
4945 case O_EXISTS:
4946 case O_SETOF:
4947 case O_BUILD:
4948 /* iterated operation */
4949 clean_domain(mpl, code->arg.loop.domain);
4950 clean_code(mpl, code->arg.loop.x);
4951 break;
4952 default:
4953 xassert(code->op != code->op);
4954 }
4955 done: return;
4956 }
4958 #if 1 /* 11/II-2008 */
4959 /**********************************************************************/
4960 /* * * DATA TABLES * * */
4961 /**********************************************************************/
4963 int mpl_tab_num_args(TABDCA *dca)
4964 { /* returns the number of arguments */
4965 return dca->na;
4966 }
4968 const char *mpl_tab_get_arg(TABDCA *dca, int k)
4969 { /* returns pointer to k-th argument */
4970 xassert(1 <= k && k <= dca->na);
4971 return dca->arg[k];
4972 }
4974 int mpl_tab_num_flds(TABDCA *dca)
4975 { /* returns the number of fields */
4976 return dca->nf;
4977 }
4979 const char *mpl_tab_get_name(TABDCA *dca, int k)
4980 { /* returns pointer to name of k-th field */
4981 xassert(1 <= k && k <= dca->nf);
4982 return dca->name[k];
4983 }
4985 int mpl_tab_get_type(TABDCA *dca, int k)
4986 { /* returns type of k-th field */
4987 xassert(1 <= k && k <= dca->nf);
4988 return dca->type[k];
4989 }
4991 double mpl_tab_get_num(TABDCA *dca, int k)
4992 { /* returns numeric value of k-th field */
4993 xassert(1 <= k && k <= dca->nf);
4994 xassert(dca->type[k] == 'N');
4995 return dca->num[k];
4996 }
4998 const char *mpl_tab_get_str(TABDCA *dca, int k)
4999 { /* returns pointer to string value of k-th field */
5000 xassert(1 <= k && k <= dca->nf);
5001 xassert(dca->type[k] == 'S');
5002 xassert(dca->str[k] != NULL);
5003 return dca->str[k];
5004 }
5006 void mpl_tab_set_num(TABDCA *dca, int k, double num)
5007 { /* assign numeric value to k-th field */
5008 xassert(1 <= k && k <= dca->nf);
5009 xassert(dca->type[k] == '?');
5010 dca->type[k] = 'N';
5011 dca->num[k] = num;
5012 return;
5013 }
5015 void mpl_tab_set_str(TABDCA *dca, int k, const char *str)
5016 { /* assign string value to k-th field */
5017 xassert(1 <= k && k <= dca->nf);
5018 xassert(dca->type[k] == '?');
5019 xassert(strlen(str) <= MAX_LENGTH);
5020 xassert(dca->str[k] != NULL);
5021 dca->type[k] = 'S';
5022 strcpy(dca->str[k], str);
5023 return;
5024 }
5026 static int write_func(MPL *mpl, void *info)
5027 { /* this is auxiliary routine to work within domain scope */
5028 TABLE *tab = info;
5029 TABDCA *dca = mpl->dca;
5030 TABOUT *out;
5031 SYMBOL *sym;
5032 int k;
5033 char buf[MAX_LENGTH+1];
5034 /* evaluate field values */
5035 k = 0;
5036 for (out = tab->u.out.list; out != NULL; out = out->next)
5037 { k++;
5038 switch (out->code->type)
5039 { case A_NUMERIC:
5040 dca->type[k] = 'N';
5041 dca->num[k] = eval_numeric(mpl, out->code);
5042 dca->str[k][0] = '\0';
5043 break;
5044 case A_SYMBOLIC:
5045 sym = eval_symbolic(mpl, out->code);
5046 if (sym->str == NULL)
5047 { dca->type[k] = 'N';
5048 dca->num[k] = sym->num;
5049 dca->str[k][0] = '\0';
5050 }
5051 else
5052 { dca->type[k] = 'S';
5053 dca->num[k] = 0.0;
5054 fetch_string(mpl, sym->str, buf);
5055 strcpy(dca->str[k], buf);
5056 }
5057 delete_symbol(mpl, sym);
5058 break;
5059 default:
5060 xassert(out != out);
5061 }
5062 }
5063 /* write record to output table */
5064 mpl_tab_drv_write(mpl);
5065 return 0;
5066 }
5068 void execute_table(MPL *mpl, TABLE *tab)
5069 { /* execute table statement */
5070 TABARG *arg;
5071 TABFLD *fld;
5072 TABIN *in;
5073 TABOUT *out;
5074 TABDCA *dca;
5075 SET *set;
5076 int k;
5077 char buf[MAX_LENGTH+1];
5078 /* allocate table driver communication area */
5079 xassert(mpl->dca == NULL);
5080 mpl->dca = dca = xmalloc(sizeof(TABDCA));
5081 dca->id = 0;
5082 dca->link = NULL;
5083 dca->na = 0;
5084 dca->arg = NULL;
5085 dca->nf = 0;
5086 dca->name = NULL;
5087 dca->type = NULL;
5088 dca->num = NULL;
5089 dca->str = NULL;
5090 /* allocate arguments */
5091 xassert(dca->na == 0);
5092 for (arg = tab->arg; arg != NULL; arg = arg->next)
5093 dca->na++;
5094 dca->arg = xcalloc(1+dca->na, sizeof(char *));
5095 #if 1 /* 28/IX-2008 */
5096 for (k = 1; k <= dca->na; k++) dca->arg[k] = NULL;
5097 #endif
5098 /* evaluate argument values */
5099 k = 0;
5100 for (arg = tab->arg; arg != NULL; arg = arg->next)
5101 { SYMBOL *sym;
5102 k++;
5103 xassert(arg->code->type == A_SYMBOLIC);
5104 sym = eval_symbolic(mpl, arg->code);
5105 if (sym->str == NULL)
5106 sprintf(buf, "%.*g", DBL_DIG, sym->num);
5107 else
5108 fetch_string(mpl, sym->str, buf);
5109 delete_symbol(mpl, sym);
5110 dca->arg[k] = xmalloc(strlen(buf)+1);
5111 strcpy(dca->arg[k], buf);
5112 }
5113 /* perform table input/output */
5114 switch (tab->type)
5115 { case A_INPUT: goto read_table;
5116 case A_OUTPUT: goto write_table;
5117 default: xassert(tab != tab);
5118 }
5119 read_table:
5120 /* read data from input table */
5121 /* add the only member to the control set and assign it empty
5122 elemental set */
5123 set = tab->u.in.set;
5124 if (set != NULL)
5125 { if (set->data)
5126 error(mpl, "%s already provided with data", set->name);
5127 xassert(set->array->head == NULL);
5128 add_member(mpl, set->array, NULL)->value.set =
5129 create_elemset(mpl, set->dimen);
5130 set->data = 1;
5131 }
5132 /* check parameters specified in the input list */
5133 for (in = tab->u.in.list; in != NULL; in = in->next)
5134 { if (in->par->data)
5135 error(mpl, "%s already provided with data", in->par->name);
5136 in->par->data = 1;
5137 }
5138 /* allocate and initialize fields */
5139 xassert(dca->nf == 0);
5140 for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
5141 dca->nf++;
5142 for (in = tab->u.in.list; in != NULL; in = in->next)
5143 dca->nf++;
5144 dca->name = xcalloc(1+dca->nf, sizeof(char *));
5145 dca->type = xcalloc(1+dca->nf, sizeof(int));
5146 dca->num = xcalloc(1+dca->nf, sizeof(double));
5147 dca->str = xcalloc(1+dca->nf, sizeof(char *));
5148 k = 0;
5149 for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
5150 { k++;
5151 dca->name[k] = fld->name;
5152 dca->type[k] = '?';
5153 dca->num[k] = 0.0;
5154 dca->str[k] = xmalloc(MAX_LENGTH+1);
5155 dca->str[k][0] = '\0';
5156 }
5157 for (in = tab->u.in.list; in != NULL; in = in->next)
5158 { k++;
5159 dca->name[k] = in->name;
5160 dca->type[k] = '?';
5161 dca->num[k] = 0.0;
5162 dca->str[k] = xmalloc(MAX_LENGTH+1);
5163 dca->str[k][0] = '\0';
5164 }
5165 /* open input table */
5166 mpl_tab_drv_open(mpl, 'R');
5167 /* read and process records */
5168 for (;;)
5169 { TUPLE *tup;
5170 /* reset field types */
5171 for (k = 1; k <= dca->nf; k++)
5172 dca->type[k] = '?';
5173 /* read next record */
5174 if (mpl_tab_drv_read(mpl)) break;
5175 /* all fields must be set by the driver */
5176 for (k = 1; k <= dca->nf; k++)
5177 { if (dca->type[k] == '?')
5178 error(mpl, "field %s missing in input table",
5179 dca->name[k]);
5180 }
5181 /* construct n-tuple */
5182 tup = create_tuple(mpl);
5183 k = 0;
5184 for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)
5185 { k++;
5186 xassert(k <= dca->nf);
5187 switch (dca->type[k])
5188 { case 'N':
5189 tup = expand_tuple(mpl, tup, create_symbol_num(mpl,
5190 dca->num[k]));
5191 break;
5192 case 'S':
5193 xassert(strlen(dca->str[k]) <= MAX_LENGTH);
5194 tup = expand_tuple(mpl, tup, create_symbol_str(mpl,
5195 create_string(mpl, dca->str[k])));
5196 break;
5197 default:
5198 xassert(dca != dca);
5199 }
5200 }
5201 /* add n-tuple just read to the control set */
5202 if (tab->u.in.set != NULL)
5203 check_then_add(mpl, tab->u.in.set->array->head->value.set,
5204 copy_tuple(mpl, tup));
5205 /* assign values to the parameters in the input list */
5206 for (in = tab->u.in.list; in != NULL; in = in->next)
5207 { MEMBER *memb;
5208 k++;
5209 xassert(k <= dca->nf);
5210 /* there must be no member with the same n-tuple */
5211 if (find_member(mpl, in->par->array, tup) != NULL)
5212 error(mpl, "%s%s already defined", in->par->name,
5213 format_tuple(mpl, '[', tup));
5214 /* create new parameter member with given n-tuple */
5215 memb = add_member(mpl, in->par->array, copy_tuple(mpl, tup))
5216 ;
5217 /* assign value to the parameter member */
5218 switch (in->par->type)
5219 { case A_NUMERIC:
5220 case A_INTEGER:
5221 case A_BINARY:
5222 if (dca->type[k] != 'N')
5223 error(mpl, "%s requires numeric data",
5224 in->par->name);
5225 memb->value.num = dca->num[k];
5226 break;
5227 case A_SYMBOLIC:
5228 switch (dca->type[k])
5229 { case 'N':
5230 memb->value.sym = create_symbol_num(mpl,
5231 dca->num[k]);
5232 break;
5233 case 'S':
5234 xassert(strlen(dca->str[k]) <= MAX_LENGTH);
5235 memb->value.sym = create_symbol_str(mpl,
5236 create_string(mpl,dca->str[k]));
5237 break;
5238 default:
5239 xassert(dca != dca);
5240 }
5241 break;
5242 default:
5243 xassert(in != in);
5244 }
5245 }
5246 /* n-tuple is no more needed */
5247 delete_tuple(mpl, tup);
5248 }
5249 /* close input table */
5250 mpl_tab_drv_close(mpl);
5251 goto done;
5252 write_table:
5253 /* write data to output table */
5254 /* allocate and initialize fields */
5255 xassert(dca->nf == 0);
5256 for (out = tab->u.out.list; out != NULL; out = out->next)
5257 dca->nf++;
5258 dca->name = xcalloc(1+dca->nf, sizeof(char *));
5259 dca->type = xcalloc(1+dca->nf, sizeof(int));
5260 dca->num = xcalloc(1+dca->nf, sizeof(double));
5261 dca->str = xcalloc(1+dca->nf, sizeof(char *));
5262 k = 0;
5263 for (out = tab->u.out.list; out != NULL; out = out->next)
5264 { k++;
5265 dca->name[k] = out->name;
5266 dca->type[k] = '?';
5267 dca->num[k] = 0.0;
5268 dca->str[k] = xmalloc(MAX_LENGTH+1);
5269 dca->str[k][0] = '\0';
5270 }
5271 /* open output table */
5272 mpl_tab_drv_open(mpl, 'W');
5273 /* evaluate fields and write records */
5274 loop_within_domain(mpl, tab->u.out.domain, tab, write_func);
5275 /* close output table */
5276 mpl_tab_drv_close(mpl);
5277 done: /* free table driver communication area */
5278 free_dca(mpl);
5279 return;
5280 }
5282 void free_dca(MPL *mpl)
5283 { /* free table driver communucation area */
5284 TABDCA *dca = mpl->dca;
5285 int k;
5286 if (dca != NULL)
5287 { if (dca->link != NULL)
5288 mpl_tab_drv_close(mpl);
5289 if (dca->arg != NULL)
5290 { for (k = 1; k <= dca->na; k++)
5291 #if 1 /* 28/IX-2008 */
5292 if (dca->arg[k] != NULL)
5293 #endif
5294 xfree(dca->arg[k]);
5295 xfree(dca->arg);
5296 }
5297 if (dca->name != NULL) xfree(dca->name);
5298 if (dca->type != NULL) xfree(dca->type);
5299 if (dca->num != NULL) xfree(dca->num);
5300 if (dca->str != NULL)
5301 { for (k = 1; k <= dca->nf; k++)
5302 xfree(dca->str[k]);
5303 xfree(dca->str);
5304 }
5305 xfree(dca), mpl->dca = NULL;
5306 }
5307 return;
5308 }
5310 void clean_table(MPL *mpl, TABLE *tab)
5311 { /* clean table statement */
5312 TABARG *arg;
5313 TABOUT *out;
5314 /* clean string list */
5315 for (arg = tab->arg; arg != NULL; arg = arg->next)
5316 clean_code(mpl, arg->code);
5317 switch (tab->type)
5318 { case A_INPUT:
5319 break;
5320 case A_OUTPUT:
5321 /* clean subscript domain */
5322 clean_domain(mpl, tab->u.out.domain);
5323 /* clean output list */
5324 for (out = tab->u.out.list; out != NULL; out = out->next)
5325 clean_code(mpl, out->code);
5326 break;
5327 default:
5328 xassert(tab != tab);
5329 }
5330 return;
5331 }
5332 #endif
5334 /**********************************************************************/
5335 /* * * MODEL STATEMENTS * * */
5336 /**********************************************************************/
5338 /*----------------------------------------------------------------------
5339 -- execute_check - execute check statement.
5340 --
5341 -- This routine executes specified check statement. */
5343 static int check_func(MPL *mpl, void *info)
5344 { /* this is auxiliary routine to work within domain scope */
5345 CHECK *chk = (CHECK *)info;
5346 if (!eval_logical(mpl, chk->code))
5347 error(mpl, "check%s failed", format_tuple(mpl, '[',
5348 get_domain_tuple(mpl, chk->domain)));
5349 return 0;
5350 }
5352 void execute_check(MPL *mpl, CHECK *chk)
5353 { loop_within_domain(mpl, chk->domain, chk, check_func);
5354 return;
5355 }
5357 /*----------------------------------------------------------------------
5358 -- clean_check - clean check statement.
5359 --
5360 -- This routine cleans specified check statement that assumes deleting
5361 -- all stuff dynamically allocated on generating/postsolving phase. */
5363 void clean_check(MPL *mpl, CHECK *chk)
5364 { /* clean subscript domain */
5365 clean_domain(mpl, chk->domain);
5366 /* clean pseudo-code for computing predicate */
5367 clean_code(mpl, chk->code);
5368 return;
5369 }
5371 /*----------------------------------------------------------------------
5372 -- execute_display - execute display statement.
5373 --
5374 -- This routine executes specified display statement. */
5376 static void display_set(MPL *mpl, SET *set, MEMBER *memb)
5377 { /* display member of model set */
5378 ELEMSET *s = memb->value.set;
5379 MEMBER *m;
5380 write_text(mpl, "%s%s%s\n", set->name,
5381 format_tuple(mpl, '[', memb->tuple),
5382 s->head == NULL ? " is empty" : ":");
5383 for (m = s->head; m != NULL; m = m->next)
5384 write_text(mpl, " %s\n", format_tuple(mpl, '(', m->tuple));
5385 return;
5386 }
5388 static void display_par(MPL *mpl, PARAMETER *par, MEMBER *memb)
5389 { /* display member of model parameter */
5390 switch (par->type)
5391 { case A_NUMERIC:
5392 case A_INTEGER:
5393 case A_BINARY:
5394 write_text(mpl, "%s%s = %.*g\n", par->name,
5395 format_tuple(mpl, '[', memb->tuple),
5396 DBL_DIG, memb->value.num);
5397 break;
5398 case A_SYMBOLIC:
5399 write_text(mpl, "%s%s = %s\n", par->name,
5400 format_tuple(mpl, '[', memb->tuple),
5401 format_symbol(mpl, memb->value.sym));
5402 break;
5403 default:
5404 xassert(par != par);
5405 }
5406 return;
5407 }
5409 #if 1 /* 15/V-2010 */
5410 static void display_var(MPL *mpl, VARIABLE *var, MEMBER *memb,
5411 int suff)
5412 { /* display member of model variable */
5413 if (suff == DOT_NONE || suff == DOT_VAL)
5414 write_text(mpl, "%s%s.val = %.*g\n", var->name,
5415 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5416 memb->value.var->prim);
5417 else if (suff == DOT_LB)
5418 write_text(mpl, "%s%s.lb = %.*g\n", var->name,
5419 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5420 memb->value.var->var->lbnd == NULL ? -DBL_MAX :
5421 memb->value.var->lbnd);
5422 else if (suff == DOT_UB)
5423 write_text(mpl, "%s%s.ub = %.*g\n", var->name,
5424 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5425 memb->value.var->var->ubnd == NULL ? +DBL_MAX :
5426 memb->value.var->ubnd);
5427 else if (suff == DOT_STATUS)
5428 write_text(mpl, "%s%s.status = %d\n", var->name, format_tuple
5429 (mpl, '[', memb->tuple), memb->value.var->stat);
5430 else if (suff == DOT_DUAL)
5431 write_text(mpl, "%s%s.dual = %.*g\n", var->name,
5432 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5433 memb->value.var->dual);
5434 else
5435 xassert(suff != suff);
5436 return;
5437 }
5438 #endif
5440 #if 1 /* 15/V-2010 */
5441 static void display_con(MPL *mpl, CONSTRAINT *con, MEMBER *memb,
5442 int suff)
5443 { /* display member of model constraint */
5444 if (suff == DOT_NONE || suff == DOT_VAL)
5445 write_text(mpl, "%s%s.val = %.*g\n", con->name,
5446 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5447 memb->value.con->prim);
5448 else if (suff == DOT_LB)
5449 write_text(mpl, "%s%s.lb = %.*g\n", con->name,
5450 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5451 memb->value.con->con->lbnd == NULL ? -DBL_MAX :
5452 memb->value.con->lbnd);
5453 else if (suff == DOT_UB)
5454 write_text(mpl, "%s%s.ub = %.*g\n", con->name,
5455 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5456 memb->value.con->con->ubnd == NULL ? +DBL_MAX :
5457 memb->value.con->ubnd);
5458 else if (suff == DOT_STATUS)
5459 write_text(mpl, "%s%s.status = %d\n", con->name, format_tuple
5460 (mpl, '[', memb->tuple), memb->value.con->stat);
5461 else if (suff == DOT_DUAL)
5462 write_text(mpl, "%s%s.dual = %.*g\n", con->name,
5463 format_tuple(mpl, '[', memb->tuple), DBL_DIG,
5464 memb->value.con->dual);
5465 else
5466 xassert(suff != suff);
5467 return;
5468 }
5469 #endif
5471 static void display_memb(MPL *mpl, CODE *code)
5472 { /* display member specified by pseudo-code */
5473 MEMBER memb;
5474 ARG_LIST *e;
5475 xassert(code->op == O_MEMNUM || code->op == O_MEMSYM
5476 || code->op == O_MEMSET || code->op == O_MEMVAR
5477 || code->op == O_MEMCON);
5478 memb.tuple = create_tuple(mpl);
5479 for (e = code->arg.par.list; e != NULL; e = e->next)
5480 memb.tuple = expand_tuple(mpl, memb.tuple, eval_symbolic(mpl,
5481 e->x));
5482 switch (code->op)
5483 { case O_MEMNUM:
5484 memb.value.num = eval_member_num(mpl, code->arg.par.par,
5485 memb.tuple);
5486 display_par(mpl, code->arg.par.par, &memb);
5487 break;
5488 case O_MEMSYM:
5489 memb.value.sym = eval_member_sym(mpl, code->arg.par.par,
5490 memb.tuple);
5491 display_par(mpl, code->arg.par.par, &memb);
5492 delete_symbol(mpl, memb.value.sym);
5493 break;
5494 case O_MEMSET:
5495 memb.value.set = eval_member_set(mpl, code->arg.set.set,
5496 memb.tuple);
5497 display_set(mpl, code->arg.set.set, &memb);
5498 break;
5499 case O_MEMVAR:
5500 memb.value.var = eval_member_var(mpl, code->arg.var.var,
5501 memb.tuple);
5502 display_var
5503 (mpl, code->arg.var.var, &memb, code->arg.var.suff);
5504 break;
5505 case O_MEMCON:
5506 memb.value.con = eval_member_con(mpl, code->arg.con.con,
5507 memb.tuple);
5508 display_con
5509 (mpl, code->arg.con.con, &memb, code->arg.con.suff);
5510 break;
5511 default:
5512 xassert(code != code);
5513 }
5514 delete_tuple(mpl, memb.tuple);
5515 return;
5516 }
5518 static void display_code(MPL *mpl, CODE *code)
5519 { /* display value of expression */
5520 switch (code->type)
5521 { case A_NUMERIC:
5522 /* numeric value */
5523 { double num;
5524 num = eval_numeric(mpl, code);
5525 write_text(mpl, "%.*g\n", DBL_DIG, num);
5526 }
5527 break;
5528 case A_SYMBOLIC:
5529 /* symbolic value */
5530 { SYMBOL *sym;
5531 sym = eval_symbolic(mpl, code);
5532 write_text(mpl, "%s\n", format_symbol(mpl, sym));
5533 delete_symbol(mpl, sym);
5534 }
5535 break;
5536 case A_LOGICAL:
5537 /* logical value */
5538 { int bit;
5539 bit = eval_logical(mpl, code);
5540 write_text(mpl, "%s\n", bit ? "true" : "false");
5541 }
5542 break;
5543 case A_TUPLE:
5544 /* n-tuple */
5545 { TUPLE *tuple;
5546 tuple = eval_tuple(mpl, code);
5547 write_text(mpl, "%s\n", format_tuple(mpl, '(', tuple));
5548 delete_tuple(mpl, tuple);
5549 }
5550 break;
5551 case A_ELEMSET:
5552 /* elemental set */
5553 { ELEMSET *set;
5554 MEMBER *memb;
5555 set = eval_elemset(mpl, code);
5556 if (set->head == 0)
5557 write_text(mpl, "set is empty\n");
5558 for (memb = set->head; memb != NULL; memb = memb->next)
5559 write_text(mpl, " %s\n", format_tuple(mpl, '(',
5560 memb->tuple));
5561 delete_elemset(mpl, set);
5562 }
5563 break;
5564 case A_FORMULA:
5565 /* linear form */
5566 { FORMULA *form, *term;
5567 form = eval_formula(mpl, code);
5568 if (form == NULL)
5569 write_text(mpl, "linear form is empty\n");
5570 for (term = form; term != NULL; term = term->next)
5571 { if (term->var == NULL)
5572 write_text(mpl, " %.*g\n", term->coef);
5573 else
5574 write_text(mpl, " %.*g %s%s\n", DBL_DIG,
5575 term->coef, term->var->var->name,
5576 format_tuple(mpl, '[', term->var->memb->tuple));
5577 }
5578 delete_formula(mpl, form);
5579 }
5580 break;
5581 default:
5582 xassert(code != code);
5583 }
5584 return;
5585 }
5587 static int display_func(MPL *mpl, void *info)
5588 { /* this is auxiliary routine to work within domain scope */
5589 DISPLAY *dpy = (DISPLAY *)info;
5590 DISPLAY1 *entry;
5591 for (entry = dpy->list; entry != NULL; entry = entry->next)
5592 { if (entry->type == A_INDEX)
5593 { /* dummy index */
5594 DOMAIN_SLOT *slot = entry->u.slot;
5595 write_text(mpl, "%s = %s\n", slot->name,
5596 format_symbol(mpl, slot->value));
5597 }
5598 else if (entry->type == A_SET)
5599 { /* model set */
5600 SET *set = entry->u.set;
5601 MEMBER *memb;
5602 if (set->assign != NULL)
5603 { /* the set has assignment expression; evaluate all its
5604 members over entire domain */
5605 eval_whole_set(mpl, set);
5606 }
5607 else
5608 { /* the set has no assignment expression; refer to its
5609 any existing member ignoring resultant value to check
5610 the data provided the data section */
5611 #if 1 /* 12/XII-2008 */
5612 if (set->gadget != NULL && set->data == 0)
5613 { /* initialize the set with data from a plain set */
5614 saturate_set(mpl, set);
5615 }
5616 #endif
5617 if (set->array->head != NULL)
5618 eval_member_set(mpl, set, set->array->head->tuple);
5619 }
5620 /* display all members of the set array */
5621 if (set->array->head == NULL)
5622 write_text(mpl, "%s has empty content\n", set->name);
5623 for (memb = set->array->head; memb != NULL; memb =
5624 memb->next) display_set(mpl, set, memb);
5625 }
5626 else if (entry->type == A_PARAMETER)
5627 { /* model parameter */
5628 PARAMETER *par = entry->u.par;
5629 MEMBER *memb;
5630 if (par->assign != NULL)
5631 { /* the parameter has an assignment expression; evaluate
5632 all its member over entire domain */
5633 eval_whole_par(mpl, par);
5634 }
5635 else
5636 { /* the parameter has no assignment expression; refer to
5637 its any existing member ignoring resultant value to
5638 check the data provided in the data section */
5639 if (par->array->head != NULL)
5640 { if (par->type != A_SYMBOLIC)
5641 eval_member_num(mpl, par, par->array->head->tuple);
5642 else
5643 delete_symbol(mpl, eval_member_sym(mpl, par,
5644 par->array->head->tuple));
5645 }
5646 }
5647 /* display all members of the parameter array */
5648 if (par->array->head == NULL)
5649 write_text(mpl, "%s has empty content\n", par->name);
5650 for (memb = par->array->head; memb != NULL; memb =
5651 memb->next) display_par(mpl, par, memb);
5652 }
5653 else if (entry->type == A_VARIABLE)
5654 { /* model variable */
5655 VARIABLE *var = entry->u.var;
5656 MEMBER *memb;
5657 xassert(mpl->flag_p);
5658 /* display all members of the variable array */
5659 if (var->array->head == NULL)
5660 write_text(mpl, "%s has empty content\n", var->name);
5661 for (memb = var->array->head; memb != NULL; memb =
5662 memb->next) display_var(mpl, var, memb, DOT_NONE);
5663 }
5664 else if (entry->type == A_CONSTRAINT)
5665 { /* model constraint */
5666 CONSTRAINT *con = entry->u.con;
5667 MEMBER *memb;
5668 xassert(mpl->flag_p);
5669 /* display all members of the constraint array */
5670 if (con->array->head == NULL)
5671 write_text(mpl, "%s has empty content\n", con->name);
5672 for (memb = con->array->head; memb != NULL; memb =
5673 memb->next) display_con(mpl, con, memb, DOT_NONE);
5674 }
5675 else if (entry->type == A_EXPRESSION)
5676 { /* expression */
5677 CODE *code = entry->u.code;
5678 if (code->op == O_MEMNUM || code->op == O_MEMSYM ||
5679 code->op == O_MEMSET || code->op == O_MEMVAR ||
5680 code->op == O_MEMCON)
5681 display_memb(mpl, code);
5682 else
5683 display_code(mpl, code);
5684 }
5685 else
5686 xassert(entry != entry);
5687 }
5688 return 0;
5689 }
5691 void execute_display(MPL *mpl, DISPLAY *dpy)
5692 { loop_within_domain(mpl, dpy->domain, dpy, display_func);
5693 return;
5694 }
5696 /*----------------------------------------------------------------------
5697 -- clean_display - clean display statement.
5698 --
5699 -- This routine cleans specified display statement that assumes deleting
5700 -- all stuff dynamically allocated on generating/postsolving phase. */
5702 void clean_display(MPL *mpl, DISPLAY *dpy)
5703 { DISPLAY1 *d;
5704 #if 0 /* 15/V-2010 */
5705 ARG_LIST *e;
5706 #endif
5707 /* clean subscript domain */
5708 clean_domain(mpl, dpy->domain);
5709 /* clean display list */
5710 for (d = dpy->list; d != NULL; d = d->next)
5711 { /* clean pseudo-code for computing expression */
5712 if (d->type == A_EXPRESSION)
5713 clean_code(mpl, d->u.code);
5714 #if 0 /* 15/V-2010 */
5715 /* clean pseudo-code for computing subscripts */
5716 for (e = d->list; e != NULL; e = e->next)
5717 clean_code(mpl, e->x);
5718 #endif
5719 }
5720 return;
5721 }
5723 /*----------------------------------------------------------------------
5724 -- execute_printf - execute printf statement.
5725 --
5726 -- This routine executes specified printf statement. */
5728 #if 1 /* 14/VII-2006 */
5729 static void print_char(MPL *mpl, int c)
5730 { if (mpl->prt_fp == NULL)
5731 write_char(mpl, c);
5732 else
5733 xfputc(c, mpl->prt_fp);
5734 return;
5735 }
5737 static void print_text(MPL *mpl, char *fmt, ...)
5738 { va_list arg;
5739 char buf[OUTBUF_SIZE], *c;
5740 va_start(arg, fmt);
5741 vsprintf(buf, fmt, arg);
5742 xassert(strlen(buf) < sizeof(buf));
5743 va_end(arg);
5744 for (c = buf; *c != '\0'; c++) print_char(mpl, *c);
5745 return;
5746 }
5747 #endif
5749 static int printf_func(MPL *mpl, void *info)
5750 { /* this is auxiliary routine to work within domain scope */
5751 PRINTF *prt = (PRINTF *)info;
5752 PRINTF1 *entry;
5753 SYMBOL *sym;
5754 char fmt[MAX_LENGTH+1], *c, *from, save;
5755 /* evaluate format control string */
5756 sym = eval_symbolic(mpl, prt->fmt);
5757 if (sym->str == NULL)
5758 sprintf(fmt, "%.*g", DBL_DIG, sym->num);
5759 else
5760 fetch_string(mpl, sym->str, fmt);
5761 delete_symbol(mpl, sym);
5762 /* scan format control string and perform formatting output */
5763 entry = prt->list;
5764 for (c = fmt; *c != '\0'; c++)
5765 { if (*c == '%')
5766 { /* scan format specifier */
5767 from = c++;
5768 if (*c == '%')
5769 { print_char(mpl, '%');
5770 continue;
5771 }
5772 if (entry == NULL) break;
5773 /* scan optional flags */
5774 while (*c == '-' || *c == '+' || *c == ' ' || *c == '#' ||
5775 *c == '0') c++;
5776 /* scan optional minimum field width */
5777 while (isdigit((unsigned char)*c)) c++;
5778 /* scan optional precision */
5779 if (*c == '.')
5780 { c++;
5781 while (isdigit((unsigned char)*c)) c++;
5782 }
5783 /* scan conversion specifier and perform formatting */
5784 save = *(c+1), *(c+1) = '\0';
5785 if (*c == 'd' || *c == 'i' || *c == 'e' || *c == 'E' ||
5786 *c == 'f' || *c == 'F' || *c == 'g' || *c == 'G')
5787 { /* the specifier requires numeric value */
5788 double value;
5789 xassert(entry != NULL);
5790 switch (entry->code->type)
5791 { case A_NUMERIC:
5792 value = eval_numeric(mpl, entry->code);
5793 break;
5794 case A_SYMBOLIC:
5795 sym = eval_symbolic(mpl, entry->code);
5796 if (sym->str != NULL)
5797 error(mpl, "cannot convert %s to floating-point"
5798 " number", format_symbol(mpl, sym));
5799 value = sym->num;
5800 delete_symbol(mpl, sym);
5801 break;
5802 case A_LOGICAL:
5803 if (eval_logical(mpl, entry->code))
5804 value = 1.0;
5805 else
5806 value = 0.0;
5807 break;
5808 default:
5809 xassert(entry != entry);
5810 }
5811 if (*c == 'd' || *c == 'i')
5812 { double int_max = (double)INT_MAX;
5813 if (!(-int_max <= value && value <= +int_max))
5814 error(mpl, "cannot convert %.*g to integer",
5815 DBL_DIG, value);
5816 print_text(mpl, from, (int)floor(value + 0.5));
5817 }
5818 else
5819 print_text(mpl, from, value);
5820 }
5821 else if (*c == 's')
5822 { /* the specifier requires symbolic value */
5823 char value[MAX_LENGTH+1];
5824 switch (entry->code->type)
5825 { case A_NUMERIC:
5826 sprintf(value, "%.*g", DBL_DIG, eval_numeric(mpl,
5827 entry->code));
5828 break;
5829 case A_LOGICAL:
5830 if (eval_logical(mpl, entry->code))
5831 strcpy(value, "T");
5832 else
5833 strcpy(value, "F");
5834 break;
5835 case A_SYMBOLIC:
5836 sym = eval_symbolic(mpl, entry->code);
5837 if (sym->str == NULL)
5838 sprintf(value, "%.*g", DBL_DIG, sym->num);
5839 else
5840 fetch_string(mpl, sym->str, value);
5841 delete_symbol(mpl, sym);
5842 break;
5843 default:
5844 xassert(entry != entry);
5845 }
5846 print_text(mpl, from, value);
5847 }
5848 else
5849 error(mpl, "format specifier missing or invalid");
5850 *(c+1) = save;
5851 entry = entry->next;
5852 }
5853 else if (*c == '\\')
5854 { /* write some control character */
5855 c++;
5856 if (*c == 't')
5857 print_char(mpl, '\t');
5858 else if (*c == 'n')
5859 print_char(mpl, '\n');
5860 #if 1 /* 28/X-2010 */
5861 else if (*c == '\0')
5862 { /* format string ends with backslash */
5863 error(mpl, "invalid use of escape character \\ in format"
5864 " control string");
5865 }
5866 #endif
5867 else
5868 print_char(mpl, *c);
5869 }
5870 else
5871 { /* write character without formatting */
5872 print_char(mpl, *c);
5873 }
5874 }
5875 return 0;
5876 }
5878 #if 0 /* 14/VII-2006 */
5879 void execute_printf(MPL *mpl, PRINTF *prt)
5880 { loop_within_domain(mpl, prt->domain, prt, printf_func);
5881 return;
5882 }
5883 #else
5884 void execute_printf(MPL *mpl, PRINTF *prt)
5885 { if (prt->fname == NULL)
5886 { /* switch to the standard output */
5887 if (mpl->prt_fp != NULL)
5888 { xfclose(mpl->prt_fp), mpl->prt_fp = NULL;
5889 xfree(mpl->prt_file), mpl->prt_file = NULL;
5890 }
5891 }
5892 else
5893 { /* evaluate file name string */
5894 SYMBOL *sym;
5895 char fname[MAX_LENGTH+1];
5896 sym = eval_symbolic(mpl, prt->fname);
5897 if (sym->str == NULL)
5898 sprintf(fname, "%.*g", DBL_DIG, sym->num);
5899 else
5900 fetch_string(mpl, sym->str, fname);
5901 delete_symbol(mpl, sym);
5902 /* close the current print file, if necessary */
5903 if (mpl->prt_fp != NULL &&
5904 (!prt->app || strcmp(mpl->prt_file, fname) != 0))
5905 { xfclose(mpl->prt_fp), mpl->prt_fp = NULL;
5906 xfree(mpl->prt_file), mpl->prt_file = NULL;
5907 }
5908 /* open the specified print file, if necessary */
5909 if (mpl->prt_fp == NULL)
5910 { mpl->prt_fp = xfopen(fname, prt->app ? "a" : "w");
5911 if (mpl->prt_fp == NULL)
5912 error(mpl, "unable to open `%s' for writing - %s",
5913 fname, xerrmsg());
5914 mpl->prt_file = xmalloc(strlen(fname)+1);
5915 strcpy(mpl->prt_file, fname);
5916 }
5917 }
5918 loop_within_domain(mpl, prt->domain, prt, printf_func);
5919 if (mpl->prt_fp != NULL)
5920 { xfflush(mpl->prt_fp);
5921 if (xferror(mpl->prt_fp))
5922 error(mpl, "writing error to `%s' - %s", mpl->prt_file,
5923 xerrmsg());
5924 }
5925 return;
5926 }
5927 #endif
5929 /*----------------------------------------------------------------------
5930 -- clean_printf - clean printf statement.
5931 --
5932 -- This routine cleans specified printf statement that assumes deleting
5933 -- all stuff dynamically allocated on generating/postsolving phase. */
5935 void clean_printf(MPL *mpl, PRINTF *prt)
5936 { PRINTF1 *p;
5937 /* clean subscript domain */
5938 clean_domain(mpl, prt->domain);
5939 /* clean pseudo-code for computing format string */
5940 clean_code(mpl, prt->fmt);
5941 /* clean printf list */
5942 for (p = prt->list; p != NULL; p = p->next)
5943 { /* clean pseudo-code for computing value to be printed */
5944 clean_code(mpl, p->code);
5945 }
5946 #if 1 /* 14/VII-2006 */
5947 /* clean pseudo-code for computing file name string */
5948 clean_code(mpl, prt->fname);
5949 #endif
5950 return;
5951 }
5953 /*----------------------------------------------------------------------
5954 -- execute_for - execute for statement.
5955 --
5956 -- This routine executes specified for statement. */
5958 static int for_func(MPL *mpl, void *info)
5959 { /* this is auxiliary routine to work within domain scope */
5960 FOR *fur = (FOR *)info;
5961 STATEMENT *stmt, *save;
5962 save = mpl->stmt;
5963 for (stmt = fur->list; stmt != NULL; stmt = stmt->next)
5964 execute_statement(mpl, stmt);
5965 mpl->stmt = save;
5966 return 0;
5967 }
5969 void execute_for(MPL *mpl, FOR *fur)
5970 { loop_within_domain(mpl, fur->domain, fur, for_func);
5971 return;
5972 }
5974 /*----------------------------------------------------------------------
5975 -- clean_for - clean for statement.
5976 --
5977 -- This routine cleans specified for statement that assumes deleting all
5978 -- stuff dynamically allocated on generating/postsolving phase. */
5980 void clean_for(MPL *mpl, FOR *fur)
5981 { STATEMENT *stmt;
5982 /* clean subscript domain */
5983 clean_domain(mpl, fur->domain);
5984 /* clean all sub-statements */
5985 for (stmt = fur->list; stmt != NULL; stmt = stmt->next)
5986 clean_statement(mpl, stmt);
5987 return;
5988 }
5990 /*----------------------------------------------------------------------
5991 -- execute_statement - execute specified model statement.
5992 --
5993 -- This routine executes specified model statement. */
5995 void execute_statement(MPL *mpl, STATEMENT *stmt)
5996 { mpl->stmt = stmt;
5997 switch (stmt->type)
5998 { case A_SET:
5999 case A_PARAMETER:
6000 case A_VARIABLE:
6001 break;
6002 case A_CONSTRAINT:
6003 xprintf("Generating %s...\n", stmt->u.con->name);
6004 eval_whole_con(mpl, stmt->u.con);
6005 break;
6006 case A_TABLE:
6007 switch (stmt->u.tab->type)
6008 { case A_INPUT:
6009 xprintf("Reading %s...\n", stmt->u.tab->name);
6010 break;
6011 case A_OUTPUT:
6012 xprintf("Writing %s...\n", stmt->u.tab->name);
6013 break;
6014 default:
6015 xassert(stmt != stmt);
6016 }
6017 execute_table(mpl, stmt->u.tab);
6018 break;
6019 case A_SOLVE:
6020 break;
6021 case A_CHECK:
6022 xprintf("Checking (line %d)...\n", stmt->line);
6023 execute_check(mpl, stmt->u.chk);
6024 break;
6025 case A_DISPLAY:
6026 write_text(mpl, "Display statement at line %d\n",
6027 stmt->line);
6028 execute_display(mpl, stmt->u.dpy);
6029 break;
6030 case A_PRINTF:
6031 execute_printf(mpl, stmt->u.prt);
6032 break;
6033 case A_FOR:
6034 execute_for(mpl, stmt->u.fur);
6035 break;
6036 default:
6037 xassert(stmt != stmt);
6038 }
6039 return;
6040 }
6042 /*----------------------------------------------------------------------
6043 -- clean_statement - clean specified model statement.
6044 --
6045 -- This routine cleans specified model statement that assumes deleting
6046 -- all stuff dynamically allocated on generating/postsolving phase. */
6048 void clean_statement(MPL *mpl, STATEMENT *stmt)
6049 { switch(stmt->type)
6050 { case A_SET:
6051 clean_set(mpl, stmt->u.set); break;
6052 case A_PARAMETER:
6053 clean_parameter(mpl, stmt->u.par); break;
6054 case A_VARIABLE:
6055 clean_variable(mpl, stmt->u.var); break;
6056 case A_CONSTRAINT:
6057 clean_constraint(mpl, stmt->u.con); break;
6058 #if 1 /* 11/II-2008 */
6059 case A_TABLE:
6060 clean_table(mpl, stmt->u.tab); break;
6061 #endif
6062 case A_SOLVE:
6063 break;
6064 case A_CHECK:
6065 clean_check(mpl, stmt->u.chk); break;
6066 case A_DISPLAY:
6067 clean_display(mpl, stmt->u.dpy); break;
6068 case A_PRINTF:
6069 clean_printf(mpl, stmt->u.prt); break;
6070 case A_FOR:
6071 clean_for(mpl, stmt->u.fur); break;
6072 default:
6073 xassert(stmt != stmt);
6074 }
6075 return;
6076 }
6078 /* eof */