glpk-cmake: src/glpmpl03.c@4c8956a7bdf4

     1 /* glpmpl03.c */

     3 /***********************************************************************

     4 *  This code is part of GLPK (GNU Linear Programming Kit).

5 *

     6 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

     7 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

     8 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

     9 *  E-mail: <mao@gnu.org>.

    10 *

    11 *  GLPK is free software: you can redistribute it and/or modify it

    12 *  under the terms of the GNU General Public License as published by

    13 *  the Free Software Foundation, either version 3 of the License, or

    14 *  (at your option) any later version.

    15 *

    16 *  GLPK is distributed in the hope that it will be useful, but WITHOUT

    17 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

    18 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

    19 *  License for more details.

    20 *

    21 *  You should have received a copy of the GNU General Public License

    22 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

    23 ***********************************************************************/

    25 #define _GLPSTD_ERRNO

    26 #define _GLPSTD_STDIO

    27 #include "glpenv.h"

    28 #include "glpmpl.h"

    30 /**********************************************************************/

    31 /* * *                   FLOATING-POINT NUMBERS                   * * */

    32 /**********************************************************************/

    34 /*----------------------------------------------------------------------

    35 -- fp_add - floating-point addition.

    36 --

    37 -- This routine computes the sum x + y. */

    39 double fp_add(MPL *mpl, double x, double y)

    40 {     if (x > 0.0 && y > 0.0 && x > + 0.999 * DBL_MAX - y ||

    41           x < 0.0 && y < 0.0 && x < - 0.999 * DBL_MAX - y)

    42          error(mpl, "%.*g + %.*g; floating-point overflow",

    43             DBL_DIG, x, DBL_DIG, y);

    44       return x + y;

    45 }

    47 /*----------------------------------------------------------------------

    48 -- fp_sub - floating-point subtraction.

    49 --

    50 -- This routine computes the difference x - y. */

    52 double fp_sub(MPL *mpl, double x, double y)

    53 {     if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y ||

    54           x < 0.0 && y > 0.0 && x < - 0.999 * DBL_MAX + y)

    55          error(mpl, "%.*g - %.*g; floating-point overflow",

    56             DBL_DIG, x, DBL_DIG, y);

    57       return x - y;

    58 }

    60 /*----------------------------------------------------------------------

    61 -- fp_less - floating-point non-negative subtraction.

    62 --

    63 -- This routine computes the non-negative difference max(0, x - y). */

    65 double fp_less(MPL *mpl, double x, double y)

    66 {     if (x < y) return 0.0;

    67       if (x > 0.0 && y < 0.0 && x > + 0.999 * DBL_MAX + y)

    68          error(mpl, "%.*g less %.*g; floating-point overflow",

    69             DBL_DIG, x, DBL_DIG, y);

    70       return x - y;

    71 }

    73 /*----------------------------------------------------------------------

    74 -- fp_mul - floating-point multiplication.

    75 --

    76 -- This routine computes the product x * y. */

    78 double fp_mul(MPL *mpl, double x, double y)

    79 {     if (fabs(y) > 1.0 && fabs(x) > (0.999 * DBL_MAX) / fabs(y))

    80          error(mpl, "%.*g * %.*g; floating-point overflow",

    81             DBL_DIG, x, DBL_DIG, y);

    82       return x * y;

    83 }

    85 /*----------------------------------------------------------------------

    86 -- fp_div - floating-point division.

    87 --

    88 -- This routine computes the quotient x / y. */

    90 double fp_div(MPL *mpl, double x, double y)

    91 {     if (fabs(y) < DBL_MIN)

    92          error(mpl, "%.*g / %.*g; floating-point zero divide",

    93             DBL_DIG, x, DBL_DIG, y);

    94       if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))

    95          error(mpl, "%.*g / %.*g; floating-point overflow",

    96             DBL_DIG, x, DBL_DIG, y);

    97       return x / y;

    98 }

   100 /*----------------------------------------------------------------------

   101 -- fp_idiv - floating-point quotient of exact division.

   102 --

   103 -- This routine computes the quotient of exact division x div y. */

   105 double fp_idiv(MPL *mpl, double x, double y)

   106 {     if (fabs(y) < DBL_MIN)

   107          error(mpl, "%.*g div %.*g; floating-point zero divide",

   108             DBL_DIG, x, DBL_DIG, y);

   109       if (fabs(y) < 1.0 && fabs(x) > (0.999 * DBL_MAX) * fabs(y))

   110          error(mpl, "%.*g div %.*g; floating-point overflow",

   111             DBL_DIG, x, DBL_DIG, y);

   112       x /= y;

   113       return x > 0.0 ? floor(x) : x < 0.0 ? ceil(x) : 0.0;

   114 }

   116 /*----------------------------------------------------------------------

   117 -- fp_mod - floating-point remainder of exact division.

   118 --

   119 -- This routine computes the remainder of exact division x mod y.

   120 --

   121 -- NOTE: By definition x mod y = x - y * floor(x / y). */

   123 double fp_mod(MPL *mpl, double x, double y)

   124 {     double r;

   125       xassert(mpl == mpl);

   126       if (x == 0.0)

   127          r = 0.0;

   128       else if (y == 0.0)

   129          r = x;

   130       else

   131       {  r = fmod(fabs(x), fabs(y));

   132          if (r != 0.0)

   133          {  if (x < 0.0) r = - r;

   134             if (x > 0.0 && y < 0.0 || x < 0.0 && y > 0.0) r += y;

   135          }

   136       }

   137       return r;

   138 }

   140 /*----------------------------------------------------------------------

   141 -- fp_power - floating-point exponentiation (raise to power).

   142 --

   143 -- This routine computes the exponentiation x ** y. */

   145 double fp_power(MPL *mpl, double x, double y)

   146 {     double r;

   147       if (x == 0.0 && y <= 0.0 || x < 0.0 && y != floor(y))

   148          error(mpl, "%.*g ** %.*g; result undefined",

   149             DBL_DIG, x, DBL_DIG, y);

   150       if (x == 0.0) goto eval;

   151       if (fabs(x) > 1.0 && y > +1.0 &&

   152             +log(fabs(x)) > (0.999 * log(DBL_MAX)) / y ||

   153           fabs(x) < 1.0 && y < -1.0 &&

   154             +log(fabs(x)) < (0.999 * log(DBL_MAX)) / y)

   155          error(mpl, "%.*g ** %.*g; floating-point overflow",

   156             DBL_DIG, x, DBL_DIG, y);

   157       if (fabs(x) > 1.0 && y < -1.0 &&

   158             -log(fabs(x)) < (0.999 * log(DBL_MAX)) / y ||

   159           fabs(x) < 1.0 && y > +1.0 &&

   160             -log(fabs(x)) > (0.999 * log(DBL_MAX)) / y)

   161          r = 0.0;

   162       else

   163 eval:    r = pow(x, y);

   164       return r;

   165 }

   167 /*----------------------------------------------------------------------

   168 -- fp_exp - floating-point base-e exponential.

   169 --

   170 -- This routine computes the base-e exponential e ** x. */

   172 double fp_exp(MPL *mpl, double x)

   173 {     if (x > 0.999 * log(DBL_MAX))

   174          error(mpl, "exp(%.*g); floating-point overflow", DBL_DIG, x);

   175       return exp(x);

   176 }

   178 /*----------------------------------------------------------------------

   179 -- fp_log - floating-point natural logarithm.

   180 --

   181 -- This routine computes the natural logarithm log x. */

   183 double fp_log(MPL *mpl, double x)

   184 {     if (x <= 0.0)

   185          error(mpl, "log(%.*g); non-positive argument", DBL_DIG, x);

   186       return log(x);

   187 }

   189 /*----------------------------------------------------------------------

   190 -- fp_log10 - floating-point common (decimal) logarithm.

   191 --

   192 -- This routine computes the common (decimal) logarithm lg x. */

   194 double fp_log10(MPL *mpl, double x)

   195 {     if (x <= 0.0)

   196          error(mpl, "log10(%.*g); non-positive argument", DBL_DIG, x);

   197       return log10(x);

   198 }

   200 /*----------------------------------------------------------------------

   201 -- fp_sqrt - floating-point square root.

   202 --

   203 -- This routine computes the square root x ** 0.5. */

   205 double fp_sqrt(MPL *mpl, double x)

   206 {     if (x < 0.0)

   207          error(mpl, "sqrt(%.*g); negative argument", DBL_DIG, x);

   208       return sqrt(x);

   209 }

   211 /*----------------------------------------------------------------------

   212 -- fp_sin - floating-point trigonometric sine.

   213 --

   214 -- This routine computes the trigonometric sine sin(x). */

   216 double fp_sin(MPL *mpl, double x)

   217 {     if (!(-1e6 <= x && x <= +1e6))

   218          error(mpl, "sin(%.*g); argument too large", DBL_DIG, x);

   219       return sin(x);

   220 }

   222 /*----------------------------------------------------------------------

   223 -- fp_cos - floating-point trigonometric cosine.

   224 --

   225 -- This routine computes the trigonometric cosine cos(x). */

   227 double fp_cos(MPL *mpl, double x)

   228 {     if (!(-1e6 <= x && x <= +1e6))

   229          error(mpl, "cos(%.*g); argument too large", DBL_DIG, x);

   230       return cos(x);

   231 }

   233 /*----------------------------------------------------------------------

   234 -- fp_atan - floating-point trigonometric arctangent.

   235 --

   236 -- This routine computes the trigonometric arctangent atan(x). */

   238 double fp_atan(MPL *mpl, double x)

   239 {     xassert(mpl == mpl);

   240       return atan(x);

   241 }

   243 /*----------------------------------------------------------------------

   244 -- fp_atan2 - floating-point trigonometric arctangent.

   245 --

   246 -- This routine computes the trigonometric arctangent atan(y / x). */

   248 double fp_atan2(MPL *mpl, double y, double x)

   249 {     xassert(mpl == mpl);

   250       return atan2(y, x);

   251 }

   253 /*----------------------------------------------------------------------

   254 -- fp_round - round floating-point value to n fractional digits.

   255 --

   256 -- This routine rounds given floating-point value x to n fractional

   257 -- digits with the formula:

   258 --

   259 --    round(x, n) = floor(x * 10^n + 0.5) / 10^n.

   260 --

   261 -- The parameter n is assumed to be integer. */

   263 double fp_round(MPL *mpl, double x, double n)

   264 {     double ten_to_n;

   265       if (n != floor(n))

   266          error(mpl, "round(%.*g, %.*g); non-integer second argument",

   267             DBL_DIG, x, DBL_DIG, n);

   268       if (n <= DBL_DIG + 2)

   269       {  ten_to_n = pow(10.0, n);

   270          if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)

   271          {  x = floor(x * ten_to_n + 0.5);

   272             if (x != 0.0) x /= ten_to_n;

   273          }

   274       }

   275       return x;

   276 }

   278 /*----------------------------------------------------------------------

   279 -- fp_trunc - truncate floating-point value to n fractional digits.

   280 --

   281 -- This routine truncates given floating-point value x to n fractional

   282 -- digits with the formula:

   283 --

   284 --                  ( floor(x * 10^n) / 10^n,  if x >= 0

   285 --    trunc(x, n) = <

   286 --                  ( ceil(x * 10^n) / 10^n,   if x < 0

   287 --

   288 -- The parameter n is assumed to be integer. */

   290 double fp_trunc(MPL *mpl, double x, double n)

   291 {     double ten_to_n;

   292       if (n != floor(n))

   293          error(mpl, "trunc(%.*g, %.*g); non-integer second argument",

   294             DBL_DIG, x, DBL_DIG, n);

   295       if (n <= DBL_DIG + 2)

   296       {  ten_to_n = pow(10.0, n);

   297          if (fabs(x) < (0.999 * DBL_MAX) / ten_to_n)

   298          {  x = (x >= 0.0 ? floor(x * ten_to_n) : ceil(x * ten_to_n));

   299             if (x != 0.0) x /= ten_to_n;

   300          }

   301       }

   302       return x;

   303 }

   305 /**********************************************************************/

   306 /* * *              PSEUDO-RANDOM NUMBER GENERATORS               * * */

   307 /**********************************************************************/

   309 /*----------------------------------------------------------------------

   310 -- fp_irand224 - pseudo-random integer in the range [0, 2^24).

   311 --

   312 -- This routine returns a next pseudo-random integer (converted to

   313 -- floating-point) which is uniformly distributed between 0 and 2^24-1,

   314 -- inclusive. */

   316 #define two_to_the_24 0x1000000

   318 double fp_irand224(MPL *mpl)

   319 {     return

   320          (double)rng_unif_rand(mpl->rand, two_to_the_24);

   321 }

   323 /*----------------------------------------------------------------------

   324 -- fp_uniform01 - pseudo-random number in the range [0, 1).

   325 --

   326 -- This routine returns a next pseudo-random number which is uniformly

   327 -- distributed in the range [0, 1). */

   329 #define two_to_the_31 ((unsigned int)0x80000000)

   331 double fp_uniform01(MPL *mpl)

   332 {     return

   333          (double)rng_next_rand(mpl->rand) / (double)two_to_the_31;

   334 }

   336 /*----------------------------------------------------------------------

   337 -- fp_uniform - pseudo-random number in the range [a, b).

   338 --

   339 -- This routine returns a next pseudo-random number which is uniformly

   340 -- distributed in the range [a, b). */

   342 double fp_uniform(MPL *mpl, double a, double b)

   343 {     double x;

   344       if (a >= b)

   345          error(mpl, "Uniform(%.*g, %.*g); invalid range",

   346             DBL_DIG, a, DBL_DIG, b);

   347       x = fp_uniform01(mpl);

   348 #if 0

   349       x = a * (1.0 - x) + b * x;

   350 #else

   351       x = fp_add(mpl, a * (1.0 - x), b * x);

   352 #endif

   353       return x;

   354 }

   356 /*----------------------------------------------------------------------

   357 -- fp_normal01 - Gaussian random variate with mu = 0 and sigma = 1.

   358 --

   359 -- This routine returns a Gaussian random variate with zero mean and

   360 -- unit standard deviation. The polar (Box-Mueller) method is used.

   361 --

   362 -- This code is a modified version of the routine gsl_ran_gaussian from

   363 -- the GNU Scientific Library Version 1.0. */

   365 double fp_normal01(MPL *mpl)

   366 {     double x, y, r2;

   367       do

   368       {  /* choose x, y in uniform square (-1,-1) to (+1,+1) */

   369          x = -1.0 + 2.0 * fp_uniform01(mpl);

   370          y = -1.0 + 2.0 * fp_uniform01(mpl);

   371          /* see if it is in the unit circle */

   372          r2 = x * x + y * y;

   373       } while (r2 > 1.0 || r2 == 0.0);

   374       /* Box-Muller transform */

   375       return y * sqrt(-2.0 * log (r2) / r2);

   376 }

   378 /*----------------------------------------------------------------------

   379 -- fp_normal - Gaussian random variate with specified mu and sigma.

   380 --

   381 -- This routine returns a Gaussian random variate with mean mu and

   382 -- standard deviation sigma. */

   384 double fp_normal(MPL *mpl, double mu, double sigma)

   385 {     double x;

   386 #if 0

   387       x = mu + sigma * fp_normal01(mpl);

   388 #else

   389       x = fp_add(mpl, mu, fp_mul(mpl, sigma, fp_normal01(mpl)));

   390 #endif

   391       return x;

   392 }

   394 /**********************************************************************/

   395 /* * *                SEGMENTED CHARACTER STRINGS                 * * */

   396 /**********************************************************************/

   398 /*----------------------------------------------------------------------

   399 -- create_string - create character string.

   400 --

   401 -- This routine creates a segmented character string, which is exactly

   402 -- equivalent to specified character string. */

   404 STRING *create_string

   405 (     MPL *mpl,

   406       char buf[MAX_LENGTH+1]  /* not changed */

   407 )

   408 #if 0

   409 {     STRING *head, *tail;

   410       int i, j;

   411       xassert(buf != NULL);

   412       xassert(strlen(buf) <= MAX_LENGTH);

   413       head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));

   414       for (i = j = 0; ; i++)

   415       {  if ((tail->seg[j++] = buf[i]) == '\0') break;

   416          if (j == STRSEG_SIZE)

   417 tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING))), j = 0;

   418       }

   419       tail->next = NULL;

   420       return head;

   421 }

   422 #else

   423 {     STRING *str;

   424       xassert(strlen(buf) <= MAX_LENGTH);

   425       str = dmp_get_atom(mpl->strings, strlen(buf)+1);

   426       strcpy(str, buf);

   427       return str;

   428 }

   429 #endif

   431 /*----------------------------------------------------------------------

   432 -- copy_string - make copy of character string.

   433 --

   434 -- This routine returns an exact copy of segmented character string. */

   436 STRING *copy_string

   437 (     MPL *mpl,

   438       STRING *str             /* not changed */

   439 )

   440 #if 0

   441 {     STRING *head, *tail;

   442       xassert(str != NULL);

   443       head = tail = dmp_get_atom(mpl->strings, sizeof(STRING));

   444       for (; str != NULL; str = str->next)

   445       {  memcpy(tail->seg, str->seg, STRSEG_SIZE);

   446          if (str->next != NULL)

   447 tail = (tail->next = dmp_get_atom(mpl->strings, sizeof(STRING)));

   448       }

   449       tail->next = NULL;

   450       return head;

   451 }

   452 #else

   453 {     xassert(mpl == mpl);

   454       return create_string(mpl, str);

   455 }

   456 #endif

   458 /*----------------------------------------------------------------------

   459 -- compare_strings - compare one character string with another.

   460 --

   461 -- This routine compares one segmented character strings with another

   462 -- and returns the result of comparison as follows:

   463 --

   464 -- = 0 - both strings are identical;

   465 -- < 0 - the first string precedes the second one;

   466 -- > 0 - the first string follows the second one. */

   468 int compare_strings

   469 (     MPL *mpl,

   470       STRING *str1,           /* not changed */

   471       STRING *str2            /* not changed */

   472 )

   473 #if 0

   474 {     int j, c1, c2;

   475       xassert(mpl == mpl);

   476       for (;; str1 = str1->next, str2 = str2->next)

   477       {  xassert(str1 != NULL);

   478          xassert(str2 != NULL);

   479          for (j = 0; j < STRSEG_SIZE; j++)

   480          {  c1 = (unsigned char)str1->seg[j];

   481             c2 = (unsigned char)str2->seg[j];

   482             if (c1 < c2) return -1;

   483             if (c1 > c2) return +1;

   484             if (c1 == '\0') goto done;

   485          }

   486       }

   487 done: return 0;

   488 }

   489 #else

   490 {     xassert(mpl == mpl);

   491       return strcmp(str1, str2);

   492 }

   493 #endif

   495 /*----------------------------------------------------------------------

   496 -- fetch_string - extract content of character string.

   497 --

   498 -- This routine returns a character string, which is exactly equivalent

   499 -- to specified segmented character string. */

   501 char *fetch_string

   502 (     MPL *mpl,

   503       STRING *str,            /* not changed */

   504       char buf[MAX_LENGTH+1]  /* modified */

   505 )

   506 #if 0

   507 {     int i, j;

   508       xassert(mpl == mpl);

   509       xassert(buf != NULL);

   510       for (i = 0; ; str = str->next)

   511       {  xassert(str != NULL);

   512          for (j = 0; j < STRSEG_SIZE; j++)

   513             if ((buf[i++] = str->seg[j]) == '\0') goto done;

   514       }

   515 done: xassert(strlen(buf) <= MAX_LENGTH);

   516       return buf;

   517 }

   518 #else

   519 {     xassert(mpl == mpl);

   520       return strcpy(buf, str);

   521 }

   522 #endif

   524 /*----------------------------------------------------------------------

   525 -- delete_string - delete character string.

   526 --

   527 -- This routine deletes specified segmented character string. */

   529 void delete_string

   530 (     MPL *mpl,

   531       STRING *str             /* destroyed */

   532 )

   533 #if 0

   534 {     STRING *temp;

   535       xassert(str != NULL);

   536       while (str != NULL)

   537       {  temp = str;

   538          str = str->next;

   539          dmp_free_atom(mpl->strings, temp, sizeof(STRING));

   540       }

   541       return;

   542 }

   543 #else

   544 {     dmp_free_atom(mpl->strings, str, strlen(str)+1);

   545       return;

   546 }

   547 #endif

   549 /**********************************************************************/

   550 /* * *                          SYMBOLS                           * * */

   551 /**********************************************************************/

   553 /*----------------------------------------------------------------------

   554 -- create_symbol_num - create symbol of numeric type.

   555 --

   556 -- This routine creates a symbol, which has a numeric value specified

   557 -- as floating-point number. */

   559 SYMBOL *create_symbol_num(MPL *mpl, double num)

   560 {     SYMBOL *sym;

   561       sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));

   562       sym->num = num;

   563       sym->str = NULL;

   564       return sym;

   565 }

   567 /*----------------------------------------------------------------------

   568 -- create_symbol_str - create symbol of abstract type.

   569 --

   570 -- This routine creates a symbol, which has an abstract value specified

   571 -- as segmented character string. */

   573 SYMBOL *create_symbol_str

   574 (     MPL *mpl,

   575       STRING *str             /* destroyed */

   576 )

   577 {     SYMBOL *sym;

   578       xassert(str != NULL);

   579       sym = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));

   580       sym->num = 0.0;

   581       sym->str = str;

   582       return sym;

   583 }

   585 /*----------------------------------------------------------------------

   586 -- copy_symbol - make copy of symbol.

   587 --

   588 -- This routine returns an exact copy of symbol. */

   590 SYMBOL *copy_symbol

   591 (     MPL *mpl,

   592       SYMBOL *sym             /* not changed */

   593 )

   594 {     SYMBOL *copy;

   595       xassert(sym != NULL);

   596       copy = dmp_get_atom(mpl->symbols, sizeof(SYMBOL));

   597       if (sym->str == NULL)

   598       {  copy->num = sym->num;

   599          copy->str = NULL;

   600       }

   601       else

   602       {  copy->num = 0.0;

   603          copy->str = copy_string(mpl, sym->str);

   604       }

   605       return copy;

   606 }

   608 /*----------------------------------------------------------------------

   609 -- compare_symbols - compare one symbol with another.

   610 --

   611 -- This routine compares one symbol with another and returns the result

   612 -- of comparison as follows:

   613 --

   614 -- = 0 - both symbols are identical;

   615 -- < 0 - the first symbol precedes the second one;

   616 -- > 0 - the first symbol follows the second one.

   617 --

   618 -- Note that the linear order, in which symbols follow each other, is

   619 -- implementation-dependent. It may be not an alphabetical order. */

   621 int compare_symbols

   622 (     MPL *mpl,

   623       SYMBOL *sym1,           /* not changed */

   624       SYMBOL *sym2            /* not changed */

   625 )

   626 {     xassert(sym1 != NULL);

   627       xassert(sym2 != NULL);

   628       /* let all numeric quantities precede all symbolic quantities */

   629       if (sym1->str == NULL && sym2->str == NULL)

   630       {  if (sym1->num < sym2->num) return -1;

   631          if (sym1->num > sym2->num) return +1;

   632          return 0;

   633       }

   634       if (sym1->str == NULL) return -1;

   635       if (sym2->str == NULL) return +1;

   636       return compare_strings(mpl, sym1->str, sym2->str);

   637 }

   639 /*----------------------------------------------------------------------

   640 -- delete_symbol - delete symbol.

   641 --

   642 -- This routine deletes specified symbol. */

   644 void delete_symbol

   645 (     MPL *mpl,

   646       SYMBOL *sym             /* destroyed */

   647 )

   648 {     xassert(sym != NULL);

   649       if (sym->str != NULL) delete_string(mpl, sym->str);

   650       dmp_free_atom(mpl->symbols, sym, sizeof(SYMBOL));

   651       return;

   652 }

   654 /*----------------------------------------------------------------------

   655 -- format_symbol - format symbol for displaying or printing.

   656 --

   657 -- This routine converts specified symbol to a charater string, which

   658 -- is suitable for displaying or printing.

   659 --

   660 -- The resultant string is never longer than 255 characters. If it gets

   661 -- longer, it is truncated from the right and appended by dots. */

   663 char *format_symbol

   664 (     MPL *mpl,

   665       SYMBOL *sym             /* not changed */

   666 )

   667 {     char *buf = mpl->sym_buf;

   668       xassert(sym != NULL);

   669       if (sym->str == NULL)

   670          sprintf(buf, "%.*g", DBL_DIG, sym->num);

   671       else

   672       {  char str[MAX_LENGTH+1];

   673          int quoted, j, len;

   674          fetch_string(mpl, sym->str, str);

   675          if (!(isalpha((unsigned char)str[0]) || str[0] == '_'))

   676             quoted = 1;

   677          else

   678          {  quoted = 0;

   679             for (j = 1; str[j] != '\0'; j++)

   680             {  if (!(isalnum((unsigned char)str[j]) ||

   681                      strchr("+-._", (unsigned char)str[j]) != NULL))

   682                {  quoted = 1;

   683                   break;

   684                }

   685             }

   686          }

   687 #        define safe_append(c) \

   688             (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)

   689          buf[0] = '\0', len = 0;

   690          if (quoted) safe_append('\'');

   691          for (j = 0; str[j] != '\0'; j++)

   692          {  if (quoted && str[j] == '\'') safe_append('\'');

   693             safe_append(str[j]);

   694          }

   695          if (quoted) safe_append('\'');

   696 #        undef safe_append

   697          buf[len] = '\0';

   698          if (len == 255) strcpy(buf+252, "...");

   699       }

   700       xassert(strlen(buf) <= 255);

   701       return buf;

   702 }

   704 /*----------------------------------------------------------------------

   705 -- concat_symbols - concatenate one symbol with another.

   706 --

   707 -- This routine concatenates values of two given symbols and assigns

   708 -- the resultant character string to a new symbol, which is returned on

   709 -- exit. Both original symbols are destroyed. */

   711 SYMBOL *concat_symbols

   712 (     MPL *mpl,

   713       SYMBOL *sym1,           /* destroyed */

   714       SYMBOL *sym2            /* destroyed */

   715 )

   716 {     char str1[MAX_LENGTH+1], str2[MAX_LENGTH+1];

   717       xassert(MAX_LENGTH >= DBL_DIG + DBL_DIG);

   718       if (sym1->str == NULL)

   719          sprintf(str1, "%.*g", DBL_DIG, sym1->num);

   720       else

   721          fetch_string(mpl, sym1->str, str1);

   722       if (sym2->str == NULL)

   723          sprintf(str2, "%.*g", DBL_DIG, sym2->num);

   724       else

   725          fetch_string(mpl, sym2->str, str2);

   726       if (strlen(str1) + strlen(str2) > MAX_LENGTH)

   727       {  char buf[255+1];

   728          strcpy(buf, format_symbol(mpl, sym1));

   729          xassert(strlen(buf) < sizeof(buf));

   730          error(mpl, "%s & %s; resultant symbol exceeds %d characters",

   731             buf, format_symbol(mpl, sym2), MAX_LENGTH);

   732       }

   733       delete_symbol(mpl, sym1);

   734       delete_symbol(mpl, sym2);

   735       return create_symbol_str(mpl, create_string(mpl, strcat(str1,

   736          str2)));

   737 }

   739 /**********************************************************************/

   740 /* * *                          N-TUPLES                          * * */

   741 /**********************************************************************/

   743 /*----------------------------------------------------------------------

   744 -- create_tuple - create n-tuple.

   745 --

   746 -- This routine creates a n-tuple, which initially has no components,

   747 -- i.e. which is 0-tuple. */

   749 TUPLE *create_tuple(MPL *mpl)

   750 {     TUPLE *tuple;

   751       xassert(mpl == mpl);

   752       tuple = NULL;

   753       return tuple;

   754 }

   756 /*----------------------------------------------------------------------

   757 -- expand_tuple - append symbol to n-tuple.

   758 --

   759 -- This routine expands n-tuple appending to it a given symbol, which

   760 -- becomes its new last component. */

   762 TUPLE *expand_tuple

   763 (     MPL *mpl,

   764       TUPLE *tuple,           /* destroyed */

   765       SYMBOL *sym             /* destroyed */

   766 )

   767 {     TUPLE *tail, *temp;

   768       xassert(sym != NULL);

   769       /* create a new component */

   770       tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));

   771       tail->sym = sym;

   772       tail->next = NULL;

   773       /* and append it to the component list */

   774       if (tuple == NULL)

   775          tuple = tail;

   776       else

   777       {  for (temp = tuple; temp->next != NULL; temp = temp->next);

   778          temp->next = tail;

   779       }

   780       return tuple;

   781 }

   783 /*----------------------------------------------------------------------

   784 -- tuple_dimen - determine dimension of n-tuple.

   785 --

   786 -- This routine returns dimension of n-tuple, i.e. number of components

   787 -- in the n-tuple. */

   789 int tuple_dimen

   790 (     MPL *mpl,

   791       TUPLE *tuple            /* not changed */

   792 )

   793 {     TUPLE *temp;

   794       int dim = 0;

   795       xassert(mpl == mpl);

   796       for (temp = tuple; temp != NULL; temp = temp->next) dim++;

   797       return dim;

   798 }

   800 /*----------------------------------------------------------------------

   801 -- copy_tuple - make copy of n-tuple.

   802 --

   803 -- This routine returns an exact copy of n-tuple. */

   805 TUPLE *copy_tuple

   806 (     MPL *mpl,

   807       TUPLE *tuple            /* not changed */

   808 )

   809 {     TUPLE *head, *tail;

   810       if (tuple == NULL)

   811          head = NULL;

   812       else

   813       {  head = tail = dmp_get_atom(mpl->tuples, sizeof(TUPLE));

   814          for (; tuple != NULL; tuple = tuple->next)

   815          {  xassert(tuple->sym != NULL);

   816             tail->sym = copy_symbol(mpl, tuple->sym);

   817             if (tuple->next != NULL)

   818 tail = (tail->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));

   819          }

   820          tail->next = NULL;

   821       }

   822       return head;

   823 }

   825 /*----------------------------------------------------------------------

   826 -- compare_tuples - compare one n-tuple with another.

   827 --

   828 -- This routine compares two given n-tuples, which must have the same

   829 -- dimension (not checked for the sake of efficiency), and returns one

   830 -- of the following codes:

   831 --

   832 -- = 0 - both n-tuples are identical;

   833 -- < 0 - the first n-tuple precedes the second one;

   834 -- > 0 - the first n-tuple follows the second one.

   835 --

   836 -- Note that the linear order, in which n-tuples follow each other, is

   837 -- implementation-dependent. It may be not an alphabetical order. */

   839 int compare_tuples

   840 (     MPL *mpl,

   841       TUPLE *tuple1,          /* not changed */

   842       TUPLE *tuple2           /* not changed */

   843 )

   844 {     TUPLE *item1, *item2;

   845       int ret;

   846       xassert(mpl == mpl);

   847       for (item1 = tuple1, item2 = tuple2; item1 != NULL;

   848            item1 = item1->next, item2 = item2->next)

   849       {  xassert(item2 != NULL);

   850          xassert(item1->sym != NULL);

   851          xassert(item2->sym != NULL);

   852          ret = compare_symbols(mpl, item1->sym, item2->sym);

   853          if (ret != 0) return ret;

   854       }

   855       xassert(item2 == NULL);

   856       return 0;

   857 }

   859 /*----------------------------------------------------------------------

   860 -- build_subtuple - build subtuple of given n-tuple.

   861 --

   862 -- This routine builds subtuple, which consists of first dim components

   863 -- of given n-tuple. */

   865 TUPLE *build_subtuple

   866 (     MPL *mpl,

   867       TUPLE *tuple,           /* not changed */

   868       int dim

   869 )

   870 {     TUPLE *head, *temp;

   871       int j;

   872       head = create_tuple(mpl);

   873       for (j = 1, temp = tuple; j <= dim; j++, temp = temp->next)

   874       {  xassert(temp != NULL);

   875          head = expand_tuple(mpl, head, copy_symbol(mpl, temp->sym));

   876       }

   877       return head;

   878 }

   880 /*----------------------------------------------------------------------

   881 -- delete_tuple - delete n-tuple.

   882 --

   883 -- This routine deletes specified n-tuple. */

   885 void delete_tuple

   886 (     MPL *mpl,

   887       TUPLE *tuple            /* destroyed */

   888 )

   889 {     TUPLE *temp;

   890       while (tuple != NULL)

   891       {  temp = tuple;

   892          tuple = temp->next;

   893          xassert(temp->sym != NULL);

   894          delete_symbol(mpl, temp->sym);

   895          dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));

   896       }

   897       return;

   898 }

   900 /*----------------------------------------------------------------------

   901 -- format_tuple - format n-tuple for displaying or printing.

   902 --

   903 -- This routine converts specified n-tuple to a character string, which

   904 -- is suitable for displaying or printing.

   905 --

   906 -- The resultant string is never longer than 255 characters. If it gets

   907 -- longer, it is truncated from the right and appended by dots. */

   909 char *format_tuple

   910 (     MPL *mpl,

   911       int c,

   912       TUPLE *tuple            /* not changed */

   913 )

   914 {     TUPLE *temp;

   915       int dim, j, len;

   916       char *buf = mpl->tup_buf, str[255+1], *save;

   917 #     define safe_append(c) \

   918          (void)(len < 255 ? (buf[len++] = (char)(c)) : 0)

   919       buf[0] = '\0', len = 0;

   920       dim = tuple_dimen(mpl, tuple);

   921       if (c == '[' && dim > 0) safe_append('[');

   922       if (c == '(' && dim > 1) safe_append('(');

   923       for (temp = tuple; temp != NULL; temp = temp->next)

   924       {  if (temp != tuple) safe_append(',');

   925          xassert(temp->sym != NULL);

   926          save = mpl->sym_buf;

   927          mpl->sym_buf = str;

   928          format_symbol(mpl, temp->sym);

   929          mpl->sym_buf = save;

   930          xassert(strlen(str) < sizeof(str));

   931          for (j = 0; str[j] != '\0'; j++) safe_append(str[j]);

   932       }

   933       if (c == '[' && dim > 0) safe_append(']');

   934       if (c == '(' && dim > 1) safe_append(')');

   935 #     undef safe_append

   936       buf[len] = '\0';

   937       if (len == 255) strcpy(buf+252, "...");

   938       xassert(strlen(buf) <= 255);

   939       return buf;

   940 }

   942 /**********************************************************************/

   943 /* * *                       ELEMENTAL SETS                       * * */

   944 /**********************************************************************/

   946 /*----------------------------------------------------------------------

   947 -- create_elemset - create elemental set.

   948 --

   949 -- This routine creates an elemental set, whose members are n-tuples of

   950 -- specified dimension. Being created the set is initially empty. */

   952 ELEMSET *create_elemset(MPL *mpl, int dim)

   953 {     ELEMSET *set;

   954       xassert(dim > 0);

   955       set = create_array(mpl, A_NONE, dim);

   956       return set;

   957 }

   959 /*----------------------------------------------------------------------

   960 -- find_tuple - check if elemental set contains given n-tuple.

   961 --

   962 -- This routine finds given n-tuple in specified elemental set in order

   963 -- to check if the set contains that n-tuple. If the n-tuple is found,

   964 -- the routine returns pointer to corresponding array member. Otherwise

   965 -- null pointer is returned. */

   967 MEMBER *find_tuple

   968 (     MPL *mpl,

   969       ELEMSET *set,           /* not changed */

   970       TUPLE *tuple            /* not changed */

   971 )

   972 {     xassert(set != NULL);

   973       xassert(set->type == A_NONE);

   974       xassert(set->dim == tuple_dimen(mpl, tuple));

   975       return find_member(mpl, set, tuple);

   976 }

   978 /*----------------------------------------------------------------------

   979 -- add_tuple - add new n-tuple to elemental set.

   980 --

   981 -- This routine adds given n-tuple to specified elemental set.

   982 --

   983 -- For the sake of efficiency this routine doesn't check whether the

   984 -- set already contains the same n-tuple or not. Therefore the calling

   985 -- program should use the routine find_tuple (if necessary) in order to

   986 -- make sure that the given n-tuple is not contained in the set, since

   987 -- duplicate n-tuples within the same set are not allowed. */

   989 MEMBER *add_tuple

   990 (     MPL *mpl,

   991       ELEMSET *set,           /* modified */

   992       TUPLE *tuple            /* destroyed */

   993 )

   994 {     MEMBER *memb;

   995       xassert(set != NULL);

   996       xassert(set->type == A_NONE);

   997       xassert(set->dim == tuple_dimen(mpl, tuple));

   998       memb = add_member(mpl, set, tuple);

   999       memb->value.none = NULL;

  1000       return memb;

  1001 }

  1003 /*----------------------------------------------------------------------

  1004 -- check_then_add - check and add new n-tuple to elemental set.

  1005 --

  1006 -- This routine is equivalent to the routine add_tuple except that it

  1007 -- does check for duplicate n-tuples. */

  1009 MEMBER *check_then_add

  1010 (     MPL *mpl,

  1011       ELEMSET *set,           /* modified */

  1012       TUPLE *tuple            /* destroyed */

  1013 )

  1014 {     if (find_tuple(mpl, set, tuple) != NULL)

  1015          error(mpl, "duplicate tuple %s detected", format_tuple(mpl,

  1016             '(', tuple));

  1017       return add_tuple(mpl, set, tuple);

  1018 }

  1020 /*----------------------------------------------------------------------

  1021 -- copy_elemset - make copy of elemental set.

  1022 --

  1023 -- This routine makes an exact copy of elemental set. */

  1025 ELEMSET *copy_elemset

  1026 (     MPL *mpl,

  1027       ELEMSET *set            /* not changed */

  1028 )

  1029 {     ELEMSET *copy;

  1030       MEMBER *memb;

  1031       xassert(set != NULL);

  1032       xassert(set->type == A_NONE);

  1033       xassert(set->dim > 0);

  1034       copy = create_elemset(mpl, set->dim);

  1035       for (memb = set->head; memb != NULL; memb = memb->next)

  1036          add_tuple(mpl, copy, copy_tuple(mpl, memb->tuple));

  1037       return copy;

  1038 }

  1040 /*----------------------------------------------------------------------

  1041 -- delete_elemset - delete elemental set.

  1042 --

  1043 -- This routine deletes specified elemental set. */

  1045 void delete_elemset

  1046 (     MPL *mpl,

  1047       ELEMSET *set            /* destroyed */

  1048 )

  1049 {     xassert(set != NULL);

  1050       xassert(set->type == A_NONE);

  1051       delete_array(mpl, set);

  1052       return;

  1053 }

  1055 /*----------------------------------------------------------------------

  1056 -- arelset_size - compute size of "arithmetic" elemental set.

  1057 --

  1058 -- This routine computes the size of "arithmetic" elemental set, which

  1059 -- is specified in the form of arithmetic progression:

  1060 --

  1061 --    { t0 .. tf by dt }.

  1062 --

  1063 -- The size is computed using the formula:

  1064 --

  1065 --    n = max(0, floor((tf - t0) / dt) + 1). */

  1067 int arelset_size(MPL *mpl, double t0, double tf, double dt)

  1068 {     double temp;

  1069       if (dt == 0.0)

  1070          error(mpl, "%.*g .. %.*g by %.*g; zero stride not allowed",

  1071             DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);

  1072       if (tf > 0.0 && t0 < 0.0 && tf > + 0.999 * DBL_MAX + t0)

  1073          temp = +DBL_MAX;

  1074       else if (tf < 0.0 && t0 > 0.0 && tf < - 0.999 * DBL_MAX + t0)

  1075          temp = -DBL_MAX;

  1076       else

  1077          temp = tf - t0;

  1078       if (fabs(dt) < 1.0 && fabs(temp) > (0.999 * DBL_MAX) * fabs(dt))

  1079       {  if (temp > 0.0 && dt > 0.0 || temp < 0.0 && dt < 0.0)

  1080             temp = +DBL_MAX;

  1081          else

  1082             temp = 0.0;

  1083       }

  1084       else

  1085       {  temp = floor(temp / dt) + 1.0;

  1086          if (temp < 0.0) temp = 0.0;

  1087       }

  1088       xassert(temp >= 0.0);

  1089       if (temp > (double)(INT_MAX - 1))

  1090          error(mpl, "%.*g .. %.*g by %.*g; set too large",

  1091             DBL_DIG, t0, DBL_DIG, tf, DBL_DIG, dt);

  1092       return (int)(temp + 0.5);

  1093 }

  1095 /*----------------------------------------------------------------------

  1096 -- arelset_member - compute member of "arithmetic" elemental set.

  1097 --

  1098 -- This routine returns a numeric value of symbol, which is equivalent

  1099 -- to j-th member of given "arithmetic" elemental set specified in the

  1100 -- form of arithmetic progression:

  1101 --

  1102 --    { t0 .. tf by dt }.

  1103 --

  1104 -- The symbol value is computed with the formula:

  1105 --

  1106 --    j-th member = t0 + (j - 1) * dt,

  1107 --

  1108 -- The number j must satisfy to the restriction 1 <= j <= n, where n is

  1109 -- the set size computed by the routine arelset_size. */

  1111 double arelset_member(MPL *mpl, double t0, double tf, double dt, int j)

  1112 {     xassert(1 <= j && j <= arelset_size(mpl, t0, tf, dt));

  1113       return t0 + (double)(j - 1) * dt;

  1114 }

  1116 /*----------------------------------------------------------------------

  1117 -- create_arelset - create "arithmetic" elemental set.

  1118 --

  1119 -- This routine creates "arithmetic" elemental set, which is specified

  1120 -- in the form of arithmetic progression:

  1121 --

  1122 --    { t0 .. tf by dt }.

  1123 --

  1124 -- Components of this set are 1-tuples. */

  1126 ELEMSET *create_arelset(MPL *mpl, double t0, double tf, double dt)

  1127 {     ELEMSET *set;

  1128       int j, n;

  1129       set = create_elemset(mpl, 1);

  1130       n = arelset_size(mpl, t0, tf, dt);

  1131       for (j = 1; j <= n; j++)

  1132       {  add_tuple

  1133          (  mpl,

  1134             set,

  1135             expand_tuple

  1136             (  mpl,

  1137                create_tuple(mpl),

  1138                create_symbol_num

  1139                (  mpl,

  1140                   arelset_member(mpl, t0, tf, dt, j)

  1141                )

  1142             )

  1143          );

  1144       }

  1145       return set;

  1146 }

  1148 /*----------------------------------------------------------------------

  1149 -- set_union - union of two elemental sets.

  1150 --

  1151 -- This routine computes the union:

  1152 --

  1153 --    X U Y = { j | (j in X) or (j in Y) },

  1154 --

  1155 -- where X and Y are given elemental sets (destroyed on exit). */

  1157 ELEMSET *set_union

  1158 (     MPL *mpl,

  1159       ELEMSET *X,             /* destroyed */

  1160       ELEMSET *Y              /* destroyed */

  1161 )

  1162 {     MEMBER *memb;

  1163       xassert(X != NULL);

  1164       xassert(X->type == A_NONE);

  1165       xassert(X->dim > 0);

  1166       xassert(Y != NULL);

  1167       xassert(Y->type == A_NONE);

  1168       xassert(Y->dim > 0);

  1169       xassert(X->dim == Y->dim);

  1170       for (memb = Y->head; memb != NULL; memb = memb->next)

  1171       {  if (find_tuple(mpl, X, memb->tuple) == NULL)

  1172             add_tuple(mpl, X, copy_tuple(mpl, memb->tuple));

  1173       }

  1174       delete_elemset(mpl, Y);

  1175       return X;

  1176 }

  1178 /*----------------------------------------------------------------------

  1179 -- set_diff - difference between two elemental sets.

  1180 --

  1181 -- This routine computes the difference:

  1182 --

  1183 --    X \ Y = { j | (j in X) and (j not in Y) },

  1184 --

  1185 -- where X and Y are given elemental sets (destroyed on exit). */

  1187 ELEMSET *set_diff

  1188 (     MPL *mpl,

  1189       ELEMSET *X,             /* destroyed */

  1190       ELEMSET *Y              /* destroyed */

  1191 )

  1192 {     ELEMSET *Z;

  1193       MEMBER *memb;

  1194       xassert(X != NULL);

  1195       xassert(X->type == A_NONE);

  1196       xassert(X->dim > 0);

  1197       xassert(Y != NULL);

  1198       xassert(Y->type == A_NONE);

  1199       xassert(Y->dim > 0);

  1200       xassert(X->dim == Y->dim);

  1201       Z = create_elemset(mpl, X->dim);

  1202       for (memb = X->head; memb != NULL; memb = memb->next)

  1203       {  if (find_tuple(mpl, Y, memb->tuple) == NULL)

  1204             add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));

  1205       }

  1206       delete_elemset(mpl, X);

  1207       delete_elemset(mpl, Y);

  1208       return Z;

  1209 }

  1211 /*----------------------------------------------------------------------

  1212 -- set_symdiff - symmetric difference between two elemental sets.

  1213 --

  1214 -- This routine computes the symmetric difference:

  1215 --

  1216 --    X (+) Y = (X \ Y) U (Y \ X),

  1217 --

  1218 -- where X and Y are given elemental sets (destroyed on exit). */

  1220 ELEMSET *set_symdiff

  1221 (     MPL *mpl,

  1222       ELEMSET *X,             /* destroyed */

  1223       ELEMSET *Y              /* destroyed */

  1224 )

  1225 {     ELEMSET *Z;

  1226       MEMBER *memb;

  1227       xassert(X != NULL);

  1228       xassert(X->type == A_NONE);

  1229       xassert(X->dim > 0);

  1230       xassert(Y != NULL);

  1231       xassert(Y->type == A_NONE);

  1232       xassert(Y->dim > 0);

  1233       xassert(X->dim == Y->dim);

  1234       /* Z := X \ Y */

  1235       Z = create_elemset(mpl, X->dim);

  1236       for (memb = X->head; memb != NULL; memb = memb->next)

  1237       {  if (find_tuple(mpl, Y, memb->tuple) == NULL)

  1238             add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));

  1239       }

  1240       /* Z := Z U (Y \ X) */

  1241       for (memb = Y->head; memb != NULL; memb = memb->next)

  1242       {  if (find_tuple(mpl, X, memb->tuple) == NULL)

  1243             add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));

  1244       }

  1245       delete_elemset(mpl, X);

  1246       delete_elemset(mpl, Y);

  1247       return Z;

  1248 }

  1250 /*----------------------------------------------------------------------

  1251 -- set_inter - intersection of two elemental sets.

  1252 --

  1253 -- This routine computes the intersection:

  1254 --

  1255 --    X ^ Y = { j | (j in X) and (j in Y) },

  1256 --

  1257 -- where X and Y are given elemental sets (destroyed on exit). */

  1259 ELEMSET *set_inter

  1260 (     MPL *mpl,

  1261       ELEMSET *X,             /* destroyed */

  1262       ELEMSET *Y              /* destroyed */

  1263 )

  1264 {     ELEMSET *Z;

  1265       MEMBER *memb;

  1266       xassert(X != NULL);

  1267       xassert(X->type == A_NONE);

  1268       xassert(X->dim > 0);

  1269       xassert(Y != NULL);

  1270       xassert(Y->type == A_NONE);

  1271       xassert(Y->dim > 0);

  1272       xassert(X->dim == Y->dim);

  1273       Z = create_elemset(mpl, X->dim);

  1274       for (memb = X->head; memb != NULL; memb = memb->next)

  1275       {  if (find_tuple(mpl, Y, memb->tuple) != NULL)

  1276             add_tuple(mpl, Z, copy_tuple(mpl, memb->tuple));

  1277       }

  1278       delete_elemset(mpl, X);

  1279       delete_elemset(mpl, Y);

  1280       return Z;

  1281 }

  1283 /*----------------------------------------------------------------------

  1284 -- set_cross - cross (Cartesian) product of two elemental sets.

  1285 --

  1286 -- This routine computes the cross (Cartesian) product:

  1287 --

  1288 --    X x Y = { (i,j) | (i in X) and (j in Y) },

  1289 --

  1290 -- where X and Y are given elemental sets (destroyed on exit). */

  1292 ELEMSET *set_cross

  1293 (     MPL *mpl,

  1294       ELEMSET *X,             /* destroyed */

  1295       ELEMSET *Y              /* destroyed */

  1296 )

  1297 {     ELEMSET *Z;

  1298       MEMBER *memx, *memy;

  1299       TUPLE *tuple, *temp;

  1300       xassert(X != NULL);

  1301       xassert(X->type == A_NONE);

  1302       xassert(X->dim > 0);

  1303       xassert(Y != NULL);

  1304       xassert(Y->type == A_NONE);

  1305       xassert(Y->dim > 0);

  1306       Z = create_elemset(mpl, X->dim + Y->dim);

  1307       for (memx = X->head; memx != NULL; memx = memx->next)

  1308       {  for (memy = Y->head; memy != NULL; memy = memy->next)

  1309          {  tuple = copy_tuple(mpl, memx->tuple);

  1310             for (temp = memy->tuple; temp != NULL; temp = temp->next)

  1311                tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,

  1312                   temp->sym));

  1313             add_tuple(mpl, Z, tuple);

  1314          }

  1315       }

  1316       delete_elemset(mpl, X);

  1317       delete_elemset(mpl, Y);

  1318       return Z;

  1319 }

  1321 /**********************************************************************/

  1322 /* * *                    ELEMENTAL VARIABLES                     * * */

  1323 /**********************************************************************/

  1325 /* (there are no specific routines for elemental variables) */

  1327 /**********************************************************************/

  1328 /* * *                        LINEAR FORMS                        * * */

  1329 /**********************************************************************/

  1331 /*----------------------------------------------------------------------

  1332 -- constant_term - create constant term.

  1333 --

  1334 -- This routine creates the linear form, which is a constant term. */

  1336 FORMULA *constant_term(MPL *mpl, double coef)

  1337 {     FORMULA *form;

  1338       if (coef == 0.0)

  1339          form = NULL;

  1340       else

  1341       {  form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));

  1342          form->coef = coef;

  1343          form->var = NULL;

  1344          form->next = NULL;

  1345       }

  1346       return form;

  1347 }

  1349 /*----------------------------------------------------------------------

  1350 -- single_variable - create single variable.

  1351 --

  1352 -- This routine creates the linear form, which is a single elemental

  1353 -- variable. */

  1355 FORMULA *single_variable

  1356 (     MPL *mpl,

  1357       ELEMVAR *var            /* referenced */

  1358 )

  1359 {     FORMULA *form;

  1360       xassert(var != NULL);

  1361       form = dmp_get_atom(mpl->formulae, sizeof(FORMULA));

  1362       form->coef = 1.0;

  1363       form->var = var;

  1364       form->next = NULL;

  1365       return form;

  1366 }

  1368 /*----------------------------------------------------------------------

  1369 -- copy_formula - make copy of linear form.

  1370 --

  1371 -- This routine returns an exact copy of linear form. */

  1373 FORMULA *copy_formula

  1374 (     MPL *mpl,

  1375       FORMULA *form           /* not changed */

  1376 )

  1377 {     FORMULA *head, *tail;

  1378       if (form == NULL)

  1379          head = NULL;

  1380       else

  1381       {  head = tail = dmp_get_atom(mpl->formulae, sizeof(FORMULA));

  1382          for (; form != NULL; form = form->next)

  1383          {  tail->coef = form->coef;

  1384             tail->var = form->var;

  1385             if (form->next != NULL)

  1386 tail = (tail->next = dmp_get_atom(mpl->formulae, sizeof(FORMULA)));

  1387          }

  1388          tail->next = NULL;

  1389       }

  1390       return head;

  1391 }

  1393 /*----------------------------------------------------------------------

  1394 -- delete_formula - delete linear form.

  1395 --

  1396 -- This routine deletes specified linear form. */

  1398 void delete_formula

  1399 (     MPL *mpl,

  1400       FORMULA *form           /* destroyed */

  1401 )

  1402 {     FORMULA *temp;

  1403       while (form != NULL)

  1404       {  temp = form;

  1405          form = form->next;

  1406          dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));

  1407       }

  1408       return;

  1409 }

  1411 /*----------------------------------------------------------------------

  1412 -- linear_comb - linear combination of two linear forms.

  1413 --

  1414 -- This routine computes the linear combination:

  1415 --

  1416 --    a * fx + b * fy,

  1417 --

  1418 -- where a and b are numeric coefficients, fx and fy are linear forms

  1419 -- (destroyed on exit). */

  1421 FORMULA *linear_comb

  1422 (     MPL *mpl,

  1423       double a, FORMULA *fx,  /* destroyed */

  1424       double b, FORMULA *fy   /* destroyed */

  1425 )

  1426 {     FORMULA *form = NULL, *term, *temp;

  1427       double c0 = 0.0;

  1428       for (term = fx; term != NULL; term = term->next)

  1429       {  if (term->var == NULL)

  1430             c0 = fp_add(mpl, c0, fp_mul(mpl, a, term->coef));

  1431          else

  1432             term->var->temp =

  1433                fp_add(mpl, term->var->temp, fp_mul(mpl, a, term->coef));

  1434       }

  1435       for (term = fy; term != NULL; term = term->next)

  1436       {  if (term->var == NULL)

  1437             c0 = fp_add(mpl, c0, fp_mul(mpl, b, term->coef));

  1438          else

  1439             term->var->temp =

  1440                fp_add(mpl, term->var->temp, fp_mul(mpl, b, term->coef));

  1441       }

  1442       for (term = fx; term != NULL; term = term->next)

  1443       {  if (term->var != NULL && term->var->temp != 0.0)

  1444          {  temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));

  1445             temp->coef = term->var->temp, temp->var = term->var;

  1446             temp->next = form, form = temp;

  1447             term->var->temp = 0.0;

  1448          }

  1449       }

  1450       for (term = fy; term != NULL; term = term->next)

  1451       {  if (term->var != NULL && term->var->temp != 0.0)

  1452          {  temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));

  1453             temp->coef = term->var->temp, temp->var = term->var;

  1454             temp->next = form, form = temp;

  1455             term->var->temp = 0.0;

  1456          }

  1457       }

  1458       if (c0 != 0.0)

  1459       {  temp = dmp_get_atom(mpl->formulae, sizeof(FORMULA));

  1460          temp->coef = c0, temp->var = NULL;

  1461          temp->next = form, form = temp;

  1462       }

  1463       delete_formula(mpl, fx);

  1464       delete_formula(mpl, fy);

  1465       return form;

  1466 }

  1468 /*----------------------------------------------------------------------

  1469 -- remove_constant - remove constant term from linear form.

  1470 --

  1471 -- This routine removes constant term from linear form and stores its

  1472 -- value to given location. */

  1474 FORMULA *remove_constant

  1475 (     MPL *mpl,

  1476       FORMULA *form,          /* destroyed */

  1477       double *coef            /* modified */

  1478 )

  1479 {     FORMULA *head = NULL, *temp;

  1480       *coef = 0.0;

  1481       while (form != NULL)

  1482       {  temp = form;

  1483          form = form->next;

  1484          if (temp->var == NULL)

  1485          {  /* constant term */

  1486             *coef = fp_add(mpl, *coef, temp->coef);

  1487             dmp_free_atom(mpl->formulae, temp, sizeof(FORMULA));

  1488          }

  1489          else

  1490          {  /* linear term */

  1491             temp->next = head;

  1492             head = temp;

  1493          }

  1494       }

  1495       return head;

  1496 }

  1498 /*----------------------------------------------------------------------

  1499 -- reduce_terms - reduce identical terms in linear form.

  1500 --

  1501 -- This routine reduces identical terms in specified linear form. */

  1503 FORMULA *reduce_terms

  1504 (     MPL *mpl,

  1505       FORMULA *form           /* destroyed */

  1506 )

  1507 {     FORMULA *term, *next_term;

  1508       double c0 = 0.0;

  1509       for (term = form; term != NULL; term = term->next)

  1510       {  if (term->var == NULL)

  1511             c0 = fp_add(mpl, c0, term->coef);

  1512          else

  1513             term->var->temp = fp_add(mpl, term->var->temp, term->coef);

  1514       }

  1515       next_term = form, form = NULL;

  1516       for (term = next_term; term != NULL; term = next_term)

  1517       {  next_term = term->next;

  1518          if (term->var == NULL && c0 != 0.0)

  1519          {  term->coef = c0, c0 = 0.0;

  1520             term->next = form, form = term;

  1521          }

  1522          else if (term->var != NULL && term->var->temp != 0.0)

  1523          {  term->coef = term->var->temp, term->var->temp = 0.0;

  1524             term->next = form, form = term;

  1525          }

  1526          else

  1527             dmp_free_atom(mpl->formulae, term, sizeof(FORMULA));

  1528       }

  1529       return form;

  1530 }

  1532 /**********************************************************************/

  1533 /* * *                   ELEMENTAL CONSTRAINTS                    * * */

  1534 /**********************************************************************/

  1536 /* (there are no specific routines for elemental constraints) */

  1538 /**********************************************************************/

  1539 /* * *                       GENERIC VALUES                       * * */

  1540 /**********************************************************************/

  1542 /*----------------------------------------------------------------------

  1543 -- delete_value - delete generic value.

  1544 --

  1545 -- This routine deletes specified generic value.

  1546 --

  1547 -- NOTE: The generic value to be deleted must be valid. */

  1549 void delete_value

  1550 (     MPL *mpl,

  1551       int type,

  1552       VALUE *value            /* content destroyed */

  1553 )

  1554 {     xassert(value != NULL);

  1555       switch (type)

  1556       {  case A_NONE:

  1557             value->none = NULL;

  1558             break;

  1559          case A_NUMERIC:

  1560             value->num = 0.0;

  1561             break;

  1562          case A_SYMBOLIC:

  1563             delete_symbol(mpl, value->sym), value->sym = NULL;

  1564             break;

  1565          case A_LOGICAL:

  1566             value->bit = 0;

  1567             break;

  1568          case A_TUPLE:

  1569             delete_tuple(mpl, value->tuple), value->tuple = NULL;

  1570             break;

  1571          case A_ELEMSET:

  1572             delete_elemset(mpl, value->set), value->set = NULL;

  1573             break;

  1574          case A_ELEMVAR:

  1575             value->var = NULL;

  1576             break;

  1577          case A_FORMULA:

  1578             delete_formula(mpl, value->form), value->form = NULL;

  1579             break;

  1580          case A_ELEMCON:

  1581             value->con = NULL;

  1582             break;

  1583          default:

  1584             xassert(type != type);

  1585       }

  1586       return;

  1587 }

  1589 /**********************************************************************/

  1590 /* * *                SYMBOLICALLY INDEXED ARRAYS                 * * */

  1591 /**********************************************************************/

  1593 /*----------------------------------------------------------------------

  1594 -- create_array - create array.

  1595 --

  1596 -- This routine creates an array of specified type and dimension. Being

  1597 -- created the array is initially empty.

  1598 --

  1599 -- The type indicator determines generic values, which can be assigned

  1600 -- to the array members:

  1601 --

  1602 -- A_NONE     - none (members have no assigned values)

  1603 -- A_NUMERIC  - floating-point numbers

  1604 -- A_SYMBOLIC - symbols

  1605 -- A_ELEMSET  - elemental sets

  1606 -- A_ELEMVAR  - elemental variables

  1607 -- A_ELEMCON  - elemental constraints

  1608 --

  1609 -- The dimension may be 0, in which case the array consists of the only

  1610 -- member (such arrays represent 0-dimensional objects). */

  1612 ARRAY *create_array(MPL *mpl, int type, int dim)

  1613 {     ARRAY *array;

  1614       xassert(type == A_NONE || type == A_NUMERIC ||

  1615              type == A_SYMBOLIC || type == A_ELEMSET ||

  1616              type == A_ELEMVAR || type == A_ELEMCON);

  1617       xassert(dim >= 0);

  1618       array = dmp_get_atom(mpl->arrays, sizeof(ARRAY));

  1619       array->type = type;

  1620       array->dim = dim;

  1621       array->size = 0;

  1622       array->head = NULL;

  1623       array->tail = NULL;

  1624       array->tree = NULL;

  1625       array->prev = NULL;

  1626       array->next = mpl->a_list;

  1627       /* include the array in the global array list */

  1628       if (array->next != NULL) array->next->prev = array;

  1629       mpl->a_list = array;

  1630       return array;

  1631 }

  1633 /*----------------------------------------------------------------------

  1634 -- find_member - find array member with given n-tuple.

  1635 --

  1636 -- This routine finds an array member, which has given n-tuple. If the

  1637 -- array is short, the linear search is used. Otherwise the routine

  1638 -- autimatically creates the search tree (i.e. the array index) to find

  1639 -- members for logarithmic time. */

  1641 static int compare_member_tuples(void *info, const void *key1,

  1642       const void *key2)

  1643 {     /* this is an auxiliary routine used to compare keys, which are

  1644          n-tuples assigned to array members */

  1645       return compare_tuples((MPL *)info, (TUPLE *)key1, (TUPLE *)key2);

  1646 }

  1648 MEMBER *find_member

  1649 (     MPL *mpl,

  1650       ARRAY *array,           /* not changed */

  1651       TUPLE *tuple            /* not changed */

  1652 )

  1653 {     MEMBER *memb;

  1654       xassert(array != NULL);

  1655       /* the n-tuple must have the same dimension as the array */

  1656       xassert(tuple_dimen(mpl, tuple) == array->dim);

  1657       /* if the array is large enough, create the search tree and index

  1658          all existing members of the array */

  1659       if (array->size > 30 && array->tree == NULL)

  1660       {  array->tree = avl_create_tree(compare_member_tuples, mpl);

  1661          for (memb = array->head; memb != NULL; memb = memb->next)

  1662 avl_set_node_link(avl_insert_node(array->tree, memb->tuple),

  1663                (void *)memb);

  1664       }

  1665       /* find a member, which has the given tuple */

  1666       if (array->tree == NULL)

  1667       {  /* the search tree doesn't exist; use the linear search */

  1668          for (memb = array->head; memb != NULL; memb = memb->next)

  1669             if (compare_tuples(mpl, memb->tuple, tuple) == 0) break;

  1670       }

  1671       else

  1672       {  /* the search tree exists; use the binary search */

  1673          AVLNODE *node;

  1674          node = avl_find_node(array->tree, tuple);

  1675 memb = (MEMBER *)(node == NULL ? NULL : avl_get_node_link(node));

  1676       }

  1677       return memb;

  1678 }

  1680 /*----------------------------------------------------------------------

  1681 -- add_member - add new member to array.

  1682 --

  1683 -- This routine creates a new member with given n-tuple and adds it to

  1684 -- specified array.

  1685 --

  1686 -- For the sake of efficiency this routine doesn't check whether the

  1687 -- array already contains a member with the given n-tuple or not. Thus,

  1688 -- if necessary, the calling program should use the routine find_member

  1689 -- in order to be sure that the array contains no member with the same

  1690 -- n-tuple, because members with duplicate n-tuples are not allowed.

  1691 --

  1692 -- This routine assigns no generic value to the new member, because the

  1693 -- calling program must do that. */

  1695 MEMBER *add_member

  1696 (     MPL *mpl,

  1697       ARRAY *array,           /* modified */

  1698       TUPLE *tuple            /* destroyed */

  1699 )

  1700 {     MEMBER *memb;

  1701       xassert(array != NULL);

  1702       /* the n-tuple must have the same dimension as the array */

  1703       xassert(tuple_dimen(mpl, tuple) == array->dim);

  1704       /* create new member */

  1705       memb = dmp_get_atom(mpl->members, sizeof(MEMBER));

  1706       memb->tuple = tuple;

  1707       memb->next = NULL;

  1708       memset(&memb->value, '?', sizeof(VALUE));

  1709       /* and append it to the member list */

  1710       array->size++;

  1711       if (array->head == NULL)

  1712          array->head = memb;

  1713       else

  1714          array->tail->next = memb;

  1715       array->tail = memb;

  1716       /* if the search tree exists, index the new member */

  1717       if (array->tree != NULL)

  1718 avl_set_node_link(avl_insert_node(array->tree, memb->tuple),

  1719             (void *)memb);

  1720       return memb;

  1721 }

  1723 /*----------------------------------------------------------------------

  1724 -- delete_array - delete array.

  1725 --

  1726 -- This routine deletes specified array.

  1727 --

  1728 -- Generic values assigned to the array members are not deleted by this

  1729 -- routine. The calling program itself must delete all assigned generic

  1730 -- values before deleting the array. */

  1732 void delete_array

  1733 (     MPL *mpl,

  1734       ARRAY *array            /* destroyed */

  1735 )

  1736 {     MEMBER *memb;

  1737       xassert(array != NULL);

  1738       /* delete all existing array members */

  1739       while (array->head != NULL)

  1740       {  memb = array->head;

  1741          array->head = memb->next;

  1742          delete_tuple(mpl, memb->tuple);

  1743          dmp_free_atom(mpl->members, memb, sizeof(MEMBER));

  1744       }

  1745       /* if the search tree exists, also delete it */

  1746       if (array->tree != NULL) avl_delete_tree(array->tree);

  1747       /* remove the array from the global array list */

  1748       if (array->prev == NULL)

  1749          mpl->a_list = array->next;

  1750       else

  1751          array->prev->next = array->next;

  1752       if (array->next == NULL)

  1753          ;

  1754       else

  1755          array->next->prev = array->prev;

  1756       /* delete the array descriptor */

  1757       dmp_free_atom(mpl->arrays, array, sizeof(ARRAY));

  1758       return;

  1759 }

  1761 /**********************************************************************/

  1762 /* * *                 DOMAINS AND DUMMY INDICES                  * * */

  1763 /**********************************************************************/

  1765 /*----------------------------------------------------------------------

  1766 -- assign_dummy_index - assign new value to dummy index.

  1767 --

  1768 -- This routine assigns new value to specified dummy index and, that is

  1769 -- important, invalidates all temporary resultant values, which depends

  1770 -- on that dummy index. */

  1772 void assign_dummy_index

  1773 (     MPL *mpl,

  1774       DOMAIN_SLOT *slot,      /* modified */

  1775       SYMBOL *value           /* not changed */

  1776 )

  1777 {     CODE *leaf, *code;

  1778       xassert(slot != NULL);

  1779       xassert(value != NULL);

  1780       /* delete the current value assigned to the dummy index */

  1781       if (slot->value != NULL)

  1782       {  /* if the current value and the new one are identical, actual

  1783             assignment is not needed */

  1784          if (compare_symbols(mpl, slot->value, value) == 0) goto done;

  1785          /* delete a symbol, which is the current value */

  1786          delete_symbol(mpl, slot->value), slot->value = NULL;

  1787       }

  1788       /* now walk through all the pseudo-codes with op = O_INDEX, which

  1789          refer to the dummy index to be changed (these pseudo-codes are

  1790          leaves in the forest of *all* expressions in the database) */

  1791       for (leaf = slot->list; leaf != NULL; leaf = leaf->arg.index.

  1792          next)

  1793       {  xassert(leaf->op == O_INDEX);

  1794          /* invalidate all resultant values, which depend on the dummy

  1795             index, walking from the current leaf toward the root of the

  1796             corresponding expression tree */

  1797          for (code = leaf; code != NULL; code = code->up)

  1798          {  if (code->valid)

  1799             {  /* invalidate and delete resultant value */

  1800                code->valid = 0;

  1801                delete_value(mpl, code->type, &code->value);

  1802             }

  1803          }

  1804       }

  1805       /* assign new value to the dummy index */

  1806       slot->value = copy_symbol(mpl, value);

  1807 done: return;

  1808 }

  1810 /*----------------------------------------------------------------------

  1811 -- update_dummy_indices - update current values of dummy indices.

  1812 --

  1813 -- This routine assigns components of "backup" n-tuple to dummy indices

  1814 -- of specified domain block. If no "backup" n-tuple is defined for the

  1815 -- domain block, values of the dummy indices remain untouched. */

  1817 void update_dummy_indices

  1818 (     MPL *mpl,

  1819       DOMAIN_BLOCK *block     /* not changed */

  1820 )

  1821 {     DOMAIN_SLOT *slot;

  1822       TUPLE *temp;

  1823       if (block->backup != NULL)

  1824       {  for (slot = block->list, temp = block->backup; slot != NULL;

  1825             slot = slot->next, temp = temp->next)

  1826          {  xassert(temp != NULL);

  1827             xassert(temp->sym != NULL);

  1828             assign_dummy_index(mpl, slot, temp->sym);

  1829          }

  1830       }

  1831       return;

  1832 }

  1834 /*----------------------------------------------------------------------

  1835 -- enter_domain_block - enter domain block.

  1836 --

  1837 -- Let specified domain block have the form:

  1838 --

  1839 --    { ..., (j1, j2, ..., jn) in J, ... }

  1840 --

  1841 -- where j1, j2, ..., jn are dummy indices, J is a basic set.

  1842 --

  1843 -- This routine does the following:

  1844 --

  1845 -- 1. Checks if the given n-tuple is a member of the basic set J. Note

  1846 --    that J being *out of the scope* of the domain block cannot depend

  1847 --    on the dummy indices in the same and inner domain blocks, so it

  1848 --    can be computed before the dummy indices are assigned new values.

  1849 --    If this check fails, the routine returns with non-zero code.

  1850 --

  1851 -- 2. Saves current values of the dummy indices j1, j2, ..., jn.

  1852 --

  1853 -- 3. Assigns new values, which are components of the given n-tuple, to

  1854 --    the dummy indices j1, j2, ..., jn. If dimension of the n-tuple is

  1855 --    larger than n, its extra components n+1, n+2, ... are not used.

  1856 --

  1857 -- 4. Calls the formal routine func which either enters the next domain

  1858 --    block or evaluates some code within the domain scope.

  1859 --

  1860 -- 5. Restores former values of the dummy indices j1, j2, ..., jn.

  1861 --

  1862 -- Since current values assigned to the dummy indices on entry to this

  1863 -- routine are restored on exit, the formal routine func is allowed to

  1864 -- call this routine recursively. */

  1866 int enter_domain_block

  1867 (     MPL *mpl,

  1868       DOMAIN_BLOCK *block,    /* not changed */

  1869       TUPLE *tuple,           /* not changed */

  1870       void *info, void (*func)(MPL *mpl, void *info)

  1871 )

  1872 {     TUPLE *backup;

  1873       int ret = 0;

  1874       /* check if the given n-tuple is a member of the basic set */

  1875       xassert(block->code != NULL);

  1876       if (!is_member(mpl, block->code, tuple))

  1877       {  ret = 1;

  1878          goto done;

  1879       }

  1880       /* save reference to "backup" n-tuple, which was used to assign

  1881          current values of the dummy indices (it is sufficient to save

  1882          reference, not value, because that n-tuple is defined in some

  1883          outer level of recursion and therefore cannot be changed on

  1884          this and deeper recursive calls) */

  1885       backup = block->backup;

  1886       /* set up new "backup" n-tuple, which defines new values of the

  1887          dummy indices */

  1888       block->backup = tuple;

  1889       /* assign new values to the dummy indices */

  1890       update_dummy_indices(mpl, block);

  1891       /* call the formal routine that does the rest part of the job */

  1892       func(mpl, info);

  1893       /* restore reference to the former "backup" n-tuple */

  1894       block->backup = backup;

  1895       /* restore former values of the dummy indices; note that if the

  1896          domain block just escaped has no other active instances which

  1897          may exist due to recursion (it is indicated by a null pointer

  1898          to the former n-tuple), former values of the dummy indices are

  1899          undefined; therefore in this case the routine keeps currently

  1900          assigned values of the dummy indices that involves keeping all

  1901          dependent temporary results and thereby, if this domain block

  1902          is not used recursively, allows improving efficiency */

  1903       update_dummy_indices(mpl, block);

  1904 done: return ret;

  1905 }

  1907 /*----------------------------------------------------------------------

  1908 -- eval_within_domain - perform evaluation within domain scope.

  1909 --

  1910 -- This routine assigns new values (symbols) to all dummy indices of

  1911 -- specified domain and calls the formal routine func, which is used to

  1912 -- evaluate some code in the domain scope. Each free dummy index in the

  1913 -- domain is assigned a value specified in the corresponding component

  1914 -- of given n-tuple. Non-free dummy indices are assigned values, which

  1915 -- are computed by this routine.

  1916 --

  1917 -- Number of components in the given n-tuple must be the same as number

  1918 -- of free indices in the domain.

  1919 --

  1920 -- If the given n-tuple is not a member of the domain set, the routine

  1921 -- func is not called, and non-zero code is returned.

  1922 --

  1923 -- For the sake of convenience it is allowed to specify domain as NULL

  1924 -- (then n-tuple also must be 0-tuple, i.e. empty), in which case this

  1925 -- routine just calls the routine func and returns zero.

  1926 --

  1927 -- This routine allows recursive calls from the routine func providing

  1928 -- correct values of dummy indices for each instance.

  1929 --

  1930 -- NOTE: The n-tuple passed to this routine must not be changed by any

  1931 --       other routines called from the formal routine func until this

  1932 --       routine has returned. */

  1934 struct eval_domain_info

  1935 {     /* working info used by the routine eval_within_domain */

  1936       DOMAIN *domain;

  1937       /* domain, which has to be entered */

  1938       DOMAIN_BLOCK *block;

  1939       /* domain block, which is currently processed */

  1940       TUPLE *tuple;

  1941       /* tail of original n-tuple, whose components have to be assigned

  1942          to free dummy indices in the current domain block */

  1943       void *info;

  1944       /* transit pointer passed to the formal routine func */

  1945       void (*func)(MPL *mpl, void *info);

  1946       /* routine, which has to be executed in the domain scope */

  1947       int failure;

  1948       /* this flag indicates that given n-tuple is not a member of the

  1949          domain set */

  1950 };

  1952 static void eval_domain_func(MPL *mpl, void *_my_info)

  1953 {     /* this routine recursively enters into the domain scope and then

  1954          calls the routine func */

  1955       struct eval_domain_info *my_info = _my_info;

  1956       if (my_info->block != NULL)

  1957       {  /* the current domain block to be entered exists */

  1958          DOMAIN_BLOCK *block;

  1959          DOMAIN_SLOT *slot;

  1960          TUPLE *tuple = NULL, *temp = NULL;

  1961          /* save pointer to the current domain block */

  1962          block = my_info->block;

  1963          /* and get ready to enter the next block (if it exists) */

  1964          my_info->block = block->next;

  1965          /* construct temporary n-tuple, whose components correspond to

  1966             dummy indices (slots) of the current domain; components of

  1967             the temporary n-tuple that correspond to free dummy indices

  1968             are assigned references (not values!) to symbols specified

  1969             in the corresponding components of the given n-tuple, while

  1970             other components that correspond to non-free dummy indices

  1971             are assigned symbolic values computed here */

  1972          for (slot = block->list; slot != NULL; slot = slot->next)

  1973          {  /* create component that corresponds to the current slot */

  1974             if (tuple == NULL)

  1975                tuple = temp = dmp_get_atom(mpl->tuples, sizeof(TUPLE));

  1976             else

  1977 temp = (temp->next = dmp_get_atom(mpl->tuples, sizeof(TUPLE)));

  1978             if (slot->code == NULL)

  1979             {  /* dummy index is free; take reference to symbol, which

  1980                   is specified in the corresponding component of given

  1981                   n-tuple */

  1982                xassert(my_info->tuple != NULL);

  1983                temp->sym = my_info->tuple->sym;

  1984                xassert(temp->sym != NULL);

  1985                my_info->tuple = my_info->tuple->next;

  1986             }

  1987             else

  1988             {  /* dummy index is non-free; compute symbolic value to be

  1989                   temporarily assigned to the dummy index */

  1990                temp->sym = eval_symbolic(mpl, slot->code);

  1991             }

  1992          }

  1993          temp->next = NULL;

  1994          /* enter the current domain block */

  1995          if (enter_domain_block(mpl, block, tuple, my_info,

  1996                eval_domain_func)) my_info->failure = 1;

  1997          /* delete temporary n-tuple as well as symbols that correspond

  1998             to non-free dummy indices (they were computed here) */

  1999          for (slot = block->list; slot != NULL; slot = slot->next)

  2000          {  xassert(tuple != NULL);

  2001             temp = tuple;

  2002             tuple = tuple->next;

  2003             if (slot->code != NULL)

  2004             {  /* dummy index is non-free; delete symbolic value */

  2005                delete_symbol(mpl, temp->sym);

  2006             }

  2007             /* delete component that corresponds to the current slot */

  2008             dmp_free_atom(mpl->tuples, temp, sizeof(TUPLE));

  2009          }

  2010       }

  2011       else

  2012       {  /* there are no more domain blocks, i.e. we have reached the

  2013             domain scope */

  2014          xassert(my_info->tuple == NULL);

  2015          /* check optional predicate specified for the domain */

  2016          if (my_info->domain->code != NULL && !eval_logical(mpl,

  2017             my_info->domain->code))

  2018          {  /* the predicate is false */

  2019             my_info->failure = 2;

  2020          }

  2021          else

  2022          {  /* the predicate is true; do the job */

  2023             my_info->func(mpl, my_info->info);

  2024          }

  2025       }

  2026       return;

  2027 }

  2029 int eval_within_domain

  2030 (     MPL *mpl,

  2031       DOMAIN *domain,         /* not changed */

  2032       TUPLE *tuple,           /* not changed */

  2033       void *info, void (*func)(MPL *mpl, void *info)

  2034 )

  2035 {     /* this routine performs evaluation within domain scope */

  2036       struct eval_domain_info _my_info, *my_info = &_my_info;

  2037       if (domain == NULL)

  2038       {  xassert(tuple == NULL);

  2039          func(mpl, info);

  2040          my_info->failure = 0;

  2041       }

  2042       else

  2043       {  xassert(tuple != NULL);

  2044          my_info->domain = domain;

  2045          my_info->block = domain->list;

  2046          my_info->tuple = tuple;

  2047          my_info->info = info;

  2048          my_info->func = func;

  2049          my_info->failure = 0;

  2050          /* enter the very first domain block */

  2051          eval_domain_func(mpl, my_info);

  2052       }

  2053       return my_info->failure;

  2054 }

  2056 /*----------------------------------------------------------------------

  2057 -- loop_within_domain - perform iterations within domain scope.

  2058 --

  2059 -- This routine iteratively assigns new values (symbols) to the dummy

  2060 -- indices of specified domain by enumerating all n-tuples, which are

  2061 -- members of the domain set, and for every n-tuple it calls the formal

  2062 -- routine func to evaluate some code within the domain scope.

  2063 --

  2064 -- If the routine func returns non-zero, enumeration within the domain

  2065 -- is prematurely terminated.

  2066 --

  2067 -- For the sake of convenience it is allowed to specify domain as NULL,

  2068 -- in which case this routine just calls the routine func only once and

  2069 -- returns zero.

  2070 --

  2071 -- This routine allows recursive calls from the routine func providing

  2072 -- correct values of dummy indices for each instance. */

  2074 struct loop_domain_info

  2075 {     /* working info used by the routine loop_within_domain */

  2076       DOMAIN *domain;

  2077       /* domain, which has to be entered */

  2078       DOMAIN_BLOCK *block;

  2079       /* domain block, which is currently processed */

  2080       int looping;

  2081       /* clearing this flag leads to terminating enumeration */

  2082       void *info;

  2083       /* transit pointer passed to the formal routine func */

  2084       int (*func)(MPL *mpl, void *info);

  2085       /* routine, which needs to be executed in the domain scope */

  2086 };

  2088 static void loop_domain_func(MPL *mpl, void *_my_info)

  2089 {     /* this routine enumerates all n-tuples in the basic set of the

  2090          current domain block, enters recursively into the domain scope

  2091          for every n-tuple, and then calls the routine func */

  2092       struct loop_domain_info *my_info = _my_info;

  2093       if (my_info->block != NULL)

  2094       {  /* the current domain block to be entered exists */

  2095          DOMAIN_BLOCK *block;

  2096          DOMAIN_SLOT *slot;

  2097          TUPLE *bound;

  2098          /* save pointer to the current domain block */

  2099          block = my_info->block;

  2100          /* and get ready to enter the next block (if it exists) */

  2101          my_info->block = block->next;

  2102          /* compute symbolic values, at which non-free dummy indices of

  2103             the current domain block are bound; since that values don't

  2104             depend on free dummy indices of the current block, they can

  2105             be computed once out of the enumeration loop */

  2106          bound = create_tuple(mpl);

  2107          for (slot = block->list; slot != NULL; slot = slot->next)

  2108          {  if (slot->code != NULL)

  2109                bound = expand_tuple(mpl, bound, eval_symbolic(mpl,

  2110                   slot->code));

  2111          }

  2112          /* start enumeration */

  2113          xassert(block->code != NULL);

  2114          if (block->code->op == O_DOTS)

  2115          {  /* the basic set is "arithmetic", in which case it doesn't

  2116                need to be computed explicitly */

  2117             TUPLE *tuple;

  2118             int n, j;

  2119             double t0, tf, dt;

  2120             /* compute "parameters" of the basic set */

  2121             t0 = eval_numeric(mpl, block->code->arg.arg.x);

  2122             tf = eval_numeric(mpl, block->code->arg.arg.y);

  2123             if (block->code->arg.arg.z == NULL)

  2124                dt = 1.0;

  2125             else

  2126                dt = eval_numeric(mpl, block->code->arg.arg.z);

  2127             /* determine cardinality of the basic set */

  2128             n = arelset_size(mpl, t0, tf, dt);

  2129             /* create dummy 1-tuple for members of the basic set */

  2130             tuple = expand_tuple(mpl, create_tuple(mpl),

  2131                create_symbol_num(mpl, 0.0));

  2132             /* in case of "arithmetic" set there is exactly one dummy

  2133                index, which cannot be non-free */

  2134             xassert(bound == NULL);

  2135             /* walk through 1-tuples of the basic set */

  2136             for (j = 1; j <= n && my_info->looping; j++)

  2137             {  /* construct dummy 1-tuple for the current member */

  2138                tuple->sym->num = arelset_member(mpl, t0, tf, dt, j);

  2139                /* enter the current domain block */

  2140                enter_domain_block(mpl, block, tuple, my_info,

  2141                   loop_domain_func);

  2142             }

  2143             /* delete dummy 1-tuple */

  2144             delete_tuple(mpl, tuple);

  2145          }

  2146          else

  2147          {  /* the basic set is of general kind, in which case it needs

  2148                to be explicitly computed */

  2149             ELEMSET *set;

  2150             MEMBER *memb;

  2151             TUPLE *temp1, *temp2;

  2152             /* compute the basic set */

  2153             set = eval_elemset(mpl, block->code);

  2154             /* walk through all n-tuples of the basic set */

  2155             for (memb = set->head; memb != NULL && my_info->looping;

  2156                memb = memb->next)

  2157             {  /* all components of the current n-tuple that correspond

  2158                   to non-free dummy indices must be feasible; otherwise

  2159                   the n-tuple is not in the basic set */

  2160                temp1 = memb->tuple;

  2161                temp2 = bound;

  2162                for (slot = block->list; slot != NULL; slot = slot->next)

  2163                {  xassert(temp1 != NULL);

  2164                   if (slot->code != NULL)

  2165                   {  /* non-free dummy index */

  2166                      xassert(temp2 != NULL);

  2167                      if (compare_symbols(mpl, temp1->sym, temp2->sym)

  2168                         != 0)

  2169                      {  /* the n-tuple is not in the basic set */

  2170                         goto skip;

  2171                      }

  2172                      temp2 = temp2->next;

  2173                   }

  2174                   temp1 = temp1->next;

  2175                }

  2176                xassert(temp1 == NULL);

  2177                xassert(temp2 == NULL);

  2178                /* enter the current domain block */

  2179                enter_domain_block(mpl, block, memb->tuple, my_info,

  2180                   loop_domain_func);

  2181 skip:          ;

  2182             }

  2183             /* delete the basic set */

  2184             delete_elemset(mpl, set);

  2185          }

  2186          /* delete symbolic values binding non-free dummy indices */

  2187          delete_tuple(mpl, bound);

  2188          /* restore pointer to the current domain block */

  2189          my_info->block = block;

  2190       }

  2191       else

  2192       {  /* there are no more domain blocks, i.e. we have reached the

  2193             domain scope */

  2194          /* check optional predicate specified for the domain */

  2195          if (my_info->domain->code != NULL && !eval_logical(mpl,

  2196             my_info->domain->code))

  2197          {  /* the predicate is false */

  2198             /* nop */;

  2199          }

  2200          else

  2201          {  /* the predicate is true; do the job */

  2202             my_info->looping = !my_info->func(mpl, my_info->info);

  2203          }

  2204       }

  2205       return;

  2206 }

  2208 void loop_within_domain

  2209 (     MPL *mpl,

  2210       DOMAIN *domain,         /* not changed */

  2211       void *info, int (*func)(MPL *mpl, void *info)

  2212 )

  2213 {     /* this routine performs iterations within domain scope */

  2214       struct loop_domain_info _my_info, *my_info = &_my_info;

  2215       if (domain == NULL)

  2216          func(mpl, info);

  2217       else

  2218       {  my_info->domain = domain;

  2219          my_info->block = domain->list;

  2220          my_info->looping = 1;

  2221          my_info->info = info;

  2222          my_info->func = func;

  2223          /* enter the very first domain block */

  2224          loop_domain_func(mpl, my_info);

  2225       }

  2226       return;

  2227 }

  2229 /*----------------------------------------------------------------------

  2230 -- out_of_domain - raise domain exception.

  2231 --

  2232 -- This routine is called when a reference is made to a member of some

  2233 -- model object, but its n-tuple is out of the object domain. */

  2235 void out_of_domain

  2236 (     MPL *mpl,

  2237       char *name,             /* not changed */

  2238       TUPLE *tuple            /* not changed */

  2239 )

  2240 {     xassert(name != NULL);

  2241       xassert(tuple != NULL);

  2242       error(mpl, "%s%s out of domain", name, format_tuple(mpl, '[',

  2243          tuple));

  2244       /* no return */

  2245 }

  2247 /*----------------------------------------------------------------------

  2248 -- get_domain_tuple - obtain current n-tuple from domain.

  2249 --

  2250 -- This routine constructs n-tuple, whose components are current values

  2251 -- assigned to *free* dummy indices of specified domain.

  2252 --

  2253 -- For the sake of convenience it is allowed to specify domain as NULL,

  2254 -- in which case this routine returns 0-tuple.

  2255 --

  2256 -- NOTE: This routine must not be called out of domain scope. */

  2258 TUPLE *get_domain_tuple

  2259 (     MPL *mpl,

  2260       DOMAIN *domain          /* not changed */

  2261 )

  2262 {     DOMAIN_BLOCK *block;

  2263       DOMAIN_SLOT *slot;

  2264       TUPLE *tuple;

  2265       tuple = create_tuple(mpl);

  2266       if (domain != NULL)

  2267       {  for (block = domain->list; block != NULL; block = block->next)

  2268          {  for (slot = block->list; slot != NULL; slot = slot->next)

  2269             {  if (slot->code == NULL)

  2270                {  xassert(slot->value != NULL);

  2271                   tuple = expand_tuple(mpl, tuple, copy_symbol(mpl,

  2272                      slot->value));

  2273                }

  2274             }

  2275          }

  2276       }

  2277       return tuple;

  2278 }

  2280 /*----------------------------------------------------------------------

  2281 -- clean_domain - clean domain.

  2282 --

  2283 -- This routine cleans specified domain that assumes deleting all stuff

  2284 -- dynamically allocated during the generation phase. */

  2286 void clean_domain(MPL *mpl, DOMAIN *domain)

  2287 {     DOMAIN_BLOCK *block;

  2288       DOMAIN_SLOT *slot;

  2289       /* if no domain is specified, do nothing */

  2290       if (domain == NULL) goto done;

  2291       /* clean all domain blocks */

  2292       for (block = domain->list; block != NULL; block = block->next)

  2293       {  /* clean all domain slots */

  2294          for (slot = block->list; slot != NULL; slot = slot->next)

  2295          {  /* clean pseudo-code for computing bound value */

  2296             clean_code(mpl, slot->code);

  2297             /* delete symbolic value assigned to dummy index */

  2298             if (slot->value != NULL)

  2299                delete_symbol(mpl, slot->value), slot->value = NULL;

  2300          }

  2301          /* clean pseudo-code for computing basic set */

  2302          clean_code(mpl, block->code);

  2303       }

  2304       /* clean pseudo-code for computing domain predicate */

  2305       clean_code(mpl, domain->code);

  2306 done: return;

  2307 }

  2309 /**********************************************************************/

  2310 /* * *                         MODEL SETS                         * * */

  2311 /**********************************************************************/

  2313 /*----------------------------------------------------------------------

  2314 -- check_elem_set - check elemental set assigned to set member.

  2315 --

  2316 -- This routine checks if given elemental set being assigned to member

  2317 -- of specified model set satisfies to all restrictions.

  2318 --

  2319 -- NOTE: This routine must not be called out of domain scope. */

  2321 void check_elem_set

  2322 (     MPL *mpl,

  2323       SET *set,               /* not changed */

  2324       TUPLE *tuple,           /* not changed */

  2325       ELEMSET *refer          /* not changed */

  2326 )

  2327 {     WITHIN *within;

  2328       MEMBER *memb;

  2329       int eqno;

  2330       /* elemental set must be within all specified supersets */

  2331       for (within = set->within, eqno = 1; within != NULL; within =

  2332          within->next, eqno++)

  2333       {  xassert(within->code != NULL);

  2334          for (memb = refer->head; memb != NULL; memb = memb->next)

  2335          {  if (!is_member(mpl, within->code, memb->tuple))

  2336             {  char buf[255+1];

  2337                strcpy(buf, format_tuple(mpl, '(', memb->tuple));

  2338                xassert(strlen(buf) < sizeof(buf));

  2339                error(mpl, "%s%s contains %s which not within specified "

  2340                   "set; see (%d)", set->name, format_tuple(mpl, '[',

  2341                      tuple), buf, eqno);

  2342             }

  2343          }

  2344       }

  2345       return;

  2346 }

  2348 /*----------------------------------------------------------------------

  2349 -- take_member_set - obtain elemental set assigned to set member.

  2350 --

  2351 -- This routine obtains a reference to elemental set assigned to given

  2352 -- member of specified model set and returns it on exit.

  2353 --

  2354 -- NOTE: This routine must not be called out of domain scope. */

  2356 ELEMSET *take_member_set      /* returns reference, not value */

  2357 (     MPL *mpl,

  2358       SET *set,               /* not changed */

  2359       TUPLE *tuple            /* not changed */

  2360 )

  2361 {     MEMBER *memb;

  2362       ELEMSET *refer;

  2363       /* find member in the set array */

  2364       memb = find_member(mpl, set->array, tuple);

  2365       if (memb != NULL)

  2366       {  /* member exists, so just take the reference */

  2367          refer = memb->value.set;

  2368       }

  2369       else if (set->assign != NULL)

  2370       {  /* compute value using assignment expression */

  2371          refer = eval_elemset(mpl, set->assign);

  2372 add:     /* check that the elemental set satisfies to all restrictions,

  2373             assign it to new member, and add the member to the array */

  2374          check_elem_set(mpl, set, tuple, refer);

  2375          memb = add_member(mpl, set->array, copy_tuple(mpl, tuple));

  2376          memb->value.set = refer;

  2377       }

  2378       else if (set->option != NULL)

  2379       {  /* compute default elemental set */

  2380          refer = eval_elemset(mpl, set->option);

  2381          goto add;

  2382       }

  2383       else

  2384       {  /* no value (elemental set) is provided */

  2385          error(mpl, "no value for %s%s", set->name, format_tuple(mpl,

  2386             '[', tuple));

  2387       }

  2388       return refer;

  2389 }

  2391 /*----------------------------------------------------------------------

  2392 -- eval_member_set - evaluate elemental set assigned to set member.

  2393 --

  2394 -- This routine evaluates a reference to elemental set assigned to given

  2395 -- member of specified model set and returns it on exit. */

  2397 struct eval_set_info

  2398 {     /* working info used by the routine eval_member_set */

  2399       SET *set;

  2400       /* model set */

  2401       TUPLE *tuple;

  2402       /* n-tuple, which defines set member */

  2403       MEMBER *memb;

  2404       /* normally this pointer is NULL; the routine uses this pointer

  2405          to check data provided in the data section, in which case it

  2406          points to a member currently checked; this check is performed

  2407          automatically only once when a reference to any member occurs

  2408          for the first time */

  2409       ELEMSET *refer;

  2410       /* evaluated reference to elemental set */

  2411 };

  2413 static void eval_set_func(MPL *mpl, void *_info)

  2414 {     /* this is auxiliary routine to work within domain scope */

  2415       struct eval_set_info *info = _info;

  2416       if (info->memb != NULL)

  2417       {  /* checking call; check elemental set being assigned */

  2418          check_elem_set(mpl, info->set, info->memb->tuple,

  2419             info->memb->value.set);

  2420       }

  2421       else

  2422       {  /* normal call; evaluate member, which has given n-tuple */

  2423          info->refer = take_member_set(mpl, info->set, info->tuple);

  2424       }

  2425       return;

  2426 }

  2428 #if 1 /* 12/XII-2008 */

  2429 static void saturate_set(MPL *mpl, SET *set)

  2430 {     GADGET *gadget = set->gadget;

  2431       ELEMSET *data;

  2432       MEMBER *elem, *memb;

  2433       TUPLE *tuple, *work[20];

  2434       int i;

  2435       xprintf("Generating %s...\n", set->name);

  2436       eval_whole_set(mpl, gadget->set);

  2437       /* gadget set must have exactly one member */

  2438       xassert(gadget->set->array != NULL);

  2439       xassert(gadget->set->array->head != NULL);

  2440       xassert(gadget->set->array->head == gadget->set->array->tail);

  2441       data = gadget->set->array->head->value.set;

  2442       xassert(data->type == A_NONE);

  2443       xassert(data->dim == gadget->set->dimen);

  2444       /* walk thru all elements of the plain set */

  2445       for (elem = data->head; elem != NULL; elem = elem->next)

  2446       {  /* create a copy of n-tuple */

  2447          tuple = copy_tuple(mpl, elem->tuple);

  2448          /* rearrange component of the n-tuple */

  2449          for (i = 0; i < gadget->set->dimen; i++)

  2450             work[i] = NULL;

  2451          for (i = 0; tuple != NULL; tuple = tuple->next)

  2452             work[gadget->ind[i++]-1] = tuple;

  2453          xassert(i == gadget->set->dimen);

  2454          for (i = 0; i < gadget->set->dimen; i++)

  2455          {  xassert(work[i] != NULL);

  2456             work[i]->next = work[i+1];

  2457          }

  2458          /* construct subscript list from first set->dim components */

  2459          if (set->dim == 0)

  2460             tuple = NULL;

  2461          else

  2462             tuple = work[0], work[set->dim-1]->next = NULL;

  2463          /* find corresponding member of the set to be initialized */

  2464          memb = find_member(mpl, set->array, tuple);

  2465          if (memb == NULL)

  2466          {  /* not found; add new member to the set and assign it empty

  2467                elemental set */

  2468             memb = add_member(mpl, set->array, tuple);

  2469             memb->value.set = create_elemset(mpl, set->dimen);

  2470          }

  2471          else

  2472          {  /* found; free subscript list */

  2473             delete_tuple(mpl, tuple);

  2474          }

  2475          /* construct new n-tuple from rest set->dimen components */

  2476          tuple = work[set->dim];

  2477          xassert(set->dim + set->dimen == gadget->set->dimen);

  2478          work[gadget->set->dimen-1]->next = NULL;

  2479          /* and add it to the elemental set assigned to the member

  2480             (no check for duplicates is needed) */

  2481          add_tuple(mpl, memb->value.set, tuple);

  2482       }

  2483       /* the set has been saturated with data */

  2484       set->data = 1;

  2485       return;

  2486 }

  2487 #endif

  2489 ELEMSET *eval_member_set      /* returns reference, not value */

  2490 (     MPL *mpl,

  2491       SET *set,               /* not changed */

  2492       TUPLE *tuple            /* not changed */

  2493 )

  2494 {     /* this routine evaluates set member */

  2495       struct eval_set_info _info, *info = &_info;

  2496       xassert(set->dim == tuple_dimen(mpl, tuple));

  2497       info->set = set;

  2498       info->tuple = tuple;

  2499 #if 1 /* 12/XII-2008 */

  2500       if (set->gadget != NULL && set->data == 0)

  2501       {  /* initialize the set with data from a plain set */

  2502          saturate_set(mpl, set);

  2503       }

  2504 #endif

  2505       if (set->data == 1)

  2506       {  /* check data, which are provided in the data section, but not

  2507             checked yet */

  2508          /* save pointer to the last array member; note that during the

  2509             check new members may be added beyond the last member due to

  2510             references to the same parameter from default expression as

  2511             well as from expressions that define restricting supersets;

  2512             however, values assigned to the new members will be checked

  2513             by other routine, so we don't need to check them here */

  2514          MEMBER *tail = set->array->tail;

  2515          /* change the data status to prevent infinite recursive loop

  2516             due to references to the same set during the check */

  2517          set->data = 2;

  2518          /* check elemental sets assigned to array members in the data

  2519             section until the marked member has been reached */

  2520          for (info->memb = set->array->head; info->memb != NULL;

  2521             info->memb = info->memb->next)

  2522          {  if (eval_within_domain(mpl, set->domain, info->memb->tuple,

  2523                info, eval_set_func))

  2524                out_of_domain(mpl, set->name, info->memb->tuple);

  2525             if (info->memb == tail) break;

  2526          }

  2527          /* the check has been finished */

  2528       }

  2529       /* evaluate member, which has given n-tuple */

  2530       info->memb = NULL;

  2531       if (eval_within_domain(mpl, info->set->domain, info->tuple, info,

  2532          eval_set_func))

  2533       out_of_domain(mpl, set->name, info->tuple);

  2534       /* bring evaluated reference to the calling program */

  2535       return info->refer;

  2536 }

  2538 /*----------------------------------------------------------------------

  2539 -- eval_whole_set - evaluate model set over entire domain.

  2540 --

  2541 -- This routine evaluates all members of specified model set over entire

  2542 -- domain. */

  2544 static int whole_set_func(MPL *mpl, void *info)

  2545 {     /* this is auxiliary routine to work within domain scope */

  2546       SET *set = (SET *)info;

  2547       TUPLE *tuple = get_domain_tuple(mpl, set->domain);

  2548       eval_member_set(mpl, set, tuple);

  2549       delete_tuple(mpl, tuple);

  2550       return 0;

  2551 }

  2553 void eval_whole_set(MPL *mpl, SET *set)

  2554 {     loop_within_domain(mpl, set->domain, set, whole_set_func);

  2555       return;

  2556 }

  2558 /*----------------------------------------------------------------------

  2559 -- clean set - clean model set.

  2560 --

  2561 -- This routine cleans specified model set that assumes deleting all

  2562 -- stuff dynamically allocated during the generation phase. */

  2564 void clean_set(MPL *mpl, SET *set)

  2565 {     WITHIN *within;

  2566       MEMBER *memb;

  2567       /* clean subscript domain */

  2568       clean_domain(mpl, set->domain);

  2569       /* clean pseudo-code for computing supersets */

  2570       for (within = set->within; within != NULL; within = within->next)

  2571          clean_code(mpl, within->code);

  2572       /* clean pseudo-code for computing assigned value */

  2573       clean_code(mpl, set->assign);

  2574       /* clean pseudo-code for computing default value */

  2575       clean_code(mpl, set->option);

  2576       /* reset data status flag */

  2577       set->data = 0;

  2578       /* delete content array */

  2579       for (memb = set->array->head; memb != NULL; memb = memb->next)

  2580          delete_value(mpl, set->array->type, &memb->value);

  2581       delete_array(mpl, set->array), set->array = NULL;

  2582       return;

  2583 }

  2585 /**********************************************************************/

  2586 /* * *                      MODEL PARAMETERS                      * * */

  2587 /**********************************************************************/

  2589 /*----------------------------------------------------------------------

  2590 -- check_value_num - check numeric value assigned to parameter member.

  2591 --

  2592 -- This routine checks if numeric value being assigned to some member

  2593 -- of specified numeric model parameter satisfies to all restrictions.

  2594 --

  2595 -- NOTE: This routine must not be called out of domain scope. */

  2597 void check_value_num

  2598 (     MPL *mpl,

  2599       PARAMETER *par,         /* not changed */

  2600       TUPLE *tuple,           /* not changed */

  2601       double value

  2602 )

  2603 {     CONDITION *cond;

  2604       WITHIN *in;

  2605       int eqno;

  2606       /* the value must satisfy to the parameter type */

  2607       switch (par->type)

  2608       {  case A_NUMERIC:

  2609             break;

  2610          case A_INTEGER:

  2611             if (value != floor(value))

  2612                error(mpl, "%s%s = %.*g not integer", par->name,

  2613                   format_tuple(mpl, '[', tuple), DBL_DIG, value);

  2614             break;

  2615          case A_BINARY:

  2616             if (!(value == 0.0 || value == 1.0))

  2617                error(mpl, "%s%s = %.*g not binary", par->name,

  2618                   format_tuple(mpl, '[', tuple), DBL_DIG, value);

  2619             break;

  2620          default:

  2621             xassert(par != par);

  2622       }

  2623       /* the value must satisfy to all specified conditions */

  2624       for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,

  2625          eqno++)

  2626       {  double bound;

  2627          char *rho;

  2628          xassert(cond->code != NULL);

  2629          bound = eval_numeric(mpl, cond->code);

  2630          switch (cond->rho)

  2631          {  case O_LT:

  2632                if (!(value < bound))

  2633                {  rho = "<";

  2634 err:              error(mpl, "%s%s = %.*g not %s %.*g; see (%d)",

  2635                      par->name, format_tuple(mpl, '[', tuple), DBL_DIG,

  2636                      value, rho, DBL_DIG, bound, eqno);

  2637                }

  2638                break;

  2639             case O_LE:

  2640                if (!(value <= bound)) { rho = "<="; goto err; }

  2641                break;

  2642             case O_EQ:

  2643                if (!(value == bound)) { rho = "="; goto err; }

  2644                break;

  2645             case O_GE:

  2646                if (!(value >= bound)) { rho = ">="; goto err; }

  2647                break;

  2648             case O_GT:

  2649                if (!(value > bound)) { rho = ">"; goto err; }

  2650                break;

  2651             case O_NE:

  2652                if (!(value != bound)) { rho = "<>"; goto err; }

  2653                break;

  2654             default:

  2655                xassert(cond != cond);

  2656          }

  2657       }

  2658       /* the value must be in all specified supersets */

  2659       for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)

  2660       {  TUPLE *dummy;

  2661          xassert(in->code != NULL);

  2662          xassert(in->code->dim == 1);

  2663          dummy = expand_tuple(mpl, create_tuple(mpl),

  2664             create_symbol_num(mpl, value));

  2665          if (!is_member(mpl, in->code, dummy))

  2666             error(mpl, "%s%s = %.*g not in specified set; see (%d)",

  2667                par->name, format_tuple(mpl, '[', tuple), DBL_DIG,

  2668                value, eqno);

  2669          delete_tuple(mpl, dummy);

  2670       }

  2671       return;

  2672 }

  2674 /*----------------------------------------------------------------------

  2675 -- take_member_num - obtain num. value assigned to parameter member.

  2676 --

  2677 -- This routine obtains a numeric value assigned to member of specified

  2678 -- numeric model parameter and returns it on exit.

  2679 --

  2680 -- NOTE: This routine must not be called out of domain scope. */

  2682 double take_member_num

  2683 (     MPL *mpl,

  2684       PARAMETER *par,         /* not changed */

  2685       TUPLE *tuple            /* not changed */

  2686 )

  2687 {     MEMBER *memb;

  2688       double value;

  2689       /* find member in the parameter array */

  2690       memb = find_member(mpl, par->array, tuple);

  2691       if (memb != NULL)

  2692       {  /* member exists, so just take its value */

  2693          value = memb->value.num;

  2694       }

  2695       else if (par->assign != NULL)

  2696       {  /* compute value using assignment expression */

  2697          value = eval_numeric(mpl, par->assign);

  2698 add:     /* check that the value satisfies to all restrictions, assign

  2699             it to new member, and add the member to the array */

  2700          check_value_num(mpl, par, tuple, value);

  2701          memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));

  2702          memb->value.num = value;

  2703       }

  2704       else if (par->option != NULL)

  2705       {  /* compute default value */

  2706          value = eval_numeric(mpl, par->option);

  2707          goto add;

  2708       }

  2709       else if (par->defval != NULL)

  2710       {  /* take default value provided in the data section */

  2711          if (par->defval->str != NULL)

  2712             error(mpl, "cannot convert %s to floating-point number",

  2713                format_symbol(mpl, par->defval));

  2714          value = par->defval->num;

  2715          goto add;

  2716       }

  2717       else

  2718       {  /* no value is provided */

  2719          error(mpl, "no value for %s%s", par->name, format_tuple(mpl,

  2720             '[', tuple));

  2721       }

  2722       return value;

  2723 }

  2725 /*----------------------------------------------------------------------

  2726 -- eval_member_num - evaluate num. value assigned to parameter member.

  2727 --

  2728 -- This routine evaluates a numeric value assigned to given member of

  2729 -- specified numeric model parameter and returns it on exit. */

  2731 struct eval_num_info

  2732 {     /* working info used by the routine eval_member_num */

  2733       PARAMETER *par;

  2734       /* model parameter  */

  2735       TUPLE *tuple;

  2736       /* n-tuple, which defines parameter member */

  2737       MEMBER *memb;

  2738       /* normally this pointer is NULL; the routine uses this pointer

  2739          to check data provided in the data section, in which case it

  2740          points to a member currently checked; this check is performed

  2741          automatically only once when a reference to any member occurs

  2742          for the first time */

  2743       double value;

  2744       /* evaluated numeric value */

  2745 };

  2747 static void eval_num_func(MPL *mpl, void *_info)

  2748 {     /* this is auxiliary routine to work within domain scope */

  2749       struct eval_num_info *info = _info;

  2750       if (info->memb != NULL)

  2751       {  /* checking call; check numeric value being assigned */

  2752          check_value_num(mpl, info->par, info->memb->tuple,

  2753             info->memb->value.num);

  2754       }

  2755       else

  2756       {  /* normal call; evaluate member, which has given n-tuple */

  2757          info->value = take_member_num(mpl, info->par, info->tuple);

  2758       }

  2759       return;

  2760 }

  2762 double eval_member_num

  2763 (     MPL *mpl,

  2764       PARAMETER *par,         /* not changed */

  2765       TUPLE *tuple            /* not changed */

  2766 )

  2767 {     /* this routine evaluates numeric parameter member */

  2768       struct eval_num_info _info, *info = &_info;

  2769       xassert(par->type == A_NUMERIC || par->type == A_INTEGER ||

  2770              par->type == A_BINARY);

  2771       xassert(par->dim == tuple_dimen(mpl, tuple));

  2772       info->par = par;

  2773       info->tuple = tuple;

  2774       if (par->data == 1)

  2775       {  /* check data, which are provided in the data section, but not

  2776             checked yet */

  2777          /* save pointer to the last array member; note that during the

  2778             check new members may be added beyond the last member due to

  2779             references to the same parameter from default expression as

  2780             well as from expressions that define restricting conditions;

  2781             however, values assigned to the new members will be checked

  2782             by other routine, so we don't need to check them here */

  2783          MEMBER *tail = par->array->tail;

  2784          /* change the data status to prevent infinite recursive loop

  2785             due to references to the same parameter during the check */

  2786          par->data = 2;

  2787          /* check values assigned to array members in the data section

  2788             until the marked member has been reached */

  2789          for (info->memb = par->array->head; info->memb != NULL;

  2790             info->memb = info->memb->next)

  2791          {  if (eval_within_domain(mpl, par->domain, info->memb->tuple,

  2792                info, eval_num_func))

  2793                out_of_domain(mpl, par->name, info->memb->tuple);

  2794             if (info->memb == tail) break;

  2795          }

  2796          /* the check has been finished */

  2797       }

  2798       /* evaluate member, which has given n-tuple */

  2799       info->memb = NULL;

  2800       if (eval_within_domain(mpl, info->par->domain, info->tuple, info,

  2801          eval_num_func))

  2802          out_of_domain(mpl, par->name, info->tuple);

  2803       /* bring evaluated value to the calling program */

  2804       return info->value;

  2805 }

  2807 /*----------------------------------------------------------------------

  2808 -- check_value_sym - check symbolic value assigned to parameter member.

  2809 --

  2810 -- This routine checks if symbolic value being assigned to some member

  2811 -- of specified symbolic model parameter satisfies to all restrictions.

  2812 --

  2813 -- NOTE: This routine must not be called out of domain scope. */

  2815 void check_value_sym

  2816 (     MPL *mpl,

  2817       PARAMETER *par,         /* not changed */

  2818       TUPLE *tuple,           /* not changed */

  2819       SYMBOL *value           /* not changed */

  2820 )

  2821 {     CONDITION *cond;

  2822       WITHIN *in;

  2823       int eqno;

  2824       /* the value must satisfy to all specified conditions */

  2825       for (cond = par->cond, eqno = 1; cond != NULL; cond = cond->next,

  2826          eqno++)

  2827       {  SYMBOL *bound;

  2828          char buf[255+1];

  2829          xassert(cond->code != NULL);

  2830          bound = eval_symbolic(mpl, cond->code);

  2831          switch (cond->rho)

  2832          {

  2833 #if 1 /* 13/VIII-2008 */

  2834             case O_LT:

  2835                if (!(compare_symbols(mpl, value, bound) < 0))

  2836                {  strcpy(buf, format_symbol(mpl, bound));

  2837                   xassert(strlen(buf) < sizeof(buf));

  2838                   error(mpl, "%s%s = %s not < %s",

  2839                      par->name, format_tuple(mpl, '[', tuple),

  2840                      format_symbol(mpl, value), buf, eqno);

  2841                }

  2842                break;

  2843             case O_LE:

  2844                if (!(compare_symbols(mpl, value, bound) <= 0))

  2845                {  strcpy(buf, format_symbol(mpl, bound));

  2846                   xassert(strlen(buf) < sizeof(buf));

  2847                   error(mpl, "%s%s = %s not <= %s",

  2848                      par->name, format_tuple(mpl, '[', tuple),

  2849                      format_symbol(mpl, value), buf, eqno);

  2850                }

  2851                break;

  2852 #endif

  2853             case O_EQ:

  2854                if (!(compare_symbols(mpl, value, bound) == 0))

  2855                {  strcpy(buf, format_symbol(mpl, bound));

  2856                   xassert(strlen(buf) < sizeof(buf));

  2857                   error(mpl, "%s%s = %s not = %s",

  2858                      par->name, format_tuple(mpl, '[', tuple),

  2859                      format_symbol(mpl, value), buf, eqno);

  2860                }

  2861                break;

  2862 #if 1 /* 13/VIII-2008 */

  2863             case O_GE:

  2864                if (!(compare_symbols(mpl, value, bound) >= 0))

  2865                {  strcpy(buf, format_symbol(mpl, bound));

  2866                   xassert(strlen(buf) < sizeof(buf));

  2867                   error(mpl, "%s%s = %s not >= %s",

  2868                      par->name, format_tuple(mpl, '[', tuple),

  2869                      format_symbol(mpl, value), buf, eqno);

  2870                }

  2871                break;

  2872             case O_GT:

  2873                if (!(compare_symbols(mpl, value, bound) > 0))

  2874                {  strcpy(buf, format_symbol(mpl, bound));

  2875                   xassert(strlen(buf) < sizeof(buf));

  2876                   error(mpl, "%s%s = %s not > %s",

  2877                      par->name, format_tuple(mpl, '[', tuple),

  2878                      format_symbol(mpl, value), buf, eqno);

  2879                }

  2880                break;

  2881 #endif

  2882             case O_NE:

  2883                if (!(compare_symbols(mpl, value, bound) != 0))

  2884                {  strcpy(buf, format_symbol(mpl, bound));

  2885                   xassert(strlen(buf) < sizeof(buf));

  2886                   error(mpl, "%s%s = %s not <> %s",

  2887                      par->name, format_tuple(mpl, '[', tuple),

  2888                      format_symbol(mpl, value), buf, eqno);

  2889                }

  2890                break;

  2891             default:

  2892                xassert(cond != cond);

  2893          }

  2894          delete_symbol(mpl, bound);

  2895       }

  2896       /* the value must be in all specified supersets */

  2897       for (in = par->in, eqno = 1; in != NULL; in = in->next, eqno++)

  2898       {  TUPLE *dummy;

  2899          xassert(in->code != NULL);

  2900          xassert(in->code->dim == 1);

  2901          dummy = expand_tuple(mpl, create_tuple(mpl), copy_symbol(mpl,

  2902             value));

  2903          if (!is_member(mpl, in->code, dummy))

  2904             error(mpl, "%s%s = %s not in specified set; see (%d)",

  2905                par->name, format_tuple(mpl, '[', tuple),

  2906                format_symbol(mpl, value), eqno);

  2907          delete_tuple(mpl, dummy);

  2908       }

  2909       return;

  2910 }

  2912 /*----------------------------------------------------------------------

  2913 -- take_member_sym - obtain symb. value assigned to parameter member.

  2914 --

  2915 -- This routine obtains a symbolic value assigned to member of specified

  2916 -- symbolic model parameter and returns it on exit.

  2917 --

  2918 -- NOTE: This routine must not be called out of domain scope. */

  2920 SYMBOL *take_member_sym       /* returns value, not reference */

  2921 (     MPL *mpl,

  2922       PARAMETER *par,         /* not changed */

  2923       TUPLE *tuple            /* not changed */

  2924 )

  2925 {     MEMBER *memb;

  2926       SYMBOL *value;

  2927       /* find member in the parameter array */

  2928       memb = find_member(mpl, par->array, tuple);

  2929       if (memb != NULL)

  2930       {  /* member exists, so just take its value */

  2931          value = copy_symbol(mpl, memb->value.sym);

  2932       }

  2933       else if (par->assign != NULL)

  2934       {  /* compute value using assignment expression */

  2935          value = eval_symbolic(mpl, par->assign);

  2936 add:     /* check that the value satisfies to all restrictions, assign

  2937             it to new member, and add the member to the array */

  2938          check_value_sym(mpl, par, tuple, value);

  2939          memb = add_member(mpl, par->array, copy_tuple(mpl, tuple));

  2940          memb->value.sym = copy_symbol(mpl, value);

  2941       }

  2942       else if (par->option != NULL)

  2943       {  /* compute default value */

  2944          value = eval_symbolic(mpl, par->option);

  2945          goto add;

  2946       }

  2947       else if (par->defval != NULL)

  2948       {  /* take default value provided in the data section */

  2949          value = copy_symbol(mpl, par->defval);

  2950          goto add;

  2951       }

  2952       else

  2953       {  /* no value is provided */

  2954          error(mpl, "no value for %s%s", par->name, format_tuple(mpl,

  2955             '[', tuple));

  2956       }

  2957       return value;

  2958 }

  2960 /*----------------------------------------------------------------------

  2961 -- eval_member_sym - evaluate symb. value assigned to parameter member.

  2962 --

  2963 -- This routine evaluates a symbolic value assigned to given member of

  2964 -- specified symbolic model parameter and returns it on exit. */

  2966 struct eval_sym_info

  2967 {     /* working info used by the routine eval_member_sym */

  2968       PARAMETER *par;

  2969       /* model parameter */

  2970       TUPLE *tuple;

  2971       /* n-tuple, which defines parameter member */

  2972       MEMBER *memb;

  2973       /* normally this pointer is NULL; the routine uses this pointer

  2974          to check data provided in the data section, in which case it

  2975          points to a member currently checked; this check is performed

  2976          automatically only once when a reference to any member occurs

  2977          for the first time */

  2978       SYMBOL *value;

  2979       /* evaluated symbolic value */

  2980 };

  2982 static void eval_sym_func(MPL *mpl, void *_info)

  2983 {     /* this is auxiliary routine to work within domain scope */

  2984       struct eval_sym_info *info = _info;

  2985       if (info->memb != NULL)

  2986       {  /* checking call; check symbolic value being assigned */

  2987          check_value_sym(mpl, info->par, info->memb->tuple,

  2988             info->memb->value.sym);

  2989       }

  2990       else

  2991       {  /* normal call; evaluate member, which has given n-tuple */

  2992          info->value = take_member_sym(mpl, info->par, info->tuple);

  2993       }

  2994       return;

  2995 }

  2997 SYMBOL *eval_member_sym       /* returns value, not reference */

  2998 (     MPL *mpl,

  2999       PARAMETER *par,         /* not changed */

  3000       TUPLE *tuple            /* not changed */

  3001 )

  3002 {     /* this routine evaluates symbolic parameter member */

  3003       struct eval_sym_info _info, *info = &_info;

  3004       xassert(par->type == A_SYMBOLIC);

  3005       xassert(par->dim == tuple_dimen(mpl, tuple));

  3006       info->par = par;

  3007       info->tuple = tuple;

  3008       if (par->data == 1)

  3009       {  /* check data, which are provided in the data section, but not

  3010             checked yet */

  3011          /* save pointer to the last array member; note that during the

  3012             check new members may be added beyond the last member due to

  3013             references to the same parameter from default expression as

  3014             well as from expressions that define restricting conditions;

  3015             however, values assigned to the new members will be checked

  3016             by other routine, so we don't need to check them here */

  3017          MEMBER *tail = par->array->tail;

  3018          /* change the data status to prevent infinite recursive loop

  3019             due to references to the same parameter during the check */

  3020          par->data = 2;

  3021          /* check values assigned to array members in the data section

  3022             until the marked member has been reached */

  3023          for (info->memb = par->array->head; info->memb != NULL;

  3024             info->memb = info->memb->next)

  3025          {  if (eval_within_domain(mpl, par->domain, info->memb->tuple,

  3026                info, eval_sym_func))

  3027                out_of_domain(mpl, par->name, info->memb->tuple);

  3028             if (info->memb == tail) break;

  3029          }

  3030          /* the check has been finished */

  3031       }

  3032       /* evaluate member, which has given n-tuple */

  3033       info->memb = NULL;

  3034       if (eval_within_domain(mpl, info->par->domain, info->tuple, info,

  3035          eval_sym_func))

  3036          out_of_domain(mpl, par->name, info->tuple);

  3037       /* bring evaluated value to the calling program */

  3038       return info->value;

  3039 }

  3041 /*----------------------------------------------------------------------

  3042 -- eval_whole_par - evaluate model parameter over entire domain.

  3043 --

  3044 -- This routine evaluates all members of specified model parameter over

  3045 -- entire domain. */

  3047 static int whole_par_func(MPL *mpl, void *info)

  3048 {     /* this is auxiliary routine to work within domain scope */

  3049       PARAMETER *par = (PARAMETER *)info;

  3050       TUPLE *tuple = get_domain_tuple(mpl, par->domain);

  3051       switch (par->type)

  3052       {  case A_NUMERIC:

  3053          case A_INTEGER:

  3054          case A_BINARY:

  3055             eval_member_num(mpl, par, tuple);

  3056             break;

  3057          case A_SYMBOLIC:

  3058             delete_symbol(mpl, eval_member_sym(mpl, par, tuple));

  3059             break;

  3060          default:

  3061             xassert(par != par);

  3062       }

  3063       delete_tuple(mpl, tuple);

  3064       return 0;

  3065 }

  3067 void eval_whole_par(MPL *mpl, PARAMETER *par)

  3068 {     loop_within_domain(mpl, par->domain, par, whole_par_func);

  3069       return;

  3070 }

  3072 /*----------------------------------------------------------------------

  3073 -- clean_parameter - clean model parameter.

  3074 --

  3075 -- This routine cleans specified model parameter that assumes deleting

  3076 -- all stuff dynamically allocated during the generation phase. */

  3078 void clean_parameter(MPL *mpl, PARAMETER *par)

  3079 {     CONDITION *cond;

  3080       WITHIN *in;

  3081       MEMBER *memb;

  3082       /* clean subscript domain */

  3083       clean_domain(mpl, par->domain);

  3084       /* clean pseudo-code for computing restricting conditions */

  3085       for (cond = par->cond; cond != NULL; cond = cond->next)

  3086          clean_code(mpl, cond->code);

  3087       /* clean pseudo-code for computing restricting supersets */

  3088       for (in = par->in; in != NULL; in = in->next)

  3089          clean_code(mpl, in->code);

  3090       /* clean pseudo-code for computing assigned value */

  3091       clean_code(mpl, par->assign);

  3092       /* clean pseudo-code for computing default value */

  3093       clean_code(mpl, par->option);

  3094       /* reset data status flag */

  3095       par->data = 0;

  3096       /* delete default symbolic value */

  3097       if (par->defval != NULL)

  3098          delete_symbol(mpl, par->defval), par->defval = NULL;

  3099       /* delete content array */

  3100       for (memb = par->array->head; memb != NULL; memb = memb->next)

  3101          delete_value(mpl, par->array->type, &memb->value);

  3102       delete_array(mpl, par->array), par->array = NULL;

  3103       return;

  3104 }

  3106 /**********************************************************************/

  3107 /* * *                      MODEL VARIABLES                       * * */

  3108 /**********************************************************************/

  3110 /*----------------------------------------------------------------------

  3111 -- take_member_var - obtain reference to elemental variable.

  3112 --

  3113 -- This routine obtains a reference to elemental variable assigned to

  3114 -- given member of specified model variable and returns it on exit. If

  3115 -- necessary, new elemental variable is created.

  3116 --

  3117 -- NOTE: This routine must not be called out of domain scope. */

  3119 ELEMVAR *take_member_var      /* returns reference */

  3120 (     MPL *mpl,

  3121       VARIABLE *var,          /* not changed */

  3122       TUPLE *tuple            /* not changed */

  3123 )

  3124 {     MEMBER *memb;

  3125       ELEMVAR *refer;

  3126       /* find member in the variable array */

  3127       memb = find_member(mpl, var->array, tuple);

  3128       if (memb != NULL)

  3129       {  /* member exists, so just take the reference */

  3130          refer = memb->value.var;

  3131       }

  3132       else

  3133       {  /* member is referenced for the first time and therefore does

  3134             not exist; create new elemental variable, assign it to new

  3135             member, and add the member to the variable array */

  3136          memb = add_member(mpl, var->array, copy_tuple(mpl, tuple));

  3137          refer = (memb->value.var =

  3138             dmp_get_atom(mpl->elemvars, sizeof(ELEMVAR)));

  3139          refer->j = 0;

  3140          refer->var = var;

  3141          refer->memb = memb;

  3142          /* compute lower bound */

  3143          if (var->lbnd == NULL)

  3144             refer->lbnd = 0.0;

  3145          else

  3146             refer->lbnd = eval_numeric(mpl, var->lbnd);

  3147          /* compute upper bound */

  3148          if (var->ubnd == NULL)

  3149             refer->ubnd = 0.0;

  3150          else if (var->ubnd == var->lbnd)

  3151             refer->ubnd = refer->lbnd;

  3152          else

  3153             refer->ubnd = eval_numeric(mpl, var->ubnd);

  3154          /* nullify working quantity */

  3155          refer->temp = 0.0;

  3156 #if 1 /* 15/V-2010 */

  3157          /* solution has not been obtained by the solver yet */

  3158          refer->stat = 0;

  3159          refer->prim = refer->dual = 0.0;

  3160 #endif

  3161       }

  3162       return refer;

  3163 }

  3165 /*----------------------------------------------------------------------

  3166 -- eval_member_var - evaluate reference to elemental variable.

  3167 --

  3168 -- This routine evaluates a reference to elemental variable assigned to

  3169 -- member of specified model variable and returns it on exit. */

  3171 struct eval_var_info

  3172 {     /* working info used by the routine eval_member_var */

  3173       VARIABLE *var;

  3174       /* model variable */

  3175       TUPLE *tuple;

  3176       /* n-tuple, which defines variable member */

  3177       ELEMVAR *refer;

  3178       /* evaluated reference to elemental variable */

  3179 };

  3181 static void eval_var_func(MPL *mpl, void *_info)

  3182 {     /* this is auxiliary routine to work within domain scope */

  3183       struct eval_var_info *info = _info;

  3184       info->refer = take_member_var(mpl, info->var, info->tuple);

  3185       return;

  3186 }

  3188 ELEMVAR *eval_member_var      /* returns reference */

  3189 (     MPL *mpl,

  3190       VARIABLE *var,          /* not changed */

  3191       TUPLE *tuple            /* not changed */

  3192 )

  3193 {     /* this routine evaluates variable member */

  3194       struct eval_var_info _info, *info = &_info;

  3195       xassert(var->dim == tuple_dimen(mpl, tuple));

  3196       info->var = var;

  3197       info->tuple = tuple;

  3198       /* evaluate member, which has given n-tuple */

  3199       if (eval_within_domain(mpl, info->var->domain, info->tuple, info,

  3200          eval_var_func))

  3201          out_of_domain(mpl, var->name, info->tuple);

  3202       /* bring evaluated reference to the calling program */

  3203       return info->refer;

  3204 }

  3206 /*----------------------------------------------------------------------

  3207 -- eval_whole_var - evaluate model variable over entire domain.

  3208 --

  3209 -- This routine evaluates all members of specified model variable over

  3210 -- entire domain. */

  3212 static int whole_var_func(MPL *mpl, void *info)

  3213 {     /* this is auxiliary routine to work within domain scope */

  3214       VARIABLE *var = (VARIABLE *)info;

  3215       TUPLE *tuple = get_domain_tuple(mpl, var->domain);

  3216       eval_member_var(mpl, var, tuple);

  3217       delete_tuple(mpl, tuple);

  3218       return 0;

  3219 }

  3221 void eval_whole_var(MPL *mpl, VARIABLE *var)

  3222 {     loop_within_domain(mpl, var->domain, var, whole_var_func);

  3223       return;

  3224 }

  3226 /*----------------------------------------------------------------------

  3227 -- clean_variable - clean model variable.

  3228 --

  3229 -- This routine cleans specified model variable that assumes deleting

  3230 -- all stuff dynamically allocated during the generation phase. */

  3232 void clean_variable(MPL *mpl, VARIABLE *var)

  3233 {     MEMBER *memb;

  3234       /* clean subscript domain */

  3235       clean_domain(mpl, var->domain);

  3236       /* clean code for computing lower bound */

  3237       clean_code(mpl, var->lbnd);

  3238       /* clean code for computing upper bound */

  3239       if (var->ubnd != var->lbnd) clean_code(mpl, var->ubnd);

  3240       /* delete content array */

  3241       for (memb = var->array->head; memb != NULL; memb = memb->next)

  3242          dmp_free_atom(mpl->elemvars, memb->value.var, sizeof(ELEMVAR));

  3243       delete_array(mpl, var->array), var->array = NULL;

  3244       return;

  3245 }

  3247 /**********************************************************************/

  3248 /* * *              MODEL CONSTRAINTS AND OBJECTIVES              * * */

  3249 /**********************************************************************/

  3251 /*----------------------------------------------------------------------

  3252 -- take_member_con - obtain reference to elemental constraint.

  3253 --

  3254 -- This routine obtains a reference to elemental constraint assigned

  3255 -- to given member of specified model constraint and returns it on exit.

  3256 -- If necessary, new elemental constraint is created.

  3257 --

  3258 -- NOTE: This routine must not be called out of domain scope. */

  3260 ELEMCON *take_member_con      /* returns reference */

  3261 (     MPL *mpl,

  3262       CONSTRAINT *con,        /* not changed */

  3263       TUPLE *tuple            /* not changed */

  3264 )

  3265 {     MEMBER *memb;

  3266       ELEMCON *refer;

  3267       /* find member in the constraint array */

  3268       memb = find_member(mpl, con->array, tuple);

  3269       if (memb != NULL)

  3270       {  /* member exists, so just take the reference */

  3271          refer = memb->value.con;

  3272       }

  3273       else

  3274       {  /* member is referenced for the first time and therefore does

  3275             not exist; create new elemental constraint, assign it to new

  3276             member, and add the member to the constraint array */

  3277          memb = add_member(mpl, con->array, copy_tuple(mpl, tuple));

  3278          refer = (memb->value.con =

  3279             dmp_get_atom(mpl->elemcons, sizeof(ELEMCON)));

  3280          refer->i = 0;

  3281          refer->con = con;

  3282          refer->memb = memb;

  3283          /* compute linear form */

  3284          xassert(con->code != NULL);

  3285          refer->form = eval_formula(mpl, con->code);

  3286          /* compute lower and upper bounds */

  3287          if (con->lbnd == NULL && con->ubnd == NULL)

  3288          {  /* objective has no bounds */

  3289             double temp;

  3290             xassert(con->type == A_MINIMIZE || con->type == A_MAXIMIZE);

  3291             /* carry the constant term to the right-hand side */

  3292             refer->form = remove_constant(mpl, refer->form, &temp);

  3293             refer->lbnd = refer->ubnd = - temp;

  3294          }

  3295          else if (con->lbnd != NULL && con->ubnd == NULL)

  3296          {  /* constraint a * x + b >= c * y + d is transformed to the

  3297                standard form a * x - c * y >= d - b */

  3298             double temp;

  3299             xassert(con->type == A_CONSTRAINT);

  3300             refer->form = linear_comb(mpl,

  3301                +1.0, refer->form,

  3302                -1.0, eval_formula(mpl, con->lbnd));

  3303             refer->form = remove_constant(mpl, refer->form, &temp);

  3304             refer->lbnd = - temp;

  3305             refer->ubnd = 0.0;

  3306          }

  3307          else if (con->lbnd == NULL && con->ubnd != NULL)

  3308          {  /* constraint a * x + b <= c * y + d is transformed to the

  3309                standard form a * x - c * y <= d - b */

  3310             double temp;

  3311             xassert(con->type == A_CONSTRAINT);

  3312             refer->form = linear_comb(mpl,

  3313                +1.0, refer->form,

  3314                -1.0, eval_formula(mpl, con->ubnd));

  3315             refer->form = remove_constant(mpl, refer->form, &temp);

  3316             refer->lbnd = 0.0;

  3317             refer->ubnd = - temp;

  3318          }

  3319          else if (con->lbnd == con->ubnd)

  3320          {  /* constraint a * x + b = c * y + d is transformed to the

  3321                standard form a * x - c * y = d - b */

  3322             double temp;

  3323             xassert(con->type == A_CONSTRAINT);

  3324             refer->form = linear_comb(mpl,

  3325                +1.0, refer->form,

  3326                -1.0, eval_formula(mpl, con->lbnd));

  3327             refer->form = remove_constant(mpl, refer->form, &temp);

  3328             refer->lbnd = refer->ubnd = - temp;

  3329          }

  3330          else

  3331          {  /* ranged constraint c <= a * x + b <= d is transformed to

  3332                the standard form c - b <= a * x <= d - b */

  3333             double temp, temp1, temp2;

  3334             xassert(con->type == A_CONSTRAINT);

  3335             refer->form = remove_constant(mpl, refer->form, &temp);

  3336             xassert(remove_constant(mpl, eval_formula(mpl, con->lbnd),

  3337                &temp1) == NULL);

  3338             xassert(remove_constant(mpl, eval_formula(mpl, con->ubnd),

  3339                &temp2) == NULL);

  3340             refer->lbnd = fp_sub(mpl, temp1, temp);

  3341             refer->ubnd = fp_sub(mpl, temp2, temp);

  3342          }

  3343 #if 1 /* 15/V-2010 */

  3344          /* solution has not been obtained by the solver yet */

  3345          refer->stat = 0;

  3346          refer->prim = refer->dual = 0.0;

  3347 #endif

  3348       }

  3349       return refer;

  3350 }

  3352 /*----------------------------------------------------------------------

  3353 -- eval_member_con - evaluate reference to elemental constraint.

  3354 --

  3355 -- This routine evaluates a reference to elemental constraint assigned

  3356 -- to member of specified model constraint and returns it on exit. */

  3358 struct eval_con_info

  3359 {     /* working info used by the routine eval_member_con */

  3360       CONSTRAINT *con;

  3361       /* model constraint */

  3362       TUPLE *tuple;

  3363       /* n-tuple, which defines constraint member */

  3364       ELEMCON *refer;

  3365       /* evaluated reference to elemental constraint */

  3366 };

  3368 static void eval_con_func(MPL *mpl, void *_info)

  3369 {     /* this is auxiliary routine to work within domain scope */

  3370       struct eval_con_info *info = _info;

  3371       info->refer = take_member_con(mpl, info->con, info->tuple);

  3372       return;

  3373 }

  3375 ELEMCON *eval_member_con      /* returns reference */

  3376 (     MPL *mpl,

  3377       CONSTRAINT *con,        /* not changed */

  3378       TUPLE *tuple            /* not changed */

  3379 )

  3380 {     /* this routine evaluates constraint member */

  3381       struct eval_con_info _info, *info = &_info;

  3382       xassert(con->dim == tuple_dimen(mpl, tuple));

  3383       info->con = con;

  3384       info->tuple = tuple;

  3385       /* evaluate member, which has given n-tuple */

  3386       if (eval_within_domain(mpl, info->con->domain, info->tuple, info,

  3387          eval_con_func))

  3388          out_of_domain(mpl, con->name, info->tuple);

  3389       /* bring evaluated reference to the calling program */

  3390       return info->refer;

  3391 }

  3393 /*----------------------------------------------------------------------

  3394 -- eval_whole_con - evaluate model constraint over entire domain.

  3395 --

  3396 -- This routine evaluates all members of specified model constraint over

  3397 -- entire domain. */

  3399 static int whole_con_func(MPL *mpl, void *info)

  3400 {     /* this is auxiliary routine to work within domain scope */

  3401       CONSTRAINT *con = (CONSTRAINT *)info;

  3402       TUPLE *tuple = get_domain_tuple(mpl, con->domain);

  3403       eval_member_con(mpl, con, tuple);

  3404       delete_tuple(mpl, tuple);

  3405       return 0;

  3406 }

  3408 void eval_whole_con(MPL *mpl, CONSTRAINT *con)

  3409 {     loop_within_domain(mpl, con->domain, con, whole_con_func);

  3410       return;

  3411 }

  3413 /*----------------------------------------------------------------------

  3414 -- clean_constraint - clean model constraint.

  3415 --

  3416 -- This routine cleans specified model constraint that assumes deleting

  3417 -- all stuff dynamically allocated during the generation phase. */

  3419 void clean_constraint(MPL *mpl, CONSTRAINT *con)

  3420 {     MEMBER *memb;

  3421       /* clean subscript domain */

  3422       clean_domain(mpl, con->domain);

  3423       /* clean code for computing main linear form */

  3424       clean_code(mpl, con->code);

  3425       /* clean code for computing lower bound */

  3426       clean_code(mpl, con->lbnd);

  3427       /* clean code for computing upper bound */

  3428       if (con->ubnd != con->lbnd) clean_code(mpl, con->ubnd);

  3429       /* delete content array */

  3430       for (memb = con->array->head; memb != NULL; memb = memb->next)

  3431       {  delete_formula(mpl, memb->value.con->form);

  3432          dmp_free_atom(mpl->elemcons, memb->value.con, sizeof(ELEMCON));

  3433       }

  3434       delete_array(mpl, con->array), con->array = NULL;

  3435       return;

  3436 }

  3438 /**********************************************************************/

  3439 /* * *                        PSEUDO-CODE                         * * */

  3440 /**********************************************************************/

  3442 /*----------------------------------------------------------------------

  3443 -- eval_numeric - evaluate pseudo-code to determine numeric value.

  3444 --

  3445 -- This routine evaluates specified pseudo-code to determine resultant

  3446 -- numeric value, which is returned on exit. */

  3448 struct iter_num_info

  3449 {     /* working info used by the routine iter_num_func */

  3450       CODE *code;

  3451       /* pseudo-code for iterated operation to be performed */

  3452       double value;

  3453       /* resultant value */

  3454 };

  3456 static int iter_num_func(MPL *mpl, void *_info)

  3457 {     /* this is auxiliary routine used to perform iterated operation

  3458          on numeric "integrand" within domain scope */

  3459       struct iter_num_info *info = _info;

  3460       double temp;

  3461       temp = eval_numeric(mpl, info->code->arg.loop.x);

  3462       switch (info->code->op)

  3463       {  case O_SUM:

  3464             /* summation over domain */

  3465             info->value = fp_add(mpl, info->value, temp);

  3466             break;

  3467          case O_PROD:

  3468             /* multiplication over domain */

  3469             info->value = fp_mul(mpl, info->value, temp);

  3470             break;

  3471          case O_MINIMUM:

  3472             /* minimum over domain */

  3473             if (info->value > temp) info->value = temp;

  3474             break;

  3475          case O_MAXIMUM:

  3476             /* maximum over domain */

  3477             if (info->value < temp) info->value = temp;

  3478             break;

  3479          default:

  3480             xassert(info != info);

  3481       }

  3482       return 0;

  3483 }

  3485 double eval_numeric(MPL *mpl, CODE *code)

  3486 {     double value;

  3487       xassert(code != NULL);

  3488       xassert(code->type == A_NUMERIC);

  3489       xassert(code->dim == 0);

  3490       /* if the operation has a side effect, invalidate and delete the

  3491          resultant value */

  3492       if (code->vflag && code->valid)

  3493       {  code->valid = 0;

  3494          delete_value(mpl, code->type, &code->value);

  3495       }

  3496       /* if resultant value is valid, no evaluation is needed */

  3497       if (code->valid)

  3498       {  value = code->value.num;

  3499          goto done;

  3500       }

  3501       /* evaluate pseudo-code recursively */

  3502       switch (code->op)

  3503       {  case O_NUMBER:

  3504             /* take floating-point number */

  3505             value = code->arg.num;

  3506             break;

  3507          case O_MEMNUM:

  3508             /* take member of numeric parameter */

  3509             {  TUPLE *tuple;

  3510                ARG_LIST *e;

  3511                tuple = create_tuple(mpl);

  3512                for (e = code->arg.par.list; e != NULL; e = e->next)

  3513                   tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,

  3514                      e->x));

  3515                value = eval_member_num(mpl, code->arg.par.par, tuple);

  3516                delete_tuple(mpl, tuple);

  3517             }

  3518             break;

  3519          case O_MEMVAR:

  3520             /* take computed value of elemental variable */

  3521             {  TUPLE *tuple;

  3522                ARG_LIST *e;

  3523 #if 1 /* 15/V-2010 */

  3524                ELEMVAR *var;

  3525 #endif

  3526                tuple = create_tuple(mpl);

  3527                for (e = code->arg.var.list; e != NULL; e = e->next)

  3528                   tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,

  3529                      e->x));

  3530 #if 0 /* 15/V-2010 */

  3531                value = eval_member_var(mpl, code->arg.var.var, tuple)

  3532                   ->value;

  3533 #else

  3534                var = eval_member_var(mpl, code->arg.var.var, tuple);

  3535                switch (code->arg.var.suff)

  3536                {  case DOT_LB:

  3537                      if (var->var->lbnd == NULL)

  3538                         value = -DBL_MAX;

  3539                      else

  3540                         value = var->lbnd;

  3541                      break;

  3542                   case DOT_UB:

  3543                      if (var->var->ubnd == NULL)

  3544                         value = +DBL_MAX;

  3545                      else

  3546                         value = var->ubnd;

  3547                      break;

  3548                   case DOT_STATUS:

  3549                      value = var->stat;

  3550                      break;

  3551                   case DOT_VAL:

  3552                      value = var->prim;

  3553                      break;

  3554                   case DOT_DUAL:

  3555                      value = var->dual;

  3556                      break;

  3557                   default:

  3558                      xassert(code != code);

  3559                }

  3560 #endif

  3561                delete_tuple(mpl, tuple);

  3562             }

  3563             break;

  3564 #if 1 /* 15/V-2010 */

  3565          case O_MEMCON:

  3566             /* take computed value of elemental constraint */

  3567             {  TUPLE *tuple;

  3568                ARG_LIST *e;

  3569                ELEMCON *con;

  3570                tuple = create_tuple(mpl);

  3571                for (e = code->arg.con.list; e != NULL; e = e->next)

  3572                   tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,

  3573                      e->x));

  3574                con = eval_member_con(mpl, code->arg.con.con, tuple);

  3575                switch (code->arg.con.suff)

  3576                {  case DOT_LB:

  3577                      if (con->con->lbnd == NULL)

  3578                         value = -DBL_MAX;

  3579                      else

  3580                         value = con->lbnd;

  3581                      break;

  3582                   case DOT_UB:

  3583                      if (con->con->ubnd == NULL)

  3584                         value = +DBL_MAX;

  3585                      else

  3586                         value = con->ubnd;

  3587                      break;

  3588                   case DOT_STATUS:

  3589                      value = con->stat;

  3590                      break;

  3591                   case DOT_VAL:

  3592                      value = con->prim;

  3593                      break;

  3594                   case DOT_DUAL:

  3595                      value = con->dual;

  3596                      break;

  3597                   default:

  3598                      xassert(code != code);

  3599                }

  3600                delete_tuple(mpl, tuple);

  3601             }

  3602             break;

  3603 #endif

  3604          case O_IRAND224:

  3605             /* pseudo-random in [0, 2^24-1] */

  3606             value = fp_irand224(mpl);

  3607             break;

  3608          case O_UNIFORM01:

  3609             /* pseudo-random in [0, 1) */

  3610             value = fp_uniform01(mpl);

  3611             break;

  3612          case O_NORMAL01:

  3613             /* gaussian random, mu = 0, sigma = 1 */

  3614             value = fp_normal01(mpl);

  3615             break;

  3616          case O_GMTIME:

  3617             /* current calendar time */

  3618             value = fn_gmtime(mpl);

  3619             break;

  3620          case O_CVTNUM:

  3621             /* conversion to numeric */

  3622             {  SYMBOL *sym;

  3623                sym = eval_symbolic(mpl, code->arg.arg.x);

  3624 #if 0 /* 23/XI-2008 */

  3625                if (sym->str != NULL)

  3626                   error(mpl, "cannot convert %s to floating-point numbe"

  3627                      "r", format_symbol(mpl, sym));

  3628                value = sym->num;

  3629 #else

  3630                if (sym->str == NULL)

  3631                   value = sym->num;

  3632                else

  3633                {  if (str2num(sym->str, &value))

  3634                      error(mpl, "cannot convert %s to floating-point nu"

  3635                         "mber", format_symbol(mpl, sym));

  3636                }

  3637 #endif

  3638                delete_symbol(mpl, sym);

  3639             }

  3640             break;

  3641          case O_PLUS:

  3642             /* unary plus */

  3643             value = + eval_numeric(mpl, code->arg.arg.x);

  3644             break;

  3645          case O_MINUS:

  3646             /* unary minus */

  3647             value = - eval_numeric(mpl, code->arg.arg.x);

  3648             break;

  3649          case O_ABS:

  3650             /* absolute value */

  3651             value = fabs(eval_numeric(mpl, code->arg.arg.x));

  3652             break;

  3653          case O_CEIL:

  3654             /* round upward ("ceiling of x") */

  3655             value = ceil(eval_numeric(mpl, code->arg.arg.x));

  3656             break;

  3657          case O_FLOOR:

  3658             /* round downward ("floor of x") */

  3659             value = floor(eval_numeric(mpl, code->arg.arg.x));

  3660             break;

  3661          case O_EXP:

  3662             /* base-e exponential */

  3663             value = fp_exp(mpl, eval_numeric(mpl, code->arg.arg.x));

  3664             break;

  3665          case O_LOG:

  3666             /* natural logarithm */

  3667             value = fp_log(mpl, eval_numeric(mpl, code->arg.arg.x));

  3668             break;

  3669          case O_LOG10:

  3670             /* common (decimal) logarithm */

  3671             value = fp_log10(mpl, eval_numeric(mpl, code->arg.arg.x));

  3672             break;

  3673          case O_SQRT:

  3674             /* square root */

  3675             value = fp_sqrt(mpl, eval_numeric(mpl, code->arg.arg.x));

  3676             break;

  3677          case O_SIN:

  3678             /* trigonometric sine */

  3679             value = fp_sin(mpl, eval_numeric(mpl, code->arg.arg.x));

  3680             break;

  3681          case O_COS:

  3682             /* trigonometric cosine */

  3683             value = fp_cos(mpl, eval_numeric(mpl, code->arg.arg.x));

  3684             break;

  3685          case O_ATAN:

  3686             /* trigonometric arctangent (one argument) */

  3687             value = fp_atan(mpl, eval_numeric(mpl, code->arg.arg.x));

  3688             break;

  3689          case O_ATAN2:

  3690             /* trigonometric arctangent (two arguments) */

  3691             value = fp_atan2(mpl,

  3692                eval_numeric(mpl, code->arg.arg.x),

  3693                eval_numeric(mpl, code->arg.arg.y));

  3694             break;

  3695          case O_ROUND:

  3696             /* round to nearest integer */

  3697             value = fp_round(mpl,

  3698                eval_numeric(mpl, code->arg.arg.x), 0.0);

  3699             break;

  3700          case O_ROUND2:

  3701             /* round to n fractional digits */

  3702             value = fp_round(mpl,

  3703                eval_numeric(mpl, code->arg.arg.x),

  3704                eval_numeric(mpl, code->arg.arg.y));

  3705             break;

  3706          case O_TRUNC:

  3707             /* truncate to nearest integer */

  3708             value = fp_trunc(mpl,

  3709                eval_numeric(mpl, code->arg.arg.x), 0.0);

  3710             break;

  3711          case O_TRUNC2:

  3712             /* truncate to n fractional digits */

  3713             value = fp_trunc(mpl,

  3714                eval_numeric(mpl, code->arg.arg.x),

  3715                eval_numeric(mpl, code->arg.arg.y));

  3716             break;

  3717          case O_ADD:

  3718             /* addition */

  3719             value = fp_add(mpl,

  3720                eval_numeric(mpl, code->arg.arg.x),

  3721                eval_numeric(mpl, code->arg.arg.y));

  3722             break;

  3723          case O_SUB:

  3724             /* subtraction */

  3725             value = fp_sub(mpl,

  3726                eval_numeric(mpl, code->arg.arg.x),

  3727                eval_numeric(mpl, code->arg.arg.y));

  3728             break;

  3729          case O_LESS:

  3730             /* non-negative subtraction */

  3731             value = fp_less(mpl,

  3732                eval_numeric(mpl, code->arg.arg.x),

  3733                eval_numeric(mpl, code->arg.arg.y));

  3734             break;

  3735          case O_MUL:

  3736             /* multiplication */

  3737             value = fp_mul(mpl,

  3738                eval_numeric(mpl, code->arg.arg.x),

  3739                eval_numeric(mpl, code->arg.arg.y));

  3740             break;

  3741          case O_DIV:

  3742             /* division */

  3743             value = fp_div(mpl,

  3744                eval_numeric(mpl, code->arg.arg.x),

  3745                eval_numeric(mpl, code->arg.arg.y));

  3746             break;

  3747          case O_IDIV:

  3748             /* quotient of exact division */

  3749             value = fp_idiv(mpl,

  3750                eval_numeric(mpl, code->arg.arg.x),

  3751                eval_numeric(mpl, code->arg.arg.y));

  3752             break;

  3753          case O_MOD:

  3754             /* remainder of exact division */

  3755             value = fp_mod(mpl,

  3756                eval_numeric(mpl, code->arg.arg.x),

  3757                eval_numeric(mpl, code->arg.arg.y));

  3758             break;

  3759          case O_POWER:

  3760             /* exponentiation (raise to power) */

  3761             value = fp_power(mpl,

  3762                eval_numeric(mpl, code->arg.arg.x),

  3763                eval_numeric(mpl, code->arg.arg.y));

  3764             break;

  3765          case O_UNIFORM:

  3766             /* pseudo-random in [a, b) */

  3767             value = fp_uniform(mpl,

  3768                eval_numeric(mpl, code->arg.arg.x),

  3769                eval_numeric(mpl, code->arg.arg.y));

  3770             break;

  3771          case O_NORMAL:

  3772             /* gaussian random, given mu and sigma */

  3773             value = fp_normal(mpl,

  3774                eval_numeric(mpl, code->arg.arg.x),

  3775                eval_numeric(mpl, code->arg.arg.y));

  3776             break;

  3777          case O_CARD:

  3778             {  ELEMSET *set;

  3779                set = eval_elemset(mpl, code->arg.arg.x);

  3780                value = set->size;

  3781                delete_array(mpl, set);

  3782             }

  3783             break;

  3784          case O_LENGTH:

  3785             {  SYMBOL *sym;

  3786                char str[MAX_LENGTH+1];

  3787                sym = eval_symbolic(mpl, code->arg.arg.x);

  3788                if (sym->str == NULL)

  3789                   sprintf(str, "%.*g", DBL_DIG, sym->num);

  3790                else

  3791                   fetch_string(mpl, sym->str, str);

  3792                delete_symbol(mpl, sym);

  3793                value = strlen(str);

  3794             }

  3795             break;

  3796          case O_STR2TIME:

  3797             {  SYMBOL *sym;

  3798                char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];

  3799                sym = eval_symbolic(mpl, code->arg.arg.x);

  3800                if (sym->str == NULL)

  3801                   sprintf(str, "%.*g", DBL_DIG, sym->num);

  3802                else

  3803                   fetch_string(mpl, sym->str, str);

  3804                delete_symbol(mpl, sym);

  3805                sym = eval_symbolic(mpl, code->arg.arg.y);

  3806                if (sym->str == NULL)

  3807                   sprintf(fmt, "%.*g", DBL_DIG, sym->num);

  3808                else

  3809                   fetch_string(mpl, sym->str, fmt);

  3810                delete_symbol(mpl, sym);

  3811                value = fn_str2time(mpl, str, fmt);

  3812             }

  3813             break;

  3814          case O_FORK:

  3815             /* if-then-else */

  3816             if (eval_logical(mpl, code->arg.arg.x))

  3817                value = eval_numeric(mpl, code->arg.arg.y);

  3818             else if (code->arg.arg.z == NULL)

  3819                value = 0.0;

  3820             else

  3821                value = eval_numeric(mpl, code->arg.arg.z);

  3822             break;

  3823          case O_MIN:

  3824             /* minimal value (n-ary) */

  3825             {  ARG_LIST *e;

  3826                double temp;

  3827                value = +DBL_MAX;

  3828                for (e = code->arg.list; e != NULL; e = e->next)

  3829                {  temp = eval_numeric(mpl, e->x);

  3830                   if (value > temp) value = temp;

  3831                }

  3832             }

  3833             break;

  3834          case O_MAX:

  3835             /* maximal value (n-ary) */

  3836             {  ARG_LIST *e;

  3837                double temp;

  3838                value = -DBL_MAX;

  3839                for (e = code->arg.list; e != NULL; e = e->next)

  3840                {  temp = eval_numeric(mpl, e->x);

  3841                   if (value < temp) value = temp;

  3842                }

  3843             }

  3844             break;

  3845          case O_SUM:

  3846             /* summation over domain */

  3847             {  struct iter_num_info _info, *info = &_info;

  3848                info->code = code;

  3849                info->value = 0.0;

  3850                loop_within_domain(mpl, code->arg.loop.domain, info,

  3851                   iter_num_func);

  3852                value = info->value;

  3853             }

  3854             break;

  3855          case O_PROD:

  3856             /* multiplication over domain */

  3857             {  struct iter_num_info _info, *info = &_info;

  3858                info->code = code;

  3859                info->value = 1.0;

  3860                loop_within_domain(mpl, code->arg.loop.domain, info,

  3861                   iter_num_func);

  3862                value = info->value;

  3863             }

  3864             break;

  3865          case O_MINIMUM:

  3866             /* minimum over domain */

  3867             {  struct iter_num_info _info, *info = &_info;

  3868                info->code = code;

  3869                info->value = +DBL_MAX;

  3870                loop_within_domain(mpl, code->arg.loop.domain, info,

  3871                   iter_num_func);

  3872                if (info->value == +DBL_MAX)

  3873                   error(mpl, "min{} over empty set; result undefined");

  3874                value = info->value;

  3875             }

  3876             break;

  3877          case O_MAXIMUM:

  3878             /* maximum over domain */

  3879             {  struct iter_num_info _info, *info = &_info;

  3880                info->code = code;

  3881                info->value = -DBL_MAX;

  3882                loop_within_domain(mpl, code->arg.loop.domain, info,

  3883                   iter_num_func);

  3884                if (info->value == -DBL_MAX)

  3885                   error(mpl, "max{} over empty set; result undefined");

  3886                value = info->value;

  3887             }

  3888             break;

  3889          default:

  3890             xassert(code != code);

  3891       }

  3892       /* save resultant value */

  3893       xassert(!code->valid);

  3894       code->valid = 1;

  3895       code->value.num = value;

  3896 done: return value;

  3897 }

  3899 /*----------------------------------------------------------------------

  3900 -- eval_symbolic - evaluate pseudo-code to determine symbolic value.

  3901 --

  3902 -- This routine evaluates specified pseudo-code to determine resultant

  3903 -- symbolic value, which is returned on exit. */

  3905 SYMBOL *eval_symbolic(MPL *mpl, CODE *code)

  3906 {     SYMBOL *value;

  3907       xassert(code != NULL);

  3908       xassert(code->type == A_SYMBOLIC);

  3909       xassert(code->dim == 0);

  3910       /* if the operation has a side effect, invalidate and delete the

  3911          resultant value */

  3912       if (code->vflag && code->valid)

  3913       {  code->valid = 0;

  3914          delete_value(mpl, code->type, &code->value);

  3915       }

  3916       /* if resultant value is valid, no evaluation is needed */

  3917       if (code->valid)

  3918       {  value = copy_symbol(mpl, code->value.sym);

  3919          goto done;

  3920       }

  3921       /* evaluate pseudo-code recursively */

  3922       switch (code->op)

  3923       {  case O_STRING:

  3924             /* take character string */

  3925             value = create_symbol_str(mpl, create_string(mpl,

  3926                code->arg.str));

  3927             break;

  3928          case O_INDEX:

  3929             /* take dummy index */

  3930             xassert(code->arg.index.slot->value != NULL);

  3931             value = copy_symbol(mpl, code->arg.index.slot->value);

  3932             break;

  3933          case O_MEMSYM:

  3934             /* take member of symbolic parameter */

  3935             {  TUPLE *tuple;

  3936                ARG_LIST *e;

  3937                tuple = create_tuple(mpl);

  3938                for (e = code->arg.par.list; e != NULL; e = e->next)

  3939                   tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,

  3940                      e->x));

  3941                value = eval_member_sym(mpl, code->arg.par.par, tuple);

  3942                delete_tuple(mpl, tuple);

  3943             }

  3944             break;

  3945          case O_CVTSYM:

  3946             /* conversion to symbolic */

  3947             value = create_symbol_num(mpl, eval_numeric(mpl,

  3948                code->arg.arg.x));

  3949             break;

  3950          case O_CONCAT:

  3951             /* concatenation */

  3952             value = concat_symbols(mpl,

  3953                eval_symbolic(mpl, code->arg.arg.x),

  3954                eval_symbolic(mpl, code->arg.arg.y));

  3955             break;

  3956          case O_FORK:

  3957             /* if-then-else */

  3958             if (eval_logical(mpl, code->arg.arg.x))

  3959                value = eval_symbolic(mpl, code->arg.arg.y);

  3960             else if (code->arg.arg.z == NULL)

  3961                value = create_symbol_num(mpl, 0.0);

  3962             else

  3963                value = eval_symbolic(mpl, code->arg.arg.z);

  3964             break;

  3965          case O_SUBSTR:

  3966          case O_SUBSTR3:

  3967             {  double pos, len;

  3968                char str[MAX_LENGTH+1];

  3969                value = eval_symbolic(mpl, code->arg.arg.x);

  3970                if (value->str == NULL)

  3971                   sprintf(str, "%.*g", DBL_DIG, value->num);

  3972                else

  3973                   fetch_string(mpl, value->str, str);

  3974                delete_symbol(mpl, value);

  3975                if (code->op == O_SUBSTR)

  3976                {  pos = eval_numeric(mpl, code->arg.arg.y);

  3977                   if (pos != floor(pos))

  3978                      error(mpl, "substr('...', %.*g); non-integer secon"

  3979                         "d argument", DBL_DIG, pos);

  3980                   if (pos < 1 || pos > strlen(str) + 1)

  3981                      error(mpl, "substr('...', %.*g); substring out of "

  3982                         "range", DBL_DIG, pos);

  3983                }

  3984                else

  3985                {  pos = eval_numeric(mpl, code->arg.arg.y);

  3986                   len = eval_numeric(mpl, code->arg.arg.z);

  3987                   if (pos != floor(pos) || len != floor(len))

  3988                      error(mpl, "substr('...', %.*g, %.*g); non-integer"

  3989                         " second and/or third argument", DBL_DIG, pos,

  3990                         DBL_DIG, len);

  3991                   if (pos < 1 || len < 0 || pos + len > strlen(str) + 1)

  3992                      error(mpl, "substr('...', %.*g, %.*g); substring o"

  3993                         "ut of range", DBL_DIG, pos, DBL_DIG, len);

  3994                   str[(int)pos + (int)len - 1] = '\0';

  3995                }

  3996                value = create_symbol_str(mpl, create_string(mpl, str +

  3997                   (int)pos - 1));

  3998             }

  3999             break;

  4000          case O_TIME2STR:

  4001             {  double num;

  4002                SYMBOL *sym;

  4003                char str[MAX_LENGTH+1], fmt[MAX_LENGTH+1];

  4004                num = eval_numeric(mpl, code->arg.arg.x);

  4005                sym = eval_symbolic(mpl, code->arg.arg.y);

  4006                if (sym->str == NULL)

  4007                   sprintf(fmt, "%.*g", DBL_DIG, sym->num);

  4008                else

  4009                   fetch_string(mpl, sym->str, fmt);

  4010                delete_symbol(mpl, sym);

  4011                fn_time2str(mpl, str, num, fmt);

  4012                value = create_symbol_str(mpl, create_string(mpl, str));

  4013             }

  4014             break;

  4015          default:

  4016             xassert(code != code);

  4017       }

  4018       /* save resultant value */

  4019       xassert(!code->valid);

  4020       code->valid = 1;

  4021       code->value.sym = copy_symbol(mpl, value);

  4022 done: return value;

  4023 }

  4025 /*----------------------------------------------------------------------

  4026 -- eval_logical - evaluate pseudo-code to determine logical value.

  4027 --

  4028 -- This routine evaluates specified pseudo-code to determine resultant

  4029 -- logical value, which is returned on exit. */

  4031 struct iter_log_info

  4032 {     /* working info used by the routine iter_log_func */

  4033       CODE *code;

  4034       /* pseudo-code for iterated operation to be performed */

  4035       int value;

  4036       /* resultant value */

  4037 };

  4039 static int iter_log_func(MPL *mpl, void *_info)

  4040 {     /* this is auxiliary routine used to perform iterated operation

  4041          on logical "integrand" within domain scope */

  4042       struct iter_log_info *info = _info;

  4043       int ret = 0;

  4044       switch (info->code->op)

  4045       {  case O_FORALL:

  4046             /* conjunction over domain */

  4047             info->value &= eval_logical(mpl, info->code->arg.loop.x);

  4048             if (!info->value) ret = 1;

  4049             break;

  4050          case O_EXISTS:

  4051             /* disjunction over domain */

  4052             info->value |= eval_logical(mpl, info->code->arg.loop.x);

  4053             if (info->value) ret = 1;

  4054             break;

  4055          default:

  4056             xassert(info != info);

  4057       }

  4058       return ret;

  4059 }

  4061 int eval_logical(MPL *mpl, CODE *code)

  4062 {     int value;

  4063       xassert(code->type == A_LOGICAL);

  4064       xassert(code->dim == 0);

  4065       /* if the operation has a side effect, invalidate and delete the

  4066          resultant value */

  4067       if (code->vflag && code->valid)

  4068       {  code->valid = 0;

  4069          delete_value(mpl, code->type, &code->value);

  4070       }

  4071       /* if resultant value is valid, no evaluation is needed */

  4072       if (code->valid)

  4073       {  value = code->value.bit;

  4074          goto done;

  4075       }

  4076       /* evaluate pseudo-code recursively */

  4077       switch (code->op)

  4078       {  case O_CVTLOG:

  4079             /* conversion to logical */

  4080             value = (eval_numeric(mpl, code->arg.arg.x) != 0.0);

  4081             break;

  4082          case O_NOT:

  4083             /* negation (logical "not") */

  4084             value = !eval_logical(mpl, code->arg.arg.x);

  4085             break;

  4086          case O_LT:

  4087             /* comparison on 'less than' */

  4088 #if 0 /* 02/VIII-2008 */

  4089             value = (eval_numeric(mpl, code->arg.arg.x) <

  4090                      eval_numeric(mpl, code->arg.arg.y));

  4091 #else

  4092             xassert(code->arg.arg.x != NULL);

  4093             if (code->arg.arg.x->type == A_NUMERIC)

  4094                value = (eval_numeric(mpl, code->arg.arg.x) <

  4095                         eval_numeric(mpl, code->arg.arg.y));

  4096             else

  4097             {  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);

  4098                SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);

  4099                value = (compare_symbols(mpl, sym1, sym2) < 0);

  4100                delete_symbol(mpl, sym1);

  4101                delete_symbol(mpl, sym2);

  4102             }

  4103 #endif

  4104             break;

  4105          case O_LE:

  4106             /* comparison on 'not greater than' */

  4107 #if 0 /* 02/VIII-2008 */

  4108             value = (eval_numeric(mpl, code->arg.arg.x) <=

  4109                      eval_numeric(mpl, code->arg.arg.y));

  4110 #else

  4111             xassert(code->arg.arg.x != NULL);

  4112             if (code->arg.arg.x->type == A_NUMERIC)

  4113                value = (eval_numeric(mpl, code->arg.arg.x) <=

  4114                         eval_numeric(mpl, code->arg.arg.y));

  4115             else

  4116             {  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);

  4117                SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);

  4118                value = (compare_symbols(mpl, sym1, sym2) <= 0);

  4119                delete_symbol(mpl, sym1);

  4120                delete_symbol(mpl, sym2);

  4121             }

  4122 #endif

  4123             break;

  4124          case O_EQ:

  4125             /* comparison on 'equal to' */

  4126             xassert(code->arg.arg.x != NULL);

  4127             if (code->arg.arg.x->type == A_NUMERIC)

  4128                value = (eval_numeric(mpl, code->arg.arg.x) ==

  4129                         eval_numeric(mpl, code->arg.arg.y));

  4130             else

  4131             {  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);

  4132                SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);

  4133                value = (compare_symbols(mpl, sym1, sym2) == 0);

  4134                delete_symbol(mpl, sym1);

  4135                delete_symbol(mpl, sym2);

  4136             }

  4137             break;

  4138          case O_GE:

  4139             /* comparison on 'not less than' */

  4140 #if 0 /* 02/VIII-2008 */

  4141             value = (eval_numeric(mpl, code->arg.arg.x) >=

  4142                      eval_numeric(mpl, code->arg.arg.y));

  4143 #else

  4144             xassert(code->arg.arg.x != NULL);

  4145             if (code->arg.arg.x->type == A_NUMERIC)

  4146                value = (eval_numeric(mpl, code->arg.arg.x) >=

  4147                         eval_numeric(mpl, code->arg.arg.y));

  4148             else

  4149             {  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);

  4150                SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);

  4151                value = (compare_symbols(mpl, sym1, sym2) >= 0);

  4152                delete_symbol(mpl, sym1);

  4153                delete_symbol(mpl, sym2);

  4154             }

  4155 #endif

  4156             break;

  4157          case O_GT:

  4158             /* comparison on 'greater than' */

  4159 #if 0 /* 02/VIII-2008 */

  4160             value = (eval_numeric(mpl, code->arg.arg.x) >

  4161                      eval_numeric(mpl, code->arg.arg.y));

  4162 #else

  4163             xassert(code->arg.arg.x != NULL);

  4164             if (code->arg.arg.x->type == A_NUMERIC)

  4165                value = (eval_numeric(mpl, code->arg.arg.x) >

  4166                         eval_numeric(mpl, code->arg.arg.y));

  4167             else

  4168             {  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);

  4169                SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);

  4170                value = (compare_symbols(mpl, sym1, sym2) > 0);

  4171                delete_symbol(mpl, sym1);

  4172                delete_symbol(mpl, sym2);

  4173             }

  4174 #endif

  4175             break;

  4176          case O_NE:

  4177             /* comparison on 'not equal to' */

  4178             xassert(code->arg.arg.x != NULL);

  4179             if (code->arg.arg.x->type == A_NUMERIC)

  4180                value = (eval_numeric(mpl, code->arg.arg.x) !=

  4181                         eval_numeric(mpl, code->arg.arg.y));

  4182             else

  4183             {  SYMBOL *sym1 = eval_symbolic(mpl, code->arg.arg.x);

  4184                SYMBOL *sym2 = eval_symbolic(mpl, code->arg.arg.y);

  4185                value = (compare_symbols(mpl, sym1, sym2) != 0);

  4186                delete_symbol(mpl, sym1);

  4187                delete_symbol(mpl, sym2);

  4188             }

  4189             break;

  4190          case O_AND:

  4191             /* conjunction (logical "and") */

  4192             value = eval_logical(mpl, code->arg.arg.x) &&

  4193                     eval_logical(mpl, code->arg.arg.y);

  4194             break;

  4195          case O_OR:

  4196             /* disjunction (logical "or") */

  4197             value = eval_logical(mpl, code->arg.arg.x) ||

  4198                     eval_logical(mpl, code->arg.arg.y);

  4199             break;

  4200          case O_IN:

  4201             /* test on 'x in Y' */

  4202             {  TUPLE *tuple;

  4203                tuple = eval_tuple(mpl, code->arg.arg.x);

  4204                value = is_member(mpl, code->arg.arg.y, tuple);

  4205                delete_tuple(mpl, tuple);

  4206             }

  4207             break;

  4208          case O_NOTIN:

  4209             /* test on 'x not in Y' */

  4210             {  TUPLE *tuple;

  4211                tuple = eval_tuple(mpl, code->arg.arg.x);

  4212                value = !is_member(mpl, code->arg.arg.y, tuple);

  4213                delete_tuple(mpl, tuple);

  4214             }

  4215             break;

  4216          case O_WITHIN:

  4217             /* test on 'X within Y' */

  4218             {  ELEMSET *set;

  4219                MEMBER *memb;

  4220                set = eval_elemset(mpl, code->arg.arg.x);

  4221                value = 1;

  4222                for (memb = set->head; memb != NULL; memb = memb->next)

  4223                {  if (!is_member(mpl, code->arg.arg.y, memb->tuple))

  4224                   {  value = 0;

  4225                      break;

  4226                   }

  4227                }

  4228                delete_elemset(mpl, set);

  4229             }

  4230             break;

  4231          case O_NOTWITHIN:

  4232             /* test on 'X not within Y' */

  4233             {  ELEMSET *set;

  4234                MEMBER *memb;

  4235                set = eval_elemset(mpl, code->arg.arg.x);

  4236                value = 1;

  4237                for (memb = set->head; memb != NULL; memb = memb->next)

  4238                {  if (is_member(mpl, code->arg.arg.y, memb->tuple))

  4239                   {  value = 0;

  4240                      break;

  4241                   }

  4242                }

  4243                delete_elemset(mpl, set);

  4244             }

  4245             break;

  4246          case O_FORALL:

  4247             /* conjunction (A-quantification) */

  4248             {  struct iter_log_info _info, *info = &_info;

  4249                info->code = code;

  4250                info->value = 1;

  4251                loop_within_domain(mpl, code->arg.loop.domain, info,

  4252                   iter_log_func);

  4253                value = info->value;

  4254             }

  4255             break;

  4256          case O_EXISTS:

  4257             /* disjunction (E-quantification) */

  4258             {  struct iter_log_info _info, *info = &_info;

  4259                info->code = code;

  4260                info->value = 0;

  4261                loop_within_domain(mpl, code->arg.loop.domain, info,

  4262                   iter_log_func);

  4263                value = info->value;

  4264             }

  4265             break;

  4266          default:

  4267             xassert(code != code);

  4268       }

  4269       /* save resultant value */

  4270       xassert(!code->valid);

  4271       code->valid = 1;

  4272       code->value.bit = value;

  4273 done: return value;

  4274 }

  4276 /*----------------------------------------------------------------------

  4277 -- eval_tuple - evaluate pseudo-code to construct n-tuple.

  4278 --

  4279 -- This routine evaluates specified pseudo-code to construct resultant

  4280 -- n-tuple, which is returned on exit. */

  4282 TUPLE *eval_tuple(MPL *mpl, CODE *code)

  4283 {     TUPLE *value;

  4284       xassert(code != NULL);

  4285       xassert(code->type == A_TUPLE);

  4286       xassert(code->dim > 0);

  4287       /* if the operation has a side effect, invalidate and delete the

  4288          resultant value */

  4289       if (code->vflag && code->valid)

  4290       {  code->valid = 0;

  4291          delete_value(mpl, code->type, &code->value);

  4292       }

  4293       /* if resultant value is valid, no evaluation is needed */

  4294       if (code->valid)

  4295       {  value = copy_tuple(mpl, code->value.tuple);

  4296          goto done;

  4297       }

  4298       /* evaluate pseudo-code recursively */

  4299       switch (code->op)

  4300       {  case O_TUPLE:

  4301             /* make n-tuple */

  4302             {  ARG_LIST *e;

  4303                value = create_tuple(mpl);

  4304                for (e = code->arg.list; e != NULL; e = e->next)

  4305                   value = expand_tuple(mpl, value, eval_symbolic(mpl,

  4306                      e->x));

  4307             }

  4308             break;

  4309          case O_CVTTUP:

  4310             /* convert to 1-tuple */

  4311             value = expand_tuple(mpl, create_tuple(mpl),

  4312                eval_symbolic(mpl, code->arg.arg.x));

  4313             break;

  4314          default:

  4315             xassert(code != code);

  4316       }

  4317       /* save resultant value */

  4318       xassert(!code->valid);

  4319       code->valid = 1;

  4320       code->value.tuple = copy_tuple(mpl, value);

  4321 done: return value;

  4322 }

  4324 /*----------------------------------------------------------------------

  4325 -- eval_elemset - evaluate pseudo-code to construct elemental set.

  4326 --

  4327 -- This routine evaluates specified pseudo-code to construct resultant

  4328 -- elemental set, which is returned on exit. */

  4330 struct iter_set_info

  4331 {     /* working info used by the routine iter_set_func */

  4332       CODE *code;

  4333       /* pseudo-code for iterated operation to be performed */

  4334       ELEMSET *value;

  4335       /* resultant value */

  4336 };

  4338 static int iter_set_func(MPL *mpl, void *_info)

  4339 {     /* this is auxiliary routine used to perform iterated operation

  4340          on n-tuple "integrand" within domain scope */

  4341       struct iter_set_info *info = _info;

  4342       TUPLE *tuple;

  4343       switch (info->code->op)

  4344       {  case O_SETOF:

  4345             /* compute next n-tuple and add it to the set; in this case

  4346                duplicate n-tuples are silently ignored */

  4347             tuple = eval_tuple(mpl, info->code->arg.loop.x);

  4348             if (find_tuple(mpl, info->value, tuple) == NULL)

  4349                add_tuple(mpl, info->value, tuple);

  4350             else

  4351                delete_tuple(mpl, tuple);

  4352             break;

  4353          case O_BUILD:

  4354             /* construct next n-tuple using current values assigned to

  4355                *free* dummy indices as its components and add it to the

  4356                set; in this case duplicate n-tuples cannot appear */

  4357             add_tuple(mpl, info->value, get_domain_tuple(mpl,

  4358                info->code->arg.loop.domain));

  4359             break;

  4360          default:

  4361             xassert(info != info);

  4362       }

  4363       return 0;

  4364 }

  4366 ELEMSET *eval_elemset(MPL *mpl, CODE *code)

  4367 {     ELEMSET *value;

  4368       xassert(code != NULL);

  4369       xassert(code->type == A_ELEMSET);

  4370       xassert(code->dim > 0);

  4371       /* if the operation has a side effect, invalidate and delete the

  4372          resultant value */

  4373       if (code->vflag && code->valid)

  4374       {  code->valid = 0;

  4375          delete_value(mpl, code->type, &code->value);

  4376       }

  4377       /* if resultant value is valid, no evaluation is needed */

  4378       if (code->valid)

  4379       {  value = copy_elemset(mpl, code->value.set);

  4380          goto done;

  4381       }

  4382       /* evaluate pseudo-code recursively */

  4383       switch (code->op)

  4384       {  case O_MEMSET:

  4385             /* take member of set */

  4386             {  TUPLE *tuple;

  4387                ARG_LIST *e;

  4388                tuple = create_tuple(mpl);

  4389                for (e = code->arg.set.list; e != NULL; e = e->next)

  4390                   tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,

  4391                      e->x));

  4392                value = copy_elemset(mpl,

  4393                   eval_member_set(mpl, code->arg.set.set, tuple));

  4394                delete_tuple(mpl, tuple);

  4395             }

  4396             break;

  4397          case O_MAKE:

  4398             /* make elemental set of n-tuples */

  4399             {  ARG_LIST *e;

  4400                value = create_elemset(mpl, code->dim);

  4401                for (e = code->arg.list; e != NULL; e = e->next)

  4402                   check_then_add(mpl, value, eval_tuple(mpl, e->x));

  4403             }

  4404             break;

  4405          case O_UNION:

  4406             /* union of two elemental sets */

  4407             value = set_union(mpl,

  4408                eval_elemset(mpl, code->arg.arg.x),

  4409                eval_elemset(mpl, code->arg.arg.y));

  4410             break;

  4411          case O_DIFF:

  4412             /* difference between two elemental sets */

  4413             value = set_diff(mpl,

  4414                eval_elemset(mpl, code->arg.arg.x),

  4415                eval_elemset(mpl, code->arg.arg.y));

  4416             break;

  4417          case O_SYMDIFF:

  4418             /* symmetric difference between two elemental sets */

  4419             value = set_symdiff(mpl,

  4420                eval_elemset(mpl, code->arg.arg.x),

  4421                eval_elemset(mpl, code->arg.arg.y));

  4422             break;

  4423          case O_INTER:

  4424             /* intersection of two elemental sets */

  4425             value = set_inter(mpl,

  4426                eval_elemset(mpl, code->arg.arg.x),

  4427                eval_elemset(mpl, code->arg.arg.y));

  4428             break;

  4429          case O_CROSS:

  4430             /* cross (Cartesian) product of two elemental sets */

  4431             value = set_cross(mpl,

  4432                eval_elemset(mpl, code->arg.arg.x),

  4433                eval_elemset(mpl, code->arg.arg.y));

  4434             break;

  4435          case O_DOTS:

  4436             /* build "arithmetic" elemental set */

  4437             value = create_arelset(mpl,

  4438                eval_numeric(mpl, code->arg.arg.x),

  4439                eval_numeric(mpl, code->arg.arg.y),

  4440                code->arg.arg.z == NULL ? 1.0 : eval_numeric(mpl,

  4441                   code->arg.arg.z));

  4442             break;

  4443          case O_FORK:

  4444             /* if-then-else */

  4445             if (eval_logical(mpl, code->arg.arg.x))

  4446                value = eval_elemset(mpl, code->arg.arg.y);

  4447             else

  4448                value = eval_elemset(mpl, code->arg.arg.z);

  4449             break;

  4450          case O_SETOF:

  4451             /* compute elemental set */

  4452             {  struct iter_set_info _info, *info = &_info;

  4453                info->code = code;

  4454                info->value = create_elemset(mpl, code->dim);

  4455                loop_within_domain(mpl, code->arg.loop.domain, info,

  4456                   iter_set_func);

  4457                value = info->value;

  4458             }

  4459             break;

  4460          case O_BUILD:

  4461             /* build elemental set identical to domain set */

  4462             {  struct iter_set_info _info, *info = &_info;

  4463                info->code = code;

  4464                info->value = create_elemset(mpl, code->dim);

  4465                loop_within_domain(mpl, code->arg.loop.domain, info,

  4466                   iter_set_func);

  4467                value = info->value;

  4468             }

  4469             break;

  4470          default:

  4471             xassert(code != code);

  4472       }

  4473       /* save resultant value */

  4474       xassert(!code->valid);

  4475       code->valid = 1;

  4476       code->value.set = copy_elemset(mpl, value);

  4477 done: return value;

  4478 }

  4480 /*----------------------------------------------------------------------

  4481 -- is_member - check if n-tuple is in set specified by pseudo-code.

  4482 --

  4483 -- This routine checks if given n-tuple is a member of elemental set

  4484 -- specified in the form of pseudo-code (i.e. by expression).

  4485 --

  4486 -- The n-tuple may have more components that dimension of the elemental

  4487 -- set, in which case the extra components are ignored. */

  4489 static void null_func(MPL *mpl, void *info)

  4490 {     /* this is dummy routine used to enter the domain scope */

  4491       xassert(mpl == mpl);

  4492       xassert(info == NULL);

  4493       return;

  4494 }

  4496 int is_member(MPL *mpl, CODE *code, TUPLE *tuple)

  4497 {     int value;

  4498       xassert(code != NULL);

  4499       xassert(code->type == A_ELEMSET);

  4500       xassert(code->dim > 0);

  4501       xassert(tuple != NULL);

  4502       switch (code->op)

  4503       {  case O_MEMSET:

  4504             /* check if given n-tuple is member of elemental set, which

  4505                is assigned to member of model set */

  4506             {  ARG_LIST *e;

  4507                TUPLE *temp;

  4508                ELEMSET *set;

  4509                /* evaluate reference to elemental set */

  4510                temp = create_tuple(mpl);

  4511                for (e = code->arg.set.list; e != NULL; e = e->next)

  4512                   temp = expand_tuple(mpl, temp, eval_symbolic(mpl,

  4513                      e->x));

  4514                set = eval_member_set(mpl, code->arg.set.set, temp);

  4515                delete_tuple(mpl, temp);

  4516                /* check if the n-tuple is contained in the set array */

  4517                temp = build_subtuple(mpl, tuple, set->dim);

  4518                value = (find_tuple(mpl, set, temp) != NULL);

  4519                delete_tuple(mpl, temp);

  4520             }

  4521             break;

  4522          case O_MAKE:

  4523             /* check if given n-tuple is member of literal set */

  4524             {  ARG_LIST *e;

  4525                TUPLE *temp, *that;

  4526                value = 0;

  4527                temp = build_subtuple(mpl, tuple, code->dim);

  4528                for (e = code->arg.list; e != NULL; e = e->next)

  4529                {  that = eval_tuple(mpl, e->x);

  4530                   value = (compare_tuples(mpl, temp, that) == 0);

  4531                   delete_tuple(mpl, that);

  4532                   if (value) break;

  4533                }

  4534                delete_tuple(mpl, temp);

  4535             }

  4536             break;

  4537          case O_UNION:

  4538             value = is_member(mpl, code->arg.arg.x, tuple) ||

  4539                     is_member(mpl, code->arg.arg.y, tuple);

  4540             break;

  4541          case O_DIFF:

  4542             value = is_member(mpl, code->arg.arg.x, tuple) &&

  4543                    !is_member(mpl, code->arg.arg.y, tuple);

  4544             break;

  4545          case O_SYMDIFF:

  4546             {  int in1 = is_member(mpl, code->arg.arg.x, tuple);

  4547                int in2 = is_member(mpl, code->arg.arg.y, tuple);

  4548                value = (in1 && !in2) || (!in1 && in2);

  4549             }

  4550             break;

  4551          case O_INTER:

  4552             value = is_member(mpl, code->arg.arg.x, tuple) &&

  4553                     is_member(mpl, code->arg.arg.y, tuple);

  4554             break;

  4555          case O_CROSS:

  4556             {  int j;

  4557                value = is_member(mpl, code->arg.arg.x, tuple);

  4558                if (value)

  4559                {  for (j = 1; j <= code->arg.arg.x->dim; j++)

  4560                   {  xassert(tuple != NULL);

  4561                      tuple = tuple->next;

  4562                   }

  4563                   value = is_member(mpl, code->arg.arg.y, tuple);

  4564                }

  4565             }

  4566             break;

  4567          case O_DOTS:

  4568             /* check if given 1-tuple is member of "arithmetic" set */

  4569             {  int j;

  4570                double x, t0, tf, dt;

  4571                xassert(code->dim == 1);

  4572                /* compute "parameters" of the "arithmetic" set */

  4573                t0 = eval_numeric(mpl, code->arg.arg.x);

  4574                tf = eval_numeric(mpl, code->arg.arg.y);

  4575                if (code->arg.arg.z == NULL)

  4576                   dt = 1.0;

  4577                else

  4578                   dt = eval_numeric(mpl, code->arg.arg.z);

  4579                /* make sure the parameters are correct */

  4580                arelset_size(mpl, t0, tf, dt);

  4581                /* if component of 1-tuple is symbolic, not numeric, the

  4582                   1-tuple cannot be member of "arithmetic" set */

  4583                xassert(tuple->sym != NULL);

  4584                if (tuple->sym->str != NULL)

  4585                {  value = 0;

  4586                   break;

  4587                }

  4588                /* determine numeric value of the component */

  4589                x = tuple->sym->num;

  4590                /* if the component value is out of the set range, the

  4591                   1-tuple is not in the set */

  4592                if (dt > 0.0 && !(t0 <= x && x <= tf) ||

  4593                    dt < 0.0 && !(tf <= x && x <= t0))

  4594                {  value = 0;

  4595                   break;

  4596                }

  4597                /* estimate ordinal number of the 1-tuple in the set */

  4598                j = (int)(((x - t0) / dt) + 0.5) + 1;

  4599                /* perform the main check */

  4600                value = (arelset_member(mpl, t0, tf, dt, j) == x);

  4601             }

  4602             break;

  4603          case O_FORK:

  4604             /* check if given n-tuple is member of conditional set */

  4605             if (eval_logical(mpl, code->arg.arg.x))

  4606                value = is_member(mpl, code->arg.arg.y, tuple);

  4607             else

  4608                value = is_member(mpl, code->arg.arg.z, tuple);

  4609             break;

  4610          case O_SETOF:

  4611             /* check if given n-tuple is member of computed set */

  4612             /* it is not clear how to efficiently perform the check not

  4613                computing the entire elemental set :+( */

  4614             error(mpl, "implementation restriction; in/within setof{} n"

  4615                "ot allowed");

  4616             break;

  4617          case O_BUILD:

  4618             /* check if given n-tuple is member of domain set */

  4619             {  TUPLE *temp;

  4620                temp = build_subtuple(mpl, tuple, code->dim);

  4621                /* try to enter the domain scope; if it is successful,

  4622                   the n-tuple is in the domain set */

  4623                value = (eval_within_domain(mpl, code->arg.loop.domain,

  4624                   temp, NULL, null_func) == 0);

  4625                delete_tuple(mpl, temp);

  4626             }

  4627             break;

  4628          default:

  4629             xassert(code != code);

  4630       }

  4631       return value;

  4632 }

  4634 /*----------------------------------------------------------------------

  4635 -- eval_formula - evaluate pseudo-code to construct linear form.

  4636 --

  4637 -- This routine evaluates specified pseudo-code to construct resultant

  4638 -- linear form, which is returned on exit. */

  4640 struct iter_form_info

  4641 {     /* working info used by the routine iter_form_func */

  4642       CODE *code;

  4643       /* pseudo-code for iterated operation to be performed */

  4644       FORMULA *value;

  4645       /* resultant value */

  4646       FORMULA *tail;

  4647       /* pointer to the last term */

  4648 };

  4650 static int iter_form_func(MPL *mpl, void *_info)

  4651 {     /* this is auxiliary routine used to perform iterated operation

  4652          on linear form "integrand" within domain scope */

  4653       struct iter_form_info *info = _info;

  4654       switch (info->code->op)

  4655       {  case O_SUM:

  4656             /* summation over domain */

  4657 #if 0

  4658             info->value =

  4659                linear_comb(mpl,

  4660                   +1.0, info->value,

  4661                   +1.0, eval_formula(mpl, info->code->arg.loop.x));

  4662 #else

  4663             /* the routine linear_comb needs to look through all terms

  4664                of both linear forms to reduce identical terms, so using

  4665                it here is not a good idea (for example, evaluation of

  4666                sum{i in 1..n} x[i] required quadratic time); the better

  4667                idea is to gather all terms of the integrand in one list

  4668                and reduce identical terms only once after all terms of

  4669                the resultant linear form have been evaluated */

  4670             {  FORMULA *form, *term;

  4671                form = eval_formula(mpl, info->code->arg.loop.x);

  4672                if (info->value == NULL)

  4673                {  xassert(info->tail == NULL);

  4674                   info->value = form;

  4675                }

  4676                else

  4677                {  xassert(info->tail != NULL);

  4678                   info->tail->next = form;

  4679                }

  4680                for (term = form; term != NULL; term = term->next)

  4681                   info->tail = term;

  4682             }

  4683 #endif

  4684             break;

  4685          default:

  4686             xassert(info != info);

  4687       }

  4688       return 0;

  4689 }

  4691 FORMULA *eval_formula(MPL *mpl, CODE *code)

  4692 {     FORMULA *value;

  4693       xassert(code != NULL);

  4694       xassert(code->type == A_FORMULA);

  4695       xassert(code->dim == 0);

  4696       /* if the operation has a side effect, invalidate and delete the

  4697          resultant value */

  4698       if (code->vflag && code->valid)

  4699       {  code->valid = 0;

  4700          delete_value(mpl, code->type, &code->value);

  4701       }

  4702       /* if resultant value is valid, no evaluation is needed */

  4703       if (code->valid)

  4704       {  value = copy_formula(mpl, code->value.form);

  4705          goto done;

  4706       }

  4707       /* evaluate pseudo-code recursively */

  4708       switch (code->op)

  4709       {  case O_MEMVAR:

  4710             /* take member of variable */

  4711             {  TUPLE *tuple;

  4712                ARG_LIST *e;

  4713                tuple = create_tuple(mpl);

  4714                for (e = code->arg.var.list; e != NULL; e = e->next)

  4715                   tuple = expand_tuple(mpl, tuple, eval_symbolic(mpl,

  4716                      e->x));

  4717 #if 1 /* 15/V-2010 */

  4718                xassert(code->arg.var.suff == DOT_NONE);

  4719 #endif

  4720                value = single_variable(mpl,

  4721                   eval_member_var(mpl, code->arg.var.var, tuple));

  4722                delete_tuple(mpl, tuple);

  4723             }

  4724             break;

  4725          case O_CVTLFM:

  4726             /* convert to linear form */

  4727             value = constant_term(mpl, eval_numeric(mpl,

  4728                code->arg.arg.x));

  4729             break;

  4730          case O_PLUS:

  4731             /* unary plus */

  4732             value = linear_comb(mpl,

  4733                 0.0, constant_term(mpl, 0.0),

  4734                +1.0, eval_formula(mpl, code->arg.arg.x));

  4735             break;

  4736          case O_MINUS:

  4737             /* unary minus */

  4738             value = linear_comb(mpl,

  4739                 0.0, constant_term(mpl, 0.0),

  4740                -1.0, eval_formula(mpl, code->arg.arg.x));

  4741             break;

  4742          case O_ADD:

  4743             /* addition */

  4744             value = linear_comb(mpl,

  4745                +1.0, eval_formula(mpl, code->arg.arg.x),

  4746                +1.0, eval_formula(mpl, code->arg.arg.y));

  4747             break;

  4748          case O_SUB:

  4749             /* subtraction */

  4750             value = linear_comb(mpl,

  4751                +1.0, eval_formula(mpl, code->arg.arg.x),

  4752                -1.0, eval_formula(mpl, code->arg.arg.y));

  4753             break;

  4754          case O_MUL:

  4755             /* multiplication */

  4756             xassert(code->arg.arg.x != NULL);

  4757             xassert(code->arg.arg.y != NULL);

  4758             if (code->arg.arg.x->type == A_NUMERIC)

  4759             {  xassert(code->arg.arg.y->type == A_FORMULA);

  4760                value = linear_comb(mpl,

  4761                   eval_numeric(mpl, code->arg.arg.x),

  4762                   eval_formula(mpl, code->arg.arg.y),

  4763                   0.0, constant_term(mpl, 0.0));

  4764             }

  4765             else

  4766             {  xassert(code->arg.arg.x->type == A_FORMULA);

  4767                xassert(code->arg.arg.y->type == A_NUMERIC);

  4768                value = linear_comb(mpl,

  4769                   eval_numeric(mpl, code->arg.arg.y),

  4770                   eval_formula(mpl, code->arg.arg.x),

  4771                   0.0, constant_term(mpl, 0.0));

  4772             }

  4773             break;

  4774          case O_DIV:

  4775             /* division */

  4776             value = linear_comb(mpl,

  4777                fp_div(mpl, 1.0, eval_numeric(mpl, code->arg.arg.y)),

  4778                eval_formula(mpl, code->arg.arg.x),

  4779                0.0, constant_term(mpl, 0.0));

  4780             break;

  4781          case O_FORK:

  4782             /* if-then-else */

  4783             if (eval_logical(mpl, code->arg.arg.x))

  4784                value = eval_formula(mpl, code->arg.arg.y);

  4785             else if (code->arg.arg.z == NULL)

  4786                value = constant_term(mpl, 0.0);

  4787             else

  4788                value = eval_formula(mpl, code->arg.arg.z);

  4789             break;

  4790          case O_SUM:

  4791             /* summation over domain */

  4792             {  struct iter_form_info _info, *info = &_info;

  4793                info->code = code;

  4794                info->value = constant_term(mpl, 0.0);

  4795                info->tail = NULL;

  4796                loop_within_domain(mpl, code->arg.loop.domain, info,

  4797                   iter_form_func);

  4798                value = reduce_terms(mpl, info->value);

  4799             }

  4800             break;

  4801          default:

  4802             xassert(code != code);

  4803       }

  4804       /* save resultant value */

  4805       xassert(!code->valid);

  4806       code->valid = 1;

  4807       code->value.form = copy_formula(mpl, value);

  4808 done: return value;

  4809 }

  4811 /*----------------------------------------------------------------------

  4812 -- clean_code - clean pseudo-code.

  4813 --

  4814 -- This routine recursively cleans specified pseudo-code that assumes

  4815 -- deleting all temporary resultant values. */

  4817 void clean_code(MPL *mpl, CODE *code)

  4818 {     ARG_LIST *e;

  4819       /* if no pseudo-code is specified, do nothing */

  4820       if (code == NULL) goto done;

  4821       /* if resultant value is valid (exists), delete it */

  4822       if (code->valid)

  4823       {  code->valid = 0;

  4824          delete_value(mpl, code->type, &code->value);

  4825       }

  4826       /* recursively clean pseudo-code for operands */

  4827       switch (code->op)

  4828       {  case O_NUMBER:

  4829          case O_STRING:

  4830          case O_INDEX:

  4831             break;

  4832          case O_MEMNUM:

  4833          case O_MEMSYM:

  4834             for (e = code->arg.par.list; e != NULL; e = e->next)

  4835                clean_code(mpl, e->x);

  4836             break;

  4837          case O_MEMSET:

  4838             for (e = code->arg.set.list; e != NULL; e = e->next)

  4839                clean_code(mpl, e->x);

  4840             break;

  4841          case O_MEMVAR:

  4842             for (e = code->arg.var.list; e != NULL; e = e->next)

  4843                clean_code(mpl, e->x);

  4844             break;

  4845 #if 1 /* 15/V-2010 */

  4846          case O_MEMCON:

  4847             for (e = code->arg.con.list; e != NULL; e = e->next)

  4848                clean_code(mpl, e->x);

  4849             break;

  4850 #endif

  4851          case O_TUPLE:

  4852          case O_MAKE:

  4853             for (e = code->arg.list; e != NULL; e = e->next)

  4854                clean_code(mpl, e->x);

  4855             break;

  4856          case O_SLICE:

  4857             xassert(code != code);

  4858          case O_IRAND224:

  4859          case O_UNIFORM01:

  4860          case O_NORMAL01:

  4861          case O_GMTIME:

  4862             break;

  4863          case O_CVTNUM:

  4864          case O_CVTSYM:

  4865          case O_CVTLOG:

  4866          case O_CVTTUP:

  4867          case O_CVTLFM:

  4868          case O_PLUS:

  4869          case O_MINUS:

  4870          case O_NOT:

  4871          case O_ABS:

  4872          case O_CEIL:

  4873          case O_FLOOR:

  4874          case O_EXP:

  4875          case O_LOG:

  4876          case O_LOG10:

  4877          case O_SQRT:

  4878          case O_SIN:

  4879          case O_COS:

  4880          case O_ATAN:

  4881          case O_ROUND:

  4882          case O_TRUNC:

  4883          case O_CARD:

  4884          case O_LENGTH:

  4885             /* unary operation */

  4886             clean_code(mpl, code->arg.arg.x);

  4887             break;

  4888          case O_ADD:

  4889          case O_SUB:

  4890          case O_LESS:

  4891          case O_MUL:

  4892          case O_DIV:

  4893          case O_IDIV:

  4894          case O_MOD:

  4895          case O_POWER:

  4896          case O_ATAN2:

  4897          case O_ROUND2:

  4898          case O_TRUNC2:

  4899          case O_UNIFORM:

  4900          case O_NORMAL:

  4901          case O_CONCAT:

  4902          case O_LT:

  4903          case O_LE:

  4904          case O_EQ:

  4905          case O_GE:

  4906          case O_GT:

  4907          case O_NE:

  4908          case O_AND:

  4909          case O_OR:

  4910          case O_UNION:

  4911          case O_DIFF:

  4912          case O_SYMDIFF:

  4913          case O_INTER:

  4914          case O_CROSS:

  4915          case O_IN:

  4916          case O_NOTIN:

  4917          case O_WITHIN:

  4918          case O_NOTWITHIN:

  4919          case O_SUBSTR:

  4920          case O_STR2TIME:

  4921          case O_TIME2STR:

  4922             /* binary operation */

  4923             clean_code(mpl, code->arg.arg.x);

  4924             clean_code(mpl, code->arg.arg.y);

  4925             break;

  4926          case O_DOTS:

  4927          case O_FORK:

  4928          case O_SUBSTR3:

  4929             /* ternary operation */

  4930             clean_code(mpl, code->arg.arg.x);

  4931             clean_code(mpl, code->arg.arg.y);

  4932             clean_code(mpl, code->arg.arg.z);

  4933             break;

  4934          case O_MIN:

  4935          case O_MAX:

  4936             /* n-ary operation */

  4937             for (e = code->arg.list; e != NULL; e = e->next)

  4938                clean_code(mpl, e->x);

  4939             break;

  4940          case O_SUM:

  4941          case O_PROD:

  4942          case O_MINIMUM:

  4943          case O_MAXIMUM:

  4944          case O_FORALL:

  4945          case O_EXISTS:

  4946          case O_SETOF:

  4947          case O_BUILD:

  4948             /* iterated operation */

  4949             clean_domain(mpl, code->arg.loop.domain);

  4950             clean_code(mpl, code->arg.loop.x);

  4951             break;

  4952          default:

  4953             xassert(code->op != code->op);

  4954       }

  4955 done: return;

  4956 }

  4958 #if 1 /* 11/II-2008 */

  4959 /**********************************************************************/

  4960 /* * *                        DATA TABLES                         * * */

  4961 /**********************************************************************/

  4963 int mpl_tab_num_args(TABDCA *dca)

  4964 {     /* returns the number of arguments */

  4965       return dca->na;

  4966 }

  4968 const char *mpl_tab_get_arg(TABDCA *dca, int k)

  4969 {     /* returns pointer to k-th argument */

  4970       xassert(1 <= k && k <= dca->na);

  4971       return dca->arg[k];

  4972 }

  4974 int mpl_tab_num_flds(TABDCA *dca)

  4975 {     /* returns the number of fields */

  4976       return dca->nf;

  4977 }

  4979 const char *mpl_tab_get_name(TABDCA *dca, int k)

  4980 {     /* returns pointer to name of k-th field */

  4981       xassert(1 <= k && k <= dca->nf);

  4982       return dca->name[k];

  4983 }

  4985 int mpl_tab_get_type(TABDCA *dca, int k)

  4986 {     /* returns type of k-th field */

  4987       xassert(1 <= k && k <= dca->nf);

  4988       return dca->type[k];

  4989 }

  4991 double mpl_tab_get_num(TABDCA *dca, int k)

  4992 {     /* returns numeric value of k-th field */

  4993       xassert(1 <= k && k <= dca->nf);

  4994       xassert(dca->type[k] == 'N');

  4995       return dca->num[k];

  4996 }

  4998 const char *mpl_tab_get_str(TABDCA *dca, int k)

  4999 {     /* returns pointer to string value of k-th field */

  5000       xassert(1 <= k && k <= dca->nf);

  5001       xassert(dca->type[k] == 'S');

  5002       xassert(dca->str[k] != NULL);

  5003       return dca->str[k];

  5004 }

  5006 void mpl_tab_set_num(TABDCA *dca, int k, double num)

  5007 {     /* assign numeric value to k-th field */

  5008       xassert(1 <= k && k <= dca->nf);

  5009       xassert(dca->type[k] == '?');

  5010       dca->type[k] = 'N';

  5011       dca->num[k] = num;

  5012       return;

  5013 }

  5015 void mpl_tab_set_str(TABDCA *dca, int k, const char *str)

  5016 {     /* assign string value to k-th field */

  5017       xassert(1 <= k && k <= dca->nf);

  5018       xassert(dca->type[k] == '?');

  5019       xassert(strlen(str) <= MAX_LENGTH);

  5020       xassert(dca->str[k] != NULL);

  5021       dca->type[k] = 'S';

  5022       strcpy(dca->str[k], str);

  5023       return;

  5024 }

  5026 static int write_func(MPL *mpl, void *info)

  5027 {     /* this is auxiliary routine to work within domain scope */

  5028       TABLE *tab = info;

  5029       TABDCA *dca = mpl->dca;

  5030       TABOUT *out;

  5031       SYMBOL *sym;

  5032       int k;

  5033       char buf[MAX_LENGTH+1];

  5034       /* evaluate field values */

  5035       k = 0;

  5036       for (out = tab->u.out.list; out != NULL; out = out->next)

  5037       {  k++;

  5038          switch (out->code->type)

  5039          {  case A_NUMERIC:

  5040                dca->type[k] = 'N';

  5041                dca->num[k] = eval_numeric(mpl, out->code);

  5042                dca->str[k][0] = '\0';

  5043                break;

  5044             case A_SYMBOLIC:

  5045                sym = eval_symbolic(mpl, out->code);

  5046                if (sym->str == NULL)

  5047                {  dca->type[k] = 'N';

  5048                   dca->num[k] = sym->num;

  5049                   dca->str[k][0] = '\0';

  5050                }

  5051                else

  5052                {  dca->type[k] = 'S';

  5053                   dca->num[k] = 0.0;

  5054                   fetch_string(mpl, sym->str, buf);

  5055                   strcpy(dca->str[k], buf);

  5056                }

  5057                delete_symbol(mpl, sym);

  5058                break;

  5059             default:

  5060                xassert(out != out);

  5061          }

  5062       }

  5063       /* write record to output table */

  5064       mpl_tab_drv_write(mpl);

  5065       return 0;

  5066 }

  5068 void execute_table(MPL *mpl, TABLE *tab)

  5069 {     /* execute table statement */

  5070       TABARG *arg;

  5071       TABFLD *fld;

  5072       TABIN *in;

  5073       TABOUT *out;

  5074       TABDCA *dca;

  5075       SET *set;

  5076       int k;

  5077       char buf[MAX_LENGTH+1];

  5078       /* allocate table driver communication area */

  5079       xassert(mpl->dca == NULL);

  5080       mpl->dca = dca = xmalloc(sizeof(TABDCA));

  5081       dca->id = 0;

  5082       dca->link = NULL;

  5083       dca->na = 0;

  5084       dca->arg = NULL;

  5085       dca->nf = 0;

  5086       dca->name = NULL;

  5087       dca->type = NULL;

  5088       dca->num = NULL;

  5089       dca->str = NULL;

  5090       /* allocate arguments */

  5091       xassert(dca->na == 0);

  5092       for (arg = tab->arg; arg != NULL; arg = arg->next)

  5093          dca->na++;

  5094       dca->arg = xcalloc(1+dca->na, sizeof(char *));

  5095 #if 1 /* 28/IX-2008 */

  5096       for (k = 1; k <= dca->na; k++) dca->arg[k] = NULL;

  5097 #endif

  5098       /* evaluate argument values */

  5099       k = 0;

  5100       for (arg = tab->arg; arg != NULL; arg = arg->next)

  5101       {  SYMBOL *sym;

  5102          k++;

  5103          xassert(arg->code->type == A_SYMBOLIC);

  5104          sym = eval_symbolic(mpl, arg->code);

  5105          if (sym->str == NULL)

  5106             sprintf(buf, "%.*g", DBL_DIG, sym->num);

  5107          else

  5108             fetch_string(mpl, sym->str, buf);

  5109          delete_symbol(mpl, sym);

  5110          dca->arg[k] = xmalloc(strlen(buf)+1);

  5111          strcpy(dca->arg[k], buf);

  5112       }

  5113       /* perform table input/output */

  5114       switch (tab->type)

  5115       {  case A_INPUT:  goto read_table;

  5116          case A_OUTPUT: goto write_table;

  5117          default:       xassert(tab != tab);

  5118       }

  5119 read_table:

  5120       /* read data from input table */

  5121       /* add the only member to the control set and assign it empty

  5122          elemental set */

  5123       set = tab->u.in.set;

  5124       if (set != NULL)

  5125       {  if (set->data)

  5126             error(mpl, "%s already provided with data", set->name);

  5127          xassert(set->array->head == NULL);

  5128          add_member(mpl, set->array, NULL)->value.set =

  5129             create_elemset(mpl, set->dimen);

  5130          set->data = 1;

  5131       }

  5132       /* check parameters specified in the input list */

  5133       for (in = tab->u.in.list; in != NULL; in = in->next)

  5134       {  if (in->par->data)

  5135             error(mpl, "%s already provided with data", in->par->name);

  5136          in->par->data = 1;

  5137       }

  5138       /* allocate and initialize fields */

  5139       xassert(dca->nf == 0);

  5140       for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)

  5141          dca->nf++;

  5142       for (in = tab->u.in.list; in != NULL; in = in->next)

  5143          dca->nf++;

  5144       dca->name = xcalloc(1+dca->nf, sizeof(char *));

  5145       dca->type = xcalloc(1+dca->nf, sizeof(int));

  5146       dca->num = xcalloc(1+dca->nf, sizeof(double));

  5147       dca->str = xcalloc(1+dca->nf, sizeof(char *));

  5148       k = 0;

  5149       for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)

  5150       {  k++;

  5151          dca->name[k] = fld->name;

  5152          dca->type[k] = '?';

  5153          dca->num[k] = 0.0;

  5154          dca->str[k] = xmalloc(MAX_LENGTH+1);

  5155          dca->str[k][0] = '\0';

  5156       }

  5157       for (in = tab->u.in.list; in != NULL; in = in->next)

  5158       {  k++;

  5159          dca->name[k] = in->name;

  5160          dca->type[k] = '?';

  5161          dca->num[k] = 0.0;

  5162          dca->str[k] = xmalloc(MAX_LENGTH+1);

  5163          dca->str[k][0] = '\0';

  5164       }

  5165       /* open input table */

  5166       mpl_tab_drv_open(mpl, 'R');

  5167       /* read and process records */

  5168       for (;;)

  5169       {  TUPLE *tup;

  5170          /* reset field types */

  5171          for (k = 1; k <= dca->nf; k++)

  5172             dca->type[k] = '?';

  5173          /* read next record */

  5174          if (mpl_tab_drv_read(mpl)) break;

  5175          /* all fields must be set by the driver */

  5176          for (k = 1; k <= dca->nf; k++)

  5177          {  if (dca->type[k] == '?')

  5178                error(mpl, "field %s missing in input table",

  5179                   dca->name[k]);

  5180          }

  5181          /* construct n-tuple */

  5182          tup = create_tuple(mpl);

  5183          k = 0;

  5184          for (fld = tab->u.in.fld; fld != NULL; fld = fld->next)

  5185          {  k++;

  5186             xassert(k <= dca->nf);

  5187             switch (dca->type[k])

  5188             {  case 'N':

  5189                   tup = expand_tuple(mpl, tup, create_symbol_num(mpl,

  5190                      dca->num[k]));

  5191                   break;

  5192                case 'S':

  5193                   xassert(strlen(dca->str[k]) <= MAX_LENGTH);

  5194                   tup = expand_tuple(mpl, tup, create_symbol_str(mpl,

  5195                      create_string(mpl, dca->str[k])));

  5196                   break;

  5197                default:

  5198                   xassert(dca != dca);

  5199             }

  5200          }

  5201          /* add n-tuple just read to the control set */

  5202          if (tab->u.in.set != NULL)

  5203             check_then_add(mpl, tab->u.in.set->array->head->value.set,

  5204                copy_tuple(mpl, tup));

  5205          /* assign values to the parameters in the input list */

  5206          for (in = tab->u.in.list; in != NULL; in = in->next)

  5207          {  MEMBER *memb;

  5208             k++;

  5209             xassert(k <= dca->nf);

  5210             /* there must be no member with the same n-tuple */

  5211             if (find_member(mpl, in->par->array, tup) != NULL)

  5212                error(mpl, "%s%s already defined", in->par->name,

  5213                format_tuple(mpl, '[', tup));

  5214             /* create new parameter member with given n-tuple */

  5215             memb = add_member(mpl, in->par->array, copy_tuple(mpl, tup))

  5216                ;

  5217             /* assign value to the parameter member */

  5218             switch (in->par->type)

  5219             {  case A_NUMERIC:

  5220                case A_INTEGER:

  5221                case A_BINARY:

  5222                   if (dca->type[k] != 'N')

  5223                      error(mpl, "%s requires numeric data",

  5224                         in->par->name);

  5225                   memb->value.num = dca->num[k];

  5226                   break;

  5227                case A_SYMBOLIC:

  5228                   switch (dca->type[k])

  5229                   {  case 'N':

  5230                         memb->value.sym = create_symbol_num(mpl,

  5231                            dca->num[k]);

  5232                         break;

  5233                      case 'S':

  5234                         xassert(strlen(dca->str[k]) <= MAX_LENGTH);

  5235                         memb->value.sym = create_symbol_str(mpl,

  5236                            create_string(mpl,dca->str[k]));

  5237                         break;

  5238                      default:

  5239                         xassert(dca != dca);

  5240                   }

  5241                   break;

  5242                default:

  5243                   xassert(in != in);

  5244             }

  5245          }

  5246          /* n-tuple is no more needed */

  5247          delete_tuple(mpl, tup);

  5248       }

  5249       /* close input table */

  5250       mpl_tab_drv_close(mpl);

  5251       goto done;

  5252 write_table:

  5253       /* write data to output table */

  5254       /* allocate and initialize fields */

  5255       xassert(dca->nf == 0);

  5256       for (out = tab->u.out.list; out != NULL; out = out->next)

  5257          dca->nf++;

  5258       dca->name = xcalloc(1+dca->nf, sizeof(char *));

  5259       dca->type = xcalloc(1+dca->nf, sizeof(int));

  5260       dca->num = xcalloc(1+dca->nf, sizeof(double));

  5261       dca->str = xcalloc(1+dca->nf, sizeof(char *));

  5262       k = 0;

  5263       for (out = tab->u.out.list; out != NULL; out = out->next)

  5264       {  k++;

  5265          dca->name[k] = out->name;

  5266          dca->type[k] = '?';

  5267          dca->num[k] = 0.0;

  5268          dca->str[k] = xmalloc(MAX_LENGTH+1);

  5269          dca->str[k][0] = '\0';

  5270       }

  5271       /* open output table */

  5272       mpl_tab_drv_open(mpl, 'W');

  5273       /* evaluate fields and write records */

  5274       loop_within_domain(mpl, tab->u.out.domain, tab, write_func);

  5275       /* close output table */

  5276       mpl_tab_drv_close(mpl);

  5277 done: /* free table driver communication area */

  5278       free_dca(mpl);

  5279       return;

  5280 }

  5282 void free_dca(MPL *mpl)

  5283 {     /* free table driver communucation area */

  5284       TABDCA *dca = mpl->dca;

  5285       int k;

  5286       if (dca != NULL)

  5287       {  if (dca->link != NULL)

  5288             mpl_tab_drv_close(mpl);

  5289          if (dca->arg != NULL)

  5290          {  for (k = 1; k <= dca->na; k++)

  5291 #if 1 /* 28/IX-2008 */

  5292                if (dca->arg[k] != NULL)

  5293 #endif

  5294                xfree(dca->arg[k]);

  5295             xfree(dca->arg);

  5296          }

  5297          if (dca->name != NULL) xfree(dca->name);

  5298          if (dca->type != NULL) xfree(dca->type);

  5299          if (dca->num != NULL) xfree(dca->num);

  5300          if (dca->str != NULL)

  5301          {  for (k = 1; k <= dca->nf; k++)

  5302                xfree(dca->str[k]);

  5303             xfree(dca->str);

  5304          }

  5305          xfree(dca), mpl->dca = NULL;

  5306       }

  5307       return;

  5308 }

  5310 void clean_table(MPL *mpl, TABLE *tab)

  5311 {     /* clean table statement */

  5312       TABARG *arg;

  5313       TABOUT *out;

  5314       /* clean string list */

  5315       for (arg = tab->arg; arg != NULL; arg = arg->next)

  5316          clean_code(mpl, arg->code);

  5317       switch (tab->type)

  5318       {  case A_INPUT:

  5319             break;

  5320          case A_OUTPUT:

  5321             /* clean subscript domain */

  5322             clean_domain(mpl, tab->u.out.domain);

  5323             /* clean output list */

  5324             for (out = tab->u.out.list; out != NULL; out = out->next)

  5325                clean_code(mpl, out->code);

  5326             break;

  5327          default:

  5328             xassert(tab != tab);

  5329       }

  5330       return;

  5331 }

  5332 #endif

  5334 /**********************************************************************/

  5335 /* * *                      MODEL STATEMENTS                      * * */

  5336 /**********************************************************************/

  5338 /*----------------------------------------------------------------------

  5339 -- execute_check - execute check statement.

  5340 --

  5341 -- This routine executes specified check statement. */

  5343 static int check_func(MPL *mpl, void *info)

  5344 {     /* this is auxiliary routine to work within domain scope */

  5345       CHECK *chk = (CHECK *)info;

  5346       if (!eval_logical(mpl, chk->code))

  5347          error(mpl, "check%s failed", format_tuple(mpl, '[',

  5348             get_domain_tuple(mpl, chk->domain)));

  5349       return 0;

  5350 }

  5352 void execute_check(MPL *mpl, CHECK *chk)

  5353 {     loop_within_domain(mpl, chk->domain, chk, check_func);

  5354       return;

  5355 }

  5357 /*----------------------------------------------------------------------

  5358 -- clean_check - clean check statement.

  5359 --

  5360 -- This routine cleans specified check statement that assumes deleting

  5361 -- all stuff dynamically allocated on generating/postsolving phase. */

  5363 void clean_check(MPL *mpl, CHECK *chk)

  5364 {     /* clean subscript domain */

  5365       clean_domain(mpl, chk->domain);

  5366       /* clean pseudo-code for computing predicate */

  5367       clean_code(mpl, chk->code);

  5368       return;

  5369 }

  5371 /*----------------------------------------------------------------------

  5372 -- execute_display - execute display statement.

  5373 --

  5374 -- This routine executes specified display statement. */

  5376 static void display_set(MPL *mpl, SET *set, MEMBER *memb)

  5377 {     /* display member of model set */

  5378       ELEMSET *s = memb->value.set;

  5379       MEMBER *m;

  5380       write_text(mpl, "%s%s%s\n", set->name,

  5381          format_tuple(mpl, '[', memb->tuple),

  5382          s->head == NULL ? " is empty" : ":");

  5383       for (m = s->head; m != NULL; m = m->next)

  5384          write_text(mpl, "   %s\n", format_tuple(mpl, '(', m->tuple));

  5385       return;

  5386 }

  5388 static void display_par(MPL *mpl, PARAMETER *par, MEMBER *memb)

  5389 {     /* display member of model parameter */

  5390       switch (par->type)

  5391       {  case A_NUMERIC:

  5392          case A_INTEGER:

  5393          case A_BINARY:

  5394             write_text(mpl, "%s%s = %.*g\n", par->name,

  5395                format_tuple(mpl, '[', memb->tuple),

  5396                DBL_DIG, memb->value.num);

  5397             break;

  5398          case A_SYMBOLIC:

  5399             write_text(mpl, "%s%s = %s\n", par->name,

  5400                format_tuple(mpl, '[', memb->tuple),

  5401                format_symbol(mpl, memb->value.sym));

  5402             break;

  5403          default:

  5404             xassert(par != par);

  5405       }

  5406       return;

  5407 }

  5409 #if 1 /* 15/V-2010 */

  5410 static void display_var(MPL *mpl, VARIABLE *var, MEMBER *memb,

  5411       int suff)

  5412 {     /* display member of model variable */

  5413       if (suff == DOT_NONE || suff == DOT_VAL)

  5414          write_text(mpl, "%s%s.val = %.*g\n", var->name,

  5415             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5416             memb->value.var->prim);

  5417       else if (suff == DOT_LB)

  5418          write_text(mpl, "%s%s.lb = %.*g\n", var->name,

  5419             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5420             memb->value.var->var->lbnd == NULL ? -DBL_MAX :

  5421             memb->value.var->lbnd);

  5422       else if (suff == DOT_UB)

  5423          write_text(mpl, "%s%s.ub = %.*g\n", var->name,

  5424             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5425             memb->value.var->var->ubnd == NULL ? +DBL_MAX :

  5426             memb->value.var->ubnd);

  5427       else if (suff == DOT_STATUS)

  5428          write_text(mpl, "%s%s.status = %d\n", var->name, format_tuple

  5429             (mpl, '[', memb->tuple), memb->value.var->stat);

  5430       else if (suff == DOT_DUAL)

  5431          write_text(mpl, "%s%s.dual = %.*g\n", var->name,

  5432             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5433             memb->value.var->dual);

  5434       else

  5435          xassert(suff != suff);

  5436       return;

  5437 }

  5438 #endif

  5440 #if 1 /* 15/V-2010 */

  5441 static void display_con(MPL *mpl, CONSTRAINT *con, MEMBER *memb,

  5442       int suff)

  5443 {     /* display member of model constraint */

  5444       if (suff == DOT_NONE || suff == DOT_VAL)

  5445          write_text(mpl, "%s%s.val = %.*g\n", con->name,

  5446             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5447             memb->value.con->prim);

  5448       else if (suff == DOT_LB)

  5449          write_text(mpl, "%s%s.lb = %.*g\n", con->name,

  5450             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5451             memb->value.con->con->lbnd == NULL ? -DBL_MAX :

  5452             memb->value.con->lbnd);

  5453       else if (suff == DOT_UB)

  5454          write_text(mpl, "%s%s.ub = %.*g\n", con->name,

  5455             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5456             memb->value.con->con->ubnd == NULL ? +DBL_MAX :

  5457             memb->value.con->ubnd);

  5458       else if (suff == DOT_STATUS)

  5459          write_text(mpl, "%s%s.status = %d\n", con->name, format_tuple

  5460             (mpl, '[', memb->tuple), memb->value.con->stat);

  5461       else if (suff == DOT_DUAL)

  5462          write_text(mpl, "%s%s.dual = %.*g\n", con->name,

  5463             format_tuple(mpl, '[', memb->tuple), DBL_DIG,

  5464             memb->value.con->dual);

  5465       else

  5466          xassert(suff != suff);

  5467       return;

  5468 }

  5469 #endif

  5471 static void display_memb(MPL *mpl, CODE *code)

  5472 {     /* display member specified by pseudo-code */

  5473       MEMBER memb;

  5474       ARG_LIST *e;

  5475       xassert(code->op == O_MEMNUM || code->op == O_MEMSYM

  5476          || code->op == O_MEMSET || code->op == O_MEMVAR

  5477          || code->op == O_MEMCON);

  5478       memb.tuple = create_tuple(mpl);

  5479       for (e = code->arg.par.list; e != NULL; e = e->next)

  5480          memb.tuple = expand_tuple(mpl, memb.tuple, eval_symbolic(mpl,

  5481             e->x));

  5482       switch (code->op)

  5483       {  case O_MEMNUM:

  5484             memb.value.num = eval_member_num(mpl, code->arg.par.par,

  5485                memb.tuple);

  5486             display_par(mpl, code->arg.par.par, &memb);

  5487             break;

  5488          case O_MEMSYM:

  5489             memb.value.sym = eval_member_sym(mpl, code->arg.par.par,

  5490                memb.tuple);

  5491             display_par(mpl, code->arg.par.par, &memb);

  5492             delete_symbol(mpl, memb.value.sym);

  5493             break;

  5494          case O_MEMSET:

  5495             memb.value.set = eval_member_set(mpl, code->arg.set.set,

  5496                memb.tuple);

  5497             display_set(mpl, code->arg.set.set, &memb);

  5498             break;

  5499          case O_MEMVAR:

  5500             memb.value.var = eval_member_var(mpl, code->arg.var.var,

  5501                memb.tuple);

  5502             display_var

  5503                (mpl, code->arg.var.var, &memb, code->arg.var.suff);

  5504             break;

  5505          case O_MEMCON:

  5506             memb.value.con = eval_member_con(mpl, code->arg.con.con,

  5507                memb.tuple);

  5508             display_con

  5509                (mpl, code->arg.con.con, &memb, code->arg.con.suff);

  5510             break;

  5511          default:

  5512             xassert(code != code);

  5513       }

  5514       delete_tuple(mpl, memb.tuple);

  5515       return;

  5516 }

  5518 static void display_code(MPL *mpl, CODE *code)

  5519 {     /* display value of expression */

  5520       switch (code->type)

  5521       {  case A_NUMERIC:

  5522             /* numeric value */

  5523             {  double num;

  5524                num = eval_numeric(mpl, code);

  5525                write_text(mpl, "%.*g\n", DBL_DIG, num);

  5526             }

  5527             break;

  5528          case A_SYMBOLIC:

  5529             /* symbolic value */

  5530             {  SYMBOL *sym;

  5531                sym = eval_symbolic(mpl, code);

  5532                write_text(mpl, "%s\n", format_symbol(mpl, sym));

  5533                delete_symbol(mpl, sym);

  5534             }

  5535             break;

  5536          case A_LOGICAL:

  5537             /* logical value */

  5538             {  int bit;

  5539                bit = eval_logical(mpl, code);

  5540                write_text(mpl, "%s\n", bit ? "true" : "false");

  5541             }

  5542             break;

  5543          case A_TUPLE:

  5544             /* n-tuple */

  5545             {  TUPLE *tuple;

  5546                tuple = eval_tuple(mpl, code);

  5547                write_text(mpl, "%s\n", format_tuple(mpl, '(', tuple));

  5548                delete_tuple(mpl, tuple);

  5549             }

  5550             break;

  5551          case A_ELEMSET:

  5552             /* elemental set */

  5553             {  ELEMSET *set;

  5554                MEMBER *memb;

  5555                set = eval_elemset(mpl, code);

  5556                if (set->head == 0)

  5557                   write_text(mpl, "set is empty\n");

  5558                for (memb = set->head; memb != NULL; memb = memb->next)

  5559                   write_text(mpl, "   %s\n", format_tuple(mpl, '(',

  5560                      memb->tuple));

  5561                delete_elemset(mpl, set);

  5562             }

  5563             break;

  5564          case A_FORMULA:

  5565             /* linear form */

  5566             {  FORMULA *form, *term;

  5567                form = eval_formula(mpl, code);

  5568                if (form == NULL)

  5569                   write_text(mpl, "linear form is empty\n");

  5570                for (term = form; term != NULL; term = term->next)

  5571                {  if (term->var == NULL)

  5572                      write_text(mpl, "   %.*g\n", term->coef);

  5573                   else

  5574                      write_text(mpl, "   %.*g %s%s\n", DBL_DIG,

  5575                         term->coef, term->var->var->name,

  5576                         format_tuple(mpl, '[', term->var->memb->tuple));

  5577                }

  5578                delete_formula(mpl, form);

  5579             }

  5580             break;

  5581          default:

  5582             xassert(code != code);

  5583       }

  5584       return;

  5585 }

  5587 static int display_func(MPL *mpl, void *info)

  5588 {     /* this is auxiliary routine to work within domain scope */

  5589       DISPLAY *dpy = (DISPLAY *)info;

  5590       DISPLAY1 *entry;

  5591       for (entry = dpy->list; entry != NULL; entry = entry->next)

  5592       {  if (entry->type == A_INDEX)

  5593          {  /* dummy index */

  5594             DOMAIN_SLOT *slot = entry->u.slot;

  5595             write_text(mpl, "%s = %s\n", slot->name,

  5596             format_symbol(mpl, slot->value));

  5597          }

  5598          else if (entry->type == A_SET)

  5599          {  /* model set */

  5600             SET *set = entry->u.set;

  5601             MEMBER *memb;

  5602             if (set->assign != NULL)

  5603             {  /* the set has assignment expression; evaluate all its

  5604                   members over entire domain */

  5605                eval_whole_set(mpl, set);

  5606             }

  5607             else

  5608             {  /* the set has no assignment expression; refer to its

  5609                   any existing member ignoring resultant value to check

  5610                   the data provided the data section */

  5611 #if 1 /* 12/XII-2008 */

  5612                if (set->gadget != NULL && set->data == 0)

  5613                {  /* initialize the set with data from a plain set */

  5614                   saturate_set(mpl, set);

  5615                }

  5616 #endif

  5617                if (set->array->head != NULL)

  5618                   eval_member_set(mpl, set, set->array->head->tuple);

  5619             }

  5620             /* display all members of the set array */

  5621             if (set->array->head == NULL)

  5622                write_text(mpl, "%s has empty content\n", set->name);

  5623             for (memb = set->array->head; memb != NULL; memb =

  5624                memb->next) display_set(mpl, set, memb);

  5625          }

  5626          else if (entry->type == A_PARAMETER)

  5627          {  /* model parameter */

  5628             PARAMETER *par = entry->u.par;

  5629             MEMBER *memb;

  5630             if (par->assign != NULL)

  5631             {  /* the parameter has an assignment expression; evaluate

  5632                   all its member over entire domain */

  5633                eval_whole_par(mpl, par);

  5634             }

  5635             else

  5636             {  /* the parameter has no assignment expression; refer to

  5637                   its any existing member ignoring resultant value to

  5638                   check the data provided in the data section */

  5639                if (par->array->head != NULL)

  5640                {  if (par->type != A_SYMBOLIC)

  5641                      eval_member_num(mpl, par, par->array->head->tuple);

  5642                   else

  5643                      delete_symbol(mpl, eval_member_sym(mpl, par,

  5644                         par->array->head->tuple));

  5645                }

  5646             }

  5647             /* display all members of the parameter array */

  5648             if (par->array->head == NULL)

  5649                write_text(mpl, "%s has empty content\n", par->name);

  5650             for (memb = par->array->head; memb != NULL; memb =

  5651                memb->next) display_par(mpl, par, memb);

  5652          }

  5653          else if (entry->type == A_VARIABLE)

  5654          {  /* model variable */

  5655             VARIABLE *var = entry->u.var;

  5656             MEMBER *memb;

  5657             xassert(mpl->flag_p);

  5658             /* display all members of the variable array */

  5659             if (var->array->head == NULL)

  5660                write_text(mpl, "%s has empty content\n", var->name);

  5661             for (memb = var->array->head; memb != NULL; memb =

  5662                memb->next) display_var(mpl, var, memb, DOT_NONE);

  5663          }

  5664          else if (entry->type == A_CONSTRAINT)

  5665          {  /* model constraint */

  5666             CONSTRAINT *con = entry->u.con;

  5667             MEMBER *memb;

  5668             xassert(mpl->flag_p);

  5669             /* display all members of the constraint array */

  5670             if (con->array->head == NULL)

  5671                write_text(mpl, "%s has empty content\n", con->name);

  5672             for (memb = con->array->head; memb != NULL; memb =

  5673                memb->next) display_con(mpl, con, memb, DOT_NONE);

  5674          }

  5675          else if (entry->type == A_EXPRESSION)

  5676          {  /* expression */

  5677             CODE *code = entry->u.code;

  5678             if (code->op == O_MEMNUM || code->op == O_MEMSYM ||

  5679                 code->op == O_MEMSET || code->op == O_MEMVAR ||

  5680                 code->op == O_MEMCON)

  5681                display_memb(mpl, code);

  5682             else

  5683                display_code(mpl, code);

  5684          }

  5685          else

  5686             xassert(entry != entry);

  5687       }

  5688       return 0;

  5689 }

  5691 void execute_display(MPL *mpl, DISPLAY *dpy)

  5692 {     loop_within_domain(mpl, dpy->domain, dpy, display_func);

  5693       return;

  5694 }

  5696 /*----------------------------------------------------------------------

  5697 -- clean_display - clean display statement.

  5698 --

  5699 -- This routine cleans specified display statement that assumes deleting

  5700 -- all stuff dynamically allocated on generating/postsolving phase. */

  5702 void clean_display(MPL *mpl, DISPLAY *dpy)

  5703 {     DISPLAY1 *d;

  5704 #if 0 /* 15/V-2010 */

  5705       ARG_LIST *e;

  5706 #endif

  5707       /* clean subscript domain */

  5708       clean_domain(mpl, dpy->domain);

  5709       /* clean display list */

  5710       for (d = dpy->list; d != NULL; d = d->next)

  5711       {  /* clean pseudo-code for computing expression */

  5712          if (d->type == A_EXPRESSION)

  5713             clean_code(mpl, d->u.code);

  5714 #if 0 /* 15/V-2010 */

  5715          /* clean pseudo-code for computing subscripts */

  5716          for (e = d->list; e != NULL; e = e->next)

  5717             clean_code(mpl, e->x);

  5718 #endif

  5719       }

  5720       return;

  5721 }

  5723 /*----------------------------------------------------------------------

  5724 -- execute_printf - execute printf statement.

  5725 --

  5726 -- This routine executes specified printf statement. */

  5728 #if 1 /* 14/VII-2006 */

  5729 static void print_char(MPL *mpl, int c)

  5730 {     if (mpl->prt_fp == NULL)

  5731          write_char(mpl, c);

  5732       else

  5733          xfputc(c, mpl->prt_fp);

  5734       return;

  5735 }

  5737 static void print_text(MPL *mpl, char *fmt, ...)

  5738 {     va_list arg;

  5739       char buf[OUTBUF_SIZE], *c;

  5740       va_start(arg, fmt);

  5741       vsprintf(buf, fmt, arg);

  5742       xassert(strlen(buf) < sizeof(buf));

  5743       va_end(arg);

  5744       for (c = buf; *c != '\0'; c++) print_char(mpl, *c);

  5745       return;

  5746 }

  5747 #endif

  5749 static int printf_func(MPL *mpl, void *info)

  5750 {     /* this is auxiliary routine to work within domain scope */

  5751       PRINTF *prt = (PRINTF *)info;

  5752       PRINTF1 *entry;

  5753       SYMBOL *sym;

  5754       char fmt[MAX_LENGTH+1], *c, *from, save;

  5755       /* evaluate format control string */

  5756       sym = eval_symbolic(mpl, prt->fmt);

  5757       if (sym->str == NULL)

  5758          sprintf(fmt, "%.*g", DBL_DIG, sym->num);

  5759       else

  5760          fetch_string(mpl, sym->str, fmt);

  5761       delete_symbol(mpl, sym);

  5762       /* scan format control string and perform formatting output */

  5763       entry = prt->list;

  5764       for (c = fmt; *c != '\0'; c++)

  5765       {  if (*c == '%')

  5766          {  /* scan format specifier */

  5767             from = c++;

  5768             if (*c == '%')

  5769             {  print_char(mpl, '%');

  5770                continue;

  5771             }

  5772             if (entry == NULL) break;

  5773             /* scan optional flags */

  5774             while (*c == '-' || *c == '+' || *c == ' ' || *c == '#' ||

  5775                    *c == '0') c++;

  5776             /* scan optional minimum field width */

  5777             while (isdigit((unsigned char)*c)) c++;

  5778             /* scan optional precision */

  5779             if (*c == '.')

  5780             {  c++;

  5781                while (isdigit((unsigned char)*c)) c++;

  5782             }

  5783             /* scan conversion specifier and perform formatting */

  5784             save = *(c+1), *(c+1) = '\0';

  5785             if (*c == 'd' || *c == 'i' || *c == 'e' || *c == 'E' ||

  5786                 *c == 'f' || *c == 'F' || *c == 'g' || *c == 'G')

  5787             {  /* the specifier requires numeric value */

  5788                double value;

  5789                xassert(entry != NULL);

  5790                switch (entry->code->type)

  5791                {  case A_NUMERIC:

  5792                      value = eval_numeric(mpl, entry->code);

  5793                      break;

  5794                   case A_SYMBOLIC:

  5795                      sym = eval_symbolic(mpl, entry->code);

  5796                      if (sym->str != NULL)

  5797                         error(mpl, "cannot convert %s to floating-point"

  5798                            " number", format_symbol(mpl, sym));

  5799                      value = sym->num;

  5800                      delete_symbol(mpl, sym);

  5801                      break;

  5802                   case A_LOGICAL:

  5803                      if (eval_logical(mpl, entry->code))

  5804                         value = 1.0;

  5805                      else

  5806                         value = 0.0;

  5807                      break;

  5808                   default:

  5809                      xassert(entry != entry);

  5810                }

  5811                if (*c == 'd' || *c == 'i')

  5812                {  double int_max = (double)INT_MAX;

  5813                   if (!(-int_max <= value && value <= +int_max))

  5814                      error(mpl, "cannot convert %.*g to integer",

  5815                         DBL_DIG, value);

  5816                   print_text(mpl, from, (int)floor(value + 0.5));

  5817                }

  5818                else

  5819                   print_text(mpl, from, value);

  5820             }

  5821             else if (*c == 's')

  5822             {  /* the specifier requires symbolic value */

  5823                char value[MAX_LENGTH+1];

  5824                switch (entry->code->type)

  5825                {  case A_NUMERIC:

  5826                      sprintf(value, "%.*g", DBL_DIG, eval_numeric(mpl,

  5827                         entry->code));

  5828                      break;

  5829                   case A_LOGICAL:

  5830                      if (eval_logical(mpl, entry->code))

  5831                         strcpy(value, "T");

  5832                      else

  5833                         strcpy(value, "F");

  5834                      break;

  5835                   case A_SYMBOLIC:

  5836                      sym = eval_symbolic(mpl, entry->code);

  5837                      if (sym->str == NULL)

  5838                         sprintf(value, "%.*g", DBL_DIG, sym->num);

  5839                      else

  5840                         fetch_string(mpl, sym->str, value);

  5841                      delete_symbol(mpl, sym);

  5842                      break;

  5843                   default:

  5844                      xassert(entry != entry);

  5845                }

  5846                print_text(mpl, from, value);

  5847             }

  5848             else

  5849                error(mpl, "format specifier missing or invalid");

  5850             *(c+1) = save;

  5851             entry = entry->next;

  5852          }

  5853          else if (*c == '\\')

  5854          {  /* write some control character */

  5855             c++;

  5856             if (*c == 't')

  5857                print_char(mpl, '\t');

  5858             else if (*c == 'n')

  5859                print_char(mpl, '\n');

  5860 #if 1 /* 28/X-2010 */

  5861             else if (*c == '\0')

  5862             {  /* format string ends with backslash */

  5863                error(mpl, "invalid use of escape character \\ in format"

  5864                   " control string");

  5865             }

  5866 #endif

  5867             else

  5868                print_char(mpl, *c);

  5869          }

  5870          else

  5871          {  /* write character without formatting */

  5872             print_char(mpl, *c);

  5873          }

  5874       }

  5875       return 0;

  5876 }

  5878 #if 0 /* 14/VII-2006 */

  5879 void execute_printf(MPL *mpl, PRINTF *prt)

  5880 {     loop_within_domain(mpl, prt->domain, prt, printf_func);

  5881       return;

  5882 }

  5883 #else

  5884 void execute_printf(MPL *mpl, PRINTF *prt)

  5885 {     if (prt->fname == NULL)

  5886       {  /* switch to the standard output */

  5887          if (mpl->prt_fp != NULL)

  5888          {  xfclose(mpl->prt_fp), mpl->prt_fp = NULL;

  5889             xfree(mpl->prt_file), mpl->prt_file = NULL;

  5890          }

  5891       }

  5892       else

  5893       {  /* evaluate file name string */

  5894          SYMBOL *sym;

  5895          char fname[MAX_LENGTH+1];

  5896          sym = eval_symbolic(mpl, prt->fname);

  5897          if (sym->str == NULL)

  5898             sprintf(fname, "%.*g", DBL_DIG, sym->num);

  5899          else

  5900             fetch_string(mpl, sym->str, fname);

  5901          delete_symbol(mpl, sym);

  5902          /* close the current print file, if necessary */

  5903          if (mpl->prt_fp != NULL &&

  5904             (!prt->app || strcmp(mpl->prt_file, fname) != 0))

  5905          {  xfclose(mpl->prt_fp), mpl->prt_fp = NULL;

  5906             xfree(mpl->prt_file), mpl->prt_file = NULL;

  5907          }

  5908          /* open the specified print file, if necessary */

  5909          if (mpl->prt_fp == NULL)

  5910          {  mpl->prt_fp = xfopen(fname, prt->app ? "a" : "w");

  5911             if (mpl->prt_fp == NULL)

  5912                error(mpl, "unable to open `%s' for writing - %s",

  5913                   fname, xerrmsg());

  5914             mpl->prt_file = xmalloc(strlen(fname)+1);

  5915             strcpy(mpl->prt_file, fname);

  5916          }

  5917       }

  5918       loop_within_domain(mpl, prt->domain, prt, printf_func);

  5919       if (mpl->prt_fp != NULL)

  5920       {  xfflush(mpl->prt_fp);

  5921          if (xferror(mpl->prt_fp))

  5922             error(mpl, "writing error to `%s' - %s", mpl->prt_file,

  5923                xerrmsg());

  5924       }

  5925       return;

  5926 }

  5927 #endif

  5929 /*----------------------------------------------------------------------

  5930 -- clean_printf - clean printf statement.

  5931 --

  5932 -- This routine cleans specified printf statement that assumes deleting

  5933 -- all stuff dynamically allocated on generating/postsolving phase. */

  5935 void clean_printf(MPL *mpl, PRINTF *prt)

  5936 {     PRINTF1 *p;

  5937       /* clean subscript domain */

  5938       clean_domain(mpl, prt->domain);

  5939       /* clean pseudo-code for computing format string */

  5940       clean_code(mpl, prt->fmt);

  5941       /* clean printf list */

  5942       for (p = prt->list; p != NULL; p = p->next)

  5943       {  /* clean pseudo-code for computing value to be printed */

  5944          clean_code(mpl, p->code);

  5945       }

  5946 #if 1 /* 14/VII-2006 */

  5947       /* clean pseudo-code for computing file name string */

  5948       clean_code(mpl, prt->fname);

  5949 #endif

  5950       return;

  5951 }

  5953 /*----------------------------------------------------------------------

  5954 -- execute_for - execute for statement.

  5955 --

  5956 -- This routine executes specified for statement. */

  5958 static int for_func(MPL *mpl, void *info)

  5959 {     /* this is auxiliary routine to work within domain scope */

  5960       FOR *fur = (FOR *)info;

  5961       STATEMENT *stmt, *save;

  5962       save = mpl->stmt;

  5963       for (stmt = fur->list; stmt != NULL; stmt = stmt->next)

  5964          execute_statement(mpl, stmt);

  5965       mpl->stmt = save;

  5966       return 0;

  5967 }

  5969 void execute_for(MPL *mpl, FOR *fur)

  5970 {     loop_within_domain(mpl, fur->domain, fur, for_func);

  5971       return;

  5972 }

  5974 /*----------------------------------------------------------------------

  5975 -- clean_for - clean for statement.

  5976 --

  5977 -- This routine cleans specified for statement that assumes deleting all

  5978 -- stuff dynamically allocated on generating/postsolving phase. */

  5980 void clean_for(MPL *mpl, FOR *fur)

  5981 {     STATEMENT *stmt;

  5982       /* clean subscript domain */

  5983       clean_domain(mpl, fur->domain);

  5984       /* clean all sub-statements */

  5985       for (stmt = fur->list; stmt != NULL; stmt = stmt->next)

  5986          clean_statement(mpl, stmt);

  5987       return;

  5988 }

  5990 /*----------------------------------------------------------------------

  5991 -- execute_statement - execute specified model statement.

  5992 --

  5993 -- This routine executes specified model statement. */

  5995 void execute_statement(MPL *mpl, STATEMENT *stmt)

  5996 {     mpl->stmt = stmt;

  5997       switch (stmt->type)

  5998       {  case A_SET:

  5999          case A_PARAMETER:

  6000          case A_VARIABLE:

  6001             break;

  6002          case A_CONSTRAINT:

  6003             xprintf("Generating %s...\n", stmt->u.con->name);

  6004             eval_whole_con(mpl, stmt->u.con);

  6005             break;

  6006          case A_TABLE:

  6007             switch (stmt->u.tab->type)

  6008             {  case A_INPUT:

  6009                   xprintf("Reading %s...\n", stmt->u.tab->name);

  6010                   break;

  6011                case A_OUTPUT:

  6012                   xprintf("Writing %s...\n", stmt->u.tab->name);

  6013                   break;

  6014                default:

  6015                   xassert(stmt != stmt);

  6016             }

  6017             execute_table(mpl, stmt->u.tab);

  6018             break;

  6019          case A_SOLVE:

  6020             break;

  6021          case A_CHECK:

  6022             xprintf("Checking (line %d)...\n", stmt->line);

  6023             execute_check(mpl, stmt->u.chk);

  6024             break;

  6025          case A_DISPLAY:

  6026             write_text(mpl, "Display statement at line %d\n",

  6027                stmt->line);

  6028             execute_display(mpl, stmt->u.dpy);

  6029             break;

  6030          case A_PRINTF:

  6031             execute_printf(mpl, stmt->u.prt);

  6032             break;

  6033          case A_FOR:

  6034             execute_for(mpl, stmt->u.fur);

  6035             break;

  6036          default:

  6037             xassert(stmt != stmt);

  6038       }

  6039       return;

  6040 }

  6042 /*----------------------------------------------------------------------

  6043 -- clean_statement - clean specified model statement.

  6044 --

  6045 -- This routine cleans specified model statement that assumes deleting

  6046 -- all stuff dynamically allocated on generating/postsolving phase. */

  6048 void clean_statement(MPL *mpl, STATEMENT *stmt)

  6049 {     switch(stmt->type)

  6050       {  case A_SET:

  6051             clean_set(mpl, stmt->u.set); break;

  6052          case A_PARAMETER:

  6053             clean_parameter(mpl, stmt->u.par); break;

  6054          case A_VARIABLE:

  6055             clean_variable(mpl, stmt->u.var); break;

  6056          case A_CONSTRAINT:

  6057             clean_constraint(mpl, stmt->u.con); break;

  6058 #if 1 /* 11/II-2008 */

  6059          case A_TABLE:

  6060             clean_table(mpl, stmt->u.tab); break;

  6061 #endif

  6062          case A_SOLVE:

  6063             break;

  6064          case A_CHECK:

  6065             clean_check(mpl, stmt->u.chk); break;

  6066          case A_DISPLAY:

  6067             clean_display(mpl, stmt->u.dpy); break;

  6068          case A_PRINTF:

  6069             clean_printf(mpl, stmt->u.prt); break;

  6070          case A_FOR:

  6071             clean_for(mpl, stmt->u.fur); break;

  6072          default:

  6073             xassert(stmt != stmt);

  6074       }

  6075       return;

  6076 }

  6078 /* eof */

author	Alpar Juttner <alpar@cs.elte.hu>
	Sun, 05 Dec 2010 17:35:23 +0100
changeset 2	4c8956a7bdf4
permissions	-rw-r--r--