Loading...
Searching...
No Matches
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 };
uResultant::denseResMat
@ denseResMat
Definition mpr_base.h:65
uResultant::none
@ none
Definition mpr_base.h:65
uResultant::sparseResMat
@ sparseResMat
Definition mpr_base.h:65

## 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::linearPoly
poly linearPoly(const resMatType rmt)
Definition mpr_base.cc:2743
uResultant::rmt
resMatType rmt
Definition mpr_base.h:91
uResultant::resMat
resMatrixBase * resMat
Definition mpr_base.h:92
uResultant::extendIdeal
ideal extendIdeal(const ideal gls, poly linPoly, const resMatType rmt)
Definition mpr_base.cc:2715
uResultant::n
int n
Definition mpr_base.h:89
uResultant::gls
ideal gls
Definition mpr_base.h:88
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
idCopy
ideal idCopy(ideal A)
Definition ideals.h:60
IDELEMS
#define IDELEMS(i)
Definition simpleideals.h:23

## ~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 );
}
i
int i
Definition cfEzgcd.cc:132
omReallocSize
#define omReallocSize(addr, o_size, size)
Definition omAllocDecl.h:220

## 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;
}
pow
Rational pow(const Rational &a, int e)
Definition GMPrat.cc:411
l
int l
Definition cfEzgcd.cc:100
p
int p
Definition cfModGcd.cc:4086
uResultant::nextPrime
int nextPrime(const int p)
Definition mpr_base.cc:3173
Print
#define Print
Definition emacs.cc:80
result
return result
Definition facAbsBiFact.cc:76
j
int j
Definition facHensel.cc:110
over
unsigned long over(const unsigned long n, const unsigned long d)
Definition mpr_base.cc:2659
exp
gmp_float exp(const gmp_float &a)
Definition mpr_complex.cc:349
mprSTICKYPROT
#define mprSTICKYPROT(msg)
Definition mpr_global.h:54
mprPROTN
#define mprPROTN(msg, nval)
Definition mpr_global.h:48
ST_BASE_EV
#define ST_BASE_EV
Definition mpr_global.h:62
mprPROTnl
#define mprPROTnl(msg)
Definition mpr_global.h:42
mprPROTNnl
#define mprPROTNnl(msg, nval)
Definition mpr_global.h:49
nDiv
#define nDiv(a, b)
Definition numbers.h:32
nDelete
#define nDelete(n)
Definition numbers.h:16
nIsZero
#define nIsZero(n)
Definition numbers.h:19
nEqual
#define nEqual(n1, n2)
Definition numbers.h:20
nPrint
#define nPrint(a)
only for debug, over any initialized currRing
Definition numbers.h:46
nNormalize
#define nNormalize(n)
Definition numbers.h:30
nInit
#define nInit(i)
Definition numbers.h:24
nTest
#define nTest(a)
Definition numbers.h:35
nPower
#define nPower(a, b, res)
Definition numbers.h:38
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
NULL
#define NULL
Definition omList.c:12
pAdd
#define pAdd(p, q)
Definition polys.h:204
pTest
#define pTest(p)
Definition polys.h:415
pSetm
#define pSetm(p)
Definition polys.h:272
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition polys.h:32
pSetExp
#define pSetExp(p, i, v)
Definition polys.h:43
pOne
#define pOne()
Definition polys.h:316
PrintS
void PrintS(const char *s)
Definition reporter.cc:288
PrintLn
void PrintLn()
Definition reporter.cc:314

## 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;
}
FALSE
#define FALSE
Definition auxiliary.h:97
rootContainer::fillContainer
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition mpr_numeric.cc:300
rootContainer::cspecial
@ cspecial
Definition mpr_numeric.h:68
rootContainer::cspecialmu
@ cspecialmu
Definition mpr_numeric.h:68
MAXEVPOINT
#define MAXEVPOINT
Definition mpr_base.cc:2653
mprPROTL
#define mprPROTL(msg, intval)
Definition mpr_global.h:46
nCopy
#define nCopy(n)
Definition numbers.h:15
omFreeSize
#define omFreeSize(addr, size)
Definition omAllocDecl.h:260
siRand
int siRand()
Definition sirandom.c:42

## 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 );
}
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
pDelete
#define pDelete(p_ptr)
Definition polys.h:187

## 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;
}
LinTree::init
void init()
Definition lintree.cc:864
IsPrime
int IsPrime(int p)
Definition prime.cc:61

## 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;
}
BOOLEAN
int BOOLEAN
Definition auxiliary.h:88
pIter
#define pIter(p)
Definition monomials.h:37
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition monomials.h:44
pTotaldegree
static long pTotaldegree(poly p)
Definition polys.h:283

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