source: git/factory/cf_linsys.cc @ 362fc67

jengelh-datetimespielwiese
Last change on this file since 362fc67 was 362fc67, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: remove $Id$
  • Property mode set to 100644
File size: 17.3 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3#include "config.h"
4
5#include "cf_assert.h"
6#include "debug.h"
7#include "timing.h"
8
9#include "cf_defs.h"
10#include "cf_primes.h"
11#include "canonicalform.h"
12#include "cf_iter.h"
13#include "cf_algorithm.h"
14#include "ffops.h"
15#include "cf_primes.h"
16
17
18TIMING_DEFINE_PRINT(det_mapping)
19TIMING_DEFINE_PRINT(det_determinant)
20TIMING_DEFINE_PRINT(det_chinese)
21TIMING_DEFINE_PRINT(det_bound)
22TIMING_DEFINE_PRINT(det_numprimes)
23
24
25static bool solve ( int **extmat, int nrows, int ncols );
26int determinant ( int **extmat, int n );
27
28static CanonicalForm bound ( const CFMatrix & M );
29CanonicalForm detbound ( const CFMatrix & M, int rows );
30
31bool
32matrix_in_Z( const CFMatrix & M, int rows )
33{
34    int i, j;
35    for ( i = 1; i <= rows; i++ )
36        for ( j = 1; j <= rows; j++ )
37            if ( ! M(i,j).inZ() )
38                return false;
39    return true;
40}
41
42bool
43matrix_in_Z( const CFMatrix & M )
44{
45    int i, j, rows = M.rows(), cols = M.columns();
46    for ( i = 1; i <= rows; i++ )
47        for ( j = 1; j <= cols; j++ )
48            if ( ! M(i,j).inZ() )
49                return false;
50    return true;
51}
52
53bool
54betterpivot ( const CanonicalForm & oldpivot, const CanonicalForm & newpivot )
55{
56    if ( newpivot.isZero() )
57        return false;
58    else  if ( oldpivot.isZero() )
59        return true;
60    else  if ( level( oldpivot ) > level( newpivot ) )
61        return true;
62    else  if ( level( oldpivot ) < level( newpivot ) )
63        return false;
64    else
65        return ( newpivot.lc() < oldpivot.lc() );
66}
67
68bool fuzzy_result;
69
70bool
71linearSystemSolve( CFMatrix & M )
72{
73    typedef int* int_ptr;
74
75    if ( ! matrix_in_Z( M ) ) {
76        int nrows = M.rows(), ncols = M.columns();
77        int i, j, k;
78        CanonicalForm rowpivot, pivotrecip;
79        // triangularization
80        for ( i = 1; i <= nrows; i++ ) {
81            //find "pivot"
82            for (j = i; j <= nrows; j++ )
83                if ( M(j,i) != 0 ) break;
84            if ( j > nrows ) return false;
85            if ( j != i )
86                M.swapRow( i, j );
87            pivotrecip = 1 / M(i,i);
88            for ( j = 1; j <= ncols; j++ )
89                M(i,j) *= pivotrecip;
90            for ( j = i+1; j <= nrows; j++ ) {
91                rowpivot = M(j,i);
92                if ( rowpivot == 0 ) continue;
93                for ( k = i; k <= ncols; k++ )
94                    M(j,k) -= M(i,k) * rowpivot;
95            }
96        }
97        // matrix is now upper triangular with 1s down the diagonal
98        // back-substitute
99        for ( i = nrows-1; i > 0; i-- ) {
100            for ( j = nrows+1; j <= ncols; j++ ) {
101                for ( k = i+1; k <= nrows; k++ )
102                    M(i,j) -= M(k,j) * M(i,k);
103            }
104        }
105        return true;
106    }
107    else {
108        int rows = M.rows(), cols = M.columns();
109        CFMatrix MM( rows, cols );
110        int ** mm = new int_ptr[rows];
111        CanonicalForm Q, Qhalf, mnew, qnew, B;
112        int i, j, p, pno;
113        bool ok;
114
115        // initialize room to hold the result and the result mod p
116        for ( i = 0; i < rows; i++ ) {
117            mm[i] = new int[cols];
118        }
119
120        // calculate the bound for the result
121        B = bound( M );
122        DEBOUTLN( cerr, "bound = " <<  B );
123
124        // find a first solution mod p
125        pno = 0;
126        do {
127            DEBOUTSL( cerr );
128            DEBOUT( cerr, "trying prime(" << pno << ") = " );
129            p = cf_getBigPrime( pno );
130            DEBOUT( cerr, p );
131            DEBOUTENDL( cerr );
132            setCharacteristic( p );
133            // map matrix into char p
134            for ( i = 1; i <= rows; i++ )
135                for ( j = 1; j <= cols; j++ )
136                    mm[i-1][j-1] = mapinto( M(i,j) ).intval();
137            // solve mod p
138            ok = solve( mm, rows, cols );
139            pno++;
140        } while ( ! ok );
141
142        // initialize the result matrix with first solution
143        setCharacteristic( 0 );
144        for ( i = 1; i <= rows; i++ )
145            for ( j = rows+1; j <= cols; j++ )
146                MM(i,j) = mm[i-1][j-1];
147
148        // Q so far
149        Q = p;
150        while ( Q < B && pno < cf_getNumBigPrimes() ) {
151            do {
152                DEBOUTSL( cerr );
153                DEBOUT( cerr, "trying prime(" << pno << ") = " );
154                p = cf_getBigPrime( pno );
155                DEBOUT( cerr, p );
156                DEBOUTENDL( cerr );
157                setCharacteristic( p );
158                for ( i = 1; i <= rows; i++ )
159                    for ( j = 1; j <= cols; j++ )
160                        mm[i-1][j-1] = mapinto( M(i,j) ).intval();
161                // solve mod p
162                ok = solve( mm, rows, cols );
163                pno++;
164            } while ( ! ok );
165            // found a solution mod p
166            // now chinese remainder it to a solution mod Q*p
167            setCharacteristic( 0 );
168            for ( i = 1; i <= rows; i++ )
169                for ( j = rows+1; j <= cols; j++ )
170                {
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        {
203            if ( M(i+1,j+1).isZero() )
204                m[i][j] = 0;
205            else
206            {
207                m[i][j] = mapinto( M(i+1,j+1) ).intval();
208//                ok = m[i][j] != 0;
209            }
210        }
211    return ok;
212}
213
214CanonicalForm
215determinant( const CFMatrix & M, int rows )
216{
217    typedef int* int_ptr;
218
219    ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
220    if ( rows == 1 )
221        return M(1,1);
222    else  if ( rows == 2 )
223        return M(1,1)*M(2,2)-M(2,1)*M(1,2);
224    else  if ( matrix_in_Z( M, rows ) )
225    {
226        int ** mm = new int_ptr[rows];
227        CanonicalForm x, q, Qhalf, B;
228        int n, i, intdet, p, pno;
229        for ( i = 0; i < rows; i++ )
230        {
231            mm[i] = new int[rows];
232        }
233        pno = 0; n = 0;
234        TIMING_START(det_bound);
235        B = detbound( M, rows );
236        TIMING_END(det_bound);
237        q = 1;
238        TIMING_START(det_numprimes);
239        while ( B > q && n < cf_getNumBigPrimes() )
240        {
241            q *= cf_getBigPrime( n );
242            n++;
243        }
244        TIMING_END(det_numprimes);
245
246        CFArray X(1,n), Q(1,n);
247
248        while ( pno < n )
249        {
250            p = cf_getBigPrime( pno );
251            setCharacteristic( p );
252            // map matrix into char p
253            TIMING_START(det_mapping);
254            fill_int_mat( M, mm, rows );
255            TIMING_END(det_mapping);
256            pno++;
257            DEBOUT( cerr, "." );
258            TIMING_START(det_determinant);
259            intdet = determinant( mm, rows );
260            TIMING_END(det_determinant);
261            setCharacteristic( 0 );
262            X[pno] = intdet;
263            Q[pno] = p;
264        }
265        TIMING_START(det_chinese);
266        chineseRemainder( X, Q, x, q );
267        TIMING_END(det_chinese);
268        Qhalf = q / 2;
269        if ( x > Qhalf )
270            x = x - q;
271        for ( i = 0; i < rows; i++ )
272            delete [] mm[i];
273        delete [] mm;
274        return x;
275    }
276    else
277    {
278        CFMatrix m( M );
279        CanonicalForm divisor = 1, pivot, mji;
280        int i, j, k, sign = 1;
281        for ( i = 1; i <= rows; i++ ) {
282            pivot = m(i,i); k = i;
283            for ( j = i+1; j <= rows; j++ ) {
284                if ( betterpivot( pivot, m(j,i) ) ) {
285                    pivot = m(j,i);
286                    k = j;
287                }
288            }
289            if ( pivot.isZero() )
290                return 0;
291            if ( i != k )
292            {
293                m.swapRow( i, k );
294                sign = -sign;
295            }
296            for ( j = i+1; j <= rows; j++ )
297            {
298                if ( ! m(j,i).isZero() )
299                {
300                    divisor *= pivot;
301                    mji = m(j,i);
302                    m(j,i) = 0;
303                    for ( k = i+1; k <= rows; k++ )
304                        m(j,k) = m(j,k) * pivot - m(i,k)*mji;
305                }
306            }
307        }
308        pivot = sign;
309        for ( i = 1; i <= rows; i++ )
310            pivot *= m(i,i);
311        return pivot / divisor;
312    }
313}
314
315CanonicalForm
316determinant2( const CFMatrix & M, int rows )
317{
318    typedef int* int_ptr;
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_ptr[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    DEBINCLEVEL( cerr, "bound" );
456    int rows = M.rows(), cols = M.columns();
457    CanonicalForm sum = 0;
458    int i, j;
459    for ( i = 1; i <= rows; i++ )
460        for ( j = 1; j <= rows; j++ )
461            sum += M(i,j) * M(i,j);
462    DEBOUTLN( cerr, "bound(matrix)^2 = " << sum );
463    CanonicalForm vmax = 0, vsum;
464    for ( j = rows+1; j <= cols; j++ ) {
465        vsum = 0;
466        for ( i = 1; i <= rows; i++ )
467            vsum += M(i,j) * M(i,j);
468        if ( vsum > vmax ) vmax = vsum;
469    }
470    DEBOUTLN( cerr, "bound(lhs)^2 = " << vmax );
471    sum += vmax;
472    DEBOUTLN( cerr, "bound(overall)^2 = " << sum );
473    DEBDECLEVEL( cerr, "bound" );
474    return sqrt( sum ) + 1;
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    DEBINCLEVEL( cerr, "solve" );
500    int i, j, k;
501    int rowpivot, pivotrecip; // all FF
502    int * rowi; // FF
503    int * rowj; // FF
504    int * swap; // FF
505    // triangularization
506    for ( i = 0; i < nrows; i++ ) {
507        //find "pivot"
508        for (j = i; j < nrows; j++ )
509            if ( extmat[j][i] != 0 ) break;
510        if ( j == nrows ) {
511            DEBOUTLN( cerr, "solve failed" );
512            DEBDECLEVEL( cerr, "solve" );
513            return false;
514        }
515        if ( j != i ) {
516            swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
517        }
518        pivotrecip = ff_inv( extmat[i][i] );
519        rowi = extmat[i];
520        for ( j = 0; j < ncols; j++ )
521            rowi[j] = ff_mul( pivotrecip, rowi[j] );
522        for ( j = i+1; j < nrows; j++ ) {
523            rowj = extmat[j];
524            rowpivot = rowj[i];
525            if ( rowpivot == 0 ) continue;
526            for ( k = i; k < ncols; k++ )
527                rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
528        }
529    }
530    // matrix is now upper triangular with 1s down the diagonal
531    // back-substitute
532    for ( i = nrows-1; i >= 0; i-- ) {
533        rowi = extmat[i];
534        for ( j = 0; j < i; j++ ) {
535            rowj = extmat[j];
536            rowpivot = rowj[i];
537            if ( rowpivot == 0 ) continue;
538            for ( k = i; k < ncols; k++ )
539                rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
540            // for (k=nrows; k<ncols; k++) rowj[k] = ff_sub(rowj[k], ff_mul(rowpivot, rowi[k]));
541        }
542    }
543    DEBOUTLN( cerr, "solve succeeded" );
544    DEBDECLEVEL( cerr, "solve" );
545    return true;
546}
547
548int
549determinant ( int **extmat, int n )
550{
551    int i, j, k;
552    int divisor, multiplier, rowii, rowji; // all FF
553    int * rowi; // FF
554    int * rowj; // FF
555    int * swap; // FF
556    // triangularization
557    multiplier = 1;
558    divisor = 1;
559
560    for ( i = 0; i < n; i++ ) {
561        //find "pivot"
562        for (j = i; j < n; j++ )
563            if ( extmat[j][i] != 0 ) break;
564        if ( j == n ) return 0;
565        if ( j != i ) {
566            multiplier = ff_neg( multiplier );
567            swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
568        }
569        rowi = extmat[i];
570        rowii = rowi[i];
571        for ( j = i+1; j < n; j++ ) {
572            rowj = extmat[j];
573            rowji = rowj[i];
574            if ( rowji == 0 ) continue;
575            divisor = ff_mul( divisor, rowii );
576            for ( k = i; k < n; k++ )
577                rowj[k] = ff_sub( ff_mul( rowj[k], rowii ), ff_mul( rowi[k], rowji ) );
578        }
579    }
580    multiplier = ff_mul( multiplier, ff_inv( divisor ) );
581    for ( i = 0; i < n; i++ )
582        multiplier = ff_mul( multiplier, extmat[i][i] );
583    return multiplier;
584}
585
586void
587solveVandermondeT ( const CFArray & a, const CFArray & w, CFArray & x, const Variable & z )
588{
589    CanonicalForm Q = 1, q, p;
590    CFIterator j;
591    int i, n = a.size();
592
593    for ( i = 1; i <= n; i++ )
594        Q *= ( z - a[i] );
595    for ( i = 1; i <= n; i++ ) {
596        q = Q / ( z - a[i] );
597        p = q / q( a[i], z );
598        x[i] = 0;
599        for ( j = p; j.hasTerms(); j++ )
600            x[i] += w[j.exp()+1] * j.coeff();
601    }
602}
Note: See TracBrowser for help on using the repository browser.