zsolve/lattice.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 "lattice.h"
  24
  25 #include <assert.h>
  26 #include <stdlib.h>
  27 #include "defs.h"
  28
  29 //                                                                            //
  30
  31 VectorArray generateLattice(LinearSystem system)
  32 {
  33         int i,j,k;
  34         int n,e;
  35         int currentindex, bestindex, currentvalue, bestvalue;
  36         int identities;
  37         int factor;
  38         bool flag;
  39         VectorArray array;
  40         Matrix H,I,C;
  41         Vector vec;
  42
  43         assert(system);
  44         assert(normVector(system->b, system->Equations)==0);
  45
  46         n = system->Variables;
  47         e = system->Equations;
  48
  49         array = createVectorArray(n);
  50         for (i=0; i<n; i++)
  51                 array->Properties[i] = system->VarProperties[i];
  52
  53         // init H
  54         H = createMatrix(n, e);
  55         for (i=0; i<n; i++)
  56                 for (j=0; j<e; j++)
  57                         H->Data[i+n*j] = system->A[i][j];
  58
  59         // search for columns with norm 1 and use them to row-echelonize to lower right
  60         identities = 0;
  61         for (i=n-1; i>=0; i--)
  62         {
  63                 flag = false;
  64                 for (j=0; j<identities; j++)
  65                         if (H->Data[i+(j+e-identities)*n]!=0)
  66                                 flag = true;
  67                 if (flag)
  68                         continue;
  69                 k = -1;
  70                 for (j=0; j<e-identities; j++)
  71                 {
  72                         if (k==-1 && abs(H->Data[i+j*n])==1)
  73                                 k = j;
  74                         else if (H->Data[i+j*n]!=0)
  75                         {
  76                                 k = -2;
  77                                 break;
  78                         }
  79                 }
  80                 if (k>=0)
  81                 {
  82                         swapMatrixRows(H, k, e-identities-1);
  83                         swapVariableProperties(array->Properties, i, n-identities-1);
  84                         swapMatrixColumns(H, i, n-identities-1);
  85                         identities++;
  86                 }
  87         }
  88
  89         // copy H for later usage
  90         C = copyMatrix(H);
  91
  92         n -= identities;
  93         e -= identities;
  94
  95         if (n>0)
  96         {
  97                 // init I
  98                 I = createIdentityMatrix(n);
  99
 100                 // hermite the matrix H[identities .. n-1][identities .. e-1]
 101
 102                 for (i=0; i<minm(n,e); i++)
 103                 {
 104                         // find gcd-minimal row
 105                         bestindex = i;
 106                         bestvalue = 0;
 107                         for (currentindex=i; currentindex<e; currentindex++)
 108                         {
 109                                 currentvalue = gcdVector(H->Data+i+H->Width*currentindex, n-i);
 110                                 if (currentvalue!=0 && (bestvalue==0 || currentvalue<bestvalue))
 111                                 {
 112                                         bestvalue = currentvalue;
 113                                         bestindex = currentindex;
 114                                 }
 115                         }
 116                         if (bestvalue==0)
 117                                 break;
 118
 119
 120                         // swap this with i-th
 121                         swapMatrixRows(H, bestindex, i);
 122
 123                         do
 124                         {
 125                                 // find reducer
 126                                 bestindex = i;
 127                                 bestvalue = 0;
 128                                 for (currentindex=i; currentindex<n; currentindex++)
 129                                 {
 130                                         currentvalue = abs(H->Data[currentindex+H->Width*i]);
 131                                         if (currentvalue!=0 && (bestvalue<=0 || currentvalue<bestvalue))
 132                                         {
 133                                                 bestindex = currentindex;
 134                                                 bestvalue = currentvalue;
 135                                         }
 136                                 }
 137                                 assert(bestvalue!=0);
 138
 139                                 // if negative, change sign of column
 140                                 if (H->Data[bestindex+H->Width*i]<0)
 141                                 {
 142                                         negateMatrixColumn(H, bestindex);
 143                                         negateMatrixColumn(I, bestindex);
 144                                 }
 145
 146                                 // reduce
 147                                 flag = false;
 148                                 for (j=0; j<n; j++)
 149                                 {
 150                                         if (j!=bestindex)
 151                                         {
 152                                                 factor = H->Data[j+H->Width*i] / -bestvalue;
 153                                                 if (H->Data[j+H->Width*i] % bestvalue != 0 && j>i)
 154                                                         flag = true;
 155                                                 combineMatrixColumns(H, j, factor, bestindex);
 156                                                 combineMatrixColumns(I, j, factor, bestindex);
 157                                         }
 158                                 }
 159                         } while (flag);
 160
 161                         if (bestindex!=i)
 162                         {
 163                                 swapMatrixColumns(H, i, bestindex);
 164                                 swapMatrixColumns(I, i, bestindex);
 165                         }
 166                 }
 167
 168                 while (i<n)
 169                 {
 170                         vec = createVector(H->Width);
 171                         // fill from I
 172                         for (j=0; j<n; j++)
 173                                 vec[j] = I->Data[i+j*n];
 174                         // fill identity part from system
 175                         for (j=0; j<identities; j++)
 176                         {
 177                                 factor = 0;
 178                                 for (k=0; k<n; k++)
 179                                         factor -= vec[k]*C->Data[C->Width*(e+j)+k];
 180                                 vec[n+j] = factor*H->Data[H->Width*(e+j)+n+j];
 181                         }
 182                         appendToVectorArray(array, vec);
 183                         i++;
 184                 }
 185                 deleteMatrix(I);
 186         }
 187
 188         deleteMatrix(H);
 189         deleteMatrix(C);
 190
 191         return array;
 192 }
 193
 194 //                                                                            //