source: git/factory/cf_linsys.cc @ fbefc9

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