source: git/factory/cf_linsys.cc @ f5d2647

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