My Project
Loading...
Searching...
No Matches
cf_algorithm.h File Reference

declarations of higher level algorithms. More...

#include "canonicalform.h"
#include "variable.h"

Go to the source code of this file.

Functions

CanonicalForm psr (const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
 CanonicalForm psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
 
CanonicalForm psq (const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
 CanonicalForm psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
 
void psqr (const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r, const Variable &x)
 void psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable & x )
 
CanonicalForm FACTORY_PUBLIC bCommonDen (const CanonicalForm &f)
 CanonicalForm bCommonDen ( const CanonicalForm & f )
 
bool fdivides (const CanonicalForm &f, const CanonicalForm &g)
 bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
 
bool fdivides (const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &quot)
 same as fdivides if true returns quotient quot of g by f otherwise quot == 0
 
bool tryFdivides (const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
 same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f
 
CanonicalForm maxNorm (const CanonicalForm &f)
 CanonicalForm maxNorm ( const CanonicalForm & f )
 
CanonicalForm euclideanNorm (const CanonicalForm &f)
 CanonicalForm euclideanNorm ( const CanonicalForm & f )
 
void FACTORY_PUBLIC chineseRemainder (const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
 void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
 
void FACTORY_PUBLIC chineseRemainder (const CFArray &x, const CFArray &q, CanonicalForm &xnew, CanonicalForm &qnew)
 void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
 
void FACTORY_PUBLIC chineseRemainderCached (const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew, CFArray &inv)
 
void FACTORY_PUBLIC chineseRemainderCached (const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
 
CanonicalForm Farey (const CanonicalForm &f, const CanonicalForm &q)
 Farey rational reconstruction.
 
bool isPurePoly (const CanonicalForm &f)
 
bool isPurePoly_m (const CanonicalForm &f)
 
CFFList FACTORY_PUBLIC factorize (const CanonicalForm &f, bool issqrfree=false)
 factorization over $ F_p $ or $ Q $
 
CFFList FACTORY_PUBLIC factorize (const CanonicalForm &f, const Variable &alpha)
 factorization over $ F_p(\alpha) $ or $ Q(\alpha) $
 
CFFList FACTORY_PUBLIC sqrFree (const CanonicalForm &f, bool sort=false)
 squarefree factorization
 
CanonicalForm homogenize (const CanonicalForm &f, const Variable &x)
 homogenize homogenizes f with Variable x
 
CanonicalForm homogenize (const CanonicalForm &f, const Variable &x, const Variable &v1, const Variable &v2)
 
Variable get_max_degree_Variable (const CanonicalForm &f)
 get_max_degree_Variable returns Variable with highest degree.
 
CFList get_Terms (const CanonicalForm &f)
 
void getTerms (const CanonicalForm &f, const CanonicalForm &t, CFList &result)
 get_Terms: Split the polynomial in the containing terms.
 
bool linearSystemSolve (CFMatrix &M)
 
CanonicalForm FACTORY_PUBLIC determinant (const CFMatrix &M, int n)
 
CFArray subResChain (const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
 CFArray subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
 
CanonicalForm FACTORY_PUBLIC resultant (const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
 CanonicalForm resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
 
CanonicalForm abs (const CanonicalForm &f)
 inline CanonicalForm abs ( const CanonicalForm & f )
 

Variables

EXTERN_VAR int singular_homog_flag
 

Detailed Description

declarations of higher level algorithms.

Header file corresponds to: cf_algorithm.cc, cf_chinese.cc, cf_factor.cc, cf_linsys.cc, cf_resultant.cc

Hierarchy: mathematical algorithms on canonical forms

Developers note:

This header file collects declarations of most of the functions in Factory which implement higher level algorithms on canonical forms (factorization, gcd, etc.) and declarations of some low level mathematical functions, too (absolute value, euclidean norm, etc.).

Definition in file cf_algorithm.h.

Function Documentation

◆ abs()

CanonicalForm abs ( const CanonicalForm & f)
inline

inline CanonicalForm abs ( const CanonicalForm & f )

abs() - return absolute value of ‘f’.

The absolute value is defined in terms of the function ‘sign()’. If it reports negative sign for ‘f’ than -‘f’ is returned, otherwise ‘f’.

This behaviour is most useful for integers and rationals. But it may be used to sign-normalize the leading coefficient of arbitrary polynomials, too.

Type info:

f: CurrentPP

Definition at line 117 of file cf_algorithm.h.

118{
119 // it is not only more general to use `sign()' instead of a
120 // direct comparison `f < 0', it is faster, too
121 if ( sign( f ) < 0 )
122 return -f;
123 else
124 return f;
125}
FILE * f
Definition checklibs.c:9
static int sign(int x)
Definition ring.cc:3503

◆ bCommonDen()

CanonicalForm bCommonDen ( const CanonicalForm & f )

bCommonDen() - calculate multivariate common denominator of coefficients of ‘f’.

The common denominator is calculated with respect to all coefficients of ‘f’ which are in a base domain. In other words, common_den( ‘f’ ) * ‘f’ is guaranteed to have integer coefficients only. The common denominator of zero is one.

Returns something non-trivial iff the current domain is Q.

Type info:

f: CurrentPP

Definition at line 293 of file cf_algorithm.cc.

294{
295 if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) )
296 {
297 // otherwise `bgcd()' returns one
298 Off( SW_RATIONAL );
300 On( SW_RATIONAL );
301 return result;
302 }
303 else
304 return CanonicalForm( 1 );
305}
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
static CanonicalForm internalBCommonDen(const CanonicalForm &f)
static CanonicalForm internalBCommonDen ( const CanonicalForm & f )
static const int SW_RATIONAL
set to 1 for computations over Q
Definition cf_defs.h:31
factory's main class
return result

◆ chineseRemainder() [1/2]

void FACTORY_PUBLIC chineseRemainder ( const CanonicalForm & x1,
const CanonicalForm & q1,
const CanonicalForm & x2,
const CanonicalForm & q2,
CanonicalForm & xnew,
CanonicalForm & qnew )

void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )

chineseRemainder - integer chinese remaindering.

Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2) and qnew = q1*q2. q1 and q2 should be positive integers, pairwise prime, x1 and x2 should be polynomials with integer coefficients. If x1 and x2 are polynomials with positive coefficients, the result is guaranteed to have positive coefficients, too.

Note: This algorithm is optimized for the case q1>>q2.

This is a standard algorithm. See, for example, Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra', par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by Homomorphic Images' in B. Buchberger - 'Computer Algebra - Symbolic and Algebraic Computation'.

Note: Be sure you are calculating in Z, and not in Q!

Definition at line 57 of file cf_chinese.cc.

58{
59 DEBINCLEVEL( cerr, "chineseRemainder" );
60
61 DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
62 DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
63
64 // We calculate xnew as follows:
65 // xnew = v1 + v2 * q1
66 // where
67 // v1 = x1 (mod q1)
68 // v2 = (x2-v1)/q1 (mod q2) (*)
69 //
70 // We do one extra test to check whether x2-v1 vanishes (mod
71 // q2) in (*) since it is not costly and may save us
72 // from calculating the inverse of q1 (mod q2).
73 //
74 // u: v1 (mod q2)
75 // d: x2-v1 (mod q2)
76 // s: 1/q1 (mod q2)
77 //
78 CanonicalForm v2, v1;
79 CanonicalForm u, d, s, dummy;
80
81 v1 = mod( x1, q1 );
82 u = mod( v1, q2 );
83 d = mod( x2-u, q2 );
84 if ( d.isZero() )
85 {
86 xnew = v1;
87 qnew = q1 * q2;
88 DEBDECLEVEL( cerr, "chineseRemainder" );
89 return;
90 }
91 (void)bextgcd( q1, q2, s, dummy );
92 v2 = mod( d*s, q2 );
93 xnew = v1 + v2*q1;
94
95 // After all, calculate new modulus. It is important that
96 // this is done at the very end of the algorithm, since q1
97 // and qnew may refer to the same object (same is true for x1
98 // and xnew).
99 qnew = q1 * q2;
100
101 DEBDECLEVEL( cerr, "chineseRemainder" );
102}
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CF_NO_INLINE bool isZero() const
int ilog2() const
int CanonicalForm::ilog2 () const
#define DEBINCLEVEL(stream, msg)
Definition debug.h:44
#define DEBOUTLN(stream, objects)
Definition debug.h:49
#define DEBDECLEVEL(stream, msg)
Definition debug.h:45
const CanonicalForm int s
Definition facAbsFact.cc:51

◆ chineseRemainder() [2/2]

void FACTORY_PUBLIC chineseRemainder ( const CFArray & x,
const CFArray & q,
CanonicalForm & xnew,
CanonicalForm & qnew )

void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )

chineseRemainder - integer chinese remaindering.

Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the product of all q[i]. q[i] should be positive integers, pairwise prime. x[i] should be polynomials with integer coefficients. If all coefficients of all x[i] are positive integers, the result is guaranteed to have positive coefficients, too.

This is a standard algorithm, too, except for the fact that we use a divide-and-conquer method instead of a linear approach to calculate the remainder.

Note: Be sure you are calculating in Z, and not in Q!

Definition at line 124 of file cf_chinese.cc.

125{
126 DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
127
128 ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
129 CFArray X(x), Q(q);
130 int i, j, n = x.size(), start = x.min();
131
132 DEBOUTLN( cerr, "array size = " << n );
133
134 while ( n != 1 )
135 {
136 i = j = start;
137 while ( i < start + n - 1 )
138 {
139 // This is a little bit dangerous: X[i] and X[j] (and
140 // Q[i] and Q[j]) may refer to the same object. But
141 // xnew and qnew in the above function are modified
142 // at the very end of the function, so we do not
143 // modify x1 and q1, resp., by accident.
144 chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
145 i += 2;
146 j++;
147 }
148
149 if ( n & 1 )
150 {
151 X[j] = X[i];
152 Q[j] = Q[i];
153 }
154 // Maybe we would get some memory back at this point if
155 // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
156 // at this point?
157
158 n = ( n + 1) / 2;
159 }
160 xnew = X[start];
161 qnew = Q[q.min()];
162
163 DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
164}
Array< CanonicalForm > CFArray
int i
Definition cfEzgcd.cc:132
Variable x
Definition cfModGcd.cc:4090
#define ASSERT(expression, message)
Definition cf_assert.h:99
void chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition cf_chinese.cc:57
int size() const
int min() const
int j
Definition facHensel.cc:110
#define Q
Definition sirandom.c:26

◆ chineseRemainderCached() [1/2]

void FACTORY_PUBLIC chineseRemainderCached ( const CanonicalForm & x1,
const CanonicalForm & q1,
const CanonicalForm & x2,
const CanonicalForm & q2,
CanonicalForm & xnew,
CanonicalForm & qnew,
CFArray & inv )

Definition at line 308 of file cf_chinese.cc.

309{
310 CFArray A(2); A[0]=a; A[1]=b;
311 CFArray Q(2); Q[0]=q1; Q[1]=q2;
312 chineseRemainderCached(A,Q,xnew,qnew,inv);
313}
CanonicalForm b
Definition cfModGcd.cc:4111
void chineseRemainderCached(const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
#define A
Definition sirandom.c:24

◆ chineseRemainderCached() [2/2]

void FACTORY_PUBLIC chineseRemainderCached ( const CFArray & a,
const CFArray & n,
CanonicalForm & xnew,
CanonicalForm & prod,
CFArray & inv )

Definition at line 290 of file cf_chinese.cc.

291{
292 CanonicalForm p, sum=0L; prod=1L;
293 int i;
294 int len=n.size();
295
296 for (i = 0; i < len; i++) prod *= n[i];
297
298 for (i = 0; i < len; i++)
299 {
300 p = prod / n[i];
301 sum += a[i] * chin_mul_inv(p, n[i], i, inv) * p;
302 }
303
304 xnew = mod(sum , prod);
305}
int p
Definition cfModGcd.cc:4086
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
fq_nmod_poly_t prod
Definition facHensel.cc:100

◆ determinant()

CanonicalForm FACTORY_PUBLIC determinant ( const CFMatrix & M,
int n )

Definition at line 222 of file cf_linsys.cc.

223{
224 typedef int* int_ptr;
225
226 ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
227 if ( rows == 1 )
228 return M(1,1);
229 else if ( rows == 2 )
230 return M(1,1)*M(2,2)-M(2,1)*M(1,2);
231 else if ( matrix_in_Z( M, rows ) )
232 {
233 int ** mm = new int_ptr[rows];
234 CanonicalForm x, q, Qhalf, B;
235 int n, i, intdet, p, pno;
236 for ( i = 0; i < rows; i++ )
237 {
238 mm[i] = new int[rows];
239 }
240 pno = 0; n = 0;
241 TIMING_START(det_bound);
242 B = detbound( M, rows );
243 TIMING_END(det_bound);
244 q = 1;
245 TIMING_START(det_numprimes);
246 while ( B > q && n < cf_getNumBigPrimes() )
247 {
248 q *= cf_getBigPrime( n );
249 n++;
250 }
251 TIMING_END(det_numprimes);
252
253 CFArray X(1,n), Q(1,n);
254
255 while ( pno < n )
256 {
257 p = cf_getBigPrime( pno );
259 // map matrix into char p
260 TIMING_START(det_mapping);
261 fill_int_mat( M, mm, rows );
262 TIMING_END(det_mapping);
263 pno++;
264 DEBOUT( cerr, "." );
265 TIMING_START(det_determinant);
266 intdet = determinant( mm, rows );
267 TIMING_END(det_determinant);
269 X[pno] = intdet;
270 Q[pno] = p;
271 }
272 TIMING_START(det_chinese);
273 chineseRemainder( X, Q, x, q );
274 TIMING_END(det_chinese);
275 Qhalf = q / 2;
276 if ( x > Qhalf )
277 x = x - q;
278 for ( i = 0; i < rows; i++ )
279 delete [] mm[i];
280 delete [] mm;
281 return x;
282 }
283 else
284 {
285 CFMatrix m( M );
286 CanonicalForm divisor = 1, pivot, mji;
287 int i, j, k, sign = 1;
288 for ( i = 1; i <= rows; i++ ) {
289 pivot = m(i,i); k = i;
290 for ( j = i+1; j <= rows; j++ ) {
291 if ( betterpivot( pivot, m(j,i) ) ) {
292 pivot = m(j,i);
293 k = j;
294 }
295 }
296 if ( pivot.isZero() )
297 return 0;
298 if ( i != k )
299 {
300 m.swapRow( i, k );
301 sign = -sign;
302 }
303 for ( j = i+1; j <= rows; j++ )
304 {
305 if ( ! m(j,i).isZero() )
306 {
307 divisor *= pivot;
308 mji = m(j,i);
309 m(j,i) = 0;
310 for ( k = i+1; k <= rows; k++ )
311 m(j,k) = m(j,k) * pivot - m(i,k)*mji;
312 }
313 }
314 }
315 pivot = sign;
316 for ( i = 1; i <= rows; i++ )
317 pivot *= m(i,i);
318 return pivot / divisor;
319 }
320}
void FACTORY_PUBLIC setCharacteristic(int c)
Definition cf_char.cc:28
Matrix< CanonicalForm > CFMatrix
int m
Definition cfEzgcd.cc:128
int k
Definition cfEzgcd.cc:99
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition cf_chinese.cc:57
bool matrix_in_Z(const CFMatrix &M, int rows)
Definition cf_linsys.cc:39
bool betterpivot(const CanonicalForm &oldpivot, const CanonicalForm &newpivot)
Definition cf_linsys.cc:61
CanonicalForm detbound(const CFMatrix &M, int rows)
Definition cf_linsys.cc:486
int determinant(int **extmat, int n)
Definition cf_linsys.cc:556
static bool fill_int_mat(const CFMatrix &M, int **m, int rows)
Definition cf_linsys.cc:203
int cf_getBigPrime(int i)
Definition cf_primes.cc:39
int cf_getNumBigPrimes()
Definition cf_primes.cc:45
#define DEBOUT(stream, objects)
Definition debug.h:47
b *CanonicalForm B
Definition facBivar.cc:52
bool isZero(const CFArray &A)
checks if entries of A are zero
bool pivot(const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring R)
This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2].
#define M
Definition sirandom.c:25
int * int_ptr
Definition structs.h:50
#define TIMING_END(t)
Definition timing.h:93
#define TIMING_START(t)
Definition timing.h:92

◆ euclideanNorm()

CanonicalForm euclideanNorm ( const CanonicalForm & f)

CanonicalForm euclideanNorm ( const CanonicalForm & f )

euclideanNorm() - return Euclidean norm of ‘f’.

Returns the largest integer smaller or equal norm(‘f’) = sqrt(sum( ‘f’[i]^2 )).

Type info:

f: UVPoly( Z )

Definition at line 565 of file cf_algorithm.cc.

566{
567 ASSERT( (f.inBaseDomain() || f.isUnivariate()) && f.LC().inZ(),
568 "type error: univariate poly over Z expected" );
569
571 for ( CFIterator i = f; i.hasTerms(); i++ ) {
572 CanonicalForm coeff = i.coeff();
573 result += coeff*coeff;
574 }
575 return sqrt( result );
576}
class to iterate through CanonicalForm's
Definition cf_iter.h:44
gmp_float sqrt(const gmp_float &a)

◆ factorize() [1/2]

CFFList FACTORY_PUBLIC factorize ( const CanonicalForm & f,
bool issqrfree = false )

factorization over $ F_p $ or $ Q $

Definition at line 409 of file cf_factor.cc.

410{
411 if ( f.inCoeffDomain() )
412 return CFFList( f );
413#ifndef NOASSERT
414 Variable a;
415 ASSERT (!hasFirstAlgVar (f, a), "f has an algebraic variable use factorize \
416 ( const CanonicalForm & f, const Variable & alpha ) instead");
417#endif
418 //out_cf("factorize:",f,"==================================\n");
419 if (! f.isUnivariate() ) // preprocess homog. polys
420 {
421 if ( singular_homog_flag && f.isHomogeneous())
422 {
424 int d_xn = degree(f,xn);
425 CFMap n;
426 CanonicalForm F = compress(f(1,xn),n);
427 CFFList Intermediatelist;
428 Intermediatelist = factorize(F);
429 CFFList Homoglist;
431 for ( j=Intermediatelist; j.hasItem(); j++ )
432 {
433 Homoglist.append(
434 CFFactor( n(j.getItem().factor()), j.getItem().exp()) );
435 }
436 CFFList Unhomoglist;
437 CanonicalForm unhomogelem;
438 for ( j=Homoglist; j.hasItem(); j++ )
439 {
440 unhomogelem= homogenize(j.getItem().factor(),xn);
441 Unhomoglist.append(CFFactor(unhomogelem,j.getItem().exp()));
442 d_xn -= (degree(unhomogelem,xn)*j.getItem().exp());
443 }
444 if ( d_xn != 0 ) // have to append xn^(d_xn)
445 Unhomoglist.append(CFFactor(CanonicalForm(xn),d_xn));
446 if(isOn(SW_USE_NTL_SORT)) Unhomoglist.sort(cmpCF);
447 return Unhomoglist;
448 }
449 }
450 CFFList F;
451 if ( getCharacteristic() > 0 )
452 {
453 if (f.isUnivariate())
454 {
455#ifdef HAVE_FLINT
456#ifdef HAVE_NTL
457 if (degree (f) < 300)
458#endif
459 {
460 // use FLINT: char p, univariate
461 nmod_poly_t f1;
463 nmod_poly_factor_t result;
464 nmod_poly_factor_init (result);
465 mp_limb_t leadingCoeff= nmod_poly_factor (result, f1);
466 F= convertFLINTnmod_poly_factor2FacCFFList (result, leadingCoeff, f.mvar());
467 nmod_poly_factor_clear (result);
468 nmod_poly_clear (f1);
470 return F;
471 }
472#endif
473#ifdef HAVE_NTL
474 { // NTL char 2, univariate
475 if (getCharacteristic()==2)
476 {
477 // Specialcase characteristic==2
478 if (fac_NTL_char != 2)
479 {
480 fac_NTL_char = 2;
481 zz_p::init(2);
482 }
483 // convert to NTL using the faster conversion routine for characteristic 2
484 GF2X f1=convertFacCF2NTLGF2X(f);
485 // no make monic necessary in GF2
486 //factorize
487 vec_pair_GF2X_long factors;
488 CanZass(factors,f1);
489
490 // convert back to factory again using the faster conversion routine for vectors over GF2X
491 F=convertNTLvec_pair_GF2X_long2FacCFFList(factors,LeadCoeff(f1),f.mvar());
493 return F;
494 }
495 }
496#endif
497#ifdef HAVE_NTL
498 {
499 // use NTL char p, univariate
501 {
503 zz_p::init(getCharacteristic());
504 }
505
506 // convert to NTL
507 zz_pX f1=convertFacCF2NTLzzpX(f);
508 zz_p leadcoeff = LeadCoeff(f1);
509
510 //make monic
511 f1=f1 / LeadCoeff(f1);
512 // factorize
513 vec_pair_zz_pX_long factors;
514 CanZass(factors,f1);
515
516 F=convertNTLvec_pair_zzpX_long2FacCFFList(factors,leadcoeff,f.mvar());
517 //test_cff(F,f);
519 return F;
520 }
521#endif
522#if !defined(HAVE_NTL) && !defined(HAVE_FLINT)
523 // Use Factory without NTL and without FLINT: char p, univariate
524 {
525 if ( isOn( SW_BERLEKAMP ) )
526 F=FpFactorizeUnivariateB( f, issqrfree );
527 else
528 F=FpFactorizeUnivariateCZ( f, issqrfree, 0, Variable(), Variable() );
529 return F;
530 }
531#endif
532 }
533 else // char p, multivariate
534 {
536 {
537 #if defined(HAVE_NTL)
538 if (issqrfree)
539 {
540 CFList factors;
542 factors= GFSqrfFactorize (f);
543 for (CFListIterator i= factors; i.hasItem(); i++)
544 F.append (CFFactor (i.getItem(), 1));
545 }
546 else
547 {
549 F= GFFactorize (f);
550 }
551 #else
552 factoryError ("multivariate factorization over GF depends on NTL(missing)");
553 return CFFList (CFFactor (f, 1));
554 #endif
555 }
556 else
557 {
558 #if defined(HAVE_FLINT) && (__FLINT_RELEASE >= 20700) && defined(HAVE_NTL)
559 if (!isOn(SW_USE_FL_FAC_P))
560 {
561 #endif
562 #if defined(HAVE_NTL)
563 if (issqrfree)
564 {
565 CFList factors;
567 factors= FpSqrfFactorize (f);
568 for (CFListIterator i= factors; i.hasItem(); i++)
569 F.append (CFFactor (i.getItem(), 1));
570 goto end_charp;
571 }
572 else
573 {
575 F= FpFactorize (f);
576 goto end_charp;
577 }
578 #endif
579 #if defined(HAVE_FLINT) && (__FLINT_RELEASE >= 20700) && defined(HAVE_NTL)
580 }
581 #endif
582 #if defined(HAVE_FLINT) && (__FLINT_RELEASE >= 20700)
583 nmod_mpoly_ctx_t ctx;
584 nmod_mpoly_ctx_init(ctx,f.level(),ORD_LEX,getCharacteristic());
585 nmod_mpoly_t Flint_f;
586 nmod_mpoly_init(Flint_f,ctx);
587 convFactoryPFlintMP(f,Flint_f,ctx,f.level());
588 nmod_mpoly_factor_t factors;
589 nmod_mpoly_factor_init(factors,ctx);
590 int okay;
591 if (issqrfree) okay=nmod_mpoly_factor_squarefree(factors,Flint_f,ctx);
592 else okay=nmod_mpoly_factor(factors,Flint_f,ctx);
593 nmod_mpoly_t fac;
594 nmod_mpoly_init(fac,ctx);
595 CanonicalForm cf_fac;
596 int cf_exp;
597 cf_fac=nmod_mpoly_factor_get_constant_ui(factors,ctx);
598 F.append(CFFactor(cf_fac,1));
599 for(int i=nmod_mpoly_factor_length(factors,ctx)-1; i>=0; i--)
600 {
601 nmod_mpoly_factor_get_base(fac,factors,i,ctx);
602 cf_fac=convFlintMPFactoryP(fac,ctx,f.level());
603 cf_exp=nmod_mpoly_factor_get_exp_si(factors,i,ctx);
604 F.append(CFFactor(cf_fac,cf_exp));
605 }
606 nmod_mpoly_factor_clear(factors,ctx);
607 nmod_mpoly_clear(Flint_f,ctx);
608 nmod_mpoly_ctx_clear(ctx);
609 if (okay==0)
610 {
613 F=factorize(f,issqrfree);
616 }
617 #endif
618 #if !defined(HAVE_FLINT) || (__FLINT_RELEASE < 20700)
619 #ifndef HAVE_NTL
620 factoryError ("multivariate factorization depends on NTL/FLINT(missing)");
621 return CFFList (CFFactor (f, 1));
622 #endif
623 #endif
624 }
625 }
626 }
627 else // char 0
628 {
629 bool on_rational = isOn(SW_RATIONAL);
632 CanonicalForm fz = f * cd;
634 if ( f.isUnivariate() )
635 {
636 CanonicalForm ic=icontent(fz);
637 fz/=ic;
638 if (fz.degree()==1)
639 {
640 F=CFFList(CFFactor(fz,1));
641 F.insert(CFFactor(ic,1));
642 }
643 else
644 #if defined(HAVE_FLINT) && (__FLINT_RELEASE>=20503) && (__FLINT_RELEASE!= 20600)
645 {
646 // FLINT 2.6.0 has a bug:
647 // factorize x^12-13*x^10-13*x^8+13*x^4+13*x^2-1 runs forever
648 // use FLINT: char 0, univariate
649 fmpz_poly_t f1;
651 fmpz_poly_factor_t result;
652 fmpz_poly_factor_init (result);
653 fmpz_poly_factor(result, f1);
655 fmpz_poly_factor_clear (result);
656 fmpz_poly_clear (f1);
657 if ( ! ic.isOne() )
658 {
659 // according to convertFLINTfmpz_polyfactor2FcaCFFlist,
660 // first entry is in CoeffDomain
661 CFFactor new_first( F.getFirst().factor() * ic );
662 F.removeFirst();
663 F.insert( new_first );
664 }
665 }
666 goto end_char0;
667 #elif defined(HAVE_NTL)
668 {
669 //use NTL
670 ZZ c;
671 vec_pair_ZZX_long factors;
672 //factorize the converted polynomial
673 factor(c,factors,convertFacCF2NTLZZX(fz));
674
675 //convert the result back to Factory
677 if ( ! ic.isOne() )
678 {
679 // according to convertNTLvec_pair_ZZX_long2FacCFFList
680 // first entry is in CoeffDomain
681 CFFactor new_first( F.getFirst().factor() * ic );
682 F.removeFirst();
683 F.insert( new_first );
684 }
685 }
686 goto end_char0;
687 #else
688 {
689 //Use Factory without NTL: char 0, univariate
690 F = ZFactorizeUnivariate( fz, issqrfree );
691 goto end_char0;
692 }
693 #endif
694 }
695 else // multivariate, char 0
696 {
697 #if defined(HAVE_FLINT) && (__FLINT_RELEASE >= 20700)
699 {
700 On (SW_RATIONAL);
701 fmpz_mpoly_ctx_t ctx;
702 fmpz_mpoly_ctx_init(ctx,f.level(),ORD_LEX);
703 fmpz_mpoly_t Flint_f;
704 fmpz_mpoly_init(Flint_f,ctx);
705 convFactoryPFlintMP(fz,Flint_f,ctx,fz.level());
706 fmpz_mpoly_factor_t factors;
707 fmpz_mpoly_factor_init(factors,ctx);
708 int rr;
709 if (issqrfree) rr=fmpz_mpoly_factor_squarefree(factors,Flint_f,ctx);
710 else rr=fmpz_mpoly_factor(factors,Flint_f,ctx);
711 if (rr==0) printf("fail\n");
712 fmpz_mpoly_t fac;
713 fmpz_mpoly_init(fac,ctx);
714 CanonicalForm cf_fac;
715 int cf_exp;
716 fmpz_t c;
717 fmpz_init(c);
718 fmpz_mpoly_factor_get_constant_fmpz(c,factors,ctx);
719 cf_fac=convertFmpz2CF(c);
720 fmpz_clear(c);
721 F.append(CFFactor(cf_fac,1));
722 for(int i=fmpz_mpoly_factor_length(factors,ctx)-1; i>=0; i--)
723 {
724 fmpz_mpoly_factor_get_base(fac,factors,i,ctx);
725 cf_fac=convFlintMPFactoryP(fac,ctx,f.level());
726 cf_exp=fmpz_mpoly_factor_get_exp_si(factors,i,ctx);
727 F.append(CFFactor(cf_fac,cf_exp));
728 }
729 fmpz_mpoly_factor_clear(factors,ctx);
730 fmpz_mpoly_clear(Flint_f,ctx);
731 fmpz_mpoly_ctx_clear(ctx);
732 goto end_char0;
733 }
734 #endif
735 #if defined(HAVE_NTL)
736 On (SW_RATIONAL);
737 if (issqrfree)
738 {
739 CFList factors= ratSqrfFactorize (fz);
740 for (CFListIterator i= factors; i.hasItem(); i++)
741 F.append (CFFactor (i.getItem(), 1));
742 }
743 else
744 {
745 F = ratFactorize (fz);
746 }
747 #endif
748 #if !defined(HAVE_FLINT) || (__FLINT_RELEASE < 20700)
749 #ifndef HAVE_NTL
750 F=ZFactorizeMultivariate(fz, issqrfree);
751 #endif
752 #endif
753 }
754
755end_char0:
756 if ( on_rational )
758 else
760 if ( ! cd.isOne() )
761 {
762 CFFactor new_first( F.getFirst().factor() / cd );
763 F.removeFirst();
764 F.insert( new_first );
765 }
766 }
767
768#if defined(HAVE_NTL)
769end_charp:
770#endif
772 return F;
773}
CanonicalForm convertFmpz2CF(const fmpz_t coefficient)
conversion of a FLINT integer to CanonicalForm
CFFList convertFLINTnmod_poly_factor2FacCFFList(const nmod_poly_factor_t fac, const mp_limb_t leadingCoeff, const Variable &x)
conversion of a FLINT factorization over Z/p (for word size p) to a CFFList
void convertFacCF2Fmpz_poly_t(fmpz_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomial over Z to a fmpz_poly_t
CFFList convertFLINTfmpz_poly_factor2FacCFFList(const fmpz_poly_factor_t fac, const Variable &x)
conversion of a FLINT factorization over Z to a CFFList
CFFList convertNTLvec_pair_GF2X_long2FacCFFList(const vec_pair_GF2X_long &e, GF2, const Variable &x)
NAME: convertNTLvec_pair_GF2X_long2FacCFFList.
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
CFFList convertNTLvec_pair_zzpX_long2FacCFFList(const vec_pair_zz_pX_long &e, const zz_p cont, const Variable &x)
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
VAR long fac_NTL_char
Definition NTLconvert.cc:46
GF2X convertFacCF2NTLGF2X(const CanonicalForm &f)
NAME: convertFacCF2NTLGF2X.
CFFList convertNTLvec_pair_ZZX_long2FacCFFList(const vec_pair_ZZX_long &e, const ZZ &cont, const Variable &x)
NAME: convertNTLvec_pair_ZZX_long2FacCFFList.
CanonicalForm FACTORY_PUBLIC icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition cf_gcd.cc:74
int degree(const CanonicalForm &f)
bool hasFirstAlgVar(const CanonicalForm &f, Variable &a)
check if poly f contains an algebraic variable a
Definition cf_ops.cc:679
ListIterator< CFFactor > CFFListIterator
ListIterator< CanonicalForm > CFListIterator
List< CFFactor > CFFList
Factor< CanonicalForm > CFFactor
List< CanonicalForm > CFList
CanonicalForm cd(bCommonDen(FF))
Definition cfModGcd.cc:4097
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
EXTERN_VAR int singular_homog_flag
static const int SW_USE_FL_GCD_P
set to 1 to use Flints gcd over F_p
Definition cf_defs.h:47
static const int SW_USE_NTL_SORT
set to 1 to sort factors in a factorization
Definition cf_defs.h:39
static const int SW_USE_FL_FAC_0
set to 1 to prefer flints multivariate factorization over Z/p
Definition cf_defs.h:57
static const int SW_USE_FL_FAC_P
set to 1 to prefer flints multivariate factorization over Z/p
Definition cf_defs.h:55
static const int SW_BERLEKAMP
set to 1 to use Factorys Berlekamp alg.
Definition cf_defs.h:51
#define GaloisFieldDomain
Definition cf_defs.h:18
Variable get_max_degree_Variable(const CanonicalForm &f)
get_max_degree_Variable returns Variable with highest degree.
Definition cf_factor.cc:264
int cmpCF(const CFFactor &f, const CFFactor &g)
Definition cf_factor.cc:398
CanonicalForm homogenize(const CanonicalForm &f, const Variable &x)
homogenize homogenizes f with Variable x
Definition cf_factor.cc:317
CFFList factorize(const CanonicalForm &f, bool issqrfree)
factorization over or
Definition cf_factor.cc:409
CanonicalForm compress(const CanonicalForm &f, CFMap &m)
CanonicalForm compress ( const CanonicalForm & f, CFMap & m )
Definition cf_map.cc:210
VAR void(* factoryError)(const char *s)
Definition cf_util.cc:80
static int gettype()
Definition cf_factory.h:28
class CFMap
Definition cf_map.h:85
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
CF_NO_INLINE bool isOne() const
int level() const
level() returns the level of CO.
T factor() const
T getFirst() const
void removeFirst()
void sort(int(*)(const T &, const T &))
void append(const T &)
void insert(const T &)
factory's class for variables
Definition factory.h:127
Variable alpha
CanonicalForm factor
Definition facAbsFact.cc:97
CFFList ratFactorize(const CanonicalForm &G, const Variable &v=Variable(1), bool substCheck=true)
factorize a multivariate polynomial over
CFList ratSqrfFactorize(const CanonicalForm &G, const Variable &v=Variable(1))
factorize a squarefree multivariate polynomial over
CFList FpSqrfFactorize(const CanonicalForm &F)
factorize a squarefree multivariate polynomial over
CFFList FpFactorize(const CanonicalForm &G, bool substCheck=true)
factorize a multivariate polynomial over
CFFList GFFactorize(const CanonicalForm &G, bool substCheck=true)
factorize a multivariate polynomial over GF
CFList GFSqrfFactorize(const CanonicalForm &F)
factorize a squarefree multivariate polynomial over GF
convertFacCF2nmod_poly_t(FLINTmipo, M)
nmod_poly_clear(FLINTmipo)
CFFList FpFactorizeUnivariateB(const CanonicalForm &f, bool issqrfree=false)
CFFList FpFactorizeUnivariateCZ(const CanonicalForm &f, bool issqrfree, int numext, const Variable alpha, const Variable beta)
CFFList ZFactorizeMultivariate(const CanonicalForm &f, bool issqrfree)
CFFList ZFactorizeUnivariate(const CanonicalForm &ff, bool issqrfree=false)
VAR int xn
Definition walk.cc:4509

◆ factorize() [2/2]

CFFList FACTORY_PUBLIC factorize ( const CanonicalForm & f,
const Variable & alpha )

factorization over $ F_p(\alpha) $ or $ Q(\alpha) $

Definition at line 778 of file cf_factor.cc.

779{
780 if ( f.inCoeffDomain() )
781 return CFFList( f );
782 //out_cf("factorize:",f,"==================================\n");
783 //out_cf("mipo:",getMipo(alpha),"\n");
784
785 CFFList F;
786 ASSERT( alpha.level() < 0 && getReduce (alpha), "not an algebraic extension" );
787#ifndef NOASSERT
789 if (hasFirstAlgVar(f, beta))
790 ASSERT (beta == alpha, "f has an algebraic variable that \
791 does not coincide with alpha");
792#endif
793 int ch=getCharacteristic();
794 if (ch>0)
795 {
796 if (f.isUnivariate())
797 {
798#ifdef HAVE_NTL
799 if (/*getCharacteristic()*/ch==2)
800 {
801 // special case : GF2
802
803 // remainder is two ==> nothing to do
804
805 // set minimal polynomial in NTL using the optimized conversion routines for characteristic 2
806 GF2X minPo=convertFacCF2NTLGF2X(getMipo(alpha,f.mvar()));
807 GF2E::init (minPo);
808
809 // convert to NTL again using the faster conversion routines
810 GF2EX f1;
811 if (isPurePoly(f))
812 {
813 GF2X f_tmp=convertFacCF2NTLGF2X(f);
814 f1=to_GF2EX(f_tmp);
815 }
816 else
817 f1=convertFacCF2NTLGF2EX(f,minPo);
818
819 // make monic (in Z/2(a))
820 GF2E f1_coef=LeadCoeff(f1);
821 MakeMonic(f1);
822
823 // factorize using NTL
824 vec_pair_GF2EX_long factors;
825 CanZass(factors,f1);
826
827 // return converted result
828 F=convertNTLvec_pair_GF2EX_long2FacCFFList(factors,f1_coef,f.mvar(),alpha);
830 return F;
831 }
832#endif
833#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
834 {
835 // use FLINT
836 nmod_poly_t FLINTmipo, leadingCoeff;
837 fq_nmod_ctx_t fq_con;
838
839 nmod_poly_init (FLINTmipo, ch);
840 nmod_poly_init (leadingCoeff, ch);
842
843 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
844 fq_nmod_poly_t FLINTF;
846 fq_nmod_poly_factor_t res;
847 fq_nmod_poly_factor_init (res, fq_con);
848 fq_nmod_poly_factor (res, leadingCoeff, FLINTF, fq_con);
850 F.insert (CFFactor (Lc (f), 1));
851
852 fq_nmod_poly_factor_clear (res, fq_con);
853 fq_nmod_poly_clear (FLINTF, fq_con);
854 nmod_poly_clear (FLINTmipo);
855 nmod_poly_clear (leadingCoeff);
858 return F;
859 }
860#endif
861#ifdef HAVE_NTL
862 {
863 // use NTL
864 if (fac_NTL_char != ch)
865 {
866 fac_NTL_char = ch;
867 zz_p::init(ch);
868 }
869
870 // set minimal polynomial in NTL
871 zz_pX minPo=convertFacCF2NTLzzpX(getMipo(alpha));
872 zz_pE::init (minPo);
873
874 // convert to NTL
875 zz_pEX f1=convertFacCF2NTLzz_pEX(f,minPo);
876 zz_pE leadcoeff= LeadCoeff(f1);
877
878 //make monic
879 f1=f1 / leadcoeff; //leadcoeff==LeadCoeff(f1);
880
881 // factorize
882 vec_pair_zz_pEX_long factors;
883 CanZass(factors,f1);
884
885 // return converted result
886 F=convertNTLvec_pair_zzpEX_long2FacCFFList(factors,leadcoeff,f.mvar(),alpha);
887 //test_cff(F,f);
889 return F;
890 }
891#endif
892#if !defined(HAVE_NTL) && !defined(HAVE_FLINT)
893 // char p, extension, univariate
894 CanonicalForm c=Lc(f);
895 CanonicalForm fc=f/c;
896 F=FpFactorizeUnivariateCZ( fc, false, 1, alpha, Variable() );
897 F.insert (CFFactor (c, 1));
898#endif
899 }
900 else // char p, multivariate
901 {
902 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
903 // use FLINT
904 nmod_poly_t FLINTmipo;
905 fq_nmod_ctx_t fq_con;
906 fq_nmod_mpoly_ctx_t ctx;
907
908 nmod_poly_init (FLINTmipo, ch);
910
911 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
912 fq_nmod_mpoly_ctx_init(ctx,f.level(),ORD_LEX,fq_con);
913
914 fq_nmod_mpoly_t FLINTF;
915 fq_nmod_mpoly_init(FLINTF,ctx);
916 convertFacCF2Fq_nmod_mpoly_t(FLINTF,f,ctx,f.level(),fq_con);
917 fq_nmod_mpoly_factor_t res;
918 fq_nmod_mpoly_factor_init (res, ctx);
919 fq_nmod_mpoly_factor (res, FLINTF, ctx);
920 F= convertFLINTFq_nmod_mpoly_factor2FacCFFList (res, ctx,f.level(),fq_con,alpha);
921 //F.insert (CFFactor (Lc (f), 1));
922
923 fq_nmod_mpoly_factor_clear (res, ctx);
924 fq_nmod_mpoly_clear (FLINTF, ctx);
925 nmod_poly_clear (FLINTmipo);
926 fq_nmod_mpoly_ctx_clear (ctx);
929 return F;
930 #elif defined(HAVE_NTL)
931 F= FqFactorize (f, alpha);
932 #else
933 factoryError ("multivariate factorization over Z/pZ(alpha) depends on NTL/Flint(missing)");
934 return CFFList (CFFactor (f, 1));
935 #endif
936 }
937 }
938 else // Q(a)[x]
939 {
940 if (f.isUnivariate())
941 {
942 F= AlgExtFactorize (f, alpha);
943 }
944 else //Q(a)[x1,...,xn]
945 {
946 #if defined(HAVE_NTL) || defined(HAVE_FLINT)
947 F= ratFactorize (f, alpha);
948 #else
949 factoryError ("multivariate factorization over Q(alpha) depends on NTL or FLINT (missing)");
950 return CFFList (CFFactor (f, 1));
951 #endif
952 }
953 }
955 return F;
956}
CFFList convertFLINTFq_nmod_poly_factor2FacCFFList(const fq_nmod_poly_factor_t fac, const Variable &x, const Variable &alpha, const fq_nmod_ctx_t fq_con)
conversion of a FLINT factorization over Fq (for word size p) to a CFFList
void convertFacCF2Fq_nmod_poly_t(fq_nmod_poly_t result, const CanonicalForm &f, const fq_nmod_ctx_t ctx)
conversion of a factory univariate poly over F_q to a FLINT fq_nmod_poly_t
zz_pEX convertFacCF2NTLzz_pEX(const CanonicalForm &f, const zz_pX &mipo)
CFFList convertNTLvec_pair_zzpEX_long2FacCFFList(const vec_pair_zz_pEX_long &e, const zz_pE &cont, const Variable &x, const Variable &alpha)
CFFList convertNTLvec_pair_GF2EX_long2FacCFFList(const vec_pair_GF2EX_long &e, const GF2E &cont, const Variable &x, const Variable &alpha)
NAME: convertNTLvec_pair_GF2EX_long2FacCFFList.
GF2EX convertFacCF2NTLGF2EX(const CanonicalForm &f, const GF2X &mipo)
CanonicalForm in Z_2(a)[X] to NTL GF2EX.
CanonicalForm Lc(const CanonicalForm &f)
bool isPurePoly(const CanonicalForm &f)
Definition cf_factor.cc:248
CanonicalForm res
Definition facAbsFact.cc:60
Variable beta
Definition facAbsFact.cc:95
CFFList AlgExtFactorize(const CanonicalForm &F, const Variable &alpha)
factorize a univariate polynomial over algebraic extension of Q
Definition facAlgExt.cc:370
CFFList FqFactorize(const CanonicalForm &G, const Variable &alpha, bool substCheck=true)
factorize a multivariate polynomial over
fq_nmod_ctx_t fq_con
Definition facHensel.cc:99
fq_nmod_ctx_clear(fq_con)
nmod_poly_init(FLINTmipo, getCharacteristic())
fq_nmod_ctx_init_modulus(fq_con, FLINTmipo, "Z")
fq_nmod_poly_clear(prod, fq_con)
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition variable.cc:207
bool getReduce(const Variable &alpha)
Definition variable.cc:232

◆ Farey()

Farey rational reconstruction.

If NTL is available it uses the fast algorithm from NTL, i.e. Encarnacion, Collins.

Definition at line 202 of file cf_chinese.cc.

203{
204 int is_rat=isOn(SW_RATIONAL);
206 Variable x = f.mvar();
210#ifdef HAVE_FLINT
211 fmpz_t FLINTq;
212 fmpz_init(FLINTq);
213 convertCF2initFmpz(FLINTq,q);
214 fmpz_t FLINTc;
215 fmpz_init(FLINTc);
216 fmpq_t FLINTres;
217 fmpq_init(FLINTres);
218#elif defined(HAVE_NTL)
219 ZZ NTLq= convertFacCF2NTLZZ (q);
220 ZZ bound;
221 SqrRoot (bound, NTLq/2);
222#else
223 factoryError("NTL/FLINT missing:Farey");
224#endif
225 for ( i = f; i.hasTerms(); i++ )
226 {
227 c = i.coeff();
228 if ( c.inCoeffDomain())
229 {
230#ifdef HAVE_FLINT
231 if (c.inZ())
232 {
233 convertCF2initFmpz(FLINTc,c);
234 fmpq_reconstruct_fmpz(FLINTres,FLINTc,FLINTq);
235 result += power (x, i.exp())*(convertFmpq2CF(FLINTres));
236 }
237#elif defined(HAVE_NTL)
238 if (c.inZ())
239 {
240 ZZ NTLc= convertFacCF2NTLZZ (c);
241 bool lessZero= (sign (NTLc) == -1);
242 if (lessZero)
243 NTL::negate (NTLc, NTLc);
244 ZZ NTLnum, NTLden;
245 if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
246 {
247 if (lessZero)
248 NTL::negate (NTLnum, NTLnum);
249 CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
250 CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
251 On (SW_RATIONAL);
252 result += power (x, i.exp())*(num/den);
254 }
255 }
256#else
257 if (c.inZ())
258 result += power (x, i.exp()) * Farey_n(c,q);
259#endif
260 else
261 result += power( x, i.exp() ) * Farey(c,q);
262 }
263 else
264 result += power( x, i.exp() ) * Farey(c,q);
265 }
266 if (is_rat) On(SW_RATIONAL);
267#ifdef HAVE_FLINT
268 fmpq_clear(FLINTres);
269 fmpz_clear(FLINTc);
270 fmpz_clear(FLINTq);
271#endif
272 return result;
273}
CanonicalForm convertFmpq2CF(const fmpq_t q)
conversion of a FLINT rational to CanonicalForm
void convertCF2initFmpz(fmpz_t result, const CanonicalForm &f)
conversion of a factory integer to fmpz_t(init.)
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
static CanonicalForm bound(const CFMatrix &M)
Definition cf_linsys.cc:460
bool inCoeffDomain() const
bool inZ() const
predicates

◆ fdivides() [1/2]

bool fdivides ( const CanonicalForm & f,
const CanonicalForm & g )

bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )

fdivides() - check whether ‘f’ divides ‘g’.

Returns true iff ‘f’ divides ‘g’. Uses some extra heuristic to avoid polynomial division. Without the heuristic, the test essentialy looks like ‘divremt(g, f, q, r) && r.isZero()’.

Type info:

f, g: Current

Elements from prime power domains (or polynomials over such domains) are admissible if ‘f’ (or lc(‘f’), resp.) is not a zero divisor. This is a slightly stronger precondition than mathematically necessary since divisibility is a well-defined notion in arbitrary rings. Hence, we decided not to declare the weaker type ‘CurrentPP’.

Developers note:

One may consider the the test ‘fdivides( f.LC(), g.LC() )’ in the main ‘if’-test superfluous since ‘divremt()’ in the ‘if’-body repeats the test. However, ‘divremt()’ does not use any heuristic to do so.

It seems not reasonable to call ‘fdivides()’ from ‘divremt()’ to check divisibility of leading coefficients. ‘fdivides()’ is on a relatively high level compared to ‘divremt()’.

Definition at line 340 of file cf_algorithm.cc.

341{
342 // trivial cases
343 if ( g.isZero() )
344 return true;
345 else if ( f.isZero() )
346 return false;
347
348 if ( (f.inCoeffDomain() || g.inCoeffDomain())
349 && ((getCharacteristic() == 0 && isOn( SW_RATIONAL ))
350 || (getCharacteristic() > 0) ))
351 {
352 // if we are in a field all elements not equal to zero are units
353 if ( f.inCoeffDomain() )
354 return true;
355 else
356 // g.inCoeffDomain()
357 return false;
358 }
359
360 // we may assume now that both levels either equal LEVELBASE
361 // or are greater zero
362 int fLevel = f.level();
363 int gLevel = g.level();
364 if ( (gLevel > 0) && (fLevel == gLevel) )
365 // f and g are polynomials in the same main variable
366 if ( degree( f ) <= degree( g )
367 && fdivides( f.tailcoeff(), g.tailcoeff() )
368 && fdivides( f.LC(), g.LC() ) )
369 {
370 CanonicalForm q, r;
371 return divremt( g, f, q, r ) && r.isZero();
372 }
373 else
374 return false;
375 else if ( gLevel < fLevel )
376 // g is a coefficient w.r.t. f
377 return false;
378 else
379 {
380 // either f is a coefficient w.r.t. polynomial g or both
381 // f and g are from a base domain (should be Z or Z/p^n,
382 // then)
383 CanonicalForm q, r;
384 return divremt( g, f, q, r ) && r.isZero();
385 }
386}
bool divremt(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
g
Definition cfModGcd.cc:4098
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )

◆ fdivides() [2/2]

bool fdivides ( const CanonicalForm & f,
const CanonicalForm & g,
CanonicalForm & quot )

same as fdivides if true returns quotient quot of g by f otherwise quot == 0

Definition at line 390 of file cf_algorithm.cc.

391{
392 quot= 0;
393 // trivial cases
394 if ( g.isZero() )
395 return true;
396 else if ( f.isZero() )
397 return false;
398
399 if ( (f.inCoeffDomain() || g.inCoeffDomain())
400 && ((getCharacteristic() == 0 && isOn( SW_RATIONAL ))
401 || (getCharacteristic() > 0) ))
402 {
403 // if we are in a field all elements not equal to zero are units
404 if ( f.inCoeffDomain() )
405 {
406 quot= g/f;
407 return true;
408 }
409 else
410 // g.inCoeffDomain()
411 return false;
412 }
413
414 // we may assume now that both levels either equal LEVELBASE
415 // or are greater zero
416 int fLevel = f.level();
417 int gLevel = g.level();
418 if ( (gLevel > 0) && (fLevel == gLevel) )
419 // f and g are polynomials in the same main variable
420 if ( degree( f ) <= degree( g )
421 && fdivides( f.tailcoeff(), g.tailcoeff() )
422 && fdivides( f.LC(), g.LC() ) )
423 {
424 CanonicalForm q, r;
425 if (divremt( g, f, q, r ) && r.isZero())
426 {
427 quot= q;
428 return true;
429 }
430 else
431 return false;
432 }
433 else
434 return false;
435 else if ( gLevel < fLevel )
436 // g is a coefficient w.r.t. f
437 return false;
438 else
439 {
440 // either f is a coefficient w.r.t. polynomial g or both
441 // f and g are from a base domain (should be Z or Z/p^n,
442 // then)
443 CanonicalForm q, r;
444 if (divremt( g, f, q, r ) && r.isZero())
445 {
446 quot= q;
447 return true;
448 }
449 else
450 return false;
451 }
452}

◆ get_max_degree_Variable()

Variable get_max_degree_Variable ( const CanonicalForm & f)

get_max_degree_Variable returns Variable with highest degree.

We assume f is not a constant!

Definition at line 264 of file cf_factor.cc.

265{
266 ASSERT( ( ! f.inCoeffDomain() ), "no constants" );
267 int max=0, maxlevel=0, n=level(f);
268 for ( int i=1; i<=n; i++ )
269 {
270 if (degree(f,Variable(i)) >= max)
271 {
272 max= degree(f,Variable(i)); maxlevel= i;
273 }
274 }
275 return Variable(maxlevel);
276}
int level(const CanonicalForm &f)
static int max(int a, int b)
Definition fast_mult.cc:264

◆ get_Terms()

CFList get_Terms ( const CanonicalForm & f)

Definition at line 293 of file cf_factor.cc.

293 {
294 CFList result,dummy,dummy2;
297
298 if ( getNumVars(f) == 0 ) result.append(f);
299 else{
300 Variable _x(level(f));
301 for ( i=f; i.hasTerms(); i++ ){
302 getTerms(i.coeff(), 1, dummy);
303 for ( j=dummy; j.hasItem(); j++ )
304 result.append(j.getItem() * power(_x, i.exp()));
305
306 dummy= dummy2; // have to initalize new
307 }
308 }
309 return result;
310}
int getNumVars(const CanonicalForm &f)
int getNumVars ( const CanonicalForm & f )
Definition cf_ops.cc:314
void getTerms(const CanonicalForm &f, const CanonicalForm &t, CFList &result)
get_Terms: Split the polynomial in the containing terms.
Definition cf_factor.cc:283

◆ getTerms()

void getTerms ( const CanonicalForm & f,
const CanonicalForm & t,
CFList & result )

get_Terms: Split the polynomial in the containing terms.

getTerms: the real work is done here.

Definition at line 283 of file cf_factor.cc.

284{
285 if ( getNumVars(f) == 0 ) result.append(f*t);
286 else{
287 Variable x(level(f));
288 for ( CFIterator i=f; i.hasTerms(); i++ )
289 getTerms( i.coeff(), t*power(x,i.exp()), result);
290 }
291}

◆ homogenize() [1/2]

CanonicalForm homogenize ( const CanonicalForm & f,
const Variable & x )

homogenize homogenizes f with Variable x

Definition at line 317 of file cf_factor.cc.

318{
319#if 0
320 int maxdeg=totaldegree(f), deg;
322 CanonicalForm elem, result(0);
323
324 for (i=f; i.hasTerms(); i++)
325 {
326 elem= i.coeff()*power(f.mvar(),i.exp());
327 deg = totaldegree(elem);
328 if ( deg < maxdeg )
329 result += elem * power(x,maxdeg-deg);
330 else
331 result+=elem;
332 }
333 return result;
334#else
335 CFList Newlist, Termlist= get_Terms(f);
336 int maxdeg=totaldegree(f), deg;
338 CanonicalForm elem, result(0);
339
340 for (i=Termlist; i.hasItem(); i++)
341 {
342 elem= i.getItem();
343 deg = totaldegree(elem);
344 if ( deg < maxdeg )
345 Newlist.append(elem * power(x,maxdeg-deg));
346 else
347 Newlist.append(elem);
348 }
349 for (i=Newlist; i.hasItem(); i++) // rebuild
350 result += i.getItem();
351
352 return result;
353#endif
354}
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition cf_ops.cc:523
CFList get_Terms(const CanonicalForm &f)
Definition cf_factor.cc:293

◆ homogenize() [2/2]

CanonicalForm homogenize ( const CanonicalForm & f,
const Variable & x,
const Variable & v1,
const Variable & v2 )

Definition at line 357 of file cf_factor.cc.

358{
359#if 0
360 int maxdeg=totaldegree(f), deg;
362 CanonicalForm elem, result(0);
363
364 for (i=f; i.hasTerms(); i++)
365 {
366 elem= i.coeff()*power(f.mvar(),i.exp());
367 deg = totaldegree(elem);
368 if ( deg < maxdeg )
369 result += elem * power(x,maxdeg-deg);
370 else
371 result+=elem;
372 }
373 return result;
374#else
375 CFList Newlist, Termlist= get_Terms(f);
376 int maxdeg=totaldegree(f), deg;
378 CanonicalForm elem, result(0);
379
380 for (i=Termlist; i.hasItem(); i++)
381 {
382 elem= i.getItem();
383 deg = totaldegree(elem,v1,v2);
384 if ( deg < maxdeg )
385 Newlist.append(elem * power(x,maxdeg-deg));
386 else
387 Newlist.append(elem);
388 }
389 for (i=Newlist; i.hasItem(); i++) // rebuild
390 result += i.getItem();
391
392 return result;
393#endif
394}

◆ isPurePoly()

bool isPurePoly ( const CanonicalForm & f)

Definition at line 248 of file cf_factor.cc.

249{
250 if (f.level()<=0) return false;
251 for (CFIterator i=f;i.hasTerms();i++)
252 {
253 if (!(i.coeff().inBaseDomain())) return false;
254 }
255 return true;
256}

◆ isPurePoly_m()

bool isPurePoly_m ( const CanonicalForm & f)

Definition at line 238 of file cf_factor.cc.

239{
240 if (f.inBaseDomain()) return true;
241 if (f.level()<0) return false;
242 for (CFIterator i=f;i.hasTerms();i++)
243 {
244 if (!isPurePoly_m(i.coeff())) return false;
245 }
246 return true;
247}
bool isPurePoly_m(const CanonicalForm &f)
Definition cf_factor.cc:238

◆ linearSystemSolve()

bool linearSystemSolve ( CFMatrix & M)

Definition at line 78 of file cf_linsys.cc.

79{
80 typedef int* int_ptr;
81
82 if ( ! matrix_in_Z( M ) ) {
83 int nrows = M.rows(), ncols = M.columns();
84 int i, j, k;
85 CanonicalForm rowpivot, pivotrecip;
86 // triangularization
87 for ( i = 1; i <= nrows; i++ ) {
88 //find "pivot"
89 for (j = i; j <= nrows; j++ )
90 if ( M(j,i) != 0 ) break;
91 if ( j > nrows ) return false;
92 if ( j != i )
93 M.swapRow( i, j );
94 pivotrecip = 1 / M(i,i);
95 for ( j = 1; j <= ncols; j++ )
96 M(i,j) *= pivotrecip;
97 for ( j = i+1; j <= nrows; j++ ) {
98 rowpivot = M(j,i);
99 if ( rowpivot == 0 ) continue;
100 for ( k = i; k <= ncols; k++ )
101 M(j,k) -= M(i,k) * rowpivot;
102 }
103 }
104 // matrix is now upper triangular with 1s down the diagonal
105 // back-substitute
106 for ( i = nrows-1; i > 0; i-- ) {
107 for ( j = nrows+1; j <= ncols; j++ ) {
108 for ( k = i+1; k <= nrows; k++ )
109 M(i,j) -= M(k,j) * M(i,k);
110 }
111 }
112 return true;
113 }
114 else {
115 int rows = M.rows(), cols = M.columns();
116 CFMatrix MM( rows, cols );
117 int ** mm = new int_ptr[rows];
118 CanonicalForm Q, Qhalf, mnew, qnew, B;
119 int i, j, p, pno;
120 bool ok;
121
122 // initialize room to hold the result and the result mod p
123 for ( i = 0; i < rows; i++ ) {
124 mm[i] = new int[cols];
125 }
126
127 // calculate the bound for the result
128 B = bound( M );
129 DEBOUTLN( cerr, "bound = " << B );
130
131 // find a first solution mod p
132 pno = 0;
133 do {
134 DEBOUTSL( cerr );
135 DEBOUT( cerr, "trying prime(" << pno << ") = " );
136 p = cf_getBigPrime( pno );
137 DEBOUT( cerr, p );
138 DEBOUTENDL( cerr );
140 // map matrix into char p
141 for ( i = 1; i <= rows; i++ )
142 for ( j = 1; j <= cols; j++ )
143 mm[i-1][j-1] = mapinto( M(i,j) ).intval();
144 // solve mod p
145 ok = solve( mm, rows, cols );
146 pno++;
147 } while ( ! ok );
148
149 // initialize the result matrix with first solution
151 for ( i = 1; i <= rows; i++ )
152 for ( j = rows+1; j <= cols; j++ )
153 MM(i,j) = mm[i-1][j-1];
154
155 // Q so far
156 Q = p;
157 while ( Q < B && pno < cf_getNumBigPrimes() ) {
158 do {
159 DEBOUTSL( cerr );
160 DEBOUT( cerr, "trying prime(" << pno << ") = " );
161 p = cf_getBigPrime( pno );
162 DEBOUT( cerr, p );
163 DEBOUTENDL( cerr );
165 for ( i = 1; i <= rows; i++ )
166 for ( j = 1; j <= cols; j++ )
167 mm[i-1][j-1] = mapinto( M(i,j) ).intval();
168 // solve mod p
169 ok = solve( mm, rows, cols );
170 pno++;
171 } while ( ! ok );
172 // found a solution mod p
173 // now chinese remainder it to a solution mod Q*p
175 for ( i = 1; i <= rows; i++ )
176 for ( j = rows+1; j <= cols; j++ )
177 {
178 chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i-1][j-1]), CanonicalForm(p), mnew, qnew );
179 MM(i, j) = mnew;
180 }
181 Q = qnew;
182 }
183 if ( pno == cf_getNumBigPrimes() )
184 fuzzy_result = true;
185 else
186 fuzzy_result = false;
187 // store the result in M
188 Qhalf = Q / 2;
189 for ( i = 1; i <= rows; i++ ) {
190 for ( j = rows+1; j <= cols; j++ )
191 if ( MM(i,j) > Qhalf )
192 M(i,j) = MM(i,j) - Q;
193 else
194 M(i,j) = MM(i,j);
195 delete [] mm[i-1];
196 }
197 delete [] mm;
198 return ! fuzzy_result;
199 }
200}
CanonicalForm mapinto(const CanonicalForm &f)
int int ncols
Definition cf_linsys.cc:32
VAR bool fuzzy_result
Definition cf_linsys.cc:75
int nrows
Definition cf_linsys.cc:32
bool solve(int **extmat, int nrows, int ncols)
Definition cf_linsys.cc:504
long intval() const
conversion functions
#define DEBOUTENDL(stream)
Definition debug.h:48
#define DEBOUTSL(stream)
Definition debug.h:46

◆ maxNorm()

CanonicalForm maxNorm ( const CanonicalForm & f )

maxNorm() - return maximum norm of ‘f’.

That is, the base coefficient of ‘f’ with the largest absolute value.

Valid for arbitrary polynomials over arbitrary domains, but most useful for multivariate polynomials over Z.

Type info:

f: CurrentPP

Definition at line 536 of file cf_algorithm.cc.

537{
538 if ( f.inBaseDomain() )
539 return abs( f );
540 else {
542 for ( CFIterator i = f; i.hasTerms(); i++ ) {
543 CanonicalForm coeffMaxNorm = maxNorm( i.coeff() );
544 if ( coeffMaxNorm > result )
545 result = coeffMaxNorm;
546 }
547 return result;
548 }
549}
Rational abs(const Rational &a)
Definition GMPrat.cc:436
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )

◆ psq()

CanonicalForm psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )

psq() - return pseudo quotient of ‘f’ and ‘g’ with respect to ‘x’.

‘g’ must not equal zero.

Type info:

f, g: Current x: Polynomial

Developers note:

This is not an optimal implementation. Better would have been an implementation in ‘InternalPoly’ avoiding the exponentiation of the leading coefficient of ‘g’. It seemed not worth to do so.

See also
psr(), psqr()

Definition at line 172 of file cf_algorithm.cc.

173{
174 ASSERT( x.level() > 0, "type error: polynomial variable expected" );
175 ASSERT( ! g.isZero(), "math error: division by zero" );
176
177 // swap variables such that x's level is larger or equal
178 // than both f's and g's levels.
179 Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
180 CanonicalForm F = swapvar( f, x, X );
181 CanonicalForm G = swapvar( g, x, X );
182
183 // now, we have to calculate the pseudo remainder of F and G
184 // w.r.t. X
185 int fDegree = degree( F, X );
186 int gDegree = degree( G, X );
187 if ( fDegree < 0 || fDegree < gDegree )
188 return 0;
189 else {
190 CanonicalForm result = (power( LC( G, X ), fDegree-gDegree+1 ) * F) / G;
191 return swapvar( result, x, X );
192 }
193}
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition cf_ops.cc:168
CanonicalForm LC(const CanonicalForm &f)
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
STATIC_VAR TreeM * G
Definition janet.cc:31

◆ psqr()

void psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable & x )

psqr() - calculate pseudo quotient and remainder of ‘f’ and ‘g’ with respect to ‘x’.

Returns the pseudo quotient of ‘f’ and ‘g’ in ‘q’, the pseudo remainder in ‘r’. ‘g’ must not equal zero.

See ‘psr()’ for more detailed information.

Type info:

f, g: Current q, r: Anything x: Polynomial

Developers note:

This is not an optimal implementation. Better would have been an implementation in ‘InternalPoly’ avoiding the exponentiation of the leading coefficient of ‘g’. It seemed not worth to do so.

See also
psr(), psq()

Definition at line 223 of file cf_algorithm.cc.

224{
225 ASSERT( x.level() > 0, "type error: polynomial variable expected" );
226 ASSERT( ! g.isZero(), "math error: division by zero" );
227
228 // swap variables such that x's level is larger or equal
229 // than both f's and g's levels.
230 Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
231 CanonicalForm F = swapvar( f, x, X );
232 CanonicalForm G = swapvar( g, x, X );
233
234 // now, we have to calculate the pseudo remainder of F and G
235 // w.r.t. X
236 int fDegree = degree( F, X );
237 int gDegree = degree( G, X );
238 if ( fDegree < 0 || fDegree < gDegree ) {
239 q = 0; r = f;
240 } else {
241 divrem( power( LC( G, X ), fDegree-gDegree+1 ) * F, G, q, r );
242 q = swapvar( q, x, X );
243 r = swapvar( r, x, X );
244 }
245}
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)

◆ psr()

CanonicalForm psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )

psr() - return pseudo remainder of ‘f’ and ‘g’ with respect to ‘x’.

‘g’ must not equal zero.

For f and g in R[x], R an arbitrary ring, g != 0, there is a representation

LC(g)^s*f = g*q + r

with r = 0 or deg(r) < deg(g) and s = 0 if f = 0 or s = max( 0, deg(f)-deg(g)+1 ) otherwise. r = psr(f, g) and q = psq(f, g) are called "pseudo remainder" and "pseudo quotient", resp. They are uniquely determined if LC(g) is not a zero divisor in R.

See H.-J. Reiffen/G. Scheja/U. Vetter - "Algebra", 2nd ed., par. 15, for a reference.

Type info:

f, g: Current x: Polynomial

Polynomials over prime power domains are admissible if lc(LC(‘g’,‘x’)) is not a zero divisor. This is a slightly stronger precondition than mathematically necessary since pseudo remainder and quotient are well-defined if LC(‘g’,‘x’) is not a zero divisor.

For example, psr(y^2, (13*x+1)*y) is well-defined in (Z/13^2[x])[y] since (13*x+1) is not a zero divisor. But calculating it with Factory would fail since 13 is a zero divisor in Z/13^2.

Due to this inconsistency with mathematical notion, we decided not to declare type ‘CurrentPP’ for ‘f’ and ‘g’.

Developers note:

This is not an optimal implementation. Better would have been an implementation in ‘InternalPoly’ avoiding the exponentiation of the leading coefficient of ‘g’. In contrast to ‘psq()’ and ‘psqr()’ it definitely seems worth to implement the pseudo remainder on the internal level.

See also
psq(), psqr()

Definition at line 117 of file cf_algorithm.cc.

118{
119 CanonicalForm r=rr, v=vv, l, test, lu, lv, t, retvalue;
120 int dr, dv, d,n=0;
121
122
123 dr = degree( r, x );
124 if (dr>0)
125 {
126 dv = degree( v, x );
127 if (dv <= dr) {l=LC(v,x); v = v -l*power(x,dv);}
128 else { l = 1; }
129 d= dr-dv+1;
130 //out_cf("psr(",rr," ");
131 //out_cf("",vv," ");
132 //printf(" var=%d\n",x.level());
133 while ( ( dv <= dr ) && ( !r.isZero()) )
134 {
135 test = power(x,dr-dv)*v*LC(r,x);
136 if ( dr == 0 ) { r= CanonicalForm(0); }
137 else { r= r - LC(r,x)*power(x,dr); }
138 r= l*r -test;
139 dr= degree(r,x);
140 n+=1;
141 }
142 r= power(l, d-n)*r;
143 }
144 return r;
145}
int l
Definition cfEzgcd.cc:100
CanonicalForm test
Definition cfModGcd.cc:4104
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
static CanonicalForm * retvalue
Definition readcf.cc:126

◆ resultant()

CanonicalForm resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )

resultant() - return resultant of f and g with respect to x.

The chain is calculated from f and g with respect to variable x which should not be an algebraic variable. If f or q equals zero, zero is returned. If f is a coefficient with respect to x, f^degree(g, x) is returned, analogously for g.

This algorithm serves as a wrapper around other resultant algorithms which do the real work. Here we use standard properties of resultants only.

Definition at line 173 of file cf_resultant.cc.

174{
175 ASSERT( x.level() > 0, "cannot calculate resultant with respect to algebraic variables" );
176
177 // some checks on triviality. We will not use degree( v )
178 // here because this may involve variable swapping.
179 if ( f.isZero() || g.isZero() )
180 return 0;
181 if ( f.mvar() < x )
182 return power( f, g.degree( x ) );
183 if ( g.mvar() < x )
184 return power( g, f.degree( x ) );
185
186 // make x main variale
187 CanonicalForm F, G;
188 Variable X;
189 if ( f.mvar() > x || g.mvar() > x ) {
190 if ( f.mvar() > g.mvar() )
191 X = f.mvar();
192 else
193 X = g.mvar();
194 F = swapvar( f, X, x );
195 G = swapvar( g, X, x );
196 }
197 else {
198 X = x;
199 F = f;
200 G = g;
201 }
202 // at this point, we have to calculate resultant( F, G, X )
203 // where X is equal to or greater than the main variables
204 // of F and G
205
206 int m = degree( F, X );
207 int n = degree( G, X );
208 // catch trivial cases
209 if ( m+n <= 2 || m == 0 || n == 0 )
210 return swapvar( trivialResultant( F, G, X ), X, x );
211
212 // exchange F and G if necessary
213 int flipFactor;
214 if ( m < n ) {
216 F = G; G = swap;
217 int degswap = m;
218 m = n; n = degswap;
219 if ( m & 1 && n & 1 )
220 flipFactor = -1;
221 else
222 flipFactor = 1;
223 } else
224 flipFactor = 1;
225
226 // this is not an effective way to calculate the resultant!
227 CanonicalForm extFactor;
228 if ( m == n ) {
229 if ( n & 1 )
230 extFactor = -LC( G, X );
231 else
232 extFactor = LC( G, X );
233 } else
234 extFactor = power( LC( F, X ), m-n-1 );
235
237 result = subResChain( F, G, X )[0] / extFactor;
238
239 return swapvar( result, X, x ) * flipFactor;
240}
#define swap(_i, _j)
CFArray subResChain(const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
CFArray subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
static CanonicalForm trivialResultant(const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
static CanonicalForm trivialResultant ( const CanonicalForm & f, const CanonicalForm & g,...

◆ sqrFree()

CFFList FACTORY_PUBLIC sqrFree ( const CanonicalForm & f,
bool sort = false )

squarefree factorization

Definition at line 961 of file cf_factor.cc.

962{
963// ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
965
966 if ( getCharacteristic() == 0 )
967 result = sqrFreeZ( f );
968 else
969 {
971 if (hasFirstAlgVar (f, alpha))
972 result = FqSqrf( f, alpha );
973 else
974 result= FpSqrf (f);
975 }
976 if (sort)
977 {
978 CFFactor buf= result.getFirst();
979 result.removeFirst();
981 result.insert (buf);
982 }
983 return result;
984}
static void sort(int **points, int sizePoints)
CFFList FqSqrf(const CanonicalForm &F, const Variable &alpha, bool sort=true)
squarefree factorization over . If input is not monic, the leading coefficient is dropped
CFFList FpSqrf(const CanonicalForm &F, bool sort=true)
squarefree factorization over . If input is not monic, the leading coefficient is dropped
CFFList sqrFreeZ(const CanonicalForm &a)
CFFList sortCFFList(CFFList &F)
int status int void * buf
Definition si_signals.h:69

◆ subResChain()

CFArray subResChain ( const CanonicalForm & f,
const CanonicalForm & g,
const Variable & x )

CFArray subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )

subResChain() - caculate extended subresultant chain.

The chain is calculated from f and g with respect to variable x which should not be an algebraic variable. If f or q equals zero, an array consisting of one zero entry is returned.

Note: this is not the standard subresultant chain but the extended chain!

This algorithm is from the article of R. Loos - 'Generalized Polynomial Remainder Sequences' in B. Buchberger - 'Computer Algebra - Symbolic and Algebraic Computation' with some necessary extensions concerning the calculation of the first step.

Definition at line 42 of file cf_resultant.cc.

43{
44 ASSERT( x.level() > 0, "cannot calculate subresultant sequence with respect to algebraic variables" );
45
46 CFArray trivialResult( 0, 0 );
48 Variable X;
49
50 // some checks on triviality
51 if ( f.isZero() || g.isZero() ) {
52 trivialResult[0] = 0;
53 return trivialResult;
54 }
55
56 // make x main variable
57 if ( f.mvar() > x || g.mvar() > x ) {
58 if ( f.mvar() > g.mvar() )
59 X = f.mvar();
60 else
61 X = g.mvar();
62 F = swapvar( f, X, x );
63 G = swapvar( g, X, x );
64 }
65 else {
66 X = x;
67 F = f;
68 G = g;
69 }
70 // at this point, we have to calculate the sequence of F and
71 // G in respect to X where X is equal to or greater than the
72 // main variables of F and G
73
74 // initialization of chain
75 int m = degree( F, X );
76 int n = degree( G, X );
77
78 int j = (m <= n) ? n : m-1;
79 int r;
80
81 CFArray S( 0, j+1 );
83 S[j+1] = F; S[j] = G;
84
85 // make sure that S[j+1] is regular and j < n
86 if ( m == n && j > 0 ) {
87 S[j-1] = LC( S[j], X ) * psr( S[j+1], S[j], X );
88 j--;
89 } else if ( m < n ) {
90 S[j-1] = LC( S[j], X ) * LC( S[j], X ) * S[j+1];
91 j--;
92 } else if ( m > n && j > 0 ) {
93 // calculate first step
94 r = degree( S[j], X );
95 R = LC( S[j+1], X );
96
97 // if there was a gap calculate similar polynomial
98 if ( j > r && r >= 0 )
99 S[r] = power( LC( S[j], X ), j - r ) * S[j] * power( R, j - r );
100
101 if ( r > 0 ) {
102 // calculate remainder
103 S[r-1] = psr( S[j+1], S[j], X ) * power( -R, j - r );
104 j = r-1;
105 }
106 }
107
108 while ( j > 0 ) {
109 // at this point, 0 < j < n and S[j+1] is regular
110 r = degree( S[j], X );
111 R = LC( S[j+1], X );
112
113 // if there was a gap calculate similar polynomial
114 if ( j > r && r >= 0 )
115 S[r] = (power( LC( S[j], X ), j - r ) * S[j]) / power( R, j - r );
116
117 if ( r <= 0 ) break;
118 // calculate remainder
119 S[r-1] = psr( S[j+1], S[j], X ) / power( -R, j - r + 2 );
120
121 j = r-1;
122 // again 0 <= j < r <= jOld and S[j+1] is regular
123 }
124
125 for ( j = 0; j <= S.max(); j++ ) {
126 // reswap variables if necessary
127 if ( X != x ) {
128 S[j] = swapvar( S[j], X, x );
129 }
130 }
131
132 return S;
133}
CanonicalForm psr(const CanonicalForm &rr, const CanonicalForm &vv, const Variable &x)
CanonicalForm psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
#define R
Definition sirandom.c:27

◆ tryFdivides()

bool tryFdivides ( const CanonicalForm & f,
const CanonicalForm & g,
const CanonicalForm & M,
bool & fail )

same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f

Definition at line 456 of file cf_algorithm.cc.

457{
458 fail= false;
459 // trivial cases
460 if ( g.isZero() )
461 return true;
462 else if ( f.isZero() )
463 return false;
464
465 if (f.inCoeffDomain() || g.inCoeffDomain())
466 {
467 // if we are in a field all elements not equal to zero are units
468 if ( f.inCoeffDomain() )
469 {
470 CanonicalForm inv;
471 tryInvert (f, M, inv, fail);
472 return !fail;
473 }
474 else
475 {
476 return false;
477 }
478 }
479
480 // we may assume now that both levels either equal LEVELBASE
481 // or are greater zero
482 int fLevel = f.level();
483 int gLevel = g.level();
484 if ( (gLevel > 0) && (fLevel == gLevel) )
485 {
486 if (degree( f ) > degree( g ))
487 return false;
488 bool dividestail= tryFdivides (f.tailcoeff(), g.tailcoeff(), M, fail);
489
490 if (fail || !dividestail)
491 return false;
492 bool dividesLC= tryFdivides (f.LC(),g.LC(), M, fail);
493 if (fail || !dividesLC)
494 return false;
495 CanonicalForm q,r;
496 bool divides= tryDivremt (g, f, q, r, M, fail);
497 if (fail || !divides)
498 return false;
499 return r.isZero();
500 }
501 else if ( gLevel < fLevel )
502 {
503 // g is a coefficient w.r.t. f
504 return false;
505 }
506 else
507 {
508 // either f is a coefficient w.r.t. polynomial g or both
509 // f and g are from a base domain (should be Z or Z/p^n,
510 // then)
511 CanonicalForm q, r;
512 bool divides= tryDivremt (g, f, q, r, M, fail);
513 if (fail || !divides)
514 return false;
515 return r.isZero();
516 }
517}
bool tryDivremt(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r, const CanonicalForm &M, bool &fail)
same as divremt but handles zero divisors in case we are in Z_p[x]/(f) where f is not irreducible
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
bool tryFdivides(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f

Variable Documentation

◆ singular_homog_flag

EXTERN_VAR int singular_homog_flag

Definition at line 65 of file cf_algorithm.h.