source: git/factory/cf_linsys.cc

spielwiese
Last change on this file was a3f0fea, checked in by Reimer Behrends <behrends@…>, 5 years ago
Modify variable declarions for pSingular.
  • Property mode set to 100644
File size: 17.4 KB
Line 
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
25TIMING_DEFINE_PRINT(det_mapping)
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
35static CanonicalForm bound ( const CFMatrix & M );
36CanonicalForm detbound ( const CFMatrix & M, int rows );
37
38bool
39matrix_in_Z( const CFMatrix & M, int rows )
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
50matrix_in_Z( const CFMatrix & M )
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
75VAR bool fuzzy_result;
76
77bool
78linearSystemSolve( CFMatrix & M )
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 );
139            setCharacteristic( p );
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
150        setCharacteristic( 0 );
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 );
164                setCharacteristic( p );
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
174            setCharacteristic( 0 );
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
221CanonicalForm
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 );
258            setCharacteristic( p );
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);
268            setCharacteristic( 0 );
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
322CanonicalForm
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 );
351            setCharacteristic( p );
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 );
360        setCharacteristic( 0 );
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 );
369                setCharacteristic( p );
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 );
378            setCharacteristic( 0 );
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 );
387                    setCharacteristic( p );
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
398                setCharacteristic( 0 );
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
485CanonicalForm
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;
597    CFIterator j;
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}
Note: See TracBrowser for help on using the repository browser.