My Project
Loading...
Searching...
No Matches
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
21void 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
33void 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
46void 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
81void 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
116int 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
160int 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
204CanonicalForm 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
213CanonicalForm 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
222typedef int * intptr;
223
224CFFList 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
276CFFList 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;
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
367CFFList 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: variable.h:33
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
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
#define Q
Definition: sirandom.c:26
int gcd(int a, int b)
Definition: walkSupport.cc:836