source: git/factory/fac_berlekamp.cc @ 160f8f7

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