uResultant Class Reference

Base class for solving 0-dim poly systems using u-resultant. More...

#include <mpr_base.h>

Public Types
enum	resMatType { none , sparseResMat , denseResMat }

Public Member Functions
	uResultant (const ideal _gls, const resMatType _rmt=sparseResMat, BOOLEAN extIdeal=true)
	~uResultant ()
poly	interpolateDense (const number subDetVal=NULL)
rootContainer **	interpolateDenseSP (BOOLEAN matchUp=false, const number subDetVal=NULL)
rootContainer **	specializeInU (BOOLEAN matchUp=false, const number subDetVal=NULL)
resMatrixBase *	accessResMat ()

Private Member Functions
	uResultant (const uResultant &)
ideal	extendIdeal (const ideal gls, poly linPoly, const resMatType rmt)
poly	linearPoly (const resMatType rmt)
int	nextPrime (const int p)

Private Attributes
ideal	gls
int	n
resMatType	rmt
resMatrixBase *	resMat

Detailed Description

Base class for solving 0-dim poly systems using u-resultant.

Definition at line 62 of file mpr_base.h.

Member Enumeration Documentation

resMatType

enum uResultant::resMatType

Enumerator
none
sparseResMat
denseResMat

Definition at line 65 of file mpr_base.h.

{ none, sparseResMat, denseResMat };

Constructor & Destructor Documentation

uResultant() [1/2]

uResultant::uResultant	(	const ideal	_gls,
		const resMatType	_rmt = sparseResMat,
		BOOLEAN	extIdeal = true
	)

Definition at line 2685 of file mpr_base.cc.

  : rmt( _rmt )
{
  if ( extIdeal )
  {
    // extend given ideal by linear poly F0=u0x0 + u1x1 +...+ unxn
    gls= extendIdeal( _gls, linearPoly( rmt ), rmt );
    n= IDELEMS( gls );
  }
  else
    gls= idCopy( _gls );
 
  switch ( rmt )
  {
  case sparseResMat:
    resMat= new resMatrixSparse( gls );
    break;
  case denseResMat:
    resMat= new resMatrixDense( gls );
    break;
  default:
    WerrorS("uResultant::uResultant: Unknown chosen resultant matrix type!");
  }
}

~uResultant()

uResultant::~uResultant

(

)

Definition at line 2710 of file mpr_base.cc.

{
  delete resMat;
}

uResultant() [2/2]

uResultant::uResultant ( const uResultant &  )

private

Member Function Documentation

accessResMat()

resMatrixBase * uResultant::accessResMat ( )

inline

Definition at line 78 of file mpr_base.h.

{ return resMat; }

extendIdeal()

ideal uResultant::extendIdeal ( const ideal  gls, poly  linPoly, const resMatType  rmt  )

private

Definition at line 2715 of file mpr_base.cc.

{
  ideal newGls= idCopy( igls );
  newGls->m= (poly *)omReallocSize( newGls->m,
                              IDELEMS(igls) * sizeof(poly),
                              (IDELEMS(igls) + 1) * sizeof(poly) );
  IDELEMS(newGls)++;
 
  switch ( rrmt )
  {
  case sparseResMat:
  case denseResMat:
    {
      int i;
      for ( i= IDELEMS(newGls)-1; i > 0; i-- )
      {
        newGls->m[i]= newGls->m[i-1];
      }
      newGls->m[0]= linPoly;
    }
    break;
  default:
    WerrorS("uResultant::extendIdeal: Unknown chosen resultant matrix type!");
  }
 
  return( newGls );
}

interpolateDense()

poly uResultant::interpolateDense

(

const number

subDetVal = NULL

)

Definition at line 2770 of file mpr_base.cc.

{
  int i,j,p;
  long tdg;
 
  // D is a Polynom homogeneous in the coeffs of F0 and of degree tdg = d0*d1*...*dn
  tdg= resMat->getDetDeg();
 
  // maximum number of terms in polynom D (homogeneous, of degree tdg)
  // long mdg= (facul(tdg+n-1) / facul( tdg )) / facul( n - 1 );
  long mdg= over( n-1, tdg );
 
  // maximal number of terms in a polynom of degree tdg
  long l=(long)pow( (double)(tdg+1), n );
 
#ifdef mprDEBUG_PROT
  Print("// total deg of D: tdg %ld\n",tdg);
  Print("// maximum number of terms in D: mdg: %ld\n",mdg);
  Print("// maximum number of terms in polynom of deg tdg: l %ld\n",l);
#endif
 
  // we need mdg results of D(p0,p1,...,pn)
  number *presults;
  presults= (number *)omAlloc( mdg * sizeof( number ) );
  for (i=0; i < mdg; i++) presults[i]= nInit(0);
 
  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
  number *pev= (number *)omAlloc( n * sizeof( number ) );
  for (i=0; i < n; i++) pev[i]= nInit(0);
 
  mprPROTnl("// initial evaluation point: ");
  // initial evaluatoin point
  p=1;
  for (i=0; i < n; i++)
  {
    // init pevpoint with primes 3,5,7,11, ...
    p= nextPrime( p );
    pevpoint[i]=nInit( p );
    nTest(pevpoint[i]);
    mprPROTNnl(" ",pevpoint[i]);
  }
 
  // evaluate the determinant in the points pev^0, pev^1, ..., pev^mdg
  mprPROTnl("// evaluating:");
  for ( i=0; i < mdg; i++ )
  {
    for (j=0; j < n; j++)
    {
      nDelete( &pev[j] );
      nPower(pevpoint[j],i,&pev[j]);
      mprPROTN(" ",pev[j]);
    }
    mprPROTnl("");
 
    nDelete( &presults[i] );
    presults[i]=resMat->getDetAt( pev );
 
    mprSTICKYPROT(ST_BASE_EV);
  }
  mprSTICKYPROT("\n");
 
  // now interpolate using vandermode interpolation
  mprPROTnl("// interpolating:");
  number *ncpoly;
  {
    vandermonde vm( mdg, n, tdg, pevpoint );
    ncpoly= vm.interpolateDense( presults );
  }
 
  if ( subDetVal != NULL )
  {   // divide by common factor
    number detdiv;
    for ( i= 0; i <= mdg; i++ )
    {
      detdiv= nDiv( ncpoly[i], subDetVal );
      nNormalize( detdiv );
      nDelete( &ncpoly[i] );
      ncpoly[i]= detdiv;
    }
  }
 
#ifdef mprDEBUG_ALL
  PrintLn();
  for ( i=0; i < mdg; i++ )
  {
    nPrint(ncpoly[i]); PrintS(" --- ");
  }
  PrintLn();
#endif
 
  // prepare ncpoly for later use
  number nn=nInit(0);
  for ( i=0; i < mdg; i++ )
  {
    if ( nEqual(ncpoly[i],nn) )
    {
      nDelete( &ncpoly[i] );
      ncpoly[i]=NULL;
    }
  }
  nDelete( &nn );
 
  // create poly presenting the determinat of the uResultant
  intvec exp( n );
  for ( i= 0; i < n; i++ ) exp[i]=0;
 
  poly result= NULL;
 
  long sum=0;
  long c=0;
 
  for ( i=0; i < l; i++ )
  {
    if ( sum == tdg )
    {
      if ( !nIsZero(ncpoly[c]) )
      {
        poly p= pOne();
        if ( rmt == denseResMat )
        {
          for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
        }
        else if ( rmt == sparseResMat )
        {
          for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
        }
        pSetCoeff( p, ncpoly[c] );
        pSetm( p );
        if (result!=NULL) result= pAdd( result, p );
        else result=  p;
      }
      c++;
    }
    sum=0;
    exp[0]++;
    for ( j= 0; j < n - 1; j++ )
    {
      if ( exp[j] > tdg )
      {
        exp[j]= 0;
        exp[j + 1]++;
      }
      sum+=exp[j];
    }
    sum+=exp[n-1];
  }
 
  pTest( result );
 
  return result;
}

interpolateDenseSP()

rootContainer ** uResultant::interpolateDenseSP	(	BOOLEAN	matchUp = false,
		const number	subDetVal = NULL
	)

Definition at line 2922 of file mpr_base.cc.

{
  int i,p,uvar;
  long tdg;
  int loops= (matchUp?n-2:n-1);
 
  mprPROTnl("uResultant::interpolateDenseSP");
 
  tdg= resMat->getDetDeg();
 
  // evaluate D in tdg+1 distinct points, so
  // we need tdg+1 results of D(p0,1,0,...,0) =
  //              c(0)*u0^tdg + c(1)*u0^tdg-1 + ... + c(tdg-1)*u0 + c(tdg)
  number *presults;
  presults= (number *)omAlloc( (tdg + 1) * sizeof( number ) );
  for ( i=0; i <= tdg; i++ ) presults[i]= nInit(0);
 
  rootContainer ** roots;
  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
 
  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
  for (i=0; i < n; i++) pevpoint[i]= nInit(0);
 
  number *pev= (number *)omAlloc( n * sizeof( number ) );
  for (i=0; i < n; i++) pev[i]= nInit(0);
 
  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
  // this gives us n-1 evaluations
  p=3;
  for ( uvar= 0; uvar < loops; uvar++ )
  {
    // generate initial evaluation point
    if ( matchUp )
    {
      for (i=0; i < n; i++)
      {
        // prime(random number) between 1 and MAXEVPOINT
        nDelete( &pevpoint[i] );
        if ( i == 0 )
        {
          //p= nextPrime( p );
          pevpoint[i]= nInit( p );
        }
        else if ( i <= uvar + 2 )
        {
          pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
          //pevpoint[i]=nInit(383);
        }
        else
          pevpoint[i]=nInit(0);
        mprPROTNnl(" ",pevpoint[i]);
      }
    }
    else
    {
      for (i=0; i < n; i++)
      {
        // init pevpoint with  prime,0,...0,1,0,...,0
        nDelete( &pevpoint[i] );
        if ( i == 0 )
        {
          //p=nextPrime( p );
          pevpoint[i]=nInit( p );
        }
        else
        {
          if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
          else pevpoint[i]= nInit(0);
        }
        mprPROTNnl(" ",pevpoint[i]);
      }
    }
 
    // prepare aktual evaluation point
    for (i=0; i < n; i++)
    {
      nDelete( &pev[i] );
      pev[i]= nCopy( pevpoint[i] );
    }
    // evaluate the determinant in the points pev^0, pev^1, ..., pev^tdg
    for ( i=0; i <= tdg; i++ )
    {
      nDelete( &pev[0] );
      nPower(pevpoint[0],i,&pev[0]);          // new evpoint
 
      nDelete( &presults[i] );
      presults[i]=resMat->getDetAt( pev );   // evaluate det at point evpoint
 
      mprPROTNnl("",presults[i]);
 
      mprSTICKYPROT(ST_BASE_EV);
      mprPROTL("",tdg-i);
    }
    mprSTICKYPROT("\n");
 
    // now interpolate
    vandermonde vm( tdg + 1, 1, tdg, pevpoint, FALSE );
    number *ncpoly= vm.interpolateDense( presults );
 
    if ( subDetVal != NULL )
    {  // divide by common factor
      number detdiv;
      for ( i= 0; i <= tdg; i++ )
      {
        detdiv= nDiv( ncpoly[i], subDetVal );
        nNormalize( detdiv );
        nDelete( &ncpoly[i] );
        ncpoly[i]= detdiv;
      }
    }
 
#ifdef mprDEBUG_ALL
    PrintLn();
    for ( i=0; i <= tdg; i++ )
    {
      nPrint(ncpoly[i]); PrintS(" --- ");
    }
    PrintLn();
#endif
 
    // save results
    roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
                                loops );
  }
 
  // free some stuff: pev, presult
  for ( i=0; i < n; i++ ) nDelete( pev + i );
  omFreeSize( (void *)pev, n * sizeof( number ) );
 
  for ( i=0; i <= tdg; i++ ) nDelete( presults+i );
  omFreeSize( (void *)presults, (tdg + 1) * sizeof( number ) );
 
  return roots;
}

linearPoly()

poly uResultant::linearPoly ( const resMatType  rmt )

private

Definition at line 2743 of file mpr_base.cc.

{
  int i;
 
  poly newlp= pOne();
  poly actlp, rootlp= newlp;
 
  for ( i= 1; i <= (currRing->N); i++ )
  {
    actlp= newlp;
    pSetExp( actlp, i, 1 );
    pSetm( actlp );
    newlp= pOne();
    actlp->next= newlp;
  }
  actlp->next= NULL;
  pDelete( &newlp );
 
  if ( rrmt == sparseResMat )
  {
    newlp= pOne();
    actlp->next= newlp;
    newlp->next= NULL;
  }
  return ( rootlp );
}

nextPrime()

int uResultant::nextPrime ( const int  p )

private

Definition at line 3173 of file mpr_base.cc.

{
  int init=i;
  int ii=i+2;
  extern int IsPrime(int p); // from Singular/ipshell.{h,cc}
  int j= IsPrime( ii );
  while ( j <= init )
  {
    ii+=2;
    j= IsPrime( ii );
  }
  return j;
}

specializeInU()

rootContainer ** uResultant::specializeInU	(	BOOLEAN	matchUp = false,
		const number	subDetVal = NULL
	)

Definition at line 3060 of file mpr_base.cc.

{
  int i,/*p,*/uvar;
  long tdg;
  poly pures,piter;
  int loops=(matchUp?n-2:n-1);
  int nn=n;
  if (loops==0) { loops=1;nn++;}
 
  mprPROTnl("uResultant::specializeInU");
 
  tdg= resMat->getDetDeg();
 
  rootContainer ** roots;
  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
 
  number *pevpoint= (number *)omAlloc( nn * sizeof( number ) );
  for (i=0; i < nn; i++) pevpoint[i]= nInit(0);
 
  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
  // p=3;
  for ( uvar= 0; uvar < loops; uvar++ )
  {
    // generate initial evaluation point
    if ( matchUp )
    {
      for (i=0; i < n; i++)
      {
        // prime(random number) between 1 and MAXEVPOINT
        nDelete( &pevpoint[i] );
        if ( i <= uvar + 2 )
        {
          pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
          //pevpoint[i]=nInit(383);
        }
        else pevpoint[i]=nInit(0);
        mprPROTNnl(" ",pevpoint[i]);
      }
    }
    else
    {
      for (i=0; i < n; i++)
      {
        // init pevpoint with  prime,0,...0,-1,0,...,0
        nDelete( &(pevpoint[i]) );
        if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
        else pevpoint[i]= nInit(0);
        mprPROTNnl(" ",pevpoint[i]);
      }
    }
 
    pures= resMat->getUDet( pevpoint );
 
    number *ncpoly= (number *)omAlloc( (tdg+1) * sizeof(number) );
 
#ifdef MPR_MASI
    BOOLEAN masi=true;
#endif
 
    piter= pures;
    for ( i= tdg; i >= 0; i-- )
    {
      if ( piter && pTotaldegree(piter) == i )
      {
        ncpoly[i]= nCopy( pGetCoeff( piter ) );
        pIter( piter );
#ifdef MPR_MASI
        masi=false;
#endif
      }
      else
      {
        ncpoly[i]= nInit(0);
      }
      mprPROTNnl("", ncpoly[i] );
    }
#ifdef MPR_MASI
    if ( masi ) mprSTICKYPROT("MASI");
#endif
 
    mprSTICKYPROT(ST_BASE_EV); // .
 
    if ( subDetVal != NULL )  // divide by common factor
    {
      number detdiv;
      for ( i= 0; i <= tdg; i++ )
      {
        detdiv= nDiv( ncpoly[i], subDetVal );
        nNormalize( detdiv );
        nDelete( &ncpoly[i] );
        ncpoly[i]= detdiv;
      }
    }
 
    pDelete( &pures );
 
    // save results
    roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
                                loops );
  }
 
  mprSTICKYPROT("\n");
 
  // free some stuff: pev, presult
  for ( i=0; i < n; i++ ) nDelete( pevpoint + i );
  omFreeSize( (void *)pevpoint, n * sizeof( number ) );
 
  return roots;
}

Field Documentation

gls

ideal uResultant::gls

private

Definition at line 88 of file mpr_base.h.

n

int uResultant::n

private

Definition at line 89 of file mpr_base.h.

resMat

resMatrixBase* uResultant::resMat

private

Definition at line 92 of file mpr_base.h.

rmt

resMatType uResultant::rmt

private

Definition at line 91 of file mpr_base.h.

---

The documentation for this class was generated from the following files:

• kernel/numeric/mpr_base.h
• kernel/numeric/mpr_base.cc