source: git/factory/cf_linsys.cc @ cc1367

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