cf_algorithm.cc File Reference

#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"

Go to the source code of this file.

Functions
void	out_cf (const char *s1, const CanonicalForm &f, const char *s2)
	cf_algorithm.cc - simple mathematical algorithms. More...
CanonicalForm	psr (const CanonicalForm &rr, const CanonicalForm &vv, const Variable &x)
	CanonicalForm psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) More...
CanonicalForm	psq (const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
	CanonicalForm psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) More...
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 ) More...
static CanonicalForm	internalBCommonDen (const CanonicalForm &f)
	static CanonicalForm internalBCommonDen ( const CanonicalForm & f ) More...
CanonicalForm	bCommonDen (const CanonicalForm &f)
	CanonicalForm bCommonDen ( const CanonicalForm & f ) More...
bool	fdivides (const CanonicalForm &f, const CanonicalForm &g)
	bool fdivides ( const CanonicalForm & f, const CanonicalForm & g ) More...
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 More...
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 More...
CanonicalForm	maxNorm (const CanonicalForm &f)
	CanonicalForm maxNorm ( const CanonicalForm & f ) More...
CanonicalForm	euclideanNorm (const CanonicalForm &f)
	CanonicalForm euclideanNorm ( const CanonicalForm & f ) More...

Function Documentation

bCommonDen()

CanonicalForm bCommonDen

(

const CanonicalForm &

f

)

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.

{
    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 );
}

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.

{
    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 );
}

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.

{
    // 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();
    }
}

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.

{
    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;
    }
}

internalBCommonDen()

static CanonicalForm internalBCommonDen ( const CanonicalForm &  f )

static

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 )

Definition at line 262 of file cf_algorithm.cc.

{
    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;
    }
}

maxNorm()

CanonicalForm maxNorm

(

const CanonicalForm &

f

)

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.

{
    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;
    }
}

out_cf()

void out_cf	(	const char *	s1,
		const CanonicalForm &	f,
		const char *	s2
	)

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.

Definition at line 99 of file cf_factor.cc.

{
  printf("%s",s1);
  if (f.isZero()) printf("+0");
  //else if (! f.inCoeffDomain() )
  else if (! f.inBaseDomain() )
  {
    int l = f.level();
    for ( CFIterator i = f; i.hasTerms(); i++ )
    {
      int e=i.exp();
      if (i.coeff().isOne())
      {
        printf("+");
        if (e==0) printf("1");
        else
        {
          printf("%c",'a'+l-1);
          if (e!=1) printf("^%d",e);
        }
      }
      else
      {
        out_cf("+(",i.coeff(),")");
        if (e!=0)
        {
          printf("*%c",'a'+l-1);
          if (e!=1) printf("^%d",e);
        }
      }
    }
  }
  else
  {
    if ( f.isImm() )
    {
      if (CFFactory::gettype()==GaloisFieldDomain)
      {
         long a= imm2int (f.getval());
         if ( a == gf_q )
           printf ("+%ld", a);
         else  if ( a == 0L )
           printf ("+1");
         else  if ( a == 1L )
           printf ("+%c",gf_name);
         else
         {
           printf ("+%c",gf_name);
           printf ("^%ld",a);
         }
      }
      else
      {
        long l=f.intval();
        if (l<0) printf("%ld",l);
        else     printf("+%ld",l);
      }
    }
    else
    {
    #ifdef NOSTREAMIO
      if (f.inZ())
      {
        mpz_t m;
        gmp_numerator(f,m);
        char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
        str = mpz_get_str( str, 10, m );
        puts(str);
        delete[] str;
        mpz_clear(m);
      }
      else if (f.inQ())
      {
        mpz_t m;
        gmp_numerator(f,m);
        char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
        str = mpz_get_str( str, 10, m );
        while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
        puts(str);putchar('/');
        delete[] str;
        mpz_clear(m);
        gmp_denominator(f,m);
        str = new char[mpz_sizeinbase( m, 10 ) + 2];
        str = mpz_get_str( str, 10, m );
        while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
        puts(str);
        delete[] str;
        mpz_clear(m);
      }
    #else
       std::cout << f;
    #endif
    }
    //if (f.inZ()) printf("(Z)");
    //else if (f.inQ()) printf("(Q)");
    //else if (f.inFF()) printf("(FF)");
    //else if (f.inPP()) printf("(PP)");
    //else if (f.inGF()) printf("(PP)");
    //else
    if (f.inExtension()) printf("E(%d)",f.level());
  }
  printf("%s",s2);
}

psq()

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

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.

{
    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 );
    }
}

psqr()

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 )

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.

{
    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 );
    }
}

psr()

CanonicalForm psr	(	const CanonicalForm &	rr,
		const CanonicalForm &	vv,
		const Variable &	x
	)

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.

{
  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;
}

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.

{
    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();
    }
}