42 for ( i = 1; i <= rows; i++ )
43 for ( j = 1; j <= rows; j++ )
53 for ( i = 1; i <= rows; i++ )
54 for ( j = 1; j <= cols; j++ )
65 else if ( oldpivot.
isZero() )
67 else if (
level( oldpivot ) >
level( newpivot ) )
69 else if (
level( oldpivot ) <
level( newpivot ) )
72 return ( newpivot.
lc() < oldpivot.
lc() );
87 for ( i = 1; i <=
nrows; i++ ) {
89 for (j = i; j <=
nrows; j++ )
90 if (
M(j,i) != 0 )
break;
91 if ( j > nrows )
return false;
94 pivotrecip = 1 /
M(i,i);
95 for ( j = 1; j <=
ncols; j++ )
97 for ( j = i+1; j <=
nrows; j++ ) {
99 if ( rowpivot == 0 )
continue;
100 for ( k = i; k <=
ncols; k++ )
101 M(j,k) -=
M(i,k) * rowpivot;
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);
117 int ** mm =
new int_ptr[rows];
123 for ( i = 0; i < rows; i++ ) {
124 mm[
i] =
new int[cols];
135 DEBOUT( cerr,
"trying prime(" << pno <<
") = " );
141 for ( i = 1; i <= rows; i++ )
142 for ( j = 1; j <= cols; j++ )
145 ok =
solve( mm, rows, cols );
151 for ( i = 1; i <= rows; i++ )
152 for ( j = rows+1; j <= cols; j++ )
153 MM(i,j) = mm[i-1][j-1];
160 DEBOUT( cerr,
"trying prime(" << pno <<
") = " );
165 for ( i = 1; i <= rows; i++ )
166 for ( j = 1; j <= cols; j++ )
169 ok =
solve( mm, rows, cols );
175 for ( i = 1; i <= rows; i++ )
176 for ( j = rows+1; j <= cols; j++ )
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;
207 for ( i = 0; i < rows && ok; i++ )
208 for ( j = 0; j < rows && ok; j++ )
226 ASSERT( rows <= M.
rows() && rows <= M.
columns() && rows > 0,
"undefined determinant" );
229 else if ( rows == 2 )
230 return M(1,1)*
M(2,2)-
M(2,1)*
M(1,2);
233 int ** mm =
new int_ptr[rows];
235 int n,
i, intdet,
p, pno;
236 for ( i = 0; i < rows; i++ )
238 mm[
i] =
new int[rows];
278 for ( i = 0; i < rows; i++ )
288 for ( i = 1; i <= rows; i++ ) {
290 for ( j = i+1; j <= rows; j++ ) {
296 if (
pivot.isZero() )
303 for ( j = i+1; j <= rows; j++ )
310 for ( k = i+1; k <= rows; k++ )
311 m(j,k) =
m(j,k) *
pivot -
m(i,k)*mji;
316 for ( i = 1; i <= rows; i++ )
318 return pivot / divisor;
327 ASSERT( rows <= M.
rows() && rows <= M.
columns() && rows > 0,
"undefined determinant" );
330 else if ( rows == 2 )
331 return M(1,1)*
M(2,2)-
M(2,1)*
M(1,2);
333 int ** mm =
new int_ptr[rows];
336 int i,
p, pcount, pno, intdet;
340 for ( i = 0; i < rows; i++ ) {
341 mm[
i] =
new int[rows];
419 for ( i = 0; i < rows; i++ )
428 for ( i = 1; i <= rows; i++ ) {
430 for ( j = i+1; j <= rows; j++ ) {
436 if (
pivot.isZero() )
442 for ( j = i+1; j <= rows; j++ ) {
447 for ( k = i+1; k <= rows; k++ )
448 m(j,k) =
m(j,k) *
pivot -
m(i,k)*mji;
453 for ( i = 1; i <= rows; i++ )
455 return pivot / divisor;
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 );
471 for ( j = rows+1; j <= cols; j++ ) {
473 for ( i = 1; i <= rows; i++ )
474 vsum +=
M(i,j) *
M(i,j);
475 if ( vsum > vmax ) vmax = vsum;
477 DEBOUTLN( cerr,
"bound(lhs)^2 = " << vmax );
479 DEBOUTLN( cerr,
"bound(overall)^2 = " << sum );
481 return sqrt( sum ) + 1;
490 for ( i = 1; i <= rows; i++ ) {
492 for ( j = 1; j <= rows; j++ )
493 sum +=
M(i,j) *
M(i,j);
508 int rowpivot, pivotrecip;
513 for ( i = 0; i <
nrows; i++ ) {
515 for (j = i; j <
nrows; j++ )
516 if ( extmat[j][i] != 0 )
break;
523 swap = extmat[
i]; extmat[
i] = extmat[
j]; extmat[
j] =
swap;
525 pivotrecip =
ff_inv( extmat[i][i] );
527 for ( j = 0; j <
ncols; j++ )
528 rowi[j] =
ff_mul( pivotrecip, rowi[j] );
529 for ( j = i+1; j <
nrows; j++ ) {
532 if ( rowpivot == 0 )
continue;
533 for ( k = i; k <
ncols; k++ )
534 rowj[k] =
ff_sub( rowj[k],
ff_mul( rowpivot, rowi[k] ) );
539 for ( i = nrows-1; i >= 0; i-- ) {
541 for ( j = 0; j <
i; j++ ) {
544 if ( rowpivot == 0 )
continue;
545 for ( k = i; k <
ncols; k++ )
546 rowj[k] =
ff_sub( rowj[k],
ff_mul( rowpivot, rowi[k] ) );
550 DEBOUTLN( cerr,
"solve succeeded" );
559 int divisor, multiplier, rowii, rowji;
567 for ( i = 0; i < n; i++ ) {
569 for (j = i; j < n; j++ )
570 if ( extmat[j][i] != 0 )
break;
571 if ( j == n )
return 0;
573 multiplier =
ff_neg( multiplier );
574 swap = extmat[
i]; extmat[
i] = extmat[
j]; extmat[
j] =
swap;
578 for ( j = i+1; j < n; j++ ) {
581 if ( rowji == 0 )
continue;
582 divisor =
ff_mul( divisor, rowii );
583 for ( k = i; k < n; k++ )
588 for ( i = 0; i < n; i++ )
589 multiplier =
ff_mul( multiplier, extmat[i][i] );
600 for ( i = 1; i <= n; i++ )
602 for ( i = 1; i <= n; i++ ) {
603 q = Q / ( z - a[
i] );
604 p = q / q( a[i], z );
TIMING_DEFINE_PRINT(det_mapping) TIMING_DEFINE_PRINT(det_determinant) TIMING_DEFINE_PRINT(det_chinese) TIMING_DEFINE_PRINT(det_bound) TIMING_DEFINE_PRINT(det_numprimes) static bool solve(int **extmat
static CanonicalForm bound(const CFMatrix &M)
functions to print debug output
CanonicalForm detbound(const CFMatrix &M, int rows)
bool betterpivot(const CanonicalForm &oldpivot, const CanonicalForm &newpivot)
#define DEBINCLEVEL(stream, msg)
bool linearSystemSolve(CFMatrix &M)
factory's class for variables
CF_NO_INLINE CanonicalForm coeff() const
get the current coefficient
#define DEBOUT(stream, objects)
void solveVandermondeT(const CFArray &a, const CFArray &w, CFArray &x, const Variable &z)
int ff_mul(const int a, const int b)
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].
bool solve(int **extmat, int nrows, int ncols)
#define DEBDECLEVEL(stream, msg)
int determinant(int **extmat, int n)
CanonicalForm determinant2(const CFMatrix &M, int rows)
#define DEBOUTENDL(stream)
static bool fill_int_mat(const CFMatrix &M, int **m, int rows)
int ff_sub(const int a, const int b)
gmp_float sqrt(const gmp_float &a)
Iterators for CanonicalForm's.
void swapRow(int i, int j)
declarations of higher level algorithms.
class to iterate through CanonicalForm's
void 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...
CF_NO_INLINE int hasTerms() const
check if iterator has reached < the end of CanonicalForm
int cf_getBigPrime(int i)
bool isZero(const CFArray &A)
checks if entries of A are zero
#define ASSERT(expression, message)
#define DEBOUTLN(stream, objects)
CF_NO_INLINE int exp() const
get the current exponent
operations in a finite prime field F_p.
bool matrix_in_Z(const CFMatrix &M, int rows)