My Project
Loading...
Searching...
No Matches
cf_linsys.cc
Go to the documentation of this file.
1/* emacs edit mode for this file is -*- C++ -*- */
2
3/**
4 * @file cf_linsys.cc
5 *
6 * solve linear systems and compute determinants of matrices
7**/
8
9#include "config.h"
10
11
12#include "cf_assert.h"
13#include "debug.h"
14#include "timing.h"
15
16#include "cf_defs.h"
17#include "cf_primes.h"
18#include "canonicalform.h"
19#include "cf_iter.h"
20#include "cf_algorithm.h"
21#include "ffops.h"
22#include "cf_primes.h"
23
24
26TIMING_DEFINE_PRINT(det_determinant)
27TIMING_DEFINE_PRINT(det_chinese)
28TIMING_DEFINE_PRINT(det_bound)
29TIMING_DEFINE_PRINT(det_numprimes)
30
31
32static bool solve ( int **extmat, int nrows, int ncols );
33int determinant ( int **extmat, int n );
34
36CanonicalForm detbound ( const CFMatrix & M, int rows );
37
38bool
40{
41 int i, j;
42 for ( i = 1; i <= rows; i++ )
43 for ( j = 1; j <= rows; j++ )
44 if ( ! M(i,j).inZ() )
45 return false;
46 return true;
47}
48
49bool
51{
52 int i, j, rows = M.rows(), cols = M.columns();
53 for ( i = 1; i <= rows; i++ )
54 for ( j = 1; j <= cols; j++ )
55 if ( ! M(i,j).inZ() )
56 return false;
57 return true;
58}
59
60bool
61betterpivot ( const CanonicalForm & oldpivot, const CanonicalForm & newpivot )
62{
63 if ( newpivot.isZero() )
64 return false;
65 else if ( oldpivot.isZero() )
66 return true;
67 else if ( level( oldpivot ) > level( newpivot ) )
68 return true;
69 else if ( level( oldpivot ) < level( newpivot ) )
70 return false;
71 else
72 return ( newpivot.lc() < oldpivot.lc() );
73}
74
76
77bool
79{
80 typedef int* int_ptr;
81
82 if ( ! matrix_in_Z( M ) ) {
83 int nrows = M.rows(), ncols = M.columns();
84 int i, j, k;
85 CanonicalForm rowpivot, pivotrecip;
86 // triangularization
87 for ( i = 1; i <= nrows; i++ ) {
88 //find "pivot"
89 for (j = i; j <= nrows; j++ )
90 if ( M(j,i) != 0 ) break;
91 if ( j > nrows ) return false;
92 if ( j != i )
93 M.swapRow( i, j );
94 pivotrecip = 1 / M(i,i);
95 for ( j = 1; j <= ncols; j++ )
96 M(i,j) *= pivotrecip;
97 for ( j = i+1; j <= nrows; j++ ) {
98 rowpivot = M(j,i);
99 if ( rowpivot == 0 ) continue;
100 for ( k = i; k <= ncols; k++ )
101 M(j,k) -= M(i,k) * rowpivot;
102 }
103 }
104 // matrix is now upper triangular with 1s down the diagonal
105 // back-substitute
106 for ( i = nrows-1; i > 0; i-- ) {
107 for ( j = nrows+1; j <= ncols; j++ ) {
108 for ( k = i+1; k <= nrows; k++ )
109 M(i,j) -= M(k,j) * M(i,k);
110 }
111 }
112 return true;
113 }
114 else {
115 int rows = M.rows(), cols = M.columns();
116 CFMatrix MM( rows, cols );
117 int ** mm = new int_ptr[rows];
118 CanonicalForm Q, Qhalf, mnew, qnew, B;
119 int i, j, p, pno;
120 bool ok;
121
122 // initialize room to hold the result and the result mod p
123 for ( i = 0; i < rows; i++ ) {
124 mm[i] = new int[cols];
125 }
126
127 // calculate the bound for the result
128 B = bound( M );
129 DEBOUTLN( cerr, "bound = " << B );
130
131 // find a first solution mod p
132 pno = 0;
133 do {
134 DEBOUTSL( cerr );
135 DEBOUT( cerr, "trying prime(" << pno << ") = " );
136 p = cf_getBigPrime( pno );
137 DEBOUT( cerr, p );
138 DEBOUTENDL( cerr );
140 // map matrix into char p
141 for ( i = 1; i <= rows; i++ )
142 for ( j = 1; j <= cols; j++ )
143 mm[i-1][j-1] = mapinto( M(i,j) ).intval();
144 // solve mod p
145 ok = solve( mm, rows, cols );
146 pno++;
147 } while ( ! ok );
148
149 // initialize the result matrix with first solution
151 for ( i = 1; i <= rows; i++ )
152 for ( j = rows+1; j <= cols; j++ )
153 MM(i,j) = mm[i-1][j-1];
154
155 // Q so far
156 Q = p;
157 while ( Q < B && pno < cf_getNumBigPrimes() ) {
158 do {
159 DEBOUTSL( cerr );
160 DEBOUT( cerr, "trying prime(" << pno << ") = " );
161 p = cf_getBigPrime( pno );
162 DEBOUT( cerr, p );
163 DEBOUTENDL( cerr );
165 for ( i = 1; i <= rows; i++ )
166 for ( j = 1; j <= cols; j++ )
167 mm[i-1][j-1] = mapinto( M(i,j) ).intval();
168 // solve mod p
169 ok = solve( mm, rows, cols );
170 pno++;
171 } while ( ! ok );
172 // found a solution mod p
173 // now chinese remainder it to a solution mod Q*p
175 for ( i = 1; i <= rows; i++ )
176 for ( j = rows+1; j <= cols; j++ )
177 {
178 chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i-1][j-1]), CanonicalForm(p), mnew, qnew );
179 MM(i, j) = mnew;
180 }
181 Q = qnew;
182 }
183 if ( pno == cf_getNumBigPrimes() )
184 fuzzy_result = true;
185 else
186 fuzzy_result = false;
187 // store the result in M
188 Qhalf = Q / 2;
189 for ( i = 1; i <= rows; i++ ) {
190 for ( j = rows+1; j <= cols; j++ )
191 if ( MM(i,j) > Qhalf )
192 M(i,j) = MM(i,j) - Q;
193 else
194 M(i,j) = MM(i,j);
195 delete [] mm[i-1];
196 }
197 delete [] mm;
198 return ! fuzzy_result;
199 }
200}
201
202static bool
203fill_int_mat( const CFMatrix & M, int ** m, int rows )
204{
205 int i, j;
206 bool ok = true;
207 for ( i = 0; i < rows && ok; i++ )
208 for ( j = 0; j < rows && ok; j++ )
209 {
210 if ( M(i+1,j+1).isZero() )
211 m[i][j] = 0;
212 else
213 {
214 m[i][j] = mapinto( M(i+1,j+1) ).intval();
215// ok = m[i][j] != 0;
216 }
217 }
218 return ok;
219}
220
222determinant( const CFMatrix & M, int rows )
223{
224 typedef int* int_ptr;
225
226 ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
227 if ( rows == 1 )
228 return M(1,1);
229 else if ( rows == 2 )
230 return M(1,1)*M(2,2)-M(2,1)*M(1,2);
231 else if ( matrix_in_Z( M, rows ) )
232 {
233 int ** mm = new int_ptr[rows];
234 CanonicalForm x, q, Qhalf, B;
235 int n, i, intdet, p, pno;
236 for ( i = 0; i < rows; i++ )
237 {
238 mm[i] = new int[rows];
239 }
240 pno = 0; n = 0;
241 TIMING_START(det_bound);
242 B = detbound( M, rows );
243 TIMING_END(det_bound);
244 q = 1;
245 TIMING_START(det_numprimes);
246 while ( B > q && n < cf_getNumBigPrimes() )
247 {
248 q *= cf_getBigPrime( n );
249 n++;
250 }
251 TIMING_END(det_numprimes);
252
253 CFArray X(1,n), Q(1,n);
254
255 while ( pno < n )
256 {
257 p = cf_getBigPrime( pno );
259 // map matrix into char p
260 TIMING_START(det_mapping);
261 fill_int_mat( M, mm, rows );
262 TIMING_END(det_mapping);
263 pno++;
264 DEBOUT( cerr, "." );
265 TIMING_START(det_determinant);
266 intdet = determinant( mm, rows );
267 TIMING_END(det_determinant);
269 X[pno] = intdet;
270 Q[pno] = p;
271 }
272 TIMING_START(det_chinese);
273 chineseRemainder( X, Q, x, q );
274 TIMING_END(det_chinese);
275 Qhalf = q / 2;
276 if ( x > Qhalf )
277 x = x - q;
278 for ( i = 0; i < rows; i++ )
279 delete [] mm[i];
280 delete [] mm;
281 return x;
282 }
283 else
284 {
285 CFMatrix m( M );
286 CanonicalForm divisor = 1, pivot, mji;
287 int i, j, k, sign = 1;
288 for ( i = 1; i <= rows; i++ ) {
289 pivot = m(i,i); k = i;
290 for ( j = i+1; j <= rows; j++ ) {
291 if ( betterpivot( pivot, m(j,i) ) ) {
292 pivot = m(j,i);
293 k = j;
294 }
295 }
296 if ( pivot.isZero() )
297 return 0;
298 if ( i != k )
299 {
300 m.swapRow( i, k );
301 sign = -sign;
302 }
303 for ( j = i+1; j <= rows; j++ )
304 {
305 if ( ! m(j,i).isZero() )
306 {
307 divisor *= pivot;
308 mji = m(j,i);
309 m(j,i) = 0;
310 for ( k = i+1; k <= rows; k++ )
311 m(j,k) = m(j,k) * pivot - m(i,k)*mji;
312 }
313 }
314 }
315 pivot = sign;
316 for ( i = 1; i <= rows; i++ )
317 pivot *= m(i,i);
318 return pivot / divisor;
319 }
320}
321
323determinant2( const CFMatrix & M, int rows )
324{
325 typedef int* int_ptr;
326
327 ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
328 if ( rows == 1 )
329 return M(1,1);
330 else if ( rows == 2 )
331 return M(1,1)*M(2,2)-M(2,1)*M(1,2);
332 else if ( matrix_in_Z( M, rows ) ) {
333 int ** mm = new int_ptr[rows];
334 CanonicalForm QQ, Q, Qhalf, mnew, q, qnew, B;
335 CanonicalForm det, detnew, qdet;
336 int i, p, pcount, pno, intdet;
337 bool ok;
338
339 // initialize room to hold the result and the result mod p
340 for ( i = 0; i < rows; i++ ) {
341 mm[i] = new int[rows];
342 }
343
344 // calculate the bound for the result
345 B = detbound( M, rows );
346
347 // find a first solution mod p
348 pno = 0;
349 do {
350 p = cf_getBigPrime( pno );
352 // map matrix into char p
353 ok = fill_int_mat( M, mm, rows );
354 pno++;
355 } while ( ! ok && pno < cf_getNumPrimes() );
356 // initialize the result matrix with first solution
357 // solve mod p
358 DEBOUT( cerr, "." );
359 intdet = determinant( mm, rows );
361 det = intdet;
362 // Q so far
363 Q = p;
364 QQ = p;
365 while ( Q < B && cf_getNumPrimes() > pno ) {
366 // find a first solution mod p
367 do {
368 p = cf_getBigPrime( pno );
370 // map matrix into char p
371 ok = fill_int_mat( M, mm, rows );
372 pno++;
373 } while ( ! ok && pno < cf_getNumPrimes() );
374 // initialize the result matrix with first solution
375 // solve mod p
376 DEBOUT( cerr, "." );
377 intdet = determinant( mm, rows );
379 qdet = intdet;
380 // Q so far
381 q = p;
382 QQ *= p;
383 pcount = 0;
384 while ( QQ < B && cf_getNumPrimes() > pno && pcount < 500 ) {
385 do {
386 p = cf_getBigPrime( pno );
388 ok = true;
389 // map matrix into char p
390 ok = fill_int_mat( M, mm, rows );
391 pno++;
392 } while ( ! ok && cf_getNumPrimes() > pno );
393 // solve mod p
394 DEBOUT( cerr, "." );
395 intdet = determinant( mm, rows );
396 // found a solution mod p
397 // now chinese remainder it to a solution mod Q*p
399 chineseRemainder( qdet, q, intdet, p, detnew, qnew );
400 qdet = detnew;
401 q = qnew;
402 QQ *= p;
403 pcount++;
404 }
405 DEBOUT( cerr, "*" );
406 chineseRemainder( det, Q, qdet, q, detnew, qnew );
407 Q = qnew;
408 QQ = Q;
409 det = detnew;
410 }
411 if ( ! ok )
412 fuzzy_result = true;
413 else
414 fuzzy_result = false;
415 // store the result in M
416 Qhalf = Q / 2;
417 if ( det > Qhalf )
418 det = det - Q;
419 for ( i = 0; i < rows; i++ )
420 delete [] mm[i];
421 delete [] mm;
422 return det;
423 }
424 else {
425 CFMatrix m( M );
426 CanonicalForm divisor = 1, pivot, mji;
427 int i, j, k, sign = 1;
428 for ( i = 1; i <= rows; i++ ) {
429 pivot = m(i,i); k = i;
430 for ( j = i+1; j <= rows; j++ ) {
431 if ( betterpivot( pivot, m(j,i) ) ) {
432 pivot = m(j,i);
433 k = j;
434 }
435 }
436 if ( pivot.isZero() )
437 return 0;
438 if ( i != k ) {
439 m.swapRow( i, k );
440 sign = -sign;
441 }
442 for ( j = i+1; j <= rows; j++ ) {
443 if ( ! m(j,i).isZero() ) {
444 divisor *= pivot;
445 mji = m(j,i);
446 m(j,i) = 0;
447 for ( k = i+1; k <= rows; k++ )
448 m(j,k) = m(j,k) * pivot - m(i,k)*mji;
449 }
450 }
451 }
452 pivot = sign;
453 for ( i = 1; i <= rows; i++ )
454 pivot *= m(i,i);
455 return pivot / divisor;
456 }
457}
458
459static CanonicalForm
460bound ( const CFMatrix & M )
461{
462 DEBINCLEVEL( cerr, "bound" );
463 int rows = M.rows(), cols = M.columns();
464 CanonicalForm sum = 0;
465 int i, j;
466 for ( i = 1; i <= rows; i++ )
467 for ( j = 1; j <= rows; j++ )
468 sum += M(i,j) * M(i,j);
469 DEBOUTLN( cerr, "bound(matrix)^2 = " << sum );
470 CanonicalForm vmax = 0, vsum;
471 for ( j = rows+1; j <= cols; j++ ) {
472 vsum = 0;
473 for ( i = 1; i <= rows; i++ )
474 vsum += M(i,j) * M(i,j);
475 if ( vsum > vmax ) vmax = vsum;
476 }
477 DEBOUTLN( cerr, "bound(lhs)^2 = " << vmax );
478 sum += vmax;
479 DEBOUTLN( cerr, "bound(overall)^2 = " << sum );
480 DEBDECLEVEL( cerr, "bound" );
481 return sqrt( sum ) + 1;
482}
483
484
486detbound ( const CFMatrix & M, int rows )
487{
488 CanonicalForm sum = 0, prod = 2;
489 int i, j;
490 for ( i = 1; i <= rows; i++ ) {
491 sum = 0;
492 for ( j = 1; j <= rows; j++ )
493 sum += M(i,j) * M(i,j);
494 prod *= 1 + sqrt(sum);
495 }
496 return prod;
497}
498
499
500// solve returns false if computation failed
501// extmat is overwritten: output is Id mat followed by solution(s)
502
503bool
504solve ( int **extmat, int nrows, int ncols )
505{
506 DEBINCLEVEL( cerr, "solve" );
507 int i, j, k;
508 int rowpivot, pivotrecip; // all FF
509 int * rowi; // FF
510 int * rowj; // FF
511 int * swap; // FF
512 // triangularization
513 for ( i = 0; i < nrows; i++ ) {
514 //find "pivot"
515 for (j = i; j < nrows; j++ )
516 if ( extmat[j][i] != 0 ) break;
517 if ( j == nrows ) {
518 DEBOUTLN( cerr, "solve failed" );
519 DEBDECLEVEL( cerr, "solve" );
520 return false;
521 }
522 if ( j != i ) {
523 swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
524 }
525 pivotrecip = ff_inv( extmat[i][i] );
526 rowi = extmat[i];
527 for ( j = 0; j < ncols; j++ )
528 rowi[j] = ff_mul( pivotrecip, rowi[j] );
529 for ( j = i+1; j < nrows; j++ ) {
530 rowj = extmat[j];
531 rowpivot = rowj[i];
532 if ( rowpivot == 0 ) continue;
533 for ( k = i; k < ncols; k++ )
534 rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
535 }
536 }
537 // matrix is now upper triangular with 1s down the diagonal
538 // back-substitute
539 for ( i = nrows-1; i >= 0; i-- ) {
540 rowi = extmat[i];
541 for ( j = 0; j < i; j++ ) {
542 rowj = extmat[j];
543 rowpivot = rowj[i];
544 if ( rowpivot == 0 ) continue;
545 for ( k = i; k < ncols; k++ )
546 rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
547 // for (k=nrows; k<ncols; k++) rowj[k] = ff_sub(rowj[k], ff_mul(rowpivot, rowi[k]));
548 }
549 }
550 DEBOUTLN( cerr, "solve succeeded" );
551 DEBDECLEVEL( cerr, "solve" );
552 return true;
553}
554
555int
556determinant ( int **extmat, int n )
557{
558 int i, j, k;
559 int divisor, multiplier, rowii, rowji; // all FF
560 int * rowi; // FF
561 int * rowj; // FF
562 int * swap; // FF
563 // triangularization
564 multiplier = 1;
565 divisor = 1;
566
567 for ( i = 0; i < n; i++ ) {
568 //find "pivot"
569 for (j = i; j < n; j++ )
570 if ( extmat[j][i] != 0 ) break;
571 if ( j == n ) return 0;
572 if ( j != i ) {
573 multiplier = ff_neg( multiplier );
574 swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
575 }
576 rowi = extmat[i];
577 rowii = rowi[i];
578 for ( j = i+1; j < n; j++ ) {
579 rowj = extmat[j];
580 rowji = rowj[i];
581 if ( rowji == 0 ) continue;
582 divisor = ff_mul( divisor, rowii );
583 for ( k = i; k < n; k++ )
584 rowj[k] = ff_sub( ff_mul( rowj[k], rowii ), ff_mul( rowi[k], rowji ) );
585 }
586 }
587 multiplier = ff_mul( multiplier, ff_inv( divisor ) );
588 for ( i = 0; i < n; i++ )
589 multiplier = ff_mul( multiplier, extmat[i][i] );
590 return multiplier;
591}
592
593void
594solveVandermondeT ( const CFArray & a, const CFArray & w, CFArray & x, const Variable & z )
595{
596 CanonicalForm Q = 1, q, p;
598 int i, n = a.size();
599
600 for ( i = 1; i <= n; i++ )
601 Q *= ( z - a[i] );
602 for ( i = 1; i <= n; i++ ) {
603 q = Q / ( z - a[i] );
604 p = q / q( a[i], z );
605 x[i] = 0;
606 for ( j = p; j.hasTerms(); j++ )
607 x[i] += w[j.exp()+1] * j.coeff();
608 }
609}
#define swap(_i, _j)
Header for factory's main class CanonicalForm.
CanonicalForm mapinto(const CanonicalForm &f)
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
int level(const CanonicalForm &f)
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
declarations of higher level algorithms.
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition: cf_chinese.cc:57
assertions for Factory
#define ASSERT(expression, message)
Definition: cf_assert.h:99
factory switches.
Iterators for CanonicalForm's.
bool matrix_in_Z(const CFMatrix &M, int rows)
Definition: cf_linsys.cc:39
bool betterpivot(const CanonicalForm &oldpivot, const CanonicalForm &newpivot)
Definition: cf_linsys.cc:61
int int ncols
Definition: cf_linsys.cc:32
VAR bool fuzzy_result
Definition: cf_linsys.cc:75
int nrows
Definition: cf_linsys.cc:32
CanonicalForm detbound(const CFMatrix &M, int rows)
Definition: cf_linsys.cc:486
void solveVandermondeT(const CFArray &a, const CFArray &w, CFArray &x, const Variable &z)
Definition: cf_linsys.cc:594
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
bool linearSystemSolve(CFMatrix &M)
Definition: cf_linsys.cc:78
int determinant(int **extmat, int n)
Definition: cf_linsys.cc:556
bool solve(int **extmat, int nrows, int ncols)
Definition: cf_linsys.cc:504
CanonicalForm determinant2(const CFMatrix &M, int rows)
Definition: cf_linsys.cc:323
static bool fill_int_mat(const CFMatrix &M, int **m, int rows)
Definition: cf_linsys.cc:203
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
int cf_getNumPrimes()
Definition: cf_primes.cc:23
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
access to prime tables
int size() const
Definition: ftmpl_array.cc:92
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
long intval() const
conversion functions
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
factory's class for variables
Definition: variable.h:33
functions to print debug output
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
#define DEBOUTENDL(stream)
Definition: debug.h:48
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
#define DEBOUTSL(stream)
Definition: debug.h:46
#define DEBOUT(stream, objects)
Definition: debug.h:47
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
const CanonicalForm & w
Definition: facAbsFact.cc:51
b *CanonicalForm B
Definition: facBivar.cc:52
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
bool isZero(const CFArray &A)
checks if entries of A are zero
#define const
Definition: fegetopt.c:39
operations in a finite prime field F_p.
int ff_inv(const int a)
Definition: ffops.h:149
int ff_mul(const int a, const int b)
Definition: ffops.h:139
int ff_neg(const int a)
Definition: ffops.h:126
int ff_sub(const int a, const int b)
Definition: ffops.h:112
#define VAR
Definition: globaldefs.h:5
bool pivot(const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring R)
This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2].
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:327
static int sign(int x)
Definition: ring.cc:3427
#define M
Definition: sirandom.c:25
#define Q
Definition: sirandom.c:26
int * int_ptr
Definition: structs.h:54
#define TIMING_DEFINE_PRINT(t)
Definition: timing.h:95
#define TIMING_END(t)
Definition: timing.h:93
#define TIMING_START(t)
Definition: timing.h:92