My Project
cf_algorithm.cc
Go to the documentation of this file.
1 /* emacs edit mode for this file is -*- C++ -*- */
2 
3 /**
4  *
5  *
6  * cf_algorithm.cc - simple mathematical algorithms.
7  *
8  * Hierarchy: mathematical algorithms on canonical forms
9  *
10  * Developers note:
11  * ----------------
12  * A "mathematical" algorithm is an algorithm which calculates
13  * some mathematical function in contrast to a "structural"
14  * algorithm which gives structural information on polynomials.
15  *
16  * Compare these functions to the functions in `cf_ops.cc', which
17  * are structural algorithms.
18  *
19 **/
20 
21 
22 #include "config.h"
23 
24 
25 #include "cf_assert.h"
26 
27 #include "cf_factory.h"
28 #include "cf_defs.h"
29 #include "canonicalform.h"
30 #include "cf_algorithm.h"
31 #include "variable.h"
32 #include "cf_iter.h"
34 #include "cfGcdAlgExt.h"
35 
36 void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
37 
38 /** CanonicalForm psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
39  *
40  *
41  * psr() - return pseudo remainder of `f' and `g' with respect
42  * to `x'.
43  *
44  * `g' must not equal zero.
45  *
46  * For f and g in R[x], R an arbitrary ring, g != 0, there is a
47  * representation
48  *
49  * LC(g)^s*f = g*q + r
50  *
51  * with r = 0 or deg(r) < deg(g) and s = 0 if f = 0 or
52  * s = max( 0, deg(f)-deg(g)+1 ) otherwise.
53  * r = psr(f, g) and q = psq(f, g) are called "pseudo remainder"
54  * and "pseudo quotient", resp. They are uniquely determined if
55  * LC(g) is not a zero divisor in R.
56  *
57  * See H.-J. Reiffen/G. Scheja/U. Vetter - "Algebra", 2nd ed.,
58  * par. 15, for a reference.
59  *
60  * Type info:
61  * ----------
62  * f, g: Current
63  * x: Polynomial
64  *
65  * Polynomials over prime power domains are admissible if
66  * lc(LC(`g',`x')) is not a zero divisor. This is a slightly
67  * stronger precondition than mathematically necessary since
68  * pseudo remainder and quotient are well-defined if LC(`g',`x')
69  * is not a zero divisor.
70  *
71  * For example, psr(y^2, (13*x+1)*y) is well-defined in
72  * (Z/13^2[x])[y] since (13*x+1) is not a zero divisor. But
73  * calculating it with Factory would fail since 13 is a zero
74  * divisor in Z/13^2.
75  *
76  * Due to this inconsistency with mathematical notion, we decided
77  * not to declare type `CurrentPP' for `f' and `g'.
78  *
79  * Developers note:
80  * ----------------
81  * This is not an optimal implementation. Better would have been
82  * an implementation in `InternalPoly' avoiding the
83  * exponentiation of the leading coefficient of `g'. In contrast
84  * to `psq()' and `psqr()' it definitely seems worth to implement
85  * the pseudo remainder on the internal level.
86  *
87  * @sa psq(), psqr()
88 **/
90 #if 0
91 psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
92 {
93 
94  ASSERT( x.level() > 0, "type error: polynomial variable expected" );
95  ASSERT( ! g.isZero(), "math error: division by zero" );
96 
97  // swap variables such that x's level is larger or equal
98  // than both f's and g's levels.
99  Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
100  CanonicalForm F = swapvar( f, x, X );
101  CanonicalForm G = swapvar( g, x, X );
102 
103  // now, we have to calculate the pseudo remainder of F and G
104  // w.r.t. X
105  int fDegree = degree( F, X );
106  int gDegree = degree( G, X );
107  if ( (fDegree < 0) || (fDegree < gDegree) )
108  return f;
109  else
110  {
111  CanonicalForm xresult = (power( LC( G, X ), fDegree-gDegree+1 ) * F) ;
112  CanonicalForm result = xresult -(xresult/G)*G;
113  return swapvar( result, x, X );
114  }
115 }
116 #else
117 psr ( const CanonicalForm &rr, const CanonicalForm &vv, const Variable & x )
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 }
146 #endif
147 
148 /** CanonicalForm psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
149  *
150  *
151  * psq() - return pseudo quotient of `f' and `g' with respect
152  * to `x'.
153  *
154  * `g' must not equal zero.
155  *
156  * Type info:
157  * ----------
158  * f, g: Current
159  * x: Polynomial
160  *
161  * Developers note:
162  * ----------------
163  * This is not an optimal implementation. Better would have been
164  * an implementation in `InternalPoly' avoiding the
165  * exponentiation of the leading coefficient of `g'. It seemed
166  * not worth to do so.
167  *
168  * @sa psr(), psqr()
169  *
170 **/
172 psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
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 }
194 
195 /** void psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable & x )
196  *
197  *
198  * psqr() - calculate pseudo quotient and remainder of `f' and
199  * `g' with respect to `x'.
200  *
201  * Returns the pseudo quotient of `f' and `g' in `q', the pseudo
202  * remainder in `r'. `g' must not equal zero.
203  *
204  * See `psr()' for more detailed information.
205  *
206  * Type info:
207  * ----------
208  * f, g: Current
209  * q, r: Anything
210  * x: Polynomial
211  *
212  * Developers note:
213  * ----------------
214  * This is not an optimal implementation. Better would have been
215  * an implementation in `InternalPoly' avoiding the
216  * exponentiation of the leading coefficient of `g'. It seemed
217  * not worth to do so.
218  *
219  * @sa psr(), psq()
220  *
221 **/
222 void
223 psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable& x )
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 }
246 
247 /** static CanonicalForm internalBCommonDen ( const CanonicalForm & f )
248  *
249  *
250  * internalBCommonDen() - recursively calculate multivariate
251  * common denominator of coefficients of `f'.
252  *
253  * Used by: bCommonDen()
254  *
255  * Type info:
256  * ----------
257  * f: Poly( Q )
258  * Switches: isOff( SW_RATIONAL )
259  *
260 **/
261 static CanonicalForm
263 {
264  if ( f.inBaseDomain() )
265  return f.den();
266  else {
267  CanonicalForm result = 1;
268  for ( CFIterator i = f; i.hasTerms(); i++ )
269  result = blcm( result, internalBCommonDen( i.coeff() ) );
270  return result;
271  }
272 }
273 
274 /** CanonicalForm bCommonDen ( const CanonicalForm & f )
275  *
276  *
277  * bCommonDen() - calculate multivariate common denominator of
278  * coefficients of `f'.
279  *
280  * The common denominator is calculated with respect to all
281  * coefficients of `f' which are in a base domain. In other
282  * words, common_den( `f' ) * `f' is guaranteed to have integer
283  * coefficients only. The common denominator of zero is one.
284  *
285  * Returns something non-trivial iff the current domain is Q.
286  *
287  * Type info:
288  * ----------
289  * f: CurrentPP
290  *
291 **/
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 }
306 
307 /** bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
308  *
309  *
310  * fdivides() - check whether `f' divides `g'.
311  *
312  * Returns true iff `f' divides `g'. Uses some extra heuristic
313  * to avoid polynomial division. Without the heuristic, the test
314  * essentialy looks like `divremt(g, f, q, r) && r.isZero()'.
315  *
316  * Type info:
317  * ----------
318  * f, g: Current
319  *
320  * Elements from prime power domains (or polynomials over such
321  * domains) are admissible if `f' (or lc(`f'), resp.) is not a
322  * zero divisor. This is a slightly stronger precondition than
323  * mathematically necessary since divisibility is a well-defined
324  * notion in arbitrary rings. Hence, we decided not to declare
325  * the weaker type `CurrentPP'.
326  *
327  * Developers note:
328  * ----------------
329  * One may consider the the test `fdivides( f.LC(), g.LC() )' in
330  * the main `if'-test superfluous since `divremt()' in the
331  * `if'-body repeats the test. However, `divremt()' does not use
332  * any heuristic to do so.
333  *
334  * It seems not reasonable to call `fdivides()' from `divremt()'
335  * to check divisibility of leading coefficients. `fdivides()' is
336  * on a relatively high level compared to `divremt()'.
337  *
338 **/
339 bool
340 fdivides ( const CanonicalForm & f, const CanonicalForm & g )
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 }
387 
388 /// same as fdivides if true returns quotient quot of g by f otherwise quot == 0
389 bool
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 }
453 
454 /// same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f
455 bool
456 tryFdivides ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm& M, bool& fail )
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 }
518 
519 /** CanonicalForm maxNorm ( const CanonicalForm & f )
520  *
521  *
522  * maxNorm() - return maximum norm of `f'.
523  *
524  * That is, the base coefficient of `f' with the largest absolute
525  * value.
526  *
527  * Valid for arbitrary polynomials over arbitrary domains, but
528  * most useful for multivariate polynomials over Z.
529  *
530  * Type info:
531  * ----------
532  * f: CurrentPP
533  *
534 **/
537 {
538  if ( f.inBaseDomain() )
539  return abs( f );
540  else {
541  CanonicalForm result = 0;
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 }
550 
551 /** CanonicalForm euclideanNorm ( const CanonicalForm & f )
552  *
553  *
554  * euclideanNorm() - return Euclidean norm of `f'.
555  *
556  * Returns the largest integer smaller or equal norm(`f') =
557  * sqrt(sum( `f'[i]^2 )).
558  *
559  * Type info:
560  * ----------
561  * f: UVPoly( Z )
562  *
563 **/
566 {
567  ASSERT( (f.inBaseDomain() || f.isUnivariate()) && f.LC().inZ(),
568  "type error: univariate poly over Z expected" );
569 
570  CanonicalForm result = 0;
571  for ( CFIterator i = f; i.hasTerms(); i++ ) {
572  CanonicalForm coeff = i.coeff();
573  result += coeff*coeff;
574  }
575  return sqrt( result );
576 }
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
bool isOn(int sw)
switches
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 On(int sw)
switches
CanonicalForm blcm(const CanonicalForm &f, const CanonicalForm &g)
bool divremt(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
void Off(int sw)
switches
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Header for factory's main class CanonicalForm.
int degree(const CanonicalForm &f)
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)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
int l
Definition: cfEzgcd.cc:100
int i
Definition: cfEzgcd.cc:132
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
Definition: cfGcdAlgExt.cc:221
GCD over Q(a)
Variable x
Definition: cfModGcd.cc:4082
g
Definition: cfModGcd.cc:4090
CanonicalForm test
Definition: cfModGcd.cc:4096
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,...
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
static CanonicalForm internalBCommonDen(const CanonicalForm &f)
static CanonicalForm internalBCommonDen ( const CanonicalForm & f )
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
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
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
cf_algorithm.cc - simple mathematical algorithms.
Definition: cf_factor.cc:99
CanonicalForm psr(const CanonicalForm &rr, const CanonicalForm &vv, 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 )
CanonicalForm euclideanNorm(const CanonicalForm &f)
CanonicalForm euclideanNorm ( const CanonicalForm & f )
declarations of higher level algorithms.
assertions for Factory
#define ASSERT(expression, message)
Definition: cf_assert.h:99
factory switches.
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:31
Interface to generate InternalCF's over various domains from intrinsic types or mpz_t's.
Iterators for CanonicalForm's.
FILE * f
Definition: checklibs.c:9
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
factory's class for variables
Definition: factory.h:127
int level() const
Definition: factory.h:143
return result
Definition: facAbsBiFact.cc:75
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
some useful template functions.
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
STATIC_VAR TreeM * G
Definition: janet.cc:31
gmp_float sqrt(const gmp_float &a)
Definition: mpr_complex.cc:327
static CanonicalForm * retvalue
Definition: readcf.cc:126
#define M
Definition: sirandom.c:25
operations on variables