Changeset b4e057 in git
 Timestamp:
 Jul 15, 1996, 10:33:18 AM (27 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 9f2b6f23a8e1f2de750bc3241e6e692d3f57c0d1
 Parents:
 7155e15ec589aec0c180c45fd0f345978d42c34b
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/cf_linsys.cc
r7155e1 rb4e057 1 1 // emacs edit mode for this file is * C++ * 2 // $Id: cf_linsys.cc,v 1. 1 19960708 08:22:51stobbe Exp $2 // $Id: cf_linsys.cc,v 1.2 19960715 08:33:18 stobbe Exp $ 3 3 4 4 /* 5 5 $Log: not supported by cvs2svn $ 6 Revision 1.1 1996/07/08 08:22:51 stobbe 7 "New function determinant. 8 " 9 6 10 Revision 1.0 1996/05/17 10:59:44 stobbe 7 11 Initial revision … … 22 26 int determinant ( int **extmat, int n ); 23 27 24 static CanonicalForm bound ( CanonicalForm ** M, int rows, int cols);28 static CanonicalForm bound ( const CFMatrix & M ); 25 29 CanonicalForm detbound ( const CFMatrix & M, int rows ); 26 30 27 bool matrix_in_Z( const CFMatrix & M, int rows ) 31 bool 32 matrix_in_Z( const CFMatrix & M, int rows ) 28 33 { 29 34 int i, j; … … 35 40 } 36 41 37 bool betterpivot ( const CanonicalForm & oldpivot, const CanonicalForm & newpivot ) 42 bool 43 matrix_in_Z( const CFMatrix & M ) 44 { 45 int i, j, rows = M.rows(), cols = M.columns(); 46 for ( i = 1; i <= rows; i++ ) 47 for ( j = 1; j <= cols; j++ ) 48 if ( ! M(i,j).inZ() ) 49 return false; 50 return true; 51 } 52 53 bool 54 betterpivot ( const CanonicalForm & oldpivot, const CanonicalForm & newpivot ) 38 55 { 39 56 if ( newpivot.isZero() ) … … 51 68 bool fuzzy_result; 52 69 53 void 54 linearSystemSolve( CanonicalForm ** M, int rows, int cols ) 55 { 56 CanonicalForm ** MM = new (CanonicalForm*)[rows]; 57 int ** mm = new (int*)[rows]; 58 CanonicalForm Q, Qhalf, mnew, qnew, B; 59 int i, j, p, pno; 60 bool ok; 61 62 // initialize room to hold the result and the result mod p 63 for ( i = 0; i < rows; i++ ) { 64 MM[i] = new CanonicalForm[cols]; 65 mm[i] = new int[cols]; 66 } 67 68 // calculate the bound for the result 69 B = bound( M, rows, cols ); 70 cout << "bound = " << B << endl; 71 72 // find a first solution mod p 73 pno = 0; 74 do { 75 cout << "trying prime(" << pno << ") = " << flush; 76 p = cf_getBigPrime( pno ); 77 cout << p << endl; 78 setCharacteristic( p ); 79 // map matrix into char p 80 for ( i = 0; i < rows; i++ ) 81 for ( j = 0; j < cols; j++ ) 82 mm[i][j] = mapinto( M[i][j] ).intval(); 83 // solve mod p 84 ok = solve( mm, rows, cols ); 85 pno++; 86 } while ( ! ok ); 87 // initialize the result matrix with first solution 88 setCharacteristic( 0 ); 89 for ( i = 0; i < rows; i++ ) 90 for ( j = rows; j < cols; j++ ) 91 MM[i][j] = mm[i][j]; 92 // Q so far 93 Q = p; 94 while ( Q < B && pno < cf_getNumBigPrimes() ) { 70 bool 71 linearSystemSolve( CFMatrix & M ) 72 { 73 if ( ! matrix_in_Z( M ) ) { 74 int nrows = M.rows(), ncols = M.columns(); 75 int i, j, k; 76 CanonicalForm rowpivot, pivotrecip; 77 // triangularization 78 for ( i = 1; i <= nrows; i++ ) { 79 //find "pivot" 80 for (j = i; j <= nrows; j++ ) 81 if ( M(j,i) != 0 ) break; 82 if ( j > nrows ) return false; 83 if ( j != i ) 84 M.swapRow( i, j ); 85 pivotrecip = 1 / M(i,i); 86 for ( j = 1; j <= ncols; j++ ) 87 M(i,j) *= pivotrecip; 88 for ( j = i+1; j <= nrows; j++ ) { 89 rowpivot = M(j,i); 90 if ( rowpivot == 0 ) continue; 91 for ( k = i; k <= ncols; k++ ) 92 M(j,k) = M(i,k) * rowpivot; 93 } 94 } 95 // matrix is now upper triangular with 1s down the diagonal 96 // backsubstitute 97 for ( i = nrows1; i > 0; i ) { 98 for ( j = nrows+1; j <= ncols; j++ ) { 99 for ( k = i+1; k <= nrows; k++ ) 100 M(i,j) = M(k,j) * M(i,k); 101 } 102 } 103 return true; 104 } 105 else { 106 int rows = M.rows(), cols = M.columns(); 107 CFMatrix MM( rows, cols ); 108 int ** mm = new (int*)[rows]; 109 CanonicalForm Q, Qhalf, mnew, qnew, B; 110 int i, j, p, pno; 111 bool ok; 112 113 // initialize room to hold the result and the result mod p 114 for ( i = 0; i < rows; i++ ) { 115 mm[i] = new int[cols]; 116 } 117 118 // calculate the bound for the result 119 B = bound( M ); 120 DEBOUTLN( cerr, "bound = ", B ); 121 122 // find a first solution mod p 123 pno = 0; 95 124 do { 96 cout << "trying prime(" << pno << ") = " << flush;125 DEBOUT( cerr, "trying prime(", pno ); DEBOUTLN( cerr, ") = ", ' ' ); 97 126 p = cf_getBigPrime( pno ); 98 127 cout << p << endl; 99 128 setCharacteristic( p ); 129 // map matrix into char p 100 130 for ( i = 0; i < rows; i++ ) 101 131 for ( j = 0; j < cols; j++ ) 102 mm[i][j] = mapinto( M [i][j]).intval();132 mm[i][j] = mapinto( M(i,j) ).intval(); 103 133 // solve mod p 104 134 ok = solve( mm, rows, cols ); 105 135 pno++; 106 136 } while ( ! ok ); 107 // found a solution mod p 108 // now chinese remainder it to a solution mod Q*p 137 // initialize the result matrix with first solution 109 138 setCharacteristic( 0 ); 110 139 for ( i = 0; i < rows; i++ ) 111 for ( j = rows; j < cols; j++ ) { 112 chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i][j]), CanonicalForm(p), mnew, qnew ); 113 MM[i][j] = mnew; 114 } 115 Q = qnew; 116 } 117 if ( pno == cf_getNumBigPrimes() ) 118 fuzzy_result = true; 119 else 120 fuzzy_result = false; 121 // store the result in M 122 Qhalf = Q / 2; 123 for ( i = 0; i < rows; i++ ) { 124 for ( j = rows; j < cols; j++ ) 125 if ( MM[i][j] > Qhalf ) 126 M[i][j] = MM[i][j]  Q; 127 else 128 M[i][j] = MM[i][j]; 129 delete [] MM[i]; 130 delete [] mm[i]; 131 } 132 delete [] MM; 133 delete [] mm; 140 for ( j = rows; j < cols; j++ ) 141 MM(i,j) = mm[i][j]; 142 // Q so far 143 Q = p; 144 while ( Q < B && pno < cf_getNumBigPrimes() ) { 145 do { 146 cout << "trying prime(" << pno << ") = " << flush; 147 p = cf_getBigPrime( pno ); 148 cout << p << endl; 149 setCharacteristic( p ); 150 for ( i = 0; i < rows; i++ ) 151 for ( j = 0; j < cols; j++ ) 152 mm[i][j] = mapinto( M(i,j) ).intval(); 153 // solve mod p 154 ok = solve( mm, rows, cols ); 155 pno++; 156 } while ( ! ok ); 157 // found a solution mod p 158 // now chinese remainder it to a solution mod Q*p 159 setCharacteristic( 0 ); 160 for ( i = 0; i < rows; i++ ) 161 for ( j = rows; j < cols; j++ ) { 162 chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i][j]), CanonicalForm(p), mnew, qnew ); 163 MM(i,j) = mnew; 164 } 165 Q = qnew; 166 } 167 if ( pno == cf_getNumBigPrimes() ) 168 fuzzy_result = true; 169 else 170 fuzzy_result = false; 171 // store the result in M 172 Qhalf = Q / 2; 173 for ( i = 0; i < rows; i++ ) { 174 for ( j = rows; j < cols; j++ ) 175 if ( MM(i,j) > Qhalf ) 176 M(i,j) = MM(i,j)  Q; 177 else 178 M(i,j) = MM(i,j); 179 delete [] mm[i]; 180 } 181 delete [] mm; 182 return ! fuzzy_result; 183 } 134 184 } 135 185 … … 257 307 258 308 static CanonicalForm 259 bound ( CanonicalForm ** M, int rows, int cols ) 260 { 309 bound ( const CFMatrix & M ) 310 { 311 int rows = M.rows(), cols = M.columns(); 261 312 CanonicalForm sum = 0; 262 313 int i, j; 263 314 for ( i = 0; i < rows; i++ ) 264 315 for ( j = 0; j < rows; j++ ) 265 sum += M [i][j] * M[i][j];316 sum += M(i,j) * M(i,j); 266 317 CanonicalForm vmax = 0, vsum; 267 318 for ( j = rows; j < cols; j++ ) { 268 319 vsum = 0; 269 320 for ( i = 0; i < rows; i++ ) 270 vsum += M [i][j] * M[i][j];321 vsum += M(i,j) * M(i,j); 271 322 if ( vsum > vmax ) vmax = vsum; 272 323 }
Note: See TracChangeset
for help on using the changeset viewer.