My Project
Loading...
Searching...
No Matches
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
36void 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
91psr ( 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
117psr ( 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**/
172psq ( 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**/
222void
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**/
261static CanonicalForm
263{
264 if ( f.inBaseDomain() )
265 return f.den();
266 else {
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**/
339bool
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
389bool
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
455bool
456tryFdivides ( 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 {
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
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: variable.h:33
int level() const
Definition: variable.h:49
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
#define M
Definition: sirandom.c:25
operations on variables