Changeset db4770 in git for factory/cf_linsys.cc
 Timestamp:
 Jul 8, 1996, 10:22:51 AM (27 years ago)
 Branches:
 (u'spielwiese', 'e7cc1ebecb61be8b9ca6c18016352af89940b21a')
 Children:
 0af18dde5af4397ff45cd495f0dbf94f4f78bb46
 Parents:
 7b4bfe63f5e8232ccce6c9f00ab9238ceabc1041
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/cf_linsys.cc
r7b4bfe6 rdb4770 1 1 // emacs edit mode for this file is * C++ * 2 // $Id: cf_linsys.cc,v 1. 0 19960517 10:59:44stobbe Exp $2 // $Id: cf_linsys.cc,v 1.1 19960708 08:22:51 stobbe Exp $ 3 3 4 4 /* 5 5 $Log: not supported by cvs2svn $ 6 Revision 1.0 1996/05/17 10:59:44 stobbe 7 Initial revision 8 6 9 */ 7 10 … … 17 20 18 21 static bool solve ( int **extmat, int nrows, int ncols ); 22 int determinant ( int **extmat, int n ); 19 23 20 24 static CanonicalForm bound ( CanonicalForm ** M, int rows, int cols ); 21 25 CanonicalForm detbound ( const CFMatrix & M, int rows ); 26 27 bool matrix_in_Z( const CFMatrix & M, int rows ) 28 { 29 int i, j; 30 for ( i = 1; i <= rows; i++ ) 31 for ( j = 1; j <= rows; j++ ) 32 if ( ! M(i,j).inZ() ) 33 return false; 34 return true; 35 } 36 37 bool betterpivot ( const CanonicalForm & oldpivot, const CanonicalForm & newpivot ) 38 { 39 if ( newpivot.isZero() ) 40 return false; 41 else if ( oldpivot.isZero() ) 42 return true; 43 else if ( level( oldpivot ) > level( newpivot ) ) 44 return true; 45 else if ( level( oldpivot ) < level( newpivot ) ) 46 return false; 47 else 48 return ( newpivot.lc() < oldpivot.lc() ); 49 } 50 51 bool fuzzy_result; 22 52 23 53 void … … 62 92 // Q so far 63 93 Q = p; 64 while ( Q < B ) {94 while ( Q < B && pno < cf_getNumBigPrimes() ) { 65 95 do { 66 96 cout << "trying prime(" << pno << ") = " << flush; … … 85 115 Q = qnew; 86 116 } 117 if ( pno == cf_getNumBigPrimes() ) 118 fuzzy_result = true; 119 else 120 fuzzy_result = false; 87 121 // store the result in M 88 122 Qhalf = Q / 2; … … 98 132 delete [] MM; 99 133 delete [] mm; 134 } 135 136 CanonicalForm 137 determinant( const CFMatrix & M, int rows ) 138 { 139 ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" ); 140 if ( rows == 1 ) 141 return M(1,1); 142 else if ( rows == 2 ) 143 return M(1,1)*M(2,2)M(2,1)*M(1,2); 144 else if ( matrix_in_Z( M, rows ) ) { 145 int ** mm = new (int*)[rows]; 146 CanonicalForm Q, Qhalf, mnew, qnew, B; 147 CanonicalForm det, detnew; 148 int i, j, p, pno, intdet; 149 bool ok; 150 151 // initialize room to hold the result and the result mod p 152 for ( i = 0; i < rows; i++ ) { 153 mm[i] = new int[rows]; 154 } 155 156 // calculate the bound for the result 157 B = detbound( M, rows ); 158 159 // find a first solution mod p 160 pno = 0; 161 do { 162 p = cf_getBigPrime( pno ); 163 setCharacteristic( p ); 164 ok = true; 165 // map matrix into char p 166 for ( i = 0; i < rows && ok; i++ ) 167 for ( j = 0; j < rows && ok; j++ ) { 168 if ( M(i+1,j+1).isZero() ) 169 mm[i][j] = 0; 170 else { 171 mm[i][j] = mapinto( M(i+1,j+1) ).intval(); 172 ok = mm[i][j] != 0; 173 } 174 } 175 pno++; 176 } while ( ! ok && pno < cf_getNumPrimes() ); 177 // initialize the result matrix with first solution 178 // solve mod p 179 intdet = determinant( mm, rows ); 180 setCharacteristic( 0 ); 181 det = intdet; 182 // Q so far 183 Q = p; 184 while ( Q < B && cf_getNumPrimes() > pno ) { 185 do { 186 p = cf_getBigPrime( pno ); 187 setCharacteristic( p ); 188 ok = true; 189 // map matrix into char p 190 for ( i = 0; i < rows && ok; i++ ) 191 for ( j = 0; j < rows && ok; j++ ) { 192 if ( M(i+1,j+1).isZero() ) 193 mm[i][j] = 0; 194 else { 195 mm[i][j] = mapinto( M(i+1,j+1) ).intval(); 196 ok = mm[i][j] != 0; 197 } 198 } 199 pno++; 200 } while ( ! ok && cf_getNumPrimes() > pno ); 201 // solve mod p 202 intdet = determinant( mm, rows ); 203 // found a solution mod p 204 // now chinese remainder it to a solution mod Q*p 205 setCharacteristic( 0 ); 206 chineseRemainder( det, Q, intdet, p, detnew, qnew ); 207 det = detnew; 208 Q = qnew; 209 } 210 if ( ! ok ) 211 fuzzy_result = true; 212 else 213 fuzzy_result = false; 214 // store the result in M 215 Qhalf = Q / 2; 216 if ( det > Qhalf ) 217 det = det  Q; 218 for ( i = 0; i < rows; i++ ) 219 delete [] mm[i]; 220 delete [] mm; 221 return det; 222 } 223 else { 224 CFMatrix m( M ); 225 CanonicalForm divisor = 1, pivot, mji; 226 int i, j, k, sign = 1; 227 for ( i = 1; i <= rows; i++ ) { 228 pivot = m(i,i); k = i; 229 for ( j = i+1; j <= rows; j++ ) { 230 if ( betterpivot( pivot, m(j,i) ) ) { 231 pivot = m(j,i); 232 k = j; 233 } 234 } 235 if ( pivot.isZero() ) 236 return 0; 237 if ( i != k ) { 238 m.swapRow( i, k ); 239 sign = sign; 240 } 241 for ( j = i+1; j <= rows; j++ ) { 242 if ( ! m(j,i).isZero() ) { 243 divisor *= pivot; 244 mji = m(j,i); 245 m(j,i) = 0; 246 for ( k = i+1; k <= rows; k++ ) 247 m(j,k) = m(j,k) * pivot  m(i,k)*mji; 248 } 249 } 250 } 251 pivot = sign; 252 for ( i = 1; i <= rows; i++ ) 253 pivot *= m(i,i); 254 return pivot / divisor; 255 } 100 256 } 101 257 … … 116 272 } 117 273 sum += vmax; 118 return sum; 274 return sqrt( sum ); 275 } 276 277 278 CanonicalForm 279 detbound ( const CFMatrix & M, int rows ) 280 { 281 CanonicalForm sum = 0, prod = 1; 282 int i, j; 283 for ( i = 1; i <= rows; i++ ) { 284 sum = 0; 285 for ( j = 1; j <= rows; j++ ) 286 sum += M(i,j) * M(i,j); 287 prod *= sum; 288 } 289 return 2 * sqrt( prod ); 119 290 } 120 291 … … 168 339 } 169 340 341 int 342 determinant ( int **extmat, int n ) 343 { 344 int i, j, k; 345 int divisor, multiplier, rowii, rowji; // all FF 346 int * rowi; // FF 347 int * rowj; // FF 348 int * swap; // FF 349 // triangularization 350 multiplier = 1; 351 divisor = 1; 352 353 for ( i = 0; i < n; i++ ) { 354 //find "pivot" 355 for (j = i; j < n; j++ ) 356 if ( extmat[j][i] != 0 ) break; 357 if ( j == n ) return 0; 358 if ( j != i ) { 359 multiplier = ff_neg( multiplier ); 360 swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap; 361 } 362 rowi = extmat[i]; 363 rowii = rowi[i]; 364 for ( j = i+1; j < n; j++ ) { 365 rowj = extmat[j]; 366 rowji = rowj[i]; 367 if ( rowji == 0 ) continue; 368 divisor = ff_mul( divisor, rowii ); 369 for ( k = i; k < n; k++ ) 370 rowj[k] = ff_sub( ff_mul( rowj[k], rowii ), ff_mul( rowi[k], rowji ) ); 371 } 372 } 373 multiplier = ff_mul( multiplier, ff_inv( divisor ) ); 374 for ( i = 0; i < n; i++ ) 375 multiplier = ff_mul( multiplier, extmat[i][i] ); 376 return multiplier; 377 } 378 170 379 void 171 380 solveVandermondeT ( const CFArray & a, const CFArray & w, CFArray & x, const Variable & z )
Note: See TracChangeset
for help on using the changeset viewer.