My Project
fac_berlekamp.cc
Go to the documentation of this file.
1 /* emacs edit mode for this file is -*- C++ -*- */
2 
3 #include <config.h>
4 
5 #include "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 #if !defined(HAVE_FLINT) && !defined(HAVE_NTL)
19 
20 #ifdef DEBUGOUTPUT
21 void QprintFF( int ** Q, int n )
22 {
23  for ( int i = 0; i < n; i++ ) {
24  for ( int j = 0; j < n; j++ )
25  std::cerr << Q[i][j] << " ";
26  std::cerr << std::endl;
27  }
28  std::cerr << std::endl;
29 }
30 #endif /* DEBUGOUTPUT */
31 
32 #ifdef DEBUGOUTPUT
33 void QprintGF( int ** Q, int n )
34 {
35  for ( int i = 0; i < n; i++ ) {
36  for ( int j = 0; j < n; j++ ) {
37  gf_print( std::cerr, Q[i][j] );
38  std::cerr << " ";
39  }
40  std::cerr << std::endl;
41  }
42  std::cerr << std::endl;
43 }
44 #endif /* DEBUGOUTPUT */
45 
46 void QmatFF ( const CanonicalForm & f, int ** Q, int p )
47 {
48  int n = degree( f ), nn = (n-1)*p + 1;
49  int i, m, rn;
50  int * a = new int [n];
51  int * r = new int [n];
52  int * q;
53 
54  q = Q[0]; *q = r[0] = 1; a[0] = 0; q++;
55 
56  for ( i = 1; i < n; i++, q++ )
57  *q = r[i] = a[i] = 0;
58  CFIterator I = f; I++;
59  while ( I.hasTerms() ) {
60  a[I.exp()] = I.coeff().intval();
61  I++;
62  }
63  for ( m = 1; m < nn; m++ ) {
64  rn = r[n-1];
65  for ( i = n-1; i > 0; i-- )
66  r[i] = ff_sub( r[i-1], ff_mul( rn, a[i] ) );
67  r[0] = ff_mul( ff_neg( rn ), a[0] );
68  if ( m % p == 0 ) {
69  q = Q[m/p];
70  for ( i = 0; i < n; i++, q++ )
71  *q = r[i];
72  }
73  }
74  for ( i = 0; i < n; i++ )
75  Q[i][i] = ff_sub( Q[i][i], 1 );
76 
77  delete [] a;
78  delete [] r;
79 }
80 
81 void QmatGF ( const CanonicalForm & f, int ** Q, int p )
82 {
83  int n = degree( f ), nn = (n-1)*p + 1;
84  int i, m, rn;
85  int * a = new int [n];
86  int * r = new int [n];
87  int * q;
88 
89  q = Q[0]; *q = r[0] = gf_one(); a[0] = gf_zero(); q++;
90 
91  for ( i = 1; i < n; i++, q++ )
92  *q = r[i] = a[i] = gf_zero();
93  CFIterator I = f; I++;
94  while ( I.hasTerms() ) {
95  a[I.exp()] = imm2int( I.coeff().getval() );
96  I++;
97  }
98  for ( m = 1; m < nn; m++ ) {
99  rn = r[n-1];
100  for ( i = n-1; i > 0; i-- )
101  r[i] = gf_sub( r[i-1], gf_mul( rn, a[i] ) );
102  r[0] = gf_mul( gf_neg( rn ), a[0] );
103  if ( m % p == 0 ) {
104  q = Q[m/p];
105  for ( i = 0; i < n; i++, q++ )
106  *q = r[i];
107  }
108  }
109  for ( i = 0; i < n; i++ )
110  Q[i][i] = gf_sub( Q[i][i], gf_one() );
111 
112  delete [] a;
113  delete [] r;
114 }
115 
116 int nullSpaceFF ( int ** Q, int ** b, int n )
117 {
118  int * c = new int[n];
119  int r, i, j, k, h, s, d;
120 
121  r = 0;
122  for ( s = 0; s < n; s++ ) c[s] = -1;
123  for ( h = 0; h < n; h++ ) {
124  j = 0;
125  while ( j < n && ! ( Q[h][j] != 0 && c[j] < 0 ) ) j++;
126  if ( j < n ) {
127  d = ff_neg( ff_inv( Q[h][j] ) );
128  for ( s = 0; s < n; s++ )
129  Q[s][j] = ff_mul( d, Q[s][j] );
130  for ( i = 0; i < n; i++ ) {
131  if ( i != j ) {
132  d = Q[h][i];
133  for ( s = 0; s < n; s++ )
134  Q[s][i] = ff_add( ff_mul( d, Q[s][j] ), Q[s][i] );
135  }
136  }
137  c[j] = h;
138  }
139  else {
140  b[r] = new int[n];
141  for ( j = 0; j < n; j++ ) {
142  if ( j == h )
143  b[r][j] = 1;
144  else {
145  k = 0;
146  while ( k < n && c[k] != j ) k++;
147  if ( k < n )
148  b[r][j] = Q[h][k];
149  else
150  b[r][j] = 0;
151  }
152  }
153  r++;
154  }
155  }
156  delete [] c;
157  return r;
158 }
159 
160 int nullSpaceGF ( int ** Q, int ** b, int n )
161 {
162  int * c = new int[n];
163  int r, i, j, k, h, s, d;
164 
165  r = 0;
166  for ( s = 0; s < n; s++ ) c[s] = -1;
167  for ( h = 0; h < n; h++ ) {
168  j = 0;
169  while ( j < n && ! ( ! gf_iszero( Q[h][j] ) && c[j] < 0 ) ) j++;
170  if ( j < n ) {
171  d = gf_neg( gf_inv( Q[h][j] ) );
172  for ( s = 0; s < n; s++ )
173  Q[s][j] = gf_mul( d, Q[s][j] );
174  for ( i = 0; i < n; i++ ) {
175  if ( i != j ) {
176  d = Q[h][i];
177  for ( s = 0; s < n; s++ )
178  Q[s][i] = gf_add( gf_mul( d, Q[s][j] ), Q[s][i] );
179  }
180  }
181  c[j] = h;
182  }
183  else {
184  b[r] = new int[n];
185  for ( j = 0; j < n; j++ ) {
186  if ( j == h )
187  b[r][j] = gf_one();
188  else {
189  k = 0;
190  while ( k < n && c[k] != j ) k++;
191  if ( k < n )
192  b[r][j] = Q[h][k];
193  else
194  b[r][j] = gf_zero();
195  }
196  }
197  r++;
198  }
199  }
200  delete [] c;
201  return r;
202 }
203 
204 CanonicalForm cfFromIntVec( int * a, int n, const Variable & x )
205 {
206  CanonicalForm result = power( x, n-1 ) * a[n-1];
207  for ( int i = n-2; i >= 0; i-- )
208  if ( a[i] != 0 )
209  result += power( x, i ) * a[i];
210  return result;
211 }
212 
213 CanonicalForm cfFromGFVec( int * a, int n, const Variable & x )
214 {
215  CanonicalForm result = power( x, n-1 ) * CanonicalForm( int2imm_gf( a[n-1] ) );
216  for ( int i = n-2; i >= 0; i-- )
217  if ( ! gf_iszero( a[i] ) )
218  result += power( x, i ) * CanonicalForm( int2imm_gf( a[i] ) );
219  return result;
220 }
221 
222 typedef int * intptr;
223 
224 CFFList BerlekampFactorFF ( const CanonicalForm & f )
225 {
226  CFFList F;
227  int p = getCharacteristic();
228  int r, s, len, i, k, n = degree( f );
229  Variable x = f.mvar();
230  CanonicalForm u, g;
231  intptr* Q = new intptr [n];
232  intptr* B = new intptr [n];
233  for ( i = 0; i < n; i++ )
234  Q[i] = new int[n];
235  QmatFF( f, Q, p );
236 #ifdef DEBUGOUTPUT
237  DEBOUTLN( cerr, "Q = " );
238  QprintFF( Q, n );
239 #endif /* DEBUGOUTPUT */
240  k = nullSpaceFF( Q, B, n );
241 #ifdef DEBUGOUTPUT
242  DEBOUTLN( cerr, "Q = " );
243  QprintFF( Q, n );
244 #endif /* DEBUGOUTPUT */
245  F.insert( CFFactor( f, 1 ) );
246  r = 1;
247  len = 1;
248  while ( len < k ) {
249  ASSERT( r < k, "fatal fatal" );
251  while ( I.hasItem() && len < k ) {
252  u = I.getItem().factor();
253  for ( s = 0; s < p && len < k; s++ ) {
254  g = gcd( cfFromIntVec( B[r], n, x ) - s, u );
255  if ( degree( g ) > 0 && g != u ) {
256  u /= g;
257  I.append( CFFactor( g, 1 ) );
258  I.append( CFFactor( u, 1 ) );
259  I.remove( 1 );
260  len++;
261  }
262  }
263  I++;
264  }
265  r++;
266  }
267  for ( i = 0; i < n; i++ )
268  delete [] Q[i];
269  for ( i = 0; i < r; i++ )
270  delete [] B[i];
271  delete [] B;
272  delete [] Q;
273  return F;
274 }
275 
276 CFFList BerlekampFactorGF ( const CanonicalForm & f )
277 {
278  CFFList F;
279  int r, len, i, k, n = degree( f );
280  Variable x = f.mvar();
281  CanonicalForm u, g;
282  intptr* Q = new intptr [n];
283  intptr* B = new intptr [n];
284  for ( i = 0; i < n; i++ )
285  Q[i] = new int[n];
286  QmatGF( f, Q, gf_q );
287 #ifdef DEBUGOUTPUT
288  DEBOUTLN( cerr, "Q = " );
289  QprintGF( Q, n );
290 #endif /* DEBUGOUTPUT */
291  k = nullSpaceGF( Q, B, n );
292 #ifdef DEBUGOUTPUT
293  DEBOUTLN( cerr, "Q = " );
294  QprintFF( Q, n );
295 #endif /* DEBUGOUTPUT */
296  F.insert( CFFactor( f, 1 ) );
297  r = 1;
298  len = 1;
299  GFGenerator s;
300  while ( len < k ) {
301  ASSERT( r < k, "fatal fatal" );
303  while ( I.hasItem() && len < k ) {
304  u = I.getItem().factor();
305  for ( s.reset(); s.hasItems() && len < k; s++ ) {
306  g = gcd( cfFromGFVec( B[r], n, x ) - s.item(), u );
307  if ( degree( g ) > 0 && g != u ) {
308  u /= g;
309  I.append( CFFactor( g, 1 ) );
310  I.append( CFFactor( u, 1 ) );
311  I.remove( 1 );
312  len++;
313  }
314  }
315  I++;
316  }
317  r++;
318  }
319  for ( i = 0; i < n; i++ )
320  delete [] Q[i];
321  for ( i = 0; i < r; i++ )
322  delete [] B[i];
323  delete [] B;
324  delete [] Q;
325  return F;
326 }
327 // {
328 // CFFList F;
329 // int p = getCharacteristic();
330 // int r, len, k, n = degree( f );
331 // Variable x = f.mvar();
332 // CanonicalForm u, g;
333 // intptr* Q = new intptr [n];
334 // for ( int i = 0; i < n; i++ )
335 // Q[i] = new int[n];
336 // QmatGF( f, Q, p );
337 // // Qprint( Q, n );
338 // k = nullSpaceGF( Q, n );
339 // // Qprint( Q, n );
340 // F.insert( CFFactor( f, 1 ) );
341 // r = 1;
342 // len = 1;
343 // GFIterator s;
344 // while ( len < k ) {
345 // ListIterator<CFFactor> I = F;
346 // while ( I.hasItem() && len < k ) {
347 // u = I.getItem().factor();
348 // for ( s.reset(); s.hasItems() && len < k; s++ ) {
349 // g = gcd( cfFromGFVec( Q[r], n, x ) - s.item(), u );
350 // if ( degree( g ) > 0 && g != u ) {
351 // u /= g;
352 // I.append( CFFactor( g, 1 ) );
353 // I.append( CFFactor( u, 1 ) );
354 // I.remove( 1 );
355 // len++;
356 // }
357 // }
358 // I++;
359 // }
360 // r++;
361 // }
362 // for ( i = 0; i < n; i++ )
363 // delete [] Q[i];
364 // return F;
365 // }
366 
367 CFFList FpFactorizeUnivariateB( const CanonicalForm& f, bool issqrfree )
368 {
369  CFFList F, G, H;
370  CanonicalForm fac;
372  int d;
373  bool galoisfield = getGFDegree() > 1;
374 
375  if ( LC( f ).isOne() )
376  if ( issqrfree )
377  F.append( CFFactor( f, 1 ) );
378  else
379  F = sqrFreeFp( f );
380  else {
381  H.append( LC( f ) );
382  if ( issqrfree )
383  F.append( CFFactor( f / LC( f ), 1 ) );
384  else
385  F = sqrFreeFp( f / LC( f ) );
386  }
387  for ( i = F; i.hasItem(); ++i ) {
388  d = i.getItem().exp();
389  fac = i.getItem().factor();
390  if ( galoisfield )
391  G = BerlekampFactorGF( fac / LC( fac ) );
392  else
393  G = BerlekampFactorFF( fac / LC( fac ) );
394  for ( k = G; k.hasItem(); ++k ) {
395  fac = k.getItem().factor();
396  H.append( CFFactor( fac / LC( fac ), d ) );
397  }
398  }
399  return H;
400 }
401 #endif
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Header for factory's main class CanonicalForm.
int degree(const CanonicalForm &f)
int getGFDegree()
Definition: cf_char.cc:75
Factor< CanonicalForm > CFFactor
CanonicalForm LC(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
#define ASSERT(expression, message)
Definition: cf_assert.h:99
factory switches.
generate integers, elements of finite fields
Iterators for CanonicalForm's.
FILE * f
Definition: checklibs.c:9
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
CF_NO_INLINE int exp() const
get the current exponent
CF_NO_INLINE CanonicalForm coeff() const
get the current coefficient
CF_NO_INLINE int hasTerms() const
check if iterator has reached the end of CanonicalForm
factory's main class
Definition: canonicalform.h:86
InternalCF * getval() const
long intval() const
conversion functions
generate all elements in GF starting from 0
Definition: cf_generator.h:75
void remove(int moveright)
Definition: ftmpl_list.cc:526
void append(const T &)
Definition: ftmpl_list.cc:509
T & getItem() const
Definition: ftmpl_list.cc:431
void append(const T &)
Definition: ftmpl_list.cc:256
void insert(const T &)
Definition: ftmpl_list.cc:193
factory's class for variables
Definition: factory.h:127
functions to print debug output
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm H
Definition: facAbsFact.cc:60
b *CanonicalForm B
Definition: facBivar.cc:52
int j
Definition: facHensel.cc:110
CFFList FpFactorizeUnivariateB(const CanonicalForm &f, bool issqrfree=false)
squarefree part and factorization over Q, Q(a)
CFFList sqrFreeFp(const CanonicalForm &f)
operations in a finite prime field F_p.
int ff_inv(const int a)
Definition: ffops.h:149
int ff_add(const int a, const int b)
Definition: ffops.h:97
int ff_mul(const int a, const int b)
Definition: ffops.h:139
int ff_neg(const int a)
Definition: ffops.h:126
int ff_sub(const int a, const int b)
Definition: ffops.h:112
VAR int gf_q
Definition: gfops.cc:47
Operations in GF, where GF is a finite field of size less than 2^16 represented by a root of Conway p...
int gf_sub(int a, int b)
Definition: gfops.h:158
void gf_print(OSTREAM &os, int a)
Definition: gfops.h:207
int gf_neg(int a)
Definition: gfops.h:123
int gf_inv(int a)
Definition: gfops.h:198
int gf_zero()
Definition: gfops.h:99
bool gf_iszero(int a)
Definition: gfops.h:43
int gf_one()
Definition: gfops.h:104
int gf_mul(int a, int b)
Definition: gfops.h:163
int gf_add(int a, int b)
Definition: gfops.h:133
operations on immediates, that is elements of F_p, GF, Z, Q that fit into intrinsic int,...
static long imm2int(const InternalCF *const imm)
Definition: imm.h:70
InternalCF * int2imm_gf(long i)
Definition: imm.h:106
STATIC_VAR TreeM * G
Definition: janet.cc:31
STATIC_VAR Poly * h
Definition: janet.cc:971
STATIC_VAR jList * Q
Definition: janet.cc:30
int gcd(int a, int b)
Definition: walkSupport.cc:836