Loading...
Searching...
No Matches
cf_algorithm.cc
Go to the documentation of this file.
/* emacs edit mode for this file is -*- C++ -*- */
/**
*
*
* cf_algorithm.cc - simple mathematical algorithms.
*
* Hierarchy: mathematical algorithms on canonical forms
*
* Developers note:
* ----------------
* A "mathematical" algorithm is an algorithm which calculates
* some mathematical function in contrast to a "structural"
* algorithm which gives structural information on polynomials.
*
* Compare these functions to the functions in `cf_ops.cc', which
* are structural algorithms.
*
**/
#include "config.h"
#include "cf_assert.h"
#include "cf_factory.h"
#include "cf_defs.h"
#include "canonicalform.h"
#include "cf_algorithm.h"
#include "variable.h"
#include "cf_iter.h"
#include "templates/ftmpl_functions.h"
#include "cfGcdAlgExt.h"
void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
/** 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.
*
* @sa psq(), psqr()
**/
CanonicalForm
#if 0
psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
{
ASSERT( x.level() > 0, "type error: polynomial variable expected" );
ASSERT( ! g.isZero(), "math error: division by zero" );
// swap variables such that x's level is larger or equal
// than both f's and g's levels.
Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
CanonicalForm F = swapvar( f, x, X );
CanonicalForm G = swapvar( g, x, X );
// now, we have to calculate the pseudo remainder of F and G
// w.r.t. X
int fDegree = degree( F, X );
int gDegree = degree( G, X );
if ( (fDegree < 0) || (fDegree < gDegree) )
return f;
else
{
CanonicalForm xresult = (power( LC( G, X ), fDegree-gDegree+1 ) * F) ;
CanonicalForm result = xresult -(xresult/G)*G;
return swapvar( result, x, X );
}
}
#else
psr ( const CanonicalForm &rr, const CanonicalForm &vv, const Variable & x )
{
CanonicalForm r=rr, v=vv, l, test, lu, lv, t, retvalue;
int dr, dv, d,n=0;
dr = degree( r, x );
if (dr>0)
{
dv = degree( v, x );
if (dv <= dr) {l=LC(v,x); v = v -l*power(x,dv);}
else { l = 1; }
d= dr-dv+1;
//out_cf("psr(",rr," ");
//out_cf("",vv," ");
//printf(" var=%d\n",x.level());
while ( ( dv <= dr ) && ( !r.isZero()) )
{
test = power(x,dr-dv)*v*LC(r,x);
if ( dr == 0 ) { r= CanonicalForm(0); }
else { r= r - LC(r,x)*power(x,dr); }
r= l*r -test;
dr= degree(r,x);
n+=1;
}
r= power(l, d-n)*r;
}
return r;
}
#endif
/** 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.
*
* @sa psr(), psqr()
*
**/
CanonicalForm
psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
{
ASSERT( x.level() > 0, "type error: polynomial variable expected" );
ASSERT( ! g.isZero(), "math error: division by zero" );
// swap variables such that x's level is larger or equal
// than both f's and g's levels.
Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
CanonicalForm F = swapvar( f, x, X );
CanonicalForm G = swapvar( g, x, X );
// now, we have to calculate the pseudo remainder of F and G
// w.r.t. X
int fDegree = degree( F, X );
int gDegree = degree( G, X );
if ( fDegree < 0 || fDegree < gDegree )
return 0;
else {
CanonicalForm result = (power( LC( G, X ), fDegree-gDegree+1 ) * F) / G;
return swapvar( result, x, X );
}
}
/** 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.
*
* @sa psr(), psq()
*
**/
void
psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable& x )
{
ASSERT( x.level() > 0, "type error: polynomial variable expected" );
ASSERT( ! g.isZero(), "math error: division by zero" );
// swap variables such that x's level is larger or equal
// than both f's and g's levels.
Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
CanonicalForm F = swapvar( f, x, X );
CanonicalForm G = swapvar( g, x, X );
// now, we have to calculate the pseudo remainder of F and G
// w.r.t. X
int fDegree = degree( F, X );
int gDegree = degree( G, X );
if ( fDegree < 0 || fDegree < gDegree ) {
q = 0; r = f;
} else {
divrem( power( LC( G, X ), fDegree-gDegree+1 ) * F, G, q, r );
q = swapvar( q, x, X );
r = swapvar( r, x, X );
}
}
/** static CanonicalForm internalBCommonDen ( const CanonicalForm & f )
*
*
* internalBCommonDen() - recursively calculate multivariate
* common denominator of coefficients of `f'.
*
* Used by: bCommonDen()
*
* Type info:
* ----------
* f: Poly( Q )
* Switches: isOff( SW_RATIONAL )
*
**/
static CanonicalForm
internalBCommonDen ( const CanonicalForm & f )
{
if ( f.inBaseDomain() )
return f.den();
else {
CanonicalForm result = 1;
for ( CFIterator i = f; i.hasTerms(); i++ )
result = blcm( result, internalBCommonDen( i.coeff() ) );
return result;
}
}
/** 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
*
**/
CanonicalForm
bCommonDen ( const CanonicalForm & f )
{
if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) )
{
// otherwise `bgcd()' returns one
Off( SW_RATIONAL );
CanonicalForm result = internalBCommonDen( f );
On( SW_RATIONAL );
return result;
}
else
return CanonicalForm( 1 );
}
/** 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()'.
*
**/
bool
fdivides ( const CanonicalForm & f, const CanonicalForm & g )
{
// trivial cases
if ( g.isZero() )
return true;
else if ( f.isZero() )
return false;
if ( (f.inCoeffDomain() || g.inCoeffDomain())
&& ((getCharacteristic() == 0 && isOn( SW_RATIONAL ))
|| (getCharacteristic() > 0) ))
{
// if we are in a field all elements not equal to zero are units
if ( f.inCoeffDomain() )
return true;
else
// g.inCoeffDomain()
return false;
}
// we may assume now that both levels either equal LEVELBASE
// or are greater zero
int fLevel = f.level();
int gLevel = g.level();
if ( (gLevel > 0) && (fLevel == gLevel) )
// f and g are polynomials in the same main variable
if ( degree( f ) <= degree( g )
&& fdivides( f.tailcoeff(), g.tailcoeff() )
&& fdivides( f.LC(), g.LC() ) )
{
CanonicalForm q, r;
return divremt( g, f, q, r ) && r.isZero();
}
else
return false;
else if ( gLevel < fLevel )
// g is a coefficient w.r.t. f
return false;
else
{
// either f is a coefficient w.r.t. polynomial g or both
// f and g are from a base domain (should be Z or Z/p^n,
// then)
CanonicalForm q, r;
return divremt( g, f, q, r ) && r.isZero();
}
}
/// same as fdivides if true returns quotient quot of g by f otherwise quot == 0
bool
fdivides ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm& quot )
{
quot= 0;
// trivial cases
if ( g.isZero() )
return true;
else if ( f.isZero() )
return false;
if ( (f.inCoeffDomain() || g.inCoeffDomain())
&& ((getCharacteristic() == 0 && isOn( SW_RATIONAL ))
|| (getCharacteristic() > 0) ))
{
// if we are in a field all elements not equal to zero are units
if ( f.inCoeffDomain() )
{
quot= g/f;
return true;
}
else
// g.inCoeffDomain()
return false;
}
// we may assume now that both levels either equal LEVELBASE
// or are greater zero
int fLevel = f.level();
int gLevel = g.level();
if ( (gLevel > 0) && (fLevel == gLevel) )
// f and g are polynomials in the same main variable
if ( degree( f ) <= degree( g )
&& fdivides( f.tailcoeff(), g.tailcoeff() )
&& fdivides( f.LC(), g.LC() ) )
{
CanonicalForm q, r;
if (divremt( g, f, q, r ) && r.isZero())
{
quot= q;
return true;
}
else
return false;
}
else
return false;
else if ( gLevel < fLevel )
// g is a coefficient w.r.t. f
return false;
else
{
// either f is a coefficient w.r.t. polynomial g or both
// f and g are from a base domain (should be Z or Z/p^n,
// then)
CanonicalForm q, r;
if (divremt( g, f, q, r ) && r.isZero())
{
quot= q;
return true;
}
else
return false;
}
}
/// same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f
bool
tryFdivides ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm& M, bool& fail )
{
fail= false;
// trivial cases
if ( g.isZero() )
return true;
else if ( f.isZero() )
return false;
if (f.inCoeffDomain() || g.inCoeffDomain())
{
// if we are in a field all elements not equal to zero are units
if ( f.inCoeffDomain() )
{
CanonicalForm inv;
tryInvert (f, M, inv, fail);
return !fail;
}
else
{
return false;
}
}
// we may assume now that both levels either equal LEVELBASE
// or are greater zero
int fLevel = f.level();
int gLevel = g.level();
if ( (gLevel > 0) && (fLevel == gLevel) )
{
if (degree( f ) > degree( g ))
return false;
bool dividestail= tryFdivides (f.tailcoeff(), g.tailcoeff(), M, fail);
if (fail || !dividestail)
return false;
bool dividesLC= tryFdivides (f.LC(),g.LC(), M, fail);
if (fail || !dividesLC)
return false;
CanonicalForm q,r;
bool divides= tryDivremt (g, f, q, r, M, fail);
if (fail || !divides)
return false;
return r.isZero();
}
else if ( gLevel < fLevel )
{
// g is a coefficient w.r.t. f
return false;
}
else
{
// either f is a coefficient w.r.t. polynomial g or both
// f and g are from a base domain (should be Z or Z/p^n,
// then)
CanonicalForm q, r;
bool divides= tryDivremt (g, f, q, r, M, fail);
if (fail || !divides)
return false;
return r.isZero();
}
}
/** 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
*
**/
CanonicalForm
maxNorm ( const CanonicalForm & f )
{
if ( f.inBaseDomain() )
return abs( f );
else {
CanonicalForm result = 0;
for ( CFIterator i = f; i.hasTerms(); i++ ) {
CanonicalForm coeffMaxNorm = maxNorm( i.coeff() );
if ( coeffMaxNorm > result )
result = coeffMaxNorm;
}
return result;
}
}
/** 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 )
*
**/
CanonicalForm
euclideanNorm ( const CanonicalForm & f )
{
ASSERT( (f.inBaseDomain() || f.isUnivariate()) && f.LC().inZ(),
"type error: univariate poly over Z expected" );
CanonicalForm result = 0;
for ( CFIterator i = f; i.hasTerms(); i++ ) {
CanonicalForm coeff = i.coeff();
result += coeff*coeff;
}
return sqrt( result );
}
abs
Rational abs(const Rational &a)
Definition GMPrat.cc:436
divrem
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
Definition canonicalform.cc:1032
isOn
bool isOn(int sw)
switches
Definition canonicalform.cc:1971
tryDivremt
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
Definition canonicalform.cc:1108
On
void On(int sw)
switches
Definition canonicalform.cc:1957
blcm
CanonicalForm blcm(const CanonicalForm &f, const CanonicalForm &g)
Definition canonicalform.cc:1816
divremt
bool divremt(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
Definition canonicalform.cc:1066
Off
void Off(int sw)
switches
Definition canonicalform.cc:1964
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition canonicalform.cc:1896
canonicalform.h
Header for factory's main class CanonicalForm.
degree
int degree(const CanonicalForm &f)
Definition canonicalform.h:316
swapvar
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition cf_ops.cc:168
LC
CanonicalForm LC(const CanonicalForm &f)
Definition canonicalform.h:310
getCharacteristic
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
l
int l
Definition cfEzgcd.cc:100
i
int i
Definition cfEzgcd.cc:132
tryInvert
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
Definition cfGcdAlgExt.cc:222
cfGcdAlgExt.h
GCD over Q(a).
x
Variable x
Definition cfModGcd.cc:4090
g
g
Definition cfModGcd.cc:4098
test
CanonicalForm test
Definition cfModGcd.cc:4104
bCommonDen
CanonicalForm bCommonDen(const CanonicalForm &f)
void psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r,...
Definition cf_algorithm.cc:293
maxNorm
CanonicalForm maxNorm(const CanonicalForm &f)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
Definition cf_algorithm.cc:536
psq
CanonicalForm psq(const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
cf_algorithm.cc - simple mathematical algorithms.
Definition cf_algorithm.cc:172
cf_algorithm.h
declarations of higher level algorithms.
psqr
void psqr(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r, const Variable &x)
fdivides
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
tryFdivides
bool tryFdivides(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
psr
CanonicalForm psr(const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
euclideanNorm
CanonicalForm euclideanNorm(const CanonicalForm &f)
cf_assert.h
assertions for Factory
ASSERT
#define ASSERT(expression, message)
Definition cf_assert.h:99
cf_defs.h
factory switches.
SW_RATIONAL
static const int SW_RATIONAL
set to 1 for computations over Q
Definition cf_defs.h:31
cf_factory.h
Interface to generate InternalCF's over various domains from intrinsic types or mpz_t's.
cf_iter.h
Iterators for CanonicalForm's.
f
FILE * f
Definition checklibs.c:9
out_cf
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
Definition cf_factor.cc:103
CFIterator
class to iterate through CanonicalForm's
Definition cf_iter.h:44
CanonicalForm
factory's main class
Definition canonicalform.h:86
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
Variable
factory's class for variables
Definition factory.h:127
result
return result
Definition facAbsBiFact.cc:76
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
ftmpl_functions.h
some useful template functions.
tmax
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
G
STATIC_VAR TreeM * G
Definition janet.cc:31
sqrt
gmp_float sqrt(const gmp_float &a)
Definition mpr_complex.cc:319
retvalue
static CanonicalForm * retvalue
Definition readcf.cc:126
M
#define M
Definition sirandom.c:25
variable.h
operations on variables