← Back to branch summary

~ubuntu-branches/ubuntu/wily/gargoyle-free/wily-proposed

« back to all changes in this revision

Viewing changes to terps/nitfol/solve.c

  • Committer: Bazaar Package Importer
  • Author(s): Sylvain Beucler
  • Date: 2009-09-11 20:09:43 UTC
  • Revision ID: james.westby@ubuntu.com-20090911200943-idgzoyupq6650zpn
Tags: upstream-2009-08-25
Import upstream version 2009-08-25

Show diffs side-by-side

added added

removed removed

Lines of Context:
terps/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., 59 Temple Place, Suite 330, Boston, MA 02111, 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