interpreters/nitfol/solve.c

   1 /*  solve.c: solve cycles for the automapper
   2     Copyright (C) 1999  Evin Robertson
   3
   4     This program is free software; you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation; either version 2 of the License, or
   7     (at your option) any later version.
   8
   9     This program is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with this program; if not, write to the Free Software
  16     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  17
  18     The author can be reached at nitfol@deja.com
  19 */
  20
  21
  22 #include "nitfol.h"
  23
  24
  25 /* Rational number stuff - doesn't handle overflow conditions */
  26
  27 typedef struct rational rational;
  28 struct rational {
  29   int num;
  30   int den;
  31 };
  32
  33 static int gcd(int a, int b)
  34 {
  35   if(a == 0)
  36     return b;
  37   return gcd(b % a, a);
  38 }
  39
  40 static int lcm(int a, int b)
  41 {
  42   return a * b / gcd(a, b);
  43 }
  44
  45 static rational rational_reduce(rational r)
  46 {
  47   int g = gcd(r.num, r.den);
  48   if(g == 0)
  49     return r;
  50   r.num /= g;
  51   r.den /= g;
  52   if(r.den < 0) {
  53     r.num = -r.num;
  54     r.den = -r.den;
  55   }
  56   return r;
  57 }
  58
  59 static BOOL rational_iszero(rational r)
  60 {
  61   return r.num == 0;
  62 }
  63
  64 static BOOL rational_isone(rational r)
  65 {
  66   return r.num == r.den;
  67 }
  68
  69 static rational rational_int(int i)
  70 {
  71   rational r;
  72   r.num = i;
  73   r.den = 1;
  74   return r;
  75 }
  76
  77 static rational rational_add(rational a, rational b)
  78 {
  79   rational r;
  80   r.num = a.num * b.den + a.den * b.num;
  81   r.den = a.den * b.den;
  82   return rational_reduce(r);
  83 }
  84
  85 static rational rational_sub(rational a, rational b)
  86 {
  87   rational r;
  88   r.num = a.num * b.den - a.den * b.num;
  89   r.den = a.den * b.den;
  90   return rational_reduce(r);
  91 }
  92
  93 static rational rational_mult(rational a, rational b)
  94 {
  95   a.num *= b.num;
  96   a.den *= b.den;
  97   return rational_reduce(a);
  98 }
  99
 100 static rational rational_div(rational a, rational b)
 101 {
 102   a.num *= b.den;
 103   a.den *= b.num;
 104   return rational_reduce(a);
 105 }
 106
 107
 108 #ifdef HEADER
 109 typedef struct cycleequation cycleequation;
 110 struct cycleequation {
 111   cycleequation *next;
 112
 113   int *var;
 114   const int *min;
 115   int xcoefficient;
 116   int ycoefficient;
 117 };
 118 #endif
 119
 120 typedef struct equation equation;
 121 struct equation {
 122   equation *next;
 123
 124   int *var;
 125   const int *min;
 126   rational coefficient;
 127   int count; /* Number of times variable is used */
 128 };
 129
 130 typedef struct equalist equalist;
 131 struct equalist {
 132   equalist *next;
 133   equation *eq;
 134 };
 135
 136 static equalist *equats = NULL;
 137
 138 void automap_add_cycle(cycleequation *cycle)
 139 {
 140   equation *newx = NULL;
 141   equation *newy = NULL;
 142   equalist newlist;
 143   cycleequation *p;
 144
 145   for(p = cycle; p; p=p->next) {
 146     equation newnode;
 147     newnode.var = p->var;
 148     newnode.min = p->min;
 149     newnode.coefficient.den = 1;
 150     newnode.coefficient.num = p->xcoefficient;
 151     LEadd(newx, newnode);
 152
 153     newnode.coefficient.num = p->ycoefficient;
 154     LEadd(newy, newnode);
 155   }
 156
 157   newlist.eq = newx;
 158   LEadd(equats, newlist);
 159   newlist.eq = newy;
 160   LEadd(equats, newlist);
 161
 162   LEdestroy(cycle);
 163 }
 164
 165 void automap_delete_cycles(void)
 166 {
 167   equalist *p;
 168   for(p = equats; p; p=p->next)
 169     LEdestroy(p->eq);
 170   LEdestroy(equats);
 171 }
 172
 173
 174 /* Do Gauss-Jordan elimination. */
 175 /* This could be done with lists instead of arrays with clever hash stuff, but
 176    this is simpler, and the elimination stage takes O(n^3) anyway */
 177 void automap_cycle_elimination(void)
 178 {
 179   equalist *p;
 180   equation *q, *v;
 181   equation *variablelist = NULL;
 182   int width = 0, height = 0;
 183   rational *array;
 184   rational *r, *s;
 185   int row, column, i, n;
 186
 187
 188   /* First, figure out how many variables we're dealing with */
 189   for(p = equats; p; p=p->next) {
 190     for(q = p->eq; q; q=q->next) {
 191       for(v = variablelist; v; v=v->next)
 192         if(v->var == q->var)
 193           break;
 194       if(!v) {
 195         LEadd(variablelist, *q);
 196         width++;
 197       }
 198     }
 199     height++;
 200   }
 201
 202   if(height == 0 || width == 0)
 203     return;
 204
 205   /* An array implementation is easier to think about than a linked-list one */
 206   array = (rational *) n_malloc(width * height * sizeof(rational));
 207
 208   /* Fill the array */
 209   r = array;
 210   for(p = equats; p; p=p->next) {
 211     for(v = variablelist; v; v=v->next) {
 212       for(q = p->eq; q; q=q->next)
 213         if(q->var == v->var)
 214           break;
 215       *r++ = q ? q->coefficient : rational_int(0);
 216     }
 217   }
 218
 219   /* Do the actual elimination */
 220   row = 0; column = 0;
 221   r = array;
 222   while(row < height && column < width) {
 223     rational c;
 224     /* If zero, swap with a non-zero row */
 225     if(rational_iszero(r[column])) {
 226       BOOL found = FALSE;
 227       for(i = row+1; i < height; i++) {
 228         if(!rational_iszero(array[i * width + column])) {
 229           n_memswap(&array[row * width], &array[i * width],
 230                     width * sizeof(*array));
 231           found = TRUE;
 232           break;
 233         }
 234       }
 235       if(!found) { /* Zeroed column; try next */
 236         column++;
 237         continue;
 238       }
 239     }
 240
 241     /* Divide row by leading coefficient */
 242     c = r[column];
 243     for(n = 0; n < width; n++)
 244       r[n] = rational_div(r[n], c);
 245
 246     /* Eliminate other entries in current column */
 247     s = array;
 248     for(i = 0; i < height; i++) {
 249       if(i != row && !rational_iszero(s[column])) {
 250         c = s[column];
 251         for(n = 0; n < width; n++)
 252           s[n] = rational_sub(s[n], rational_mult(r[n], c));
 253       }
 254       s += width;
 255     }
 256
 257     r += width;
 258     row++; column++;
 259   }
 260
 261   /* Delete the old lists */
 262   automap_delete_cycles();
 263
 264   /* Count how many times each variable is used */
 265   for(v = variablelist; v; v=v->next)
 266     v->count = 0;
 267   r = array;
 268   for(i = 0; i < height; i++) {
 269     for(v = variablelist; v; v=v->next) {
 270       if(!rational_iszero(*r))
 271         v->count++;
 272       r++;
 273     }
 274   }
 275
 276   /* Make new lists from the array */
 277   r = array;
 278   for(i = 0; i < height; i++) {
 279     equation *neweq = NULL;
 280     for(v = variablelist; v; v=v->next) {
 281       if(!rational_iszero(*r)) {
 282         equation newnode;
 283         newnode = *v;
 284         newnode.coefficient = *r;
 285         LEadd(neweq, newnode);
 286       }
 287       r++;
 288     }
 289     if(neweq) {
 290       equalist newlist;
 291       newlist.eq = neweq;
 292       LEadd(equats, newlist);
 293     }
 294   }
 295
 296   n_free(array);
 297 }
 298
 299
 300 /* Find an edge to lengthen which would cause the least amount of lengthening
 301    to edges in other cycles */
 302 static equation *select_edge(equation *cycle, int *need_recalc)
 303 {
 304   equation *help = NULL;     /* Increasing its length will help other cycles */
 305   equation *solitary = NULL; /* Only in one cycle */
 306   equation *nonharm = NULL;  /* Increasing its length won't destroy things */
 307   BOOL is_harmful_past = FALSE;
 308
 309   equation *p;
 310
 311   for(p = cycle; p; p=p->next) {
 312     if(p->next && p->coefficient.num < 0) {
 313       equalist *t;
 314       BOOL pastthis = FALSE;
 315       BOOL is_found = FALSE;
 316       BOOL is_harmful = FALSE;
 317       BOOL is_past = FALSE;
 318       BOOL is_help = FALSE;
 319       BOOL is_compensator = FALSE;
 320       for(t = equats; t; t=t->next) {
 321         if(t->eq == cycle)
 322           pastthis = TRUE;
 323         else {
 324           rational sum = rational_int(0);
 325           equation *r, *foundme = NULL;
 326           BOOL first_find = TRUE;
 327           for(r = t->eq; r; r=r->next) {
 328             if(r->next) {
 329               int value = *(r->var) ? *(r->var) : *(r->min);
 330               sum = rational_add(sum, rational_mult(r->coefficient,
 331                                                     rational_int(value)));
 332               if(r->count == 1 && r->coefficient.num < 0)
 333                 is_compensator = TRUE;
 334               if(r->var == p->var) {
 335                 if(foundme)
 336                   first_find = FALSE;
 337                 foundme = r;
 338                 is_past = pastthis && (is_past || !is_found);
 339                 is_found = TRUE;
 340                 if(r->coefficient.num > 0)
 341                   is_harmful = TRUE;
 342               }
 343             } else if(pastthis && foundme && -sum.num < *(r->min) && first_find
 344                       && foundme->coefficient.num < 0)
 345               is_help = TRUE;
 346           }
 347         }
 348       }
 349       if(is_help && !is_harmful && !help)
 350         help = p;
 351       if(!is_found) {
 352         solitary = p;
 353       } else if(!is_harmful || is_compensator) {
 354         if(!nonharm || is_past) {
 355           is_harmful_past = !is_past;
 356           nonharm = p;
 357         }
 358       }
 359     }
 360   }
 361
 362   if(help)     return help;
 363   if(solitary) return solitary;
 364   if(nonharm)  { if(is_harmful_past) *need_recalc = 2; return nonharm; }
 365
 366   return NULL;
 367 }
 368
 369
 370
 371 /* Fill variables with valid numbers. Assumes Gauss-Jordan elimination has
 372    already happened. */
 373 void automap_cycles_fill_values(void)
 374 {
 375   equalist *p;
 376   equation *q;
 377   int calccount;
 378   int recalculate = 0;
 379
 380   for(p = equats; p; p=p->next)
 381     for(q = p->eq; q; q=q->next)
 382       *(q->var) = 0;
 383
 384   for(calccount = 0; calccount <= recalculate; calccount++) {
 385     recalculate = 0;
 386
 387   /* Last variable in each list is the dependant; all others are independant */
 388   /* Fill independant variables with their minimums, then calculate the
 389      dependant one; if it's less than its minimum, play with independants */
 390     for(p = equats; p; p=p->next) {
 391       rational sum = rational_int(0);
 392       for(q = p->eq; q; q=q->next) {
 393         if(q->next) { /* Independant */
 394           int value = *(q->var) ? *(q->var) : *(q->min);
 395           sum = rational_add(sum, rational_mult(q->coefficient,
 396                                                 rational_int(value)));
 397           *(q->var) = value;
 398         } else { /* Dependant */
 399           int value = -sum.num;
 400           if(!rational_isone(q->coefficient))
 401             n_show_error(E_SYSTEM, "last coefficient not 1", q->coefficient.num);
 402           if(sum.den != 1)
 403             n_show_error(E_SYSTEM, "unimplemented case denominator != 1", sum.den);
 404           else if(value < *(q->min)) {
 405             /* Edge is not long enough - try increasing lengths of another edge
 406                in cycle to lengthen it */
 407             equation *m = select_edge(p->eq, &recalculate);
 408
 409             if(m) {
 410               rational oldval = rational_mult(m->coefficient,
 411                                               rational_int(*(m->var)));
 412               rational newval;
 413               int diff = value - *(q->min);
 414               sum = rational_sub(sum, oldval);
 415               if(oldval.den != 1)
 416                 n_show_error(E_SYSTEM, "unimplemented case denominator != 1", oldval.den);
 417               diff += oldval.num;
 418               newval = rational_div(rational_int(diff), m->coefficient);
 419               *(m->var) = newval.num;
 420               sum = rational_add(sum, rational_mult(m->coefficient, newval));
 421               value = -sum.num;
 422             }
 423             if(value > *(q->min))
 424               n_show_error(E_SYSTEM, "met more than needed", sum.num);
 425           }
 426           if(value < *(q->min))
 427             n_show_error(E_SYSTEM, "failed to meet needs", sum.num);
 428           *(q->var) = value;
 429           sum = rational_add(sum, rational_mult(q->coefficient,
 430                                                 rational_int(*(q->var))));
 431           if(!rational_iszero(sum))
 432             n_show_error(E_SYSTEM, "doesn't add up", sum.num);
 433
 434         }
 435       }
 436     }
 437 #if 0
 438     {
 439       rational checksum = rational_int(0);
 440       equation *cq;
 441       for(cq = p->eq; cq; cq=cq->next) {
 442         checksum = rational_add(checksum,rational_mult(cq->coefficient,
 443                                                     rational_int(*(cq->var))));
 444       }
 445       if(checksum.num != sum.num || checksum.den != sum.den) {
 446         n_show_error(E_SYSTEM, "correction for correction incorrect", checksum.num);
 447         sum = checksum;
 448       }
 449     }
 450 #endif
 451   }
 452
 453
 454 #if 0
 455   for(p = equats; p; p=p->next) {
 456     rational sum = rational_int(0);
 457     for(q = p->eq; q; q=q->next) {
 458       if(*(q->var) < *(q->min))
 459         n_show_error(E_SYSTEM, "variable less than minimum", *(q->var));
 460       sum = rational_add(sum, rational_mult(q->coefficient,
 461                                             rational_int(*(q->var))));
 462     }
 463     if(!rational_iszero(sum))
 464       n_show_error(E_SYSTEM, "equation not equal", sum.num / sum.den);
 465   }
 466 #endif
 467
 468
 469 }
 470
 471