zsolve/linearsystem.c

   1 /*
   2 4ti2 -- A software package for algebraic, geometric and combinatorial
   3 problems on linear spaces.
   4
   5 Copyright (C) 2006 4ti2 team.
   6 Main author(s): Matthias Walter.
   7
   8 This program is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public License
  10 as published by the Free Software Foundation; either version 2
  11 of the License, or (at your option) any later version.
  12
  13 This program is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 GNU General Public License for more details.
  17
  18 You should have received a copy of the GNU General Public License
  19 along with this program; if not, write to the Free Software
  20 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  21 */
  22
  23 #include "linearsystem.h"
  24
  25 #include <assert.h>
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28 #include "defs.h"
  29
  30 //                                                                            //
  31
  32 LinearSystem createLinearSystem()
  33 {
  34         LinearSystem system = (LinearSystem)malloc(sizeof(linearsystem_t));
  35
  36         if (system==NULL)
  37         {
  38                 fprintf(stderr, "Fatal Error (%s/%d): Could not allocate memory for LinearSystem!\n", __FILE__, __LINE__);
  39                 exit(1);
  40         }
  41
  42         system->Equations = 0;
  43         system->Variables = 0;
  44         system->A = NULL;
  45         system->b = NULL;
  46         system->EqProperties = NULL;
  47         system->VarProperties = NULL;
  48
  49         return system;
  50 }
  51
  52 //                                                                            //
  53
  54 void deleteLinearSystem(LinearSystem system)
  55 {
  56         int i;
  57
  58         if (system)
  59         {
  60                 for (i=0; i<system->Variables; i++)
  61                         deleteVector(system->A[i]);
  62                 free(system->A);
  63                 deleteVector(system->b);
  64                 deleteVariableProperties(system->VarProperties);
  65                 if (system->EqProperties)
  66                         free(system->EqProperties);
  67                 free(system);
  68         }
  69 }
  70
  71 //                                                                            //
  72
  73 void setLinearSystemSize(LinearSystem system, int variables, int equations)
  74 {
  75         int i;
  76
  77         assert(system);
  78         assert(variables>=0);
  79         assert(equations>=0);
  80
  81         if (variables!=system->Variables || equations!=system->Equations)
  82         {
  83                 if (system->A!=NULL)
  84                 {
  85                         for (i=0; i<system->Variables; i++)
  86                                 deleteVector(system->A[i]);
  87                         free(system->A);
  88                         system->A = NULL;
  89                 }
  90
  91                 if (variables>0)
  92                 {
  93                         system->A = (Vector *)malloc(variables*sizeof(Vector));
  94                         if (system->A==NULL)
  95                         {
  96                                 fprintf(stderr, "Fatal Error (%s/%d): Could not allocate memory for LinearSystem->A!\n", __FILE__, __LINE__);
  97                                 exit(1);
  98                         }
  99                         for (i=0; i<variables; i++)
 100                                 system->A[i] = createVector(equations);
 101                 }
 102         }
 103
 104         if (equations!=system->Equations)
 105         {
 106                 if (system->b!=NULL)
 107                 {
 108                         deleteVector(system->b);
 109                         system->b = NULL;
 110                         free(system->EqProperties);
 111                         system->EqProperties = NULL;
 112                 }
 113                 if (equations>0)
 114                 {
 115                         system->b = createVector(equations);
 116                         system->EqProperties = (EquationProperties)malloc(equations*sizeof(equationproperty_t));
 117                         if (system->EqProperties==NULL)
 118                         {
 119                                 fprintf(stderr, "Fatal Error (%s/%d): Could not allocate memory for LinearSystem->EqProperties!\n", __FILE__, __LINE__);
 120                                 exit(1);
 121                         }
 122                         for (i=0; i<equations; i++)
 123                         {
 124                                 system->EqProperties[i].Type = EQUATION_EQUAL;
 125                                 system->EqProperties[i].Modulus = 0;
 126                         }
 127                 }
 128         }
 129
 130         if (variables!=system->Variables)
 131         {
 132                 if (system->VarProperties!=NULL)
 133                 {
 134                         deleteVariableProperties(system->VarProperties);
 135                         system->VarProperties = NULL;
 136                 }
 137                 if (variables>0)
 138                         system->VarProperties = createVariableProperties(variables);
 139         }
 140
 141         system->Variables = variables;
 142         system->Equations = equations;
 143 }
 144
 145 //                                                                            //
 146
 147 void setLinearSystemMatrix(LinearSystem system, Matrix matrix)
 148 {
 149         int i,j;
 150
 151         assert(system);
 152         assert(matrix);
 153
 154         setLinearSystemSize(system, matrix->Width, matrix->Height);
 155
 156         for (i=0; i<system->Variables; i++)
 157                 for (j=0; j<system->Equations; j++)
 158                         system->A[i][j] = matrix->Data[i+j*system->Variables];
 159 }
 160
 161 //                                                                            //
 162
 163 void setLinearSystemRHS(LinearSystem system, Vector vector)
 164 {
 165         assert(system);
 166         assert(vector);
 167
 168         if (system->b!=NULL)
 169         {
 170                 deleteVector(system->b);
 171                 system->b = NULL;
 172         }
 173
 174         if (system->Equations)
 175                 system->b = copyVector(vector, system->Equations);
 176 }
 177
 178 //                                                                            //
 179
 180 void setLinearSystemLimit(LinearSystem system, int id, int lower, int upper, bool free)
 181 {
 182         assert(system);
 183         assert(lower<=0);
 184         assert(upper>=0);
 185
 186         if (id>=0 && id<system->Variables)
 187         {
 188                 system->VarProperties[id].Free = free;
 189                 system->VarProperties[id].Lower = lower;
 190                 system->VarProperties[id].Upper = upper;
 191         }
 192         else if (id==-1)
 193         {
 194                 for (id=0; id<system->Variables; id++)
 195                 {
 196                         system->VarProperties[id].Free = free;
 197                         system->VarProperties[id].Lower = lower;
 198                         system->VarProperties[id].Upper = upper;
 199                 }
 200         }
 201 }
 202
 203 //                                                                            //
 204
 205 void setLinearSystemBound(LinearSystem system, int id, char type, int value)
 206 {
 207         assert(system);
 208         assert(id>=0);
 209         assert(id<system->Variables);
 210         assert(type == 'u' || type == 'l');
 211
 212         if (type == 'u')
 213                 system->VarProperties[id].Upper = value;
 214         else if (type == 'l')
 215                 system->VarProperties[id].Lower = value;
 216 }
 217
 218 //                                                                            //
 219
 220 void setLinearSystemEquationType(LinearSystem system, int id, EquationType type, int modulus)
 221 {
 222         assert(system);
 223         assert(type<=5);
 224         assert(type || EQUATION_MODULO || modulus!=0);
 225
 226         if (id>=0)
 227         {
 228                 system->EqProperties[id].Type = type;
 229                 system->EqProperties[id].Modulus = modulus;
 230         }
 231         else if (id==-1)
 232         {
 233                 for (id=0; id<system->Equations; id++)
 234                 {
 235                         system->EqProperties[id].Type = type;
 236                         system->EqProperties[id].Modulus = modulus;
 237                 }
 238         }
 239 }
 240
 241 //                                                                            //
 242
 243 int imax(int a, int b)
 244 {
 245         return a > b ? a : b;
 246 }
 247
 248 //                                                                            //
 249
 250 int numberSize(int value)
 251 {
 252         int result = 1;
 253
 254         if (value<0)
 255         {
 256                 value *= -1;
 257                 result++;
 258         }
 259
 260         while (value>=10)
 261         {
 262                 value /= 10;
 263                 result++;
 264         }
 265
 266         return result;
 267 }
 268
 269 //                                                                            //
 270
 271 int equationSize(EquationType type)
 272 {
 273         if (type<4)
 274                 return 1;
 275         else
 276                 return 2;
 277 }
 278
 279 //                                                                            //
 280
 281 int propertySize(variableproperty_t property)
 282 {
 283         return (property.Free ? 1 : imax(
 284                         (property.Lower==-MAXINT ? 1 : numberSize(property.Lower)),
 285                         (property.Upper==MAXINT ? 1 : numberSize(property.Upper))
 286                         ));
 287 }
 288
 289 //                                                                            //
 290
 291 void fprintLinearSystem(FILE *stream, LinearSystem system)
 292 {
 293         Vector space;
 294         int i,j;
 295         bool allfree = true;
 296
 297         if (!stream || !system)
 298         {
 299                 fprintf(stderr, "Fatal Error (%s/%d): printLinearSystem was called with wrong arguments!\n", __FILE__, __LINE__);
 300                 exit(1);
 301         }
 302
 303         space = createVector(system->Variables+2);
 304
 305         // Collect space need by all columns
 306         for (i=0; i<system->Variables; i++)
 307         {
 308                 if (!checkVariableFree(system->VarProperties, i))
 309                         allfree = false;
 310                 space[i] = propertySize(system->VarProperties[i]);
 311                 for (j=0; j<system->Equations; j++)
 312                         space[i] = imax(space[i], numberSize(system->A[i][j]));
 313         }
 314         space[system->Variables] = 0;
 315         space[system->Variables+1] = 0;
 316         for (i=0; i < system->Equations; i++)
 317         {
 318                 space[system->Variables] = imax(space[system->Variables], equationSize(system->EqProperties[i].Type));
 319                 space[system->Variables+1] = imax(space[system->Variables+1], numberSize(system->b[i]));
 320         }
 321
 322         // Print header
 323         if (!allfree)
 324         {
 325                 for (i=0; i<system->Variables; i++)
 326                 {
 327                         if (checkVariableFree(system->VarProperties, i))
 328                                 fprintf(stream, "%*s ", space[i], "");
 329                         else if (system->VarProperties[i].Upper==MAXINT)
 330                                 fprintf(stream, "%*s ", space[i], "+");
 331                         else
 332                                 fprintf(stream, "%*d ", space[i], system->VarProperties[i].Upper);
 333                 }
 334                 fprintf(stream, "\n");
 335                 for (i=0; i<system->Variables; i++)
 336                 {
 337                         if (checkVariableFree(system->VarProperties, i))
 338                                 fprintf(stream, "%*s ", space[i], "");
 339                         else if (system->VarProperties[i].Lower==-MAXINT)
 340                                 fprintf(stream, "%*s ", space[i], "-");
 341                         else
 342                                 fprintf(stream, "%*d ", space[i], system->VarProperties[i].Lower);
 343                 }
 344                 fprintf(stream, "\n\n");
 345         }
 346
 347         for (i=0; i<system->Equations; i++)
 348         {
 349                 for (j=0; j<system->Variables; j++)
 350                         fprintf(stream, "%*d ", space[j], system->A[j][i]);
 351                 switch (system->EqProperties[i].Type)
 352                 {
 353                         case EQUATION_EQUAL:
 354                         case EQUATION_MODULO:
 355                                 fprintf(stream, "%*s ", space[system->Variables], "=");
 356                         break;
 357                         case EQUATION_LESSER:
 358                                 fprintf(stream, "%*s ", space[system->Variables], "<");
 359                         break;
 360                         case EQUATION_LESSEREQUAL:
 361                                 fprintf(stream, "%*s ", space[system->Variables], "<=");
 362                         break;
 363                         case EQUATION_GREATER:
 364                                 fprintf(stream, "%*s ", space[system->Variables], ">");
 365                         break;
 366                         case EQUATION_GREATEREQUAL:
 367                                 fprintf(stream, "%*s ", space[system->Variables], ">=");
 368                         break;
 369                 }
 370                 fprintf(stream, "%*d", space[system->Variables+1], system->b[i]);
 371                 if (system->EqProperties[i].Type==EQUATION_MODULO)
 372                         fprintf(stream, " (mod %d)", system->EqProperties[i].Modulus);
 373                 fprintf(stream, "\n");
 374         }
 375
 376         deleteVector(space);
 377 }
 378
 379 //                                                                            //
 380
 381 void printLinearSystem(LinearSystem system)
 382 {
 383         fprintLinearSystem(stdout, system);
 384 }
 385
 386 //                                                                            //
 387
 388 LinearSystem homogenizeLinearSystem(LinearSystem old)
 389 {
 390         LinearSystem new;
 391         int i,j,k;
 392
 393         if (old==NULL)
 394                 return NULL;
 395
 396         new = createLinearSystem();
 397
 398         new->Equations = old->Equations;
 399         new->Variables = old->Variables;
 400
 401         // RHS stays temporarily the same
 402         new->b = copyVector(old->b, old->Equations);
 403
 404         // Original matrix with properties
 405         new->A = (Vector *)malloc(new->Variables*sizeof(Vector));
 406         new->VarProperties = createVariableProperties(new->Variables);
 407         for (i=0; i<new->Variables; i++)
 408         {
 409                 new->A[i] = copyVector(old->A[i], old->Equations);
 410                 new->VarProperties[i] = old->VarProperties[i];
 411         }
 412
 413         // transformation
 414         for (i=0; i<old->Equations; i++)
 415         {
 416                 if (old->EqProperties[i].Type==EQUATION_LESSER)
 417                 {
 418                         new->b[i]--;
 419                 }
 420                 if (old->EqProperties[i].Type==EQUATION_GREATER)
 421                         new->b[i]++;
 422                 switch (old->EqProperties[i].Type)
 423                 {
 424                         case EQUATION_LESSER:
 425                         case EQUATION_LESSEREQUAL:
 426                                 k = 1;
 427                         break;
 428                         case EQUATION_GREATER:
 429                         case EQUATION_GREATEREQUAL:
 430                                 k = -1;
 431                         break;
 432                         case EQUATION_MODULO:
 433                                 k = old->EqProperties[i].Modulus;
 434                         break;
 435                         default:
 436                                 k = 0;
 437                         break;
 438                 }
 439                 if (k!=0)
 440                 {
 441                         new->Variables++;
 442                         new->A = (Vector *)realloc(new->A, new->Variables*sizeof(Vector));
 443                         new->VarProperties = (VariableProperties)realloc(new->VarProperties, new->Variables*sizeof(variableproperty_t));
 444                         new->A[new->Variables-1] = createVector(new->Equations);
 445                         for (j=0; j<new->Equations; j++)
 446                                 new->A[new->Variables-1][j] = (j==i ? k : 0);
 447                         new->VarProperties[new->Variables-1].Free = (old->EqProperties[i].Type==EQUATION_MODULO);
 448                         new->VarProperties[new->Variables-1].Upper = MAXINT;
 449                         new->VarProperties[new->Variables-1].Lower = (old->EqProperties[i].Type==EQUATION_MODULO ? -MAXINT : 0);
 450                         new->VarProperties[new->Variables-1].Column = -1;
 451                 }
 452         }
 453
 454         // -rhs as a new row in [0;1]
 455         if (normVector(new->b, old->Equations)>0)
 456         {
 457                 new->Variables++;
 458                 new->A = (Vector *)realloc(new->A, new->Variables*sizeof(Vector));
 459                 new->VarProperties = (VariableProperties)realloc(new->VarProperties, new->Variables*sizeof(variableproperty_t));
 460
 461                 new->A[new->Variables-1] = new->b;
 462                 for (i=0; i<old->Equations; i++)
 463                         new->A[new->Variables-1][i] *= -1;
 464                 new->VarProperties[new->Variables-1].Free = false;
 465                 new->VarProperties[new->Variables-1].Upper = 1;
 466                 new->VarProperties[new->Variables-1].Lower = 0;
 467                 new->VarProperties[new->Variables-1].Column = -2;
 468
 469                 // recreate rhs as a zero-vector
 470                 new->b = createVector(new->Equations);
 471                 for (i=0; i<new->Equations; i++)
 472                         new->b[i] = 0;
 473         }
 474
 475         // all are now equaitons
 476         new->EqProperties = (EquationProperties)malloc(new->Equations*sizeof(equationproperty_t));
 477         for (i=0; i<new->Equations; i++)
 478         {
 479                 new->EqProperties[i].Type = EQUATION_EQUAL;
 480                 new->EqProperties[i].Modulus = 0;
 481         }
 482
 483         return new;
 484 }
 485
 486 //                                                                            //