source: git/factory/cf_linsys.cc @ e4fe2b

jengelh-datetimespielwiese
Last change on this file since e4fe2b was e4fe2b, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
FIX: Fixed huge BUG in cf_gmp.h CHG: starting to cleanup factory
  • Property mode set to 100644
File size: 17.3 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include "config.h"
5
6#include "cf_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                {
172                    chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i-1][j-1]), CanonicalForm(p), mnew, qnew );
173                    MM(i, j) = mnew;
174                }
175            Q = qnew;
176        }
177        if ( pno == cf_getNumBigPrimes() )
178            fuzzy_result = true;
179        else
180            fuzzy_result = false;
181        // store the result in M
182        Qhalf = Q / 2;
183        for ( i = 1; i <= rows; i++ ) {
184            for ( j = rows+1; j <= cols; j++ )
185                if ( MM(i,j) > Qhalf )
186                    M(i,j) = MM(i,j) - Q;
187                else
188                    M(i,j) = MM(i,j);
189            delete [] mm[i-1];
190        }
191        delete [] mm;
192        return ! fuzzy_result;
193    }
194}
195
196static bool
197fill_int_mat( const CFMatrix & M, int ** m, int rows )
198{
199    int i, j;
200    bool ok = true;
201    for ( i = 0; i < rows && ok; i++ )
202        for ( j = 0; j < rows && ok; j++ )
203        {
204            if ( M(i+1,j+1).isZero() )
205                m[i][j] = 0;
206            else
207            {
208                m[i][j] = mapinto( M(i+1,j+1) ).intval();
209//                ok = m[i][j] != 0;
210            }
211        }
212    return ok;
213}
214
215CanonicalForm
216determinant( const CFMatrix & M, int rows )
217{
218    typedef int* int_ptr;
219
220    ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
221    if ( rows == 1 )
222        return M(1,1);
223    else  if ( rows == 2 )
224        return M(1,1)*M(2,2)-M(2,1)*M(1,2);
225    else  if ( matrix_in_Z( M, rows ) )
226    {
227        int ** mm = new int_ptr[rows];
228        CanonicalForm x, q, Qhalf, B;
229        int n, i, intdet, p, pno;
230        for ( i = 0; i < rows; i++ )
231        {
232            mm[i] = new int[rows];
233        }
234        pno = 0; n = 0;
235        TIMING_START(det_bound);
236        B = detbound( M, rows );
237        TIMING_END(det_bound);
238        q = 1;
239        TIMING_START(det_numprimes);
240        while ( B > q && n < cf_getNumBigPrimes() )
241        {
242            q *= cf_getBigPrime( n );
243            n++;
244        }
245        TIMING_END(det_numprimes);
246
247        CFArray X(1,n), Q(1,n);
248
249        while ( pno < n )
250        {
251            p = cf_getBigPrime( pno );
252            setCharacteristic( p );
253            // map matrix into char p
254            TIMING_START(det_mapping);
255            fill_int_mat( M, mm, rows );
256            TIMING_END(det_mapping);
257            pno++;
258            DEBOUT( cerr, "." );
259            TIMING_START(det_determinant);
260            intdet = determinant( mm, rows );
261            TIMING_END(det_determinant);
262            setCharacteristic( 0 );
263            X[pno] = intdet;
264            Q[pno] = p;
265        }
266        TIMING_START(det_chinese);
267        chineseRemainder( X, Q, x, q );
268        TIMING_END(det_chinese);
269        Qhalf = q / 2;
270        if ( x > Qhalf )
271            x = x - q;
272        for ( i = 0; i < rows; i++ )
273            delete [] mm[i];
274        delete [] mm;
275        return x;
276    }
277    else
278    {
279        CFMatrix m( M );
280        CanonicalForm divisor = 1, pivot, mji;
281        int i, j, k, sign = 1;
282        for ( i = 1; i <= rows; i++ ) {
283            pivot = m(i,i); k = i;
284            for ( j = i+1; j <= rows; j++ ) {
285                if ( betterpivot( pivot, m(j,i) ) ) {
286                    pivot = m(j,i);
287                    k = j;
288                }
289            }
290            if ( pivot.isZero() )
291                return 0;
292            if ( i != k )
293            {
294                m.swapRow( i, k );
295                sign = -sign;
296            }
297            for ( j = i+1; j <= rows; j++ )
298            {
299                if ( ! m(j,i).isZero() )
300                {
301                    divisor *= pivot;
302                    mji = m(j,i);
303                    m(j,i) = 0;
304                    for ( k = i+1; k <= rows; k++ )
305                        m(j,k) = m(j,k) * pivot - m(i,k)*mji;
306                }
307            }
308        }
309        pivot = sign;
310        for ( i = 1; i <= rows; i++ )
311            pivot *= m(i,i);
312        return pivot / divisor;
313    }
314}
315
316CanonicalForm
317determinant2( const CFMatrix & M, int rows )
318{
319    typedef int* int_ptr;
320
321    ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
322    if ( rows == 1 )
323        return M(1,1);
324    else  if ( rows == 2 )
325        return M(1,1)*M(2,2)-M(2,1)*M(1,2);
326    else  if ( matrix_in_Z( M, rows ) ) {
327        int ** mm = new int_ptr[rows];
328        CanonicalForm QQ, Q, Qhalf, mnew, q, qnew, B;
329        CanonicalForm det, detnew, qdet;
330        int i, p, pcount, pno, intdet;
331        bool ok;
332
333        // initialize room to hold the result and the result mod p
334        for ( i = 0; i < rows; i++ ) {
335            mm[i] = new int[rows];
336        }
337
338        // calculate the bound for the result
339        B = detbound( M, rows );
340
341        // find a first solution mod p
342        pno = 0;
343        do {
344            p = cf_getBigPrime( pno );
345            setCharacteristic( p );
346            // map matrix into char p
347            ok = fill_int_mat( M, mm, rows );
348            pno++;
349        } while ( ! ok && pno < cf_getNumPrimes() );
350        // initialize the result matrix with first solution
351        // solve mod p
352        DEBOUT( cerr, "." );
353        intdet = determinant( mm, rows );
354        setCharacteristic( 0 );
355        det = intdet;
356        // Q so far
357        Q = p;
358        QQ = p;
359        while ( Q < B && cf_getNumPrimes() > pno ) {
360            // find a first solution mod p
361            do {
362                p = cf_getBigPrime( pno );
363                setCharacteristic( p );
364                // map matrix into char p
365                ok = fill_int_mat( M, mm, rows );
366                pno++;
367            } while ( ! ok && pno < cf_getNumPrimes() );
368            // initialize the result matrix with first solution
369            // solve mod p
370            DEBOUT( cerr, "." );
371            intdet = determinant( mm, rows );
372            setCharacteristic( 0 );
373            qdet = intdet;
374            // Q so far
375            q = p;
376            QQ *= p;
377            pcount = 0;
378            while ( QQ < B && cf_getNumPrimes() > pno && pcount < 500 ) {
379                do {
380                    p = cf_getBigPrime( pno );
381                    setCharacteristic( p );
382                    ok = true;
383                    // map matrix into char p
384                    ok = fill_int_mat( M, mm, rows );
385                    pno++;
386                } while ( ! ok && cf_getNumPrimes() > pno );
387                // solve mod p
388                DEBOUT( cerr, "." );
389                intdet = determinant( mm, rows );
390                // found a solution mod p
391                // now chinese remainder it to a solution mod Q*p
392                setCharacteristic( 0 );
393                chineseRemainder( qdet, q, intdet, p, detnew, qnew );
394                qdet = detnew;
395                q = qnew;
396                QQ *= p;
397                pcount++;
398            }
399            DEBOUT( cerr, "*" );
400            chineseRemainder( det, Q, qdet, q, detnew, qnew );
401            Q = qnew;
402            QQ = Q;
403            det = detnew;
404        }
405        if ( ! ok )
406            fuzzy_result = true;
407        else
408            fuzzy_result = false;
409        // store the result in M
410        Qhalf = Q / 2;
411        if ( det > Qhalf )
412            det = det - Q;
413        for ( i = 0; i < rows; i++ )
414            delete [] mm[i];
415        delete [] mm;
416        return det;
417    }
418    else {
419        CFMatrix m( M );
420        CanonicalForm divisor = 1, pivot, mji;
421        int i, j, k, sign = 1;
422        for ( i = 1; i <= rows; i++ ) {
423            pivot = m(i,i); k = i;
424            for ( j = i+1; j <= rows; j++ ) {
425                if ( betterpivot( pivot, m(j,i) ) ) {
426                    pivot = m(j,i);
427                    k = j;
428                }
429            }
430            if ( pivot.isZero() )
431                return 0;
432            if ( i != k ) {
433                m.swapRow( i, k );
434                sign = -sign;
435            }
436            for ( j = i+1; j <= rows; j++ ) {
437                if ( ! m(j,i).isZero() ) {
438                    divisor *= pivot;
439                    mji = m(j,i);
440                    m(j,i) = 0;
441                    for ( k = i+1; k <= rows; k++ )
442                        m(j,k) = m(j,k) * pivot - m(i,k)*mji;
443                }
444            }
445        }
446        pivot = sign;
447        for ( i = 1; i <= rows; i++ )
448            pivot *= m(i,i);
449        return pivot / divisor;
450    }
451}
452
453static CanonicalForm
454bound ( const CFMatrix & M )
455{
456    DEBINCLEVEL( cerr, "bound" );
457    int rows = M.rows(), cols = M.columns();
458    CanonicalForm sum = 0;
459    int i, j;
460    for ( i = 1; i <= rows; i++ )
461        for ( j = 1; j <= rows; j++ )
462            sum += M(i,j) * M(i,j);
463    DEBOUTLN( cerr, "bound(matrix)^2 = " << sum );
464    CanonicalForm vmax = 0, vsum;
465    for ( j = rows+1; j <= cols; j++ ) {
466        vsum = 0;
467        for ( i = 1; i <= rows; i++ )
468            vsum += M(i,j) * M(i,j);
469        if ( vsum > vmax ) vmax = vsum;
470    }
471    DEBOUTLN( cerr, "bound(lhs)^2 = " << vmax );
472    sum += vmax;
473    DEBOUTLN( cerr, "bound(overall)^2 = " << sum );
474    DEBDECLEVEL( cerr, "bound" );
475    return sqrt( sum ) + 1;
476}
477
478
479CanonicalForm
480detbound ( const CFMatrix & M, int rows )
481{
482    CanonicalForm sum = 0, prod = 2;
483    int i, j;
484    for ( i = 1; i <= rows; i++ ) {
485        sum = 0;
486        for ( j = 1; j <= rows; j++ )
487            sum += M(i,j) * M(i,j);
488        prod *= 1 + sqrt(sum);
489    }
490    return prod;
491}
492
493
494// solve returns false if computation failed
495// extmat is overwritten: output is Id mat followed by solution(s)
496
497bool
498solve ( int **extmat, int nrows, int ncols )
499{
500    DEBINCLEVEL( cerr, "solve" );
501    int i, j, k;
502    int rowpivot, pivotrecip; // all FF
503    int * rowi; // FF
504    int * rowj; // FF
505    int * swap; // FF
506    // triangularization
507    for ( i = 0; i < nrows; i++ ) {
508        //find "pivot"
509        for (j = i; j < nrows; j++ )
510            if ( extmat[j][i] != 0 ) break;
511        if ( j == nrows ) {
512            DEBOUTLN( cerr, "solve failed" );
513            DEBDECLEVEL( cerr, "solve" );
514            return false;
515        }
516        if ( j != i ) {
517            swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
518        }
519        pivotrecip = ff_inv( extmat[i][i] );
520        rowi = extmat[i];
521        for ( j = 0; j < ncols; j++ )
522            rowi[j] = ff_mul( pivotrecip, rowi[j] );
523        for ( j = i+1; j < nrows; j++ ) {
524            rowj = extmat[j];
525            rowpivot = rowj[i];
526            if ( rowpivot == 0 ) continue;
527            for ( k = i; k < ncols; k++ )
528                rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
529        }
530    }
531    // matrix is now upper triangular with 1s down the diagonal
532    // back-substitute
533    for ( i = nrows-1; i >= 0; i-- ) {
534        rowi = extmat[i];
535        for ( j = 0; j < i; j++ ) {
536            rowj = extmat[j];
537            rowpivot = rowj[i];
538            if ( rowpivot == 0 ) continue;
539            for ( k = i; k < ncols; k++ )
540                rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
541            // for (k=nrows; k<ncols; k++) rowj[k] = ff_sub(rowj[k], ff_mul(rowpivot, rowi[k]));
542        }
543    }
544    DEBOUTLN( cerr, "solve succeeded" );
545    DEBDECLEVEL( cerr, "solve" );
546    return true;
547}
548
549int
550determinant ( int **extmat, int n )
551{
552    int i, j, k;
553    int divisor, multiplier, rowii, rowji; // all FF
554    int * rowi; // FF
555    int * rowj; // FF
556    int * swap; // FF
557    // triangularization
558    multiplier = 1;
559    divisor = 1;
560
561    for ( i = 0; i < n; i++ ) {
562        //find "pivot"
563        for (j = i; j < n; j++ )
564            if ( extmat[j][i] != 0 ) break;
565        if ( j == n ) return 0;
566        if ( j != i ) {
567            multiplier = ff_neg( multiplier );
568            swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
569        }
570        rowi = extmat[i];
571        rowii = rowi[i];
572        for ( j = i+1; j < n; j++ ) {
573            rowj = extmat[j];
574            rowji = rowj[i];
575            if ( rowji == 0 ) continue;
576            divisor = ff_mul( divisor, rowii );
577            for ( k = i; k < n; k++ )
578                rowj[k] = ff_sub( ff_mul( rowj[k], rowii ), ff_mul( rowi[k], rowji ) );
579        }
580    }
581    multiplier = ff_mul( multiplier, ff_inv( divisor ) );
582    for ( i = 0; i < n; i++ )
583        multiplier = ff_mul( multiplier, extmat[i][i] );
584    return multiplier;
585}
586
587void
588solveVandermondeT ( const CFArray & a, const CFArray & w, CFArray & x, const Variable & z )
589{
590    CanonicalForm Q = 1, q, p;
591    CFIterator j;
592    int i, n = a.size();
593
594    for ( i = 1; i <= n; i++ )
595        Q *= ( z - a[i] );
596    for ( i = 1; i <= n; i++ ) {
597        q = Q / ( z - a[i] );
598        p = q / q( a[i], z );
599        x[i] = 0;
600        for ( j = p; j.hasTerms(); j++ )
601            x[i] += w[j.exp()+1] * j.coeff();
602    }
603}
Note: See TracBrowser for help on using the repository browser.