source: git/factory/fac_berlekamp.cc @ 5ab7d6

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