source: git/factory/fac_berlekamp.cc @ 346edc8

jengelh-datetimespielwiese
Last change on this file since 346edc8 was 346edc8, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: need std::
  • Property mode set to 100644
File size: 10.3 KB
RevLine 
[493c477]1/* emacs edit mode for this file is -*- C++ -*- */
[341696]2/* $Id$ */
[b973c0]3
[e4fe2b]4#include "config.h"
[b973c0]5
[650f2d8]6#include "cf_assert.h"
[38bb4c]7#include "debug.h"
8
[2dd068]9#include "cf_defs.h"
10#include "fac_berlekamp.h"
11#include "ffops.h"
12#include "gfops.h"
13#include "imm.h"
14#include "canonicalform.h"
15#include "cf_iter.h"
16#include "cf_generator.h"
17#include "fac_sqrfree.h"
18
[38bb4c]19#ifdef DEBUGOUTPUT
20void QprintFF( int ** Q, int n )
21{
22    for ( int i = 0; i < n; i++ ) {
[806c18]23        for ( int j = 0; j < n; j++ )
[346edc8]24            std::cerr << Q[i][j] << "  ";
25        std::cerr << std::endl;
[38bb4c]26    }
[346edc8]27    std::cerr << std::endl;
[38bb4c]28}
29#endif /* DEBUGOUTPUT */
30
31#ifdef DEBUGOUTPUT
32void QprintGF( int ** Q, int n )
33{
34    for ( int i = 0; i < n; i++ ) {
[806c18]35        for ( int j = 0; j < n; j++ ) {
[346edc8]36            gf_print( std::cerr, Q[i][j] );
37            std::cerr << "  ";
[806c18]38        }
[346edc8]39        std::cerr << std::endl;
[38bb4c]40    }
[346edc8]41    std::cerr << std::endl;
[38bb4c]42}
43#endif /* DEBUGOUTPUT */
44
[2dd068]45void QmatFF ( const CanonicalForm & f, int ** Q, int p )
46{
47    int n = degree( f ), nn = (n-1)*p + 1;
48    int i, m, rn;
49    int * a = new int [n];
50    int * r = new int [n];
51    int * q;
52
53    q = Q[0]; *q = r[0] = 1; a[0] = 0; q++;
54
55    for ( i = 1; i < n; i++, q++ )
[806c18]56        *q = r[i] = a[i] = 0;
[2dd068]57    CFIterator I = f; I++;
58    while ( I.hasTerms() ) {
[806c18]59        a[I.exp()] = I.coeff().intval();
60        I++;
[2dd068]61    }
62    for ( m = 1; m < nn; m++ ) {
[806c18]63        rn = r[n-1];
64        for ( i = n-1; i > 0; i-- )
65            r[i] = ff_sub( r[i-1], ff_mul( rn, a[i] ) );
66        r[0] = ff_mul( ff_neg( rn ), a[0] );
67        if ( m % p == 0 ) {
68            q = Q[m/p];
69            for ( i = 0; i < n; i++, q++ )
70                *q = r[i];
71        }
[2dd068]72    }
73    for ( i = 0; i < n; i++ )
[806c18]74        Q[i][i] = ff_sub( Q[i][i], 1 );
[2dd068]75
76    delete [] a;
77    delete [] r;
78}
79
80void QmatGF ( const CanonicalForm & f, int ** Q, int p )
81{
82    int n = degree( f ), nn = (n-1)*p + 1;
83    int i, m, rn;
84    int * a = new int [n];
85    int * r = new int [n];
86    int * q;
87
88    q = Q[0]; *q = r[0] = gf_one(); a[0] = gf_zero(); q++;
89
90    for ( i = 1; i < n; i++, q++ )
[806c18]91        *q = r[i] = a[i] = gf_zero();
[2dd068]92    CFIterator I = f; I++;
93    while ( I.hasTerms() ) {
[806c18]94        a[I.exp()] = imm2int( I.coeff().getval() );
95        I++;
[2dd068]96    }
97    for ( m = 1; m < nn; m++ ) {
[806c18]98        rn = r[n-1];
99        for ( i = n-1; i > 0; i-- )
100            r[i] = gf_sub( r[i-1], gf_mul( rn, a[i] ) );
101        r[0] = gf_mul( gf_neg( rn ), a[0] );
102        if ( m % p == 0 ) {
103            q = Q[m/p];
104            for ( i = 0; i < n; i++, q++ )
105                *q = r[i];
106        }
[2dd068]107    }
108    for ( i = 0; i < n; i++ )
[806c18]109        Q[i][i] = gf_sub( Q[i][i], gf_one() );
[2dd068]110
111    delete [] a;
112    delete [] r;
113}
114
115int nullSpaceFF ( int ** Q, int ** b, int n )
116{
117    int * c = new int[n];
118    int r, i, j, k, h, s, d;
119
120    r = 0;
121    for ( s = 0; s < n; s++ ) c[s] = -1;
122    for ( h = 0; h < n; h++ ) {
[806c18]123        j = 0;
124        while ( j < n && ! ( Q[h][j] != 0 && c[j] < 0 ) ) j++;
125        if ( j < n ) {
126            d = ff_neg( ff_inv( Q[h][j] ) );
127            for ( s = 0; s < n; s++ )
128                Q[s][j] = ff_mul( d, Q[s][j] );
129            for ( i = 0; i < n; i++ ) {
130                if ( i != j ) {
131                    d = Q[h][i];
132                    for ( s = 0; s < n; s++ )
133                        Q[s][i] = ff_add( ff_mul( d, Q[s][j] ), Q[s][i] );
134                }
135            }
136            c[j] = h;
137        }
138        else {
139            b[r] = new int[n];
140            for ( j = 0; j < n; j++ ) {
141                if ( j == h )
142                    b[r][j] = 1;
143                else {
144                    k = 0;
145                    while ( k < n && c[k] != j ) k++;
146                    if ( k < n )
147                        b[r][j] = Q[h][k];
148                    else
149                        b[r][j] = 0;
150                }
151            }
152            r++;
153        }
[2dd068]154    }
155    delete [] c;
156    return r;
157}
158
159int nullSpaceGF ( int ** Q, int ** b, int n )
160{
161    int * c = new int[n];
162    int r, i, j, k, h, s, d;
163
164    r = 0;
165    for ( s = 0; s < n; s++ ) c[s] = -1;
166    for ( h = 0; h < n; h++ ) {
[806c18]167        j = 0;
168        while ( j < n && ! ( ! gf_iszero( Q[h][j] ) && c[j] < 0 ) ) j++;
169        if ( j < n ) {
170            d = gf_neg( gf_inv( Q[h][j] ) );
171            for ( s = 0; s < n; s++ )
172                Q[s][j] = gf_mul( d, Q[s][j] );
173            for ( i = 0; i < n; i++ ) {
174                if ( i != j ) {
175                    d = Q[h][i];
176                    for ( s = 0; s < n; s++ )
177                        Q[s][i] = gf_add( gf_mul( d, Q[s][j] ), Q[s][i] );
178                }
179            }
180            c[j] = h;
181        }
182        else {
183            b[r] = new int[n];
184            for ( j = 0; j < n; j++ ) {
185                if ( j == h )
186                    b[r][j] = gf_one();
187                else {
188                    k = 0;
189                    while ( k < n && c[k] != j ) k++;
190                    if ( k < n )
191                        b[r][j] = Q[h][k];
192                    else
193                        b[r][j] = gf_zero();
194                }
195            }
196            r++;
197        }
[2dd068]198    }
199    delete [] c;
200    return r;
201}
202
203CanonicalForm cfFromIntVec( int * a, int n, const Variable & x )
204{
205    CanonicalForm result = power( x, n-1 ) * a[n-1];
206    for ( int i = n-2; i >= 0; i-- )
[806c18]207        if ( a[i] != 0 )
208            result += power( x, i ) * a[i];
[2dd068]209    return result;
210}
211
212CanonicalForm cfFromGFVec( int * a, int n, const Variable & x )
213{
214    CanonicalForm result = power( x, n-1 ) * CanonicalForm( int2imm_gf( a[n-1] ) );
215    for ( int i = n-2; i >= 0; i-- )
[806c18]216        if ( ! gf_iszero( a[i] ) )
217            result += power( x, i ) * CanonicalForm( int2imm_gf( a[i] ) );
[2dd068]218    return result;
219}
220
221typedef int * intptr;
222
223CFFList BerlekampFactorFF ( const CanonicalForm & f )
224{
225    CFFList F;
226    int p = getCharacteristic();
227    int r, s, len, i, k, n = degree( f );
228    Variable x = f.mvar();
229    CanonicalForm u, g;
230    intptr* Q = new intptr [n];
231    intptr* B = new intptr [n];
232    for ( i = 0; i < n; i++ )
[806c18]233        Q[i] = new int[n];
[2dd068]234    QmatFF( f, Q, p );
[38bb4c]235#ifdef DEBUGOUTPUT
[160f8f7]236    DEBOUTLN( cerr, "Q = " );
[38bb4c]237    QprintFF( Q, n );
238#endif /* DEBUGOUTPUT */
[2dd068]239    k = nullSpaceFF( Q, B, n );
[38bb4c]240#ifdef DEBUGOUTPUT
[160f8f7]241    DEBOUTLN( cerr, "Q = " );
[38bb4c]242    QprintFF( Q, n );
243#endif /* DEBUGOUTPUT */
[2dd068]244    F.insert( CFFactor( f, 1 ) );
245    r = 1;
246    len = 1;
247    while ( len < k ) {
[806c18]248        ASSERT( r < k, "fatal fatal" );
249        ListIterator<CFFactor> I = F;
250        while ( I.hasItem() && len < k ) {
251            u = I.getItem().factor();
252            for ( s = 0; s < p && len < k; s++ ) {
253                g = gcd( cfFromIntVec( B[r], n, x ) - s, u );
254                if ( degree( g ) > 0 && g != u ) {
255                    u /= g;
256                    I.append( CFFactor( g, 1 ) );
257                    I.append( CFFactor( u, 1 ) );
258                    I.remove( 1 );
259                    len++;
260                }
261            }
262            I++;
263        }
264        r++;
[2dd068]265    }
266    for ( i = 0; i < n; i++ )
[806c18]267        delete [] Q[i];
[2dd068]268    for ( i = 0; i < r; i++ )
[806c18]269        delete [] B[i];
[2dd068]270    delete [] B;
271    delete [] Q;
272    return F;
273}
274
275CFFList BerlekampFactorGF ( const CanonicalForm & f )
276{
277    CFFList F;
278    int r, len, i, k, n = degree( f );
279    Variable x = f.mvar();
280    CanonicalForm u, g;
281    intptr* Q = new intptr [n];
282    intptr* B = new intptr [n];
283    for ( i = 0; i < n; i++ )
[806c18]284        Q[i] = new int[n];
[2dd068]285    QmatGF( f, Q, gf_q );
[38bb4c]286#ifdef DEBUGOUTPUT
[160f8f7]287    DEBOUTLN( cerr, "Q = " );
[38bb4c]288    QprintGF( Q, n );
289#endif /* DEBUGOUTPUT */
[2dd068]290    k = nullSpaceGF( Q, B, n );
[38bb4c]291#ifdef DEBUGOUTPUT
[160f8f7]292    DEBOUTLN( cerr, "Q = " );
[38bb4c]293    QprintFF( Q, n );
294#endif /* DEBUGOUTPUT */
[2dd068]295    F.insert( CFFactor( f, 1 ) );
296    r = 1;
297    len = 1;
298    GFGenerator s;
299    while ( len < k ) {
[806c18]300        ASSERT( r < k, "fatal fatal" );
301        ListIterator<CFFactor> I = F;
302        while ( I.hasItem() && len < k ) {
303            u = I.getItem().factor();
304            for ( s.reset(); s.hasItems() && len < k; s++ ) {
305                g = gcd( cfFromGFVec( B[r], n, x ) - s.item(), u );
306                if ( degree( g ) > 0 && g != u ) {
307                    u /= g;
308                    I.append( CFFactor( g, 1 ) );
309                    I.append( CFFactor( u, 1 ) );
310                    I.remove( 1 );
311                    len++;
312                }
313            }
314            I++;
315        }
316        r++;
[2dd068]317    }
318    for ( i = 0; i < n; i++ )
[806c18]319        delete [] Q[i];
[2dd068]320    for ( i = 0; i < r; i++ )
[806c18]321        delete [] B[i];
[2dd068]322    delete [] B;
323    delete [] Q;
324    return F;
325}
326// {
327//     CFFList F;
328//     int p = getCharacteristic();
329//     int r, len, k, n = degree( f );
330//     Variable x = f.mvar();
331//     CanonicalForm u, g;
332//     intptr* Q = new intptr [n];
333//     for ( int i = 0; i < n; i++ )
[806c18]334//         Q[i] = new int[n];
[2dd068]335//     QmatGF( f, Q, p );
336// //  Qprint( Q, n );
337//     k = nullSpaceGF( Q, n );
338// //  Qprint( Q, n );
339//     F.insert( CFFactor( f, 1 ) );
340//     r = 1;
341//     len = 1;
342//     GFIterator s;
343//     while ( len < k ) {
[806c18]344//         ListIterator<CFFactor> I = F;
345//         while ( I.hasItem() && len < k ) {
346//             u = I.getItem().factor();
347//             for ( s.reset(); s.hasItems() && len < k; s++ ) {
348//                 g = gcd( cfFromGFVec( Q[r], n, x ) - s.item(), u );
349//                 if ( degree( g ) > 0 && g != u ) {
350//                     u /= g;
351//                     I.append( CFFactor( g, 1 ) );
352//                     I.append( CFFactor( u, 1 ) );
353//                     I.remove( 1 );
354//                     len++;
355//                 }
356//             }
357//             I++;
358//         }
359//         r++;
[2dd068]360//     }
361//     for ( i = 0; i < n; i++ )
[806c18]362//         delete [] Q[i];
[2dd068]363//     return F;
364// }
365
366CFFList FpFactorizeUnivariateB( const CanonicalForm& f, bool issqrfree )
367{
368    CFFList F, G, H;
369    CanonicalForm fac;
370    ListIterator<CFFactor> i, k;
371    int d;
372    bool galoisfield = getGFDegree() > 1;
373
374    if ( LC( f ).isOne() )
[806c18]375        if ( issqrfree )
376            F.append( CFFactor( f, 1 ) );
377        else
378            F = sqrFreeFp( f );
[2dd068]379    else {
[806c18]380        H.append( LC( f ) );
381        if ( issqrfree )
382            F.append( CFFactor( f / LC( f ), 1 ) );
383        else
384            F = sqrFreeFp( f / LC( f ) );
[2dd068]385    }
386    for ( i = F; i.hasItem(); ++i ) {
[806c18]387        d = i.getItem().exp();
388        fac = i.getItem().factor();
389        if ( galoisfield )
390            G = BerlekampFactorGF( fac / LC( fac ) );
391        else
392            G = BerlekampFactorFF( fac / LC( fac ) );
393        for ( k = G; k.hasItem(); ++k ) {
394            fac = k.getItem().factor();
395            H.append( CFFactor( fac / LC( fac ), d ) );
396        }
[2dd068]397    }
398    return H;
399}
Note: See TracBrowser for help on using the repository browser.