COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpmpl03.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 210.8 KB
Line 
1/* glpmpl03.c */
2
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 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***********************************************************************/
24
25#define _GLPSTD_ERRNO
26#define _GLPSTD_STDIO
27#include "glpenv.h"
28#include "glpmpl.h"
29
30/**********************************************************************/
31/* * *                   FLOATING-POINT NUMBERS                   * * */
32/**********************************************************************/
33
34/*----------------------------------------------------------------------
35-- fp_add - floating-point addition.
36--
37-- This routine computes the sum x + y. */
38
39double 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}
46
47/*----------------------------------------------------------------------
48-- fp_sub - floating-point subtraction.
49--
50-- This routine computes the difference x - y. */
51
52double 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}
59
60/*----------------------------------------------------------------------
61-- fp_less - floating-point non-negative subtraction.
62--
63-- This routine computes the non-negative difference max(0, x - y). */
64
65double 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}
72
73/*----------------------------------------------------------------------
74-- fp_mul - floating-point multiplication.
75--
76-- This routine computes the product x * y. */
77
78double 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}
84
85/*----------------------------------------------------------------------
86-- fp_div - floating-point division.
87--
88-- This routine computes the quotient x / y. */
89
90double 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}
99
100/*----------------------------------------------------------------------
101-- fp_idiv - floating-point quotient of exact division.
102--
103-- This routine computes the quotient of exact division x div y. */
104
105double 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}
115
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). */
122
123double 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}
139
140/*----------------------------------------------------------------------
141-- fp_power - floating-point exponentiation (raise to power).
142--
143-- This routine computes the exponentiation x ** y. */
144
145double 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
163eval:    r = pow(x, y);
164      return r;
165}
166
167/*----------------------------------------------------------------------
168-- fp_exp - floating-point base-e exponential.
169--
170-- This routine computes the base-e exponential e ** x. */
171
172double 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}
177
178/*----------------------------------------------------------------------
179-- fp_log - floating-point natural logarithm.
180--
181-- This routine computes the natural logarithm log x. */
182
183double 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}
188
189/*----------------------------------------------------------------------
190-- fp_log10 - floating-point common (decimal) logarithm.
191--
192-- This routine computes the common (decimal) logarithm lg x. */
193
194double 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}
199
200/*----------------------------------------------------------------------
201-- fp_sqrt - floating-point square root.
202--
203-- This routine computes the square root x ** 0.5. */
204
205double 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}
210
211/*----------------------------------------------------------------------
212-- fp_sin - floating-point trigonometric sine.
213--
214-- This routine computes the trigonometric sine sin(x). */
215
216double 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}
221
222/*----------------------------------------------------------------------
223-- fp_cos - floating-point trigonometric cosine.
224--
225-- This routine computes the trigonometric cosine cos(x). */
226
227double 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}
232
233/*----------------------------------------------------------------------
234-- fp_atan - floating-point trigonometric arctangent.
235--
236-- This routine computes the trigonometric arctangent atan(x). */
237
238double fp_atan(MPL *mpl, double x)
239{     xassert(mpl == mpl);
240      return atan(x);
241}
242
243/*----------------------------------------------------------------------
244-- fp_atan2 - floating-point trigonometric arctangent.
245--
246-- This routine computes the trigonometric arctangent atan(y / x). */
247
248double fp_atan2(MPL *mpl, double y, double x)
249{     xassert(mpl == mpl);
250      return atan2(y, x);
251}
252
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. */
262
263double 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}
277
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. */
289
290double 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}
304
305/**********************************************************************/
306/* * *              PSEUDO-RANDOM NUMBER GENERATORS               * * */
307/**********************************************************************/
308
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. */
315
316#define two_to_the_24 0x1000000
317
318double fp_irand224(MPL *mpl)
319{     return
320         (double)rng_unif_rand(mpl->rand, two_to_the_24);
321}
322
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). */
328
329#define two_to_the_31 ((unsigned int)0x80000000)
330
331double fp_uniform01(MPL *mpl)
332{     return
333         (double)rng_next_rand(mpl->rand) / (double)two_to_the_31;
334}
335
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). */
341
342double 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}
355
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. */
364
365double 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}
377
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. */
383
384double 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}
393
394/**********************************************************************/
395/* * *                SEGMENTED CHARACTER STRINGS                 * * */
396/**********************************************************************/
397
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. */
403
404STRING *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)
417tail = (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
430
431/*----------------------------------------------------------------------
432-- copy_string - make copy of character string.
433--
434-- This routine returns an exact copy of segmented character string. */
435
436STRING *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)
447tail = (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
457
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. */
467
468int 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      }
487done: return 0;
488}
489#else
490{     xassert(mpl == mpl);
491      return strcmp(str1, str2);
492}
493#endif
494
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. */
500
501char *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      }
515done: xassert(strlen(buf) <= MAX_LENGTH);
516      return buf;
517}
518#else
519{     xassert(mpl == mpl);
520      return strcpy(buf, str);
521}
522#endif
523
524/*----------------------------------------------------------------------
525-- delete_string - delete character string.
526--
527-- This routine deletes specified segmented character string. */
528
529void 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
548
549/**********************************************************************/
550/* * *                          SYMBOLS                           * * */
551/**********************************************************************/
552
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. */
558
559SYMBOL *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}
566
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. */
572
573SYMBOL *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}
584
585/*----------------------------------------------------------------------
586-- copy_symbol - make copy of symbol.
587--
588-- This routine returns an exact copy of symbol. */
589
590SYMBOL *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}
607
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. */
620
621int 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}
638
639/*----------------------------------------------------------------------
640-- delete_symbol - delete symbol.
641--
642-- This routine deletes specified symbol. */
643
644void 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}
653
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. */
662
663char *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}
703
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. */
710
711SYMBOL *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}
738
739/**********************************************************************/
740/* * *                          N-TUPLES                          * * */
741/**********************************************************************/
742
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. */
748
749TUPLE *create_tuple(MPL *mpl)
750{     TUPLE *tuple;
751      xassert(mpl == mpl);
752      tuple = NULL;
753      return tuple;
754}
755
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. */
761
762TUPLE *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}
782
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. */
788
789int 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}
799
800/*----------------------------------------------------------------------
801-- copy_tuple - make copy of n-tuple.
802--
803-- This routine returns an exact copy of n-tuple. */
804
805TUPLE *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)
818tail = (tail->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));
819         }
820         tail->next = NULL;
821      }
822      return head;
823}
824
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. */
838
839int 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}
858
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. */
864
865TUPLE *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}
879
880/*----------------------------------------------------------------------
881-- delete_tuple - delete n-tuple.
882--
883-- This routine deletes specified n-tuple. */
884
885void 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}
899
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. */
908
909char *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}
941
942/**********************************************************************/
943/* * *                       ELEMENTAL SETS                       * * */
944/**********************************************************************/
945
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. */
951
952ELEMSET *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}
958
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. */
966
967MEMBER *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}
977
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. */
988
989MEMBER *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}
1002
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. */
1008
1009MEMBER *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}
1019
1020/*----------------------------------------------------------------------
1021-- copy_elemset - make copy of elemental set.
1022--
1023-- This routine makes an exact copy of elemental set. */
1024
1025ELEMSET *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}
1039
1040/*----------------------------------------------------------------------
1041-- delete_elemset - delete elemental set.
1042--
1043-- This routine deletes specified elemental set. */
1044
1045void 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}
1054
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). */
1066
1067int 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}
1094
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. */
1110
1111double 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}
1115
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. */
1125
1126ELEMSET *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}
1147
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). */
1156
1157ELEMSET *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}
1177
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). */
1186
1187ELEMSET *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}
1210
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). */
1219
1220ELEMSET *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}
1249
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). */
1258
1259ELEMSET *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}
1282
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). */
1291
1292ELEMSET *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}
1320
1321/**********************************************************************/
1322/* * *                    ELEMENTAL VARIABLES                     * * */
1323/**********************************************************************/
1324
1325/* (there are no specific routines for elemental variables) */
1326
1327/**********************************************************************/
1328/* * *                        LINEAR FORMS                        * * */
1329/**********************************************************************/
1330
1331/*----------------------------------------------------------------------
1332-- constant_term - create constant term.
1333--
1334-- This routine creates the linear form, which is a constant term. */
1335
1336FORMULA *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}
1348
1349/*----------------------------------------------------------------------
1350-- single_variable - create single variable.
1351--
1352-- This routine creates the linear form, which is a single elemental
1353-- variable. */
1354
1355FORMULA *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}
1367
1368/*----------------------------------------------------------------------
1369-- copy_formula - make copy of linear form.
1370--
1371-- This routine returns an exact copy of linear form. */
1372
1373FORMULA *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)
1386tail = (tail->next = dmp_get_atom(mpl->formulae, sizeof(FORMULA)));
1387         }
1388         tail->next = NULL;
1389      }
1390      return head;
1391}
1392
1393/*----------------------------------------------------------------------
1394-- delete_formula - delete linear form.
1395--
1396-- This routine deletes specified linear form. */
1397
1398void 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}
1410
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). */
1420
1421FORMULA *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}
1467
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. */
1473
1474FORMULA *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}
1497
1498/*----------------------------------------------------------------------
1499-- reduce_terms - reduce identical terms in linear form.
1500--
1501-- This routine reduces identical terms in specified linear form. */
1502
1503FORMULA *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}
1531
1532/**********************************************************************/
1533/* * *                   ELEMENTAL CONSTRAINTS                    * * */
1534/**********************************************************************/
1535
1536/* (there are no specific routines for elemental constraints) */
1537
1538/**********************************************************************/
1539/* * *                       GENERIC VALUES                       * * */
1540/**********************************************************************/
1541
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. */
1548
1549void 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}
1588
1589/**********************************************************************/
1590/* * *                SYMBOLICALLY INDEXED ARRAYS                 * * */
1591/**********************************************************************/
1592
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). */
1611
1612ARRAY *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}
1632
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. */
1640
1641static 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}
1647
1648MEMBER *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)
1662avl_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);
1675memb = (MEMBER *)(node == NULL ? NULL : avl_get_node_link(node));
1676      }
1677      return memb;
1678}
1679
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. */
1694
1695MEMBER *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)
1718avl_set_node_link(avl_insert_node(array->tree, memb->tuple),
1719            (void *)memb);
1720      return memb;
1721}
1722
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. */
1731
1732void 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}
1760
1761/**********************************************************************/
1762/* * *                 DOMAINS AND DUMMY INDICES                  * * */
1763/**********************************************************************/
1764
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. */
1771
1772void 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);
1807done: return;
1808}
1809
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. */
1816
1817void 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}
1833
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. */
1865
1866int 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);
1904done: return ret;
1905}
1906
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. */
1933
1934struct 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};
1951
1952static 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
1977temp = (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}
2028
2029int 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}
2055
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. */
2073
2074struct 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};
2087
2088static 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);
2181skip:          ;
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}
2207
2208void 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}
2228
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. */
2234
2235void 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}
2246
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. */
2257
2258TUPLE *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}
2279
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. */
2285
2286void 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);
2306done: return;
2307}
2308
2309/**********************************************************************/
2310/* * *                         MODEL SETS                         * * */
2311/**********************************************************************/
2312
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. */
2320
2321void 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}
2347
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. */
2355
2356ELEMSET *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);
2372add:     /* 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}
2390
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. */
2396
2397struct 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};
2412
2413static 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}
2427
2428#if 1 /* 12/XII-2008 */
2429static 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
2488
2489ELEMSET *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}
2537
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. */
2543
2544static 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}
2552
2553void eval_whole_set(MPL *mpl, SET *set)
2554{     loop_within_domain(mpl, set->domain, set, whole_set_func);
2555      return;
2556}
2557
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. */
2563
2564void 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}
2584
2585/**********************************************************************/
2586/* * *                      MODEL PARAMETERS                      * * */
2587/**********************************************************************/
2588
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. */
2596
2597void 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 = "<";
2634err:              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}
2673
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. */
2681
2682double 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);
2698add:     /* 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}
2724
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. */
2730
2731struct 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};
2746
2747static 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}
2761
2762double 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}
2806
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. */
2814
2815void 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}
2911
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. */
2919
2920SYMBOL *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);
2936add:     /* 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}
2959
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. */
2965
2966struct 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};
2981
2982static 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}
2996
2997SYMBOL *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}
3040
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. */
3046
3047static 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}
3066
3067void eval_whole_par(MPL *mpl, PARAMETER *par)
3068{     loop_within_domain(mpl, par->domain, par, whole_par_func);
3069      return;
3070}
3071
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. */
3077
3078void 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}
3105
3106/**********************************************************************/
3107/* * *                      MODEL VARIABLES                       * * */
3108/**********************************************************************/
3109
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. */
3118
3119ELEMVAR *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}
3164
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. */
3170
3171struct 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};
3180
3181static 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}
3187
3188ELEMVAR *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}
3205
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. */
3211
3212static 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}
3220
3221void eval_whole_var(MPL *mpl, VARIABLE *var)
3222{     loop_within_domain(mpl, var->domain, var, whole_var_func);
3223      return;
3224}
3225
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. */
3231
3232void 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}
3246
3247/**********************************************************************/
3248/* * *              MODEL CONSTRAINTS AND OBJECTIVES              * * */
3249/**********************************************************************/
3250
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. */
3259
3260ELEMCON *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}
3351
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. */
3357
3358struct 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};
3367
3368static 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}
3374
3375ELEMCON *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}
3392
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. */
3398
3399static 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}
3407
3408void eval_whole_con(MPL *mpl, CONSTRAINT *con)
3409{     loop_within_domain(mpl, con->domain, con, whole_con_func);
3410      return;
3411}
3412
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. */
3418
3419void 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}
3437
3438/**********************************************************************/
3439/* * *                        PSEUDO-CODE                         * * */
3440/**********************************************************************/
3441
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. */
3447
3448struct 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};
3455
3456static 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}
3484
3485double 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;
3896done: return value;
3897}
3898
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. */
3904
3905SYMBOL *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);
4022done: return value;
4023}
4024
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. */
4030
4031struct 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};
4038
4039static 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}
4060
4061int 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;
4273done: return value;
4274}
4275
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. */
4281
4282TUPLE *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);
4321done: return value;
4322}
4323
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. */
4329
4330struct 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};
4337
4338static 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}
4365
4366ELEMSET *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);
4477done: return value;
4478}
4479
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. */
4488
4489static 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}
4495
4496int 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}
4633
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. */
4639
4640struct 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};
4649
4650static 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}
4690
4691FORMULA *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);
4808done: return value;
4809}
4810
4811/*----------------------------------------------------------------------
4812-- clean_code - clean pseudo-code.
4813--
4814-- This routine recursively cleans specified pseudo-code that assumes
4815-- deleting all temporary resultant values. */
4816
4817void 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      }
4955done: return;
4956}
4957
4958#if 1 /* 11/II-2008 */
4959/**********************************************************************/
4960/* * *                        DATA TABLES                         * * */
4961/**********************************************************************/
4962
4963int mpl_tab_num_args(TABDCA *dca)
4964{     /* returns the number of arguments */
4965      return dca->na;
4966}
4967
4968const 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}
4973
4974int mpl_tab_num_flds(TABDCA *dca)
4975{     /* returns the number of fields */
4976      return dca->nf;
4977}
4978
4979const 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}
4984
4985int 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}
4990
4991double 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}
4997
4998const 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}
5005
5006void 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}
5014
5015void 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}
5025
5026static 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}
5067
5068void 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      }
5119read_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;
5252write_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);
5277done: /* free table driver communication area */
5278      free_dca(mpl);
5279      return;
5280}
5281
5282void 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}
5309
5310void 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
5333
5334/**********************************************************************/
5335/* * *                      MODEL STATEMENTS                      * * */
5336/**********************************************************************/
5337
5338/*----------------------------------------------------------------------
5339-- execute_check - execute check statement.
5340--
5341-- This routine executes specified check statement. */
5342
5343static 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}
5351
5352void execute_check(MPL *mpl, CHECK *chk)
5353{     loop_within_domain(mpl, chk->domain, chk, check_func);
5354      return;
5355}
5356
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. */
5362
5363void 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}
5370
5371/*----------------------------------------------------------------------
5372-- execute_display - execute display statement.
5373--
5374-- This routine executes specified display statement. */
5375
5376static 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}
5387
5388static 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}
5408
5409#if 1 /* 15/V-2010 */
5410static 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
5439
5440#if 1 /* 15/V-2010 */
5441static 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
5470
5471static 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}
5517
5518static 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}
5586
5587static 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}
5690
5691void execute_display(MPL *mpl, DISPLAY *dpy)
5692{     loop_within_domain(mpl, dpy->domain, dpy, display_func);
5693      return;
5694}
5695
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. */
5701
5702void 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}
5722
5723/*----------------------------------------------------------------------
5724-- execute_printf - execute printf statement.
5725--
5726-- This routine executes specified printf statement. */
5727
5728#if 1 /* 14/VII-2006 */
5729static 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}
5736
5737static 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
5748
5749static 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}
5877
5878#if 0 /* 14/VII-2006 */
5879void execute_printf(MPL *mpl, PRINTF *prt)
5880{     loop_within_domain(mpl, prt->domain, prt, printf_func);
5881      return;
5882}
5883#else
5884void 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
5928
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. */
5934
5935void 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}
5952
5953/*----------------------------------------------------------------------
5954-- execute_for - execute for statement.
5955--
5956-- This routine executes specified for statement. */
5957
5958static 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}
5968
5969void execute_for(MPL *mpl, FOR *fur)
5970{     loop_within_domain(mpl, fur->domain, fur, for_func);
5971      return;
5972}
5973
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. */
5979
5980void 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}
5989
5990/*----------------------------------------------------------------------
5991-- execute_statement - execute specified model statement.
5992--
5993-- This routine executes specified model statement. */
5994
5995void 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}
6041
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. */
6047
6048void 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}
6077
6078/* eof */
Note: See TracBrowser for help on using the repository browser.