source: git/factory/cf_linsys.cc @ b973c0

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